Asymmetric hybridization in Cordulegaster (Odonata: Cordulegastridae): Secondary postglacial contact and the possible role of mechanical constraints

Abstract Two Cordulegaster dragonflies present in Italy, the Palaearctic and northern distributed Cordulegaster boltonii and the endemic to the south of the peninsula Cordulegaster trinacriae, meet in central Italy and give rise to individuals of intermediate morphology. By means of mitochondrial and nuclear markers and of Geometric Morphometrics applied to sexual appendages, we defined i) the geographical boundaries between the two species in Italy and ii) we determined the presence, the extent, and the genetic characteristics of the hybridization. Genetic data evidenced asymmetric hybridization with the males of C. trinacriae able to mate both interspecifically and intraspecifically. The results contrast with expectations under neutral gene introgression and sexual selection. This data, along with the morphological evidence of significant differences in size and shape of sexual appendages between the males of the two species, seem indicative of the role of mechanical constraints in intraspecific matings. The origin of the two species is dated about to 1.32 Mya and the hybridization resulted related to range expansion of the two species after Last Glacial Maximum and this led to the secondary contact between the two taxa in central Italy. At last, our results indicate that the range of C. trinacriae, a threatened and protected species, has been moving northward probably driven by climate changes. As a result, the latter species is currently intruding into the range of C. boltonii. The hybrid area is quite extended and the hybrids seem well adapted to the environment. From a conservation point of view, even if C. trinacriae has a strong genetic identity, the discovery of hybridization between the two species should be considered in a future species management.


| INTRODUC TI ON
The evolutionary history of European biota is closely linked to postglacial expansion, which has structured the current genetic distribution of species (Hewitt, 2000;Schmitt, 2007). The Last Glacial Maximum (LGM), the latest event of global cooling, shaped the phylogeny and phylogeography of many extant populations, and the effects of the postglacial processes have been investigated thoroughly. The first reviews by Hewitt (1999Hewitt ( , 2000Hewitt ( , 2001 highlighted that most postglacial colonization events in animals and plants followed a similar pattern of expansion from refugia.
A common phenomenon in the model of expansion after the LGM contraction of the Northern Hemisphere is the formation of hybrid zones after secondary contact between populations (Hewitt, 2001;Schmitt, 2007). Although it is often difficult to determine whether secondary contact occurred in an area, most cases reported in the literature proved to be secondary contacts related to climatic fluctuations in the Quaternary (37%, reviewed by Barton & Hewitt, 1985; but see also Taberlet, Fumagalli, Wust-Saucy, & Cosson, 1998;Davis & Shaw, 2001;Currat, Ruedi, Petit, & Excoffier, 2008). The comparison between expansion patterns has revealed the presence in Europe of hybridization "hotspots" for thermophilous taxa (called suture zones, Remington, 1968;Hewitt, 2000 for a review; Stewart & Lister, 2001) and one of these is located in northern Italy (Schmitt, 2007;Taberlet et al., 1998). The Alps were a barrier for species dispersal (Schmitt, 2007;Slechtová, Bohlen, Freyhof, Persat, & Delmastro, 2004) and the south of Italy has been indicated as a glacial refugium for several taxa (Canestrelli, Sacco, & Nascetti, 2012;Schmitt, 2007;Taberlet et al., 1998). It is interesting that the existence of refugia in sub-Alpine xerothermic areas has also been reported (Bhagwat & Willis, 2008;Kaltenrieder et al., 2009).
When two allopatric populations come into secondary contact and the reproductive barriers fail, the formation of hybrid zones takes place, which can be more or less stable and extended (Hayashi, Dobata, & Futahashi, 2005;Sánchez-Guillén et al., 2014;Wolf, Takebayashi, & Rieseberg, 2001). Interspecific hybridization is recognized as one of the potential drivers of speciation (Abbott et al., 2013;Barton, 2001;Jiggins & Mallet, 2000;Mallet, 2007). In fact, introgression, that is, acquisition of new genes through hybridization and backcrossing (Andersson, 1949), could be an important source of new variability and subsequent evolution. Hybridization can be symmetric (bidirectional), if males and females of both species are involved in heterospecific mating, or asymmetric (unidirectional), if one sex is predominantly involved in heterospecific mating. This is reflected in differential transfer of genes (hybrids contain a discordant mitochondrial/nuclear pattern) and of specific matrilinear and patrilinear markers.
When a species colonizes a new area where a sister species is already present, a massive introgression of neutral genes between the intrusive and the local species may take place. A number of factors can promote unidirectional hybridization (Canestrelli et al., 2016;Currat et al., 2008;Grant & Grant, 1997;Wirtz, 1999) and these are different in case of selective/not selective processes. If interbreeding is not prevented or if selective processes are not present (including asymmetric mating preference), introgression occurs almost exclusively from the local to the intrusive species, irrespective of the relative densities of the two species (Currat et al., 2008).
In the presence of selective processes, mechanical and ethological barriers are among the main prezygotic causes of asymmetric hybridization. Mechanical incompatibility of structures involved in mating is an efficient prezygotic isolation mechanism during tandem formation and copulation. This is an important constraint to interspecific reproduction in Odonata (Monetti, Sánchez-Guillén, & Cordero Rivera, 2002;Paulson, 1974;Sánchez-Guillén, Wellenreuther, Cordero-Rivera, & Hansson, 2011). Under sexual selection, when females of one species are locally rare and cannot find suitable (i.e., conspecific) males, they eventually mate with males of the locally common species (Turgeon, Stoks, Thum, Brown, & McPeek, 2005).
Thus, if an immigrant species invades an area of a closely related species, the rare females mate with both abundant and rare males to enhance their fitness (Wirtz, 1999). These circumstances are in contrast to what occurs in the absence of selective processes (Currat et al., 2008). At last, Canestrelli et al. (2016) highlighted that, during range expansion, asymmetric hybridization can be influenced by sexual preference, as a character linked to dispersal propensity, and spatial sorting of traits could affect the direction and extent of gene exchange. In some cases, the tendency to mate assortatively could increase at the expansion front and differential introgression could be favoured (examples in Canestrelli et al., 2016).
From a morphological point of view, in extant hybrid zones relying on gene flow, individuals may show admixture of morphological characters or a spatial transition between hybridizing forms of the two parental species. In other cases, character displacement evolves as an adaptive response of one or both species to the hybridization and when the speciation has occurred the hybrid zone is expected to disappear over time (Barton & Hewitt, 1985).
The genus Cordulegaster (Odonata: Cordulegastridae) has an Holarctic distribution, comprises 29 species and four of these are present in Italy: C. heros Theischinger, 1979, C. bidentata Sélys 1843 and the two sister species C. boltonii (Donovan, 1807) and C. trinacriae Waterston, 1976. Cordulegaster boltonii, a species with Palaearctic distribution, is present in north-central Italy while C. trinacriae is endemic to Sicily and to the southern Italian regions Basilicata, Calabria, Campania, and Molise . The two species can easily be distinguished by morphological features (Boudot, 2001;Dijkstra & Lewington, 2006). Froufe, Ferreira, Boudot, Alves, & Harris (2013) confirmed the validity of C. trinacriae as a distinct species but investigated only a limited number of populations from Italy, and they did not include samples from the central part of the peninsula where the two sister species are likely to be in contact.  presumed that C. boltonii and C. trinacriae can be sympatric in central Italy, and both species have been caught in Gerano (Rome) (Galletti & Pavesi, 1985).
However, no detailed data are available from the contact area.
In the genus Cordulegaster, one putative hybrid between C. bilineata and the sister species C. diastatops has recently been identified with molecular tools (Pilgrim, Roush, & Krane, 2002). Therefore, it would be important to investigate the contact zone between C. boltonii and C. trinacriae, as the latter species has been listed as "near threatened" in the Italian and the European IUCN Red List (Kalkman et al., 2010;. It is also listed in Annexes II and IV of the Habitats Directive (Cardoso, 2012) along with 16 other species of Odonata, making this taxonomic group the primary invertebrates in freshwater conservation (Dijkstra & Kalkman, 2012).
The aim of this study was to define the geographical boundaries of both species in Italy and to determine if they hybridized in the contact area. Based on the results, we discuss (a) the dynamics of postglacial expansion of the two species after the LGM, that led to the secondary contact between the taxa in central Italy; (b) the hybridization in relation to morphology and sexual selection; (c) the nature and the extension of the hybrid area within a conservation perspective.

| MATERIAL S AND ME THODS
In total, 93 male dragonflies belonging to the putative species C. boltonii and C. trinacriae were collected between 2008 and 2012 at 26 localities in the Italian peninsula (Table 1, Supporting Information   Table S1 and Figure 1a). The specimens were preliminarily attributed to one of the above species according to the morphological characters given by Boudot (2001) and Dijkstra & Lewington (2006).  Table 1). The left superior appendage and the inferior appendage of all individuals were photographed with a Hitachi TM1000 environmental scanning electron microscope (ESEM). Leg tissue samples of all individuals were collected for the DNA analysis and preserved in 100° ethanol.

| Molecular data analysis (i) Phylogeography, phylogeny and hybridization
The DNA was extracted from dissected metafemoral muscles of all 93 individuals and of two individuals of C. bidentata (used as outgroups in the phylogenetic analysis) according to the salting out procedure (Aljanabi & Martinez, 1997). Two markers were amplified: a fragment of the mitochondrial (mt) cytochrome c oxidase subunit I gene (coI) and the elongation factor 1α (ef1α) nuclear gene (nu). The mt DNA is a non recombining marker, easily tracked across nuclear backgrounds (Currat et al., 2008;Excoffier, Foll, & Petit, 2009;) and therefore, useful to detect asymmetric introgression. The use of two markers with distinct modes of inheritance provides insight into the degree of reproductive isolation between lineages and lead into particular isolating mechanisms (Harrison & Bogdanowicz, 1997;Jiggins & Mallet, 2000;Wirtz, 1999), especially in terms of mate choice-based introgression asymmetries. A 570 bp fragment at the 5′ terminus of the coI gene was amplified with the universal barcode primers LCO-1490 and HCO-2198 (Folmer, Black, Hoeh, Lutz, & Vrijenhoek, 1994) under standard PCR conditions. Both forward and reverse primers were used to sequence the fragment in double strand. For the ef1α, we amplified 550 bp using the primers EF1Fa and EF1Ra described in Pilgrim & Von Dohlen (2008).
The Bayesian Inference (BI) approach was used to construct the coI phylogenetic trees. A Generalized Time-Reversible model with a proportion of invariable sites and heterogeneous substitution rates following a gamma distribution (GTR + I + G, Rodríguez, Oliver, Marín, & Medina, 1990) was selected as the best substitution model under the AIC criterion as implemented in JModelTest (Posada, 2008). The Bayesian phylogenetic analysis was performed under the selected substitution model and was carried out with MrBayes v. 3.2.1 (Huelsenbeck & Ronquist, 2001) by running 1,000,000 generations, with Markov chains sampled every 1000 generations. The convergence of the parameters was determined with the software TRACER 1.6 (Rambaut, Suchard, Xie, & Drummond, 2014). A burn-in of 10% was applied and the remaining trees were used to compute a 50% majority rule consensus tree and posterior probabilities. To reconstruct coI phylogeography within the two species we built a parsimony network on the entire coI dataset with TCS v. 1.21 (Clement, Posada, & Crandall, 2000), in default settings, calculating connection limits at 30 steps. The time of the most recent common ancestor (TMRCA) was estimated using BEAST v. 1.8 (Drummond, Suchard, Xie, & Rambaut, 2012).
Given the lack of fossil records for Cordulegastridae and its related families Chlorogomphidae and Neopetalidae, as calibration points we used an average value of the coI substitution rate in a range between 1.5% and 2.3%, values estimated for mitochondrial DNA in insects (Brower, 1994;Farrell, 2001;Quek, Itino, & Pierce, 2004).
Thus, we applied a lognormal distribution (μ = 0.019; SD = 0.15) to include the maximum and minimum substitution rate values. The analysis was performed three times independently, with 100 million generations and sampling of trees every 10,000 steps using a Yule process as tree model. TRACER 1.6 (Rambaut et al., 2014) was used in all the BEAST analyses to check for parameter convergences and stationarity.
TA B L E 1 Sample details. Number of the locality on the map (map); sampling localities (Region, Province); number of individuals per locality (n); mitochondrial lineages present at each locality (coI clade); nuclear (ef1-α) genotype for individuals at each locality and number of individual assigned to each category (n) when not all the individuals are assigned to the same one. In the column, the categories stand for HOMb: homozygous C. boltonii; HOMt: homozygous C. trinacriae; Hz: heterozygous. The columns Inferior appendage and Superior appendage show typical appendages for C. boltonii (localities 1-6); a hybrid (localities 7-11) and C. trinacriae (localities 12-26) To evaluate the relations among the nuclear haplotypes, we used the phased haplotypes to construct a parsimony network with TCS v. 1.21 (Clement et al., 2000), in default settings, calculating connection limits at 95%. To estimate gene flow patterns between C. boltonii and C. trinacriae, we used the Isolation-with-Migration model as implemented in IMa2 (Hey, 2010). The two-populations analysis was conducted constraining the nuclear DNA dataset according to the mtDNA information and therefore, treating the hybrids on the basis of their respective mtDNA assignment. Preliminary short runs were conducted to determine the most appropriate upper bound in order to incorporate the posterior distribution for each parameter. After obtaining the fitting scalars (-t15, -q30, -m3), we performed a final M mode analysis using the HKY model for both genes, with 20 independent Markov-coupled chains with a geometric heating scheme (-ha0.96 and -hb0.9) as suggested for small to medium-sized data sets (Hey, 2011). In total, 100,000 genealogies were retained and mixing properties among chains were checked by inspection of low parameter autocorrelations, high swap rates and no trend in the splitting time plots (Hey, 2010).

| Molecular data analysis (ii) Demography and spatial analysis
The mode and time of population size changes for both species were inferred using the multilocus Extended Bayesian Skyline Plot (EBSP, Drummond et al., 2012). To provide the most reliable demographic reconstruction for each species, the analysis was carried out only on individuals identified as C. boltonii or C. trinacriae by both molecular markers, therefore, we excluded the hybrid specimens from the analysis. The final run was carried out with substitution models, clock rates and tree models for the two genes kept unlinked. At last, we F I G U R E 1 (a) Study area. The map reports the sampling localities (black dots, the numbers refer to the localities reported in Table 1) (Lemey, Rambaut, Welch, & Suchard, 2010) to identify the ancestral areas and from them the process of spatial diffusion. After fine-tuning the prior parameters, we ran a final analysis using a relaxed random walk diffusion model (RRWs) with Cauchy distribution (Lemey et al., 2010). Each species was processed separately with application of a strict clock model and a Gaussian Markov Random Field (GMRF) Bayesian Skyride (Minin, Bloomquis, & Suchard, 2008) as a coalescent tree prior. After computing the MCC with TreeAnnotator 1.7.5, we visualized the processes of spatial diffusion from ancestral areas using the Time Slicer option in SPREAD v. 1.0.5 (Bielejec, Rambaut, Suchard, & Lemey, 2011). We performed both the demographic inference and the spatial diffusion reconstruction in BEAST v. 1.7 (Drummond et al., 2012), running 300 million generations, sampling parameters every 30,000 generations and providing the coI substitution rate setting as previously described for the TMRCA analysis.

| Geometric morphometrics
The two diagnostic characters, the inferior appendage and the left superior appendage, were photographed with a Hitachi TM1000 environmental scanning electron microscope (ESEM) for all 93 individuals (example figures in Figure 2 and in Table 1). The form of these appendages is used to distinguish the two species (Boudot, 2001;Dijkstra & Lewington, 2006). In male Odonata they interlock with parts of the female body during tandem formation (Corbet, 1999) and can be mechanically incompatible with sister species resulting in and efficient prezygotic isolation mechanism (Monetti et al., 2002). When the left superior appendage was damaged or missing, a mirror image of the right superior appendage was used. The geometric morphometrics approach, based on the Procrustes method (Bookstein, 1991), was used to analyze the morphological inter-intraspecific variation. To define the position of the landmarks on curved structures in an unambiguous way, we superimposed a "comb" on each image. The extremities of the comb were positioned at homologous points (included as landmarks in the analysis) and the semi-landmarks were selected at the intersection between the teeth of the comb and the margin of the structure ( Figure 2). Fifteen landmarks were digitized on the inferior appendage and 19 landmarks on the superior appendage using tpsDig2 (Rohlf, 2006). We had planned to investigate the form of the entire left superior appendage, but the inner margin proved to be highly variable in response to slight differences in orientation, and therefore only the outer margin was digitized. A Generalized Procrustes Analysis (GPA, Rohlf & Slice, 1990) was performed to convert the landmarks into separate sets of variables for size (centroid size, CS) and shape (Procrustes coordinates).
The overall shape variability was analyzed by Principal Component Analysis (PCA) and graphically evaluated through deformation grids.
A second PCA analysis was performed on a subset of individuals identified as C. boltonii and C. trinacriae by both molecular markers (see Table 1, Supporting Information Table S1). Canonical Variates Analysis (CVA) with a permutation test was applied to the same data set to test the presence of significant shape differences between the two species [evaluated by Multivariate Analysis of Variance (MANOVA)].
Size differences among C. boltonii, C. trinacriae, and hybrid individuals were investigated by ANOVA (Tukey test) performed on the CS and graphically represented by a box-plot. A regression between the shape variables and CS was carried out to investigate the influence of size on shape. At last, we carried out a regression of Procrustes coordinates on the latitude values to test if the morphological variability depended on geographical/ecological factors. In fact, the distribution of both species covers a large latitudinal range and we found that the differences in the shape of appendages were driven by the differences in size (see Results). Therefore, it was important to test whether this variation in shape depends on ecological factors, as size is known to be influenced by latitude in Odonata (Johansson, 2003;Sniegula, Golab, & Johansson, 2016). The analyses were carried out on the inferior and superior appendages separately with the software MorphoJ 1.06b (Klingenberg, 2011), with the exception of ANOVA performed with PAST 3.11 (Hammer, Harper, & Ryanet, 2001) and MANOVA performed in R (R Development Core Team, 2011).

| RE SULTS
The 111 new sequences (21 coI haplotypes and 90 ef1α genotypes) of Italian Cordulegaster are deposited in GenBank, Accession Numbers MH304646 -MH304756 (see Supporting Information Table S1 for details).
The Bayesian consensus tree for coI is reported in Figure 1b (only posterior probability (p.p.) values exceeding 70% are shown).
The estimated ESS values were greater than 200 for all parameters.

Glu Gly
(1) (1) (  These hybrids also include the two individuals from Macerata which were morphologically attributed to C. boltonii (tree in Figures 1 and   3; Table 1).
The migration parameters estimated using the IMa2 analysis

| Geometric morphometrics
The PCA of the superior appendage is reported in Figure 2a The analysis of size of both the superior and inferior appendages (Figure 2a,b)

| Preliminary overview
Emerging evidence reveals that postglacial routes of animals and plants in the Italian Peninsula followed complex species-specific patterns (Chiocchio, Bisconti, Zampiglia, Nascetti, & Canestrelli, 2017;Mouton et al., 2012;Provan & Bennett, 2008;Schmitt, 2007). The common trend is that northward range shifts initiated from southern refugia (Hewitt, 1999(Hewitt, , 2000(Hewitt, , 2001Taberlet et al., 1998). In contrast, for Cordulegaster we found that the two species survived in distinct refugia, two located in the south (C. trinacriae) and one in the north (C. boltonii), and then expanded until they reached a secondary contact zone in central Italy where they hybridized.
Cordulegaster boltonii in the north and C. trinacriae in the south occupy a large part of Italy .
However, the limits of their distributions were unknown, as was the area where the two species came into contact (Galletti & Pavesi, 1985). Individuals from the north are morphologically distinguishable from those from the south (Boudot, 2001;Dijkstra & Lewington, 2006), but we found that several specimens from central Italy show intermediate characters and this was the first empirical evidence of a possible hybridization between the two taxa.
The mitochondrial molecular analysis (Figure 1)

| Demography and postglacial expansion
Cordulegaster trinacriae and C. boltonii have an early Pleistocene origin, in the Calabrian age. Their separation is dated at ca. 1.32 Mya.
Considering C. boltonii having a Palaearctic distribution, this datum indicates a split from C. trinacriae in the Pleistocene and may reflect the geomorphological situation of the Italian Peninsula during this period. In fact, during the early Pleistocene the Italian Peninsula was characterized by a series of environmental discontinuities due to the presence of the Apennines and primed by repeated glacioeustatic sea-level fluctuations which isolated the populations along the west-east axis (Cremaschi, 2003a,b;Di Giovanni, Vlach, Giangiuliani, Goretti, & Torricelli, 1998;Giraudi, 2004;Stefani, Galli, Crosa, Zaccara, & Calamari, 2004).
The demographic analysis showed exponential growth of the populations of both species at the end of the LGM (Figure 4b,c).
The beginning of population expansion (dated between 50,000 and 25,000 ya) coincides with the decline of the Ice Age and thus is a typical expansion following the postglacial ice retreat (Hewitt, 1999(Hewitt, , 2000(Hewitt, , 2001. Over the same time span, C. trinacriae showed a more rapid increase in population size, twice that of C. boltonii (Figure 4b,c). This recent expansion may also have been responsible for the slightly higher population variability ( Figure 1a) detected with the mitochondrial marker.
The spatial diffusion analysis revealed three main ancestral areas From its refugium in the north-western sub-Alpine region,

C. boltonii expanded first to the center and later to north-eastern
Italy, accordingly to the presence of other xerothermic micro-refugia located south of the Alps (Bhagwat & Willis, 2008;Kaltenrieder et al., 2009). For C. trinacriae, the spatial diffusion analysis (Figure 5a-c) supported the presence of two ancestral distributions, the first located between Campania, Basilicata and Calabria and the second in southern Calabria (Figure 5a). From the two refugia, the colonization proceeded first southward and then toward the center. The spatial data suggest that C. boltonii was present in central Italy earlier than C. trinacriae, indicating a recent contact between the two species.

| Asymmetric hybridization and morphological constraints
The Ima2 analysis (Figure 4a) revealed the presence of asymmetric gene flow from C. trinacriae to C. boltonii. Unidirectional hybridization has also been observed in the Odonata genera Enallagma and Ischnura (Sánchez-Guillén et al., 2011;Turgeon et al., 2005), and for Enallagma it was shown that male sexual appendages are primarily used to identify potential mates (McPeek, Shen, Torrey, & Farid, 2008). The reproductive isolation in damselflies is complex because it is caused by multiple mechanisms, including both premating and postmating barriers (Wellenreuther & Sánchez-Guillén, 2016). Often, in Odonata, the structures involved in mating can act as a premating barrier, resulting in morphological incompatibility during tandem formation and copulation (Monetti et al., 2002;Paulson, 1974;Sánchez-Guillén, Wellenreuther, & Cordero Rivera, 2012;Wellenreuther & Sánchez-Guillén, 2016). In the present study, we found evidence that the two species are morphologically well characterized in regard to superior and inferior appendages (Figure 2 and MANOVA). It is interesting that the hybrid individuals fall exclusively within the range of variability of C. boltonii. The morphological differences between the two species are not confined to shape, indeed the appendages of C. boltonii and those of the hybrid individuals are significantly smaller than those of C. trinacriae (box plots in Figure 2a,b). The analysis of allometry suggests that the change in appendages shape was driven by differences in their size. In keeping with this hypothesis, no correlation of shape with latitude was found. This excludes geographical/ ecological causes of the differences in shape and confirms the genetic bases of the appendages morphology.  (Figure 5a,b,c).
Furthermore, the data show that C. trinacriae was the species with faster population growth, invading the other one in the contact area ( Figure 3). In the perspective of neutral gene introgression during expansion, the direction of hybridization is supposed to be the opposite (Currat et al., 2008), that is, from the local (C. boltonii) to the invading species (C. trinacriae). On the contrary, in the perspective of sexual selection (Wirtz, 1999) the invading females can mate with both local and immigrant males to enhance their fitness, that is, from the invading females (C. trinacriae) to the local species males (C. boltonii). In our case, we found introgression from the invading males (C. trinacriae) to the local females (C. boltonii). In fact, C. trinacriae is the advancing species but it has longer sexual appendages.
Hence, those findings may suggest a role of sexual appendages size in constrain the direction of intraspecific matings (Wirtz, 1999), as seen also in Sánchez- Guillén et al. (2011). Experimental work on these Cordulegaster species would be necessary to ultimately confirm the effectiveness of this reproductive barrier.

| Hybrid area extension and conservation implications
The geographical distribution of nuclear haplotypes (Figures 1b and   3) indicates that the hybridization covers an extensive area stretching the entire width of the central Italian Peninsula. Such a wideranging introgression has also been observed in another species of Odonata (Hayashi et al., 2005). The hybrid zone extends from Marche to southern Latium (localities 7-11). We can set Marche as the northern boundary because of the presence of the two individuals (MC.SaG.3 and MC.SaG.4, locality 7) which are mitochondrially and morphologically C. boltonii whereas their nuclear genome was found to be hybrid, probably of a generation following the F1. In the case of secondary contact, it is important if one of the species involved is endangered, in which case it is necessary to establish contingent consequences for its conservation and to formulate specific management policies (Genovart, 2009).
Our results clearly indicate that C. trinacriae has been moving northward in the last few thousand years (Figure 5a,b,c), long before the onset of climate change, and this southern species is intruding into C. boltonii's territory. The direction of the range shift is essential from a conservation point of view, as C. trinacriae, the expanding taxon, is an endangered and protected species (see above).
When a rapidly-expanding species invades the territory of a sister species creating well-adapted hybrids, the hybridization area may be destined to be stable and might even extend further. A plausible outcome could be disappearance of one or both of the parental species, with a tendency toward the mixing of the genomes and the spreading of intermediate morphological characters along the entire peninsula (Sánchez-Guillén et al., 2014;Wolf et al., 2001). From a conservation point of view, even if C. trinacriae currently has a strong genetic identity, the presence of extended hybridization between the two species, should be considered relevant for future species management.

| CON CLUS IONS
This study assessed the presence of asymmetric hybridization after postglacial secondary contact between the two sister species and is in line with the results of recent studies on various animal species showing that the postglacial patterns can show numerous facets (i.e., Alpine xerothermic refuge), even though they follow the general dynamics described by Hewitt (2000).
The use of two loci over the entire genome imposed clear limit to our study. The use of several unlinked genes or even better of genome-wide loci is highly desirable to have the most accurate inference about hybridization and demography. However, the present study raises interesting questions for future and experimental works. In fact, studies on local dispersal in Odonata indicate that the animal's wing size and load are related to mobility, and these aspects are predictive of their flight ability (Conrad et al., 2002;Grabow & Ruppell, 1995). Therefore, if asymmetric hybridization due to sexual preference, is linked to dispersal propensity (i.e., dispersal syndrome, Travis & Dytham, 2002;Shine, Brown, & Phillips, 2011), it can be hypothesized and tested in future, that the mating of C. trinacriae males-with a strong propensity for dispersal-with females of both species could be the outcome of spatial sorting of multiple traits evolved during the range expansion.
At last, as the larvae of the Cordulegaster species live in vulnerable habitats, they can be an excellent model to study the effects of climate change on biodiversity and on the fate of these particularly stenotopic species. Dr. Gloria Antonini contributed to the data acquisition, analysis, and interpretation and coordinated the work of all co-authors.

DATA ACCE SS I B I LIT Y
The coI and ef1α sequences are deposited in GenBank, Accession Numbers MH304646 -MH304756.