From refugia to contact: Pine processionary moth hybrid zone in a complex biogeographic setting

Abstract Contact zones occur at the crossroad between specific dispersal routes and are facilitated by biogeographic discontinuities. Here, we focused on two Lepidoptera sister species that come in contact near the Turkish Straits System (TSS). We aimed to infer their phylogeographic histories in the Eastern Mediterranean and finely analyze their co‐occurrence and hybridization patterns in this biogeographic context. We used molecular mitochondrial and nuclear markers to study 224 individuals from 42 localities. We used discordances between markers and complementary assignment methods to identify and map hybrids and parental individuals. We confirmed the parapatric distribution of Thaumetopoea pityocampa (Lepidoptera: Notodontidae) in the west and Thaumetopoea wilkinsoni in the east and identified a narrow contact zone. We identified several glacial refugia of T. wilkinsoni in southern Turkey with a strong east–west differentiation in this species. Unexpectedly, T. pityocampa crossed the TSS and occur in northern Aegean Turkey and some eastern Greek islands. We found robust evidence of introgression between the two species in a restricted zone in northwestern Turkey, but we did not identify any F1 individuals. The identified hybrid zone was mostly bimodal. The distributions and genetic patterns of the studied species were strongly influenced both by the Quaternary climatic oscillations and the complex geological history of the Aegean region. T. pityocampa and T. wilkinsoni survived the last glacial maximum in disjoint refugia and met in western Turkey at the edge of the recolonization routes. Expanding population of T. wilkinsoni constrained T. pityocampa to the western Turkish shore. Additionally, we found evidence of recurrent introgression by T. wilkinsoni males in several T. pityocampa populations. Our results suggest that some prezygotic isolation mechanisms, such as differences in timing of the adult emergences, might be a driver of the isolation between the sister species.

We used molecular mitochondrial and nuclear markers to study 224 individuals from 42 localities. We used discordances between markers and complementary assignment methods to identify and map hybrids and parental individuals.
We confirmed the parapatric distribution of Thaumetopoea pityocampa (Lepidoptera: Notodontidae) in the west and Thaumetopoea wilkinsoni in the east and identified a narrow contact zone. We identified several glacial refugia of T. wilkinsoni in southern Turkey with a strong east-west differentiation in this species. Unexpectedly, T. pityocampa crossed the TSS and occur in northern Aegean Turkey and some eastern Greek islands. We found robust evidence of introgression between the two species in a restricted zone in northwestern Turkey, but we did not identify any F 1 individuals.
The identified hybrid zone was mostly bimodal.
The distributions and genetic patterns of the studied species were strongly influenced both by the Quaternary climatic oscillations and the complex geological history of the Aegean region. T. pityocampa and T. wilkinsoni survived the last glacial maximum in disjoint refugia and met in western Turkey at the edge of the recolonization routes. Expanding population of T. wilkinsoni constrained T. pityocampa to the western Turkish shore. Additionally, we found evidence of recurrent introgression by T. wilkinsoni males in several T. pityocampa populations. Our results suggest that some prezygotic isolation mechanisms, such as differences in timing of the adult emergences, might be a driver of the isolation between the sister species.

K E Y W O R D S
Aegean Sea, asymmetric introgression, natural hybridization, secondary contact, Thaumetopoea pityocampa, Thaumetopoea wilkinsoni

| INTRODUC TI ON
Climate and habitat changes can facilitate range movements, which can cause secondary contacts between species or lineages (Taylor, Larson, & Harrison, 2015). The Quaternary glacial cycles have strongly influenced the current distributions of Mediterranean species and their genetic diversity (Hewitt, 1999;Schmitt, 2007). Glacial periods occurred ca. every 100,000 years and forced species to contract their ranges into restricted refugia, while interglacial periods allowed them to expand northward (the Expansion-Contraction model) (Taberlet, Fumagalli, Wust-Saucy, & Cosson, 1998). In regions strongly affected by glaciations, most of the genetic diversity has accumulated in lower latitudes (Petit et al., 2003), while northern populations show decreased allelic richness and signs of demographic expansions, a pattern known as "southern richness and northern purity" (Hewitt, 1999). Recolonization routes of the species during interglacials were shaped by major biogeographic factors such as barriers and corridors, and by biotic interactions in some cases. The glacial cycles have also affected sea levels and land configurations, thereby modifying landmass connectivity, and thus possible dispersal routes, over time (Hewitt, 2011). The genetic footprints of species' responses to these successions of climate changes have been extensively studied for many species in Europe and North America, while studies in the Near East remain scarce.
The biogeographic dynamics of the Mediterranean region were affected by two main geologic events: (a) the opening of the Mid-Aegean Trench (MAT) (~12 Mya), which separated Greek and Turkish peninsulas, and (b) the "Messinian Salinity Crisis" (5.9-5.3 Mya), a major drop in sea level that allowed several landmasses to emerge and connect previously isolated islands (Krijgsman et al., 1999;Poulakakis et al., 2015). The Aegean region is located on a crossroad between continents where contacts between divergent lineages moving from Europe to Anatolia or vice versa are likely (Dubey et al., 2007). Various biogeographic barriers such as the Turkish Straits System (TSS) (Dardanelles-Sea of Marmara-Bosphorus) and the Anatolian Diagonal (Figure 1a), a mountain range extending from southern to northeastern Turkey, make the Eastern Mediterranean a promising area for studying suture zones. Indeed, an increasing number of studies in Turkey have identified secondary contact zones (Bilgin, 2011) and showed that Anatolia was a major glacial refugium for many organisms (e.g., Biltekin et al., 2015;Korkmaz, Lunt, Çıplak, Değerli, & Başıbüyük, 2014;Mutun, 2010). Therefore, biogeographic consequences of the particularly complex geological history and dynamics of landmass configuration in the Aegean region (e.g., Poulakakis et al., 2015) render it a good candidate to study possible past and current hybridization events between lineages or sister species, which is the general frame of the present work.
Most hybrid zones arise from secondary contact in so-called "Suture Zones" (Hewitt, 1999(Hewitt, , 2011, which usually occur at the crossroad between specific dispersal routes and are facilitated by biogeographic discontinuities. They correspond to relatively narrow regions, where gene flow between related taxa leads to recombination of parental species' alleles or to discordances between nuclear and mitochondrial genomes. This reveals the porous nature of the genome and semipermeability of species boundaries (Harrison & Larson, 2014). An informative way of considering hybrid zones is to determine where they lie in the continuum from uni-to bimodality (Harrison & Bogdanowicz, 1997;Jiggins & Mallet, 2000) and to characterize their width, which brings information about the distance reached by introgressed genes (i.e., genes of one species included in the genome of the other species).
Unimodal hybrid zones correspond to regions where most individuals have intermediate genotypes compared to the two parental species in contact; these intermediate genotypes have also been called "hybrid swarms" (e.g., Scriber & Ording, 2005). On the other hand, bimodal hybrid zones consist largely of genotypes resembling the parental forms with few intermediates (Jiggins & Mallet, 2000). They usually indicate that speciation of parental forms is nearly complete and is strongly associated with assortative mating or fertilization failure, rather than postzygotic isolation (Jiggins & Mallet, 2000). Both uni-and bimodality can occur when two allopatric sister species come in secondary contact, or at the border of the distributions of in situ parapatric species.
A major postglacial demographic expansion occurred throughout Europe after the last glacial maximum, but its history in the eastern part of its range remains poorly documented, except in Greece (Korsch et al., 2015). Previous studies revealed four differentiated mitochondrial lineages within T. wilkinsoni: Crete, Cyprus, western Turkey, and eastern Turkey-Israel Simonato et al., 2007) (the population in Crete might be a separate species, see Petsopoulos et al., 2018). However, biogeographic patterns of T.
wilkinsoni in Turkey and western limits of its distribution have been largely unexplored.
A contact zone between T. pityocampa and T. wilkinsoni was recently discovered in western Turkey with some evidence of introgression (İpekdal, Burban, Kerdelhué, & Çağlar, 2015). However, this result was based on sparse sampling and a limited number of mitochondrial and nuclear markers. Petrucco-Toffolo et al. (2018) further demonstrated the viability and fecundity of F 1 individuals without any sign of outbreeding depression in no-choice laboratory hybridization experiments. Apart from these two studies, there is no other study on PPM hybridization, details of which, therefore, have remained unknown so far.
The objectives of the present study were (a) to infer their local phylogeographic history and delineate their co-occurrence patterns, (b) to identify natural hybrids and explore the genetic legacy of interspecific gene flow between the two species, and (c) to determine how the biogeographic context of the Aegean region influenced species range shifts and formation of a hybrid zone.

| Sampling
We sampled 174 larvae from 27 localities in Greece, Bulgaria, Turkey, Cyprus, and Lebanon (major regions related to the suspected hybrid zone and regions that were ignored or insufficiently studied previously) between 2002 and 2012. Fifty specimens from İpekdal et al. (2015) were also included; thus, we obtained 224 individuals from 42 localities in total ( Figure 1; Table 1). Larvae (first-fifth instars) were collected from different host trees to avoid sampling siblings and stored at −20°C in 70% ethanol.

| Laboratory protocols
We extracted larval genomic DNA (full body or head) using DNeasy tissue kit (Qiagen) and sequenced the samples for an 810 bp fragment of the mitochondrial Cytochrome c Oxidase subunit I gene (COI) using the primers Jerry/Pat (Rousselet et al., 2010), and a 660 bp fragment of the nuclear Photolyase gene (Pho) using the primers from Simonato et al. (2013). MWG Company performed sequencing on an ABI PRISM 3730 genetic analyzer using BigDye Terminator chemistry (Applied Biosystems).
Whenever we obtained heterozygous sequences (16 individuals) for Pho, we cloned PCR products using pGEM T Easy Vector (Promega Corp.). We sequenced six to eight clones per individual to obtain the phased sequence of both alleles. Additionally, we characterized the Internal Transcribed Spacer 1 (ITS-1) by RFLP-PCR for all 224 individuals. We amplified the ITS-1 fragment using the primers from Vogler and DeSalle (1994) and digested it using the restriction enzyme Hga-1, which produces different banding patterns for T. pityocampa (6 bands of 152, 131, 97, 70, 33, and 17 bp based on the GenBank accession EF189684), and T. wilkinsoni (5 bands of 232, 153, 66, 33, and 17 bp, accession EF189687).
To check the consistency of results, we cloned and sequenced 13 individuals showing typical banding pattern for the species and 5 individuals with an ambiguous or hybrid-like banding pattern.
We checked sequences for stop codons and double peaks to avoid pseudogenes. We also used 12 microsatellite markers F I G U R E 1 (a) Distribution of the two processionary moth species in the Mediterranean basin (green:

| Analyses of sequence and RFLP data
We built COI haplotype and Pho allelic networks using TCS 1.21 (Clement, Posada, & Crandall, 2000), including T. pityocampa and T.
By examining the respective networks and corresponding p-distances, each haplotype/allele could be considered as characteristic of one or the other species. Taking also into account their ITS RFLP patterns, hybrid individuals were identified based on two criteria: (a) the co-occurrence of allele's characteristic of both species for at least one nuclear marker, or (b) the incongruence (mtnuc or nuc-nuc) between the specific diagnosis of the different markers.

| Population genetic structure
We performed a Principal Component Analysis (PCA) using the adegenet 1.4-2 package (Jombart, 2008) implemented in R (R Core Team, 2018). We assigned individuals to clusters according to their membership coefficients (q) using a Bayesian inference method implemented in STRUCTURE 2.3 (Pritchard, Stephens, & Donnelly, 2000) with Admixture Model and correlative allelic frequency option. We used 100,000 burn-in cycles followed by 100,000 MCMC simulations and performed 10 iterations for each K (number of clusters) from 1 to 6. We selected the K value that best fit the data by examining the curve of ln P(X|K) and calculating ΔK (Evanno, Regnaut, & Goudet, 2005). After checking consistency of results over all iterations, we used the greedy algorithm of CLUMPP 1.1.2 (Jakobsson & Rosenberg, 2007) to average results of all the runs for each K. Results were plotted using DISTRUCT 1.1 (Rosenberg, 2004

| Hybrid detection from microsatellites
To identify each individual as pure T. pityocampa, pure T. wilkinsoni, or hybrid (F 1 , F 2 , or backcross) from their microsatellite genotypes, we followed the strategy developed in Burban et al. (2016) using the field dataset together with simulated genotypes of known ancestry (see below). Briefly, identification of hybrid individuals is based on complementary methods whose respective performance is dependent both on the power of the data set available (type and number of markers, sampling size) and of the level of differentiation and reproductive isolation. Therefore, using a combination of methods and a case-specific evaluation of their performances through F I G U R E 2 Networks and geographic distribution of (a) COI haplotypes analyses of simulated genotypes is highly recommended (Marie, Bernatchez, & Garant, 2011;Vähä & Primmer, 2006). STRUCTURE is first used to delimit the clusters that are expected to hybridize (parental clusters). This first step then allows to infer the proportion of alleles inherited from each parental species trough h-index estimation using INTROGRESS (Buerkle, 2005) and is also useful to generate simulated genotypes of known ancestry. A complementary approach is provided by NewHybrids (Anderson & Thompson, 2002) that potentially discriminates specifically first (F 1 ) and second (F 2 and backcrosses) hybrid generations from parental individuals. Note that based on the first STRUCTURE analyses and PCA, we excluded individuals from southeastern Turkey, Cyprus, and Lebanon (localities 37-42) to develop this procedure of hybrid identification, because these sites are distant from the potential contact zone and correspond to genetically differentiated lineages of T. wilkinsoni.

Simulated dataset
We selected all individuals having STRUCTURE q-value (hereafter, and NewHybrids to support the interpretation of the field dataset.

STRUCTURE
We analyzed the field dataset together with the 6,000 simulated individuals and determined the range of individual q-values (q STR ) corresponding to sim Pit and sim Wil . Any field individual having a q STR value outside this range was considered as "nonparental," that is, hybrid.

INTROGRESS
We calculated the hybrid index (h-index) using the est.h function of the INTROGRESS package (Gompert & Buerkle, 2010). We used sim Pit and sim Wil as reference parental individuals to calculate h-index for field and simulated individuals, the values ranging from 0 (pure T. pityocampa ancestry) to 1 (pure T. wilkinsoni ancestry). Finally, hybrids were identified among the field individuals when their h-index fell out of the range of simulated parents.

NewHybrids
This Bayesian software estimates posterior probabilities (q NH ) of each individual to fall into one of the six genotypic classes described above. We first ran NewHybrids on the simulated dataset alone using Jeffreys-like priors for allele frequencies and mixing proportions. We used 5 iterations of this analysis to check consistency of results and obtain averaged q NH values across runs. We used the "majority assignment" method to assign each individual to the class corresponding to the highest q NH , as described in . We compared the known genetic class of each simulated individual to the corresponding assignment result and used Vähä and Primmer's (2006) measures of efficiency, accuracy, and overall performance to estimate the power of assignment of the method with our dataset. Calculations of these measures are as follows: Additionally, to characterize the potential resulting errors, we examined the type of misassignments obtained for each class of simulated genotype. We then ran NewHybrids as described above, using the field and simulated individuals together to avoid any bias due to unbalanced occurrence of each genotypic class.

| Mitochondrial and nuclear sequences and RFLP data COI
We found 51 COI haplotypes that corresponded to two distinct clades (between-group mean distance: 0.0971). When compared to the reference sequences, these two clades corresponded to the

Pho
10 Pho alleles were found among 217 successfully sequenced individuals. Three alleles were closer to the pityocampa reference efficiency = number of simulated F 1 correctly assigned to F 1 class for the given q∕ number of individuals in the simulated F 1 class accuracy = number of simulated F 1 correctly assigned to F 1 class for the given q∕ number of individuals assigned to F 1 class either correctly or incorrectly for the given q overall performance = efficiency∕accuracy (alleles 1-3) and corresponded to individuals from Europe and west-  (15), Tekirdağ (16), and İstanbul (17) (Figure 2b).
Apart from individuals that had alleles from two different species for at least one marker, we also found individuals having "taxonomical discordance" between their mitochondrial and nuclear sequences and/or microsatellite markers (mt-nuc), or among the nuclear sequences (nuc-nuc) studied (Table 2). All the individuals having mt-nuc discordance had a pityocampa mitochondrial haplotype (Table 2). Geographic distributions of the two species and their hybrids according to each marker are shown in Figure 5.

PCA
The first two axes of the PCA explained 7.5% and 3.4% of the total inertia, respectively. PC1 separated the two species.

STRUCTURE
From the analysis of the simulated data, we found that pure parental T. pityocampa and T. wilkinsoni individuals had q-values below 0.207 and above 0.809, respectively. Applying these thresholds to the field individuals, 11 individuals were assigned as hybrids (Table 2).

INTROGRESS
The h-index ranged from 0.000-0.183 for sim Pit and 0.838-1.000 for sim Wil . When applying these ranges for the field individuals analyzed together with the simulated individuals, 105 individuals were assigned as T. pityocampa, and 85 as T. wilkinsoni. Eleven individuals had an h-index outside of this range and were thus considered as hybrids.

NewHybrids
Efficiencies, accuracies, and overall performances were generally high for the simulated dataset, except for F 2 ( Table 3). The levels of misassignment for each class of simulated individuals are shown in Table 3. Among field individuals, 17 were assigned as hybrid (6 F 2 , 8 Bck Pit , and 3 Bck Wil , Table 2).
To sum up, we found a high but not strict consistency of results between the three methods. The 11 individuals identified as hybrids from INTROGRESS were also identified as such by NewHybrids, and 10 of them by STRUCTURE (Table 2). However, the methods used did not allow to test more complex ancestries, such as recurrent backcrosses.
All individuals considered as hybrid from these analyses were collected in a restricted region comprising the Turkish Strait System ( wilkinsoni range in the eastern part of the zone (Figure 5e).  (Poulakakis et al., 2015) or, simply, current genetic exchange facilitated by the short distance between the island and the continent . This is true for most of the islands included in the present study except for Chios where we found T. pityocampa, whereas T. wilkinsoni was on the nearest Turkish coast. Moreover, all the hybrid individuals found (14%; note that no hybrid F 1 were identified) were only in eight sampling sites restricted to western Turkey, delimiting a relatively narrow contact zone where the two species can hybridize.

| Contrasting phylogeographic patterns of the two sister species
Consistent with previous studies, we showed that T. wilkinsoni has three highly differentiated mitochondrial lineages (Cypriot, eastern, TA B L E 2 Genetic characteristic of individuals considered as hybrid either from discordance between species assignment of the Pho and ITS alleles and COI haplotypes (p: pityocampa, w: wilkinsoni), or from one of the hybrid detection methods based on microsatellite data  Simonato et al. (2007), is probably due to a founder effect and local drift of the maternal lineage    (Schmitt, 2007). Accordingly, our results suggest that the Anatolian Diagonal (Figure 1a), the north-eastern extension of the Taurus chain, has acted as a major geographic barrier separating the eastern and the western wilkinsoni lineages (Figure 2), as is documented for many taxa in Turkey (Gür, 2016). Within the western wilkinsoni lineage, the haplotype network revealed two different patterns: a star-like and a branched pattern. The haplotypes Following Bilgin (2011), we suggest that Anatolia should be considered as a collection of multiple small refugia (refugia within refugia, see Gómez & Lunt, 2006) rather than a single refugium. Our results exemplify how the diverse geography of Anatolia has influenced the response of T. wilkinsoni to Quaternary glacial cycles, with strong barriers (the Anatolian Diagonal), a diversity of southern local refugia (the Taurus Mountains), and genetically depauperate northern sites.
The recolonization by T. wilkinsoni of the northwestern part of its distribution was probably impeded near the Bosphorus by the presence of its sister species T. pityocampa. We found one single group of haplotypes for T. pityocampa, with a clear star-like network. This suggests that northeastern Aegean region would not be a long-lasting refugium for T. pityocampa but the edge of its expanded range. The main haplotype found here (haplotype-1, see Figure 2) is also the major one found all along from eastern France to Greece (Korsch et al., 2015;Rousselet et al., 2010). This confirms the extensive eastward postglacial range expansion of T. pityocampa throughout Europe. Previous studies showed that divergent lineages occur for this species in the Iberian Peninsula, Corsica and North Africa Rousselet et al., 2010). Our results confirm the low genetic diversity found elsewhere in Europe, and no new stable glacial refugium was identified in the easternmost part of its distribution. Nonetheless, we found some closely related but private haplotypes in south of the Dardanelles and in Lesvos island (as some already mentioned by Korsch et al., 2015), suggesting that T. pityocampa survived the last glacial period locally in this region. This scenario suggests that T. pityocampa is a recent invader in Anatolia sensu Poulakakis et al. (2015). It probably expanded its range in the Eastern Mediterranean in the Pleistocene, long after the MAT had formed.
As documented for a diversity of taxa (e.g., Poulakakis et al., 2015), T. pityocampa dispersed eastwards from its glacial refugia, circumventing the MAT using a terrestrial corridor, crossed the TSS and spread southward along the coast, finally reaching as south as Chios. Enlarging the study area over the Balkans might allow to determine whether rare haplotypes, such as those found in this study, would occur elsewhere in Europe.
In conclusion, we showed that T. wilkinsoni and T. pityocampa display a clear parapatric pattern and meet near the TSS in northwestern Turkey. This region corresponds to a hybrid zone between these two organisms, where European T. pityocampa recently extended its range and met the northern limit of the western wilkinsoni lineage after the last postglacial recolonization. Due to complex paleogeography of the region, where land connections have constantly moved, the TSS could successively act as a geographic barrier between Europe and Anatolia, as for several other terrestrial organisms (e.g., Arntzen & Wielstra, 2010). It could also act as a land bridge allowing species to circumvent the MAT and expand their ranges in either direction (Dubey et al., 2007;Gündüz et al., 2007), facilitating the formation of contact zones (Antoniou, Magoulas, Platis, & Kotoulas, 2013;Bilgin, 2011;Hewitt, 2011).

| The hybrid zone is restricted and mostly bimodal, suggesting limited gene exchange
Locations of contact zones are strongly influenced by biogeographic barriers, which eventually constrain hybrid individuals in certain territories (Barton & Gale, 1993;Kawakami, Butlin, Adams, Paull, & Cooper, 2008). Considering all markers, we here detected a total of 31 hybrid individuals (Table 2), all restricted to a narrow corridor among the Sea of Marmara in the north, Aegean Sea in the west, climatically unsuitable Central Turkey in the east, and İzmir TA B L E 3 Efficiency, accuracy, overall performance of NewHybrids to characterize parental and hybrid categories from simulated genotypes (lines) and their proportion to each category (columns) by using the majority assignment approach in the south ( Figure 5). We did not identify any F 1 , and only found individuals assigned as F 2 or backcross categories. A similar pattern was found in an Iris hybrid zone, in which F 1 individuals were extremely rare (Sung, Bell, Nice, & Martin, 2018). It should be noted that the analyses of the microsatellite data are not meant to explicitly identify hybridization events that would have occurred more than two generations ago. The individuals identified as backcross or F 2 by NewHybrids might therefore have a more complex ancestry.  Table 2 for details)  (Jiggins & Mallet, 2000). In the present study, we found mostly pure parents and backcrosses in sites where the corresponding parental species occurred; F 2 hybrids were only found in Burhaniye (26). Although a more extensive sampling could reveal other patterns, these features, taken together, show that the identified hybrid zone is bimodal sensu Jiggins and Mallet (2000).
Overall, the restricted width and general bimodality of the hybrid zone suggest that effective dispersal of genes between the two species is relatively limited. Effective barriers to gene flow, thus, might be at play in the hybrid zone and could be linked to either prezygotic (e.g., assortative mating) or postzygotic (e.g., selection against hybrids) isolation. A similar pattern of limited gene flow was found for the greater mouse-eared bat in the Thrace region (Furman, Emek, & Çoraman, 2018). In order to understand the nature of this limitation, hybrid fitness (see Bimova et al., 2011) and phenological characteristics of the two species (see Chunco, 2014 andOrding, Mercader, Aardema, &Scriber, 2010) should be specifically investigated in the contact zone.

| Discordance between markers reveals direction of gene flow and possible past geographic distributions
Discordance between markers corresponds to situations where a single individual is assigned to one species by some of the loci, whereas it bears allele(s) corresponding to the other species at least for one other marker. Discordances can be found either between nuclear markers (nuc-nuc discordance; i.e., traces of the nuclear genome of one species in a genomic background of the other), or between nuclear and mitochondrial markers (mt-nuc discordance; i.e., occurrence of a mitochondrial haplotype of one species in an individual assigned to the other species by all or a majority of the nuclear loci used. This genomic mosaic of introgression might emerge because of unequal exchange of genes and density-dependant processes (Harrison & Larson, 2014), in which genomic traces of the rarest species eroding with time through backcrosses (Currat, Ruedi, Petit, & Excoffier, 2008). Conflicting geographic patterns between mitochondrial and nuclear markers (mt-nuc discordance) in secondary contact zones can correspond to divergent patterns of gene flow between the two genomes, and thus to sex-biased dispersal as mtDNA is inherited from the maternal lineage in most animals. They can also sign past hybrid zone movement, the majority of nuclear markers shifting their range while a wake of mtDNA is left behind (Toews & Breslford, 2012), or introgression by the native species in an expanding one (Currat et al., 2008). In the present study, most hybrids were characterized by a general T. pityocampa genomic background, and some traces of T. wilkinsoni alleles for a minority of nuclear loci.  (Wielstra et al., 2017). Complex geographic history of this region, where ancient gulfs and lakes were once present in today's Menderes and Gediz Deltas (Hakyemez, Erkal, & Göktaş, 1999) and made the landscape more fragmented for forest insects, might have influenced past hybridization events and favored T. wilkinsoni on the mainland.

| CON CLUS I ON AND PER S PEC TIVE S
Our results allow to go far beyond previous knowledge about the evolution of the PPM sister species. They suggest that the two species diverged in allopatry, survived the successive glacial maxima in disjoint refugial areas, and met at the extreme edge of their recolonization routes during the last interglacial. Fine-scale genetic diversity patterns nonetheless suggest that T. pityocampa refugia were certainly present in Anatolia during last glacial period. Even if reproductive isolation probably due to premating barriers is found, hybridization seems to have played a major role in species displacement, in favor of T. wilkinsoni that constrained T. pityocampa to western Anatolia. The respective roles of differential ecological traits, competition, and genetic interference in shaping the present patterns remain to be elucidated. Hybrid zones were defined as natural laboratories of evolutionary studies by Hewitt (1988), or "windows on evolutionary process" by Harrison (1990). In the case study presented here, hybrid characterization will offer the possibility to explore the mosaic of introgression along the genome concerning selected adaptive genes and to detect genes involved in reproductive isolation, owing to the recent development of genomics resources in the PPM (Gschloessl et al., 2018;Leblois et al., 2018).
Monitoring the hybrid zone using systematic sampling will help us understand whether this pattern is stable over time and whether F 1 can ever be identified in nature. Factors that shape the dynamics of introgression and limit interspecific gene flow in the field (e.g., differential phenology) would be a complementary research topic.
Finally, deciphering the ecological niches of both species through species distribution modeling will shed light on the history of their past, current, and future ranges. . We finally thank the anonymous reviewers for their insightful comments and suggestions that allowed us to significantly improve the quality of the manuscript.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

AUTH O R CO NTR I B UTI O N
Author contributions: K.İ., C.B., and C.K. conceived the study; K.İ.
conducted the fieldwork in Turkey and Cyprus, and A.B. and C.K.
provided samples from other countries and collaborators; K.İ., C.B., and L.S. performed laboratory experiments; K.İ., C.B., and C.K. analyzed the data; K.İ. led the writing with assistance from C.K., C.B., and A.B. All authors read and approved the final manuscript.