Christian Burban
email: christian.burban@inra.fr

Susana Rocha

Raphaël Leblois

Jean-Pierre Rossi

Laure Saune

Manuela Branco

Carole Kerdelhue

Susana Roc

Raphaël Leblo

Laure Sauné

Carole Kerdelhué

From sympatry to parapatry: a rapid change in the spatial context of incipient allochronic speciation

Keywords: Allochronic speciation, Phenology, Sympatry, Parapatry, Hybridization, Thaumetopoea pityocampa

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

3

Introduction Automatic magic traits refer to traits subject to divergent selection, which at the same time contribute to non-random mating [START_REF] Servedio | Magic traits in speciation: 'magic' but not rare?[END_REF][START_REF] Nosil | Ecological speciation. Oxford series in ecology and evolution[END_REF]. Classical examples of such traits include phenology when assortative mating occurs via temporal isolation or habitat choice when mating occurs in the preferred habitat [START_REF] Nosil | Ecological speciation. Oxford series in ecology and evolution[END_REF]. They are involved in sympatric or parapatric differentiation processes possibly leading to speciation, because they cause both immediate reduction of gene flow and differential adaptation (Thibert-Plante and [START_REF] Thibert-Plante | Evolution of mate choice and the so-called magic traits in ecological speciation[END_REF]. While speciation studies have originally focused on the spatial context of the differentiation, this formal categorisation into sympatric versus allopatric differentiation was subsequently questioned (e.g. [START_REF] Butlin | Sympatric, parapatric or allopatric: the most important way to classify speciation?[END_REF]). Yet, studying how ecological and spatial patterns driving the differentiation interact all along the speciation process remains highly relevant (Dieckman et al. 2004;[START_REF] Rundle | Ecological speciation[END_REF][START_REF] Kisel | Speciation has a spatial scale that depends on levels of gene flow[END_REF][START_REF] Aguilée | Ecological speciation in dynamic landscapes[END_REF]. Some key information are however hardly accessible when considering past speciation events because only the final outcomes, potentially modulated by biogeographic changes, are available [START_REF] Gavrilets | Patterns of parapatric speciation[END_REF][START_REF] Bolnik | Sympatric speciation: models and empirical evidence[END_REF]. We here document a case study of incipient speciation showing a rapid evolution of spatiality with time, the diverging populations moving from sympatry to parapatry.

Allochrony refers to a situation where reproductive isolation is due to a shift in reproductive time and is recognized as a salient magic trait [START_REF] Taylor | The role of allochrony in speciation[END_REF]. In case of seasonal shift, the diverging populations can experience a significant shift in their environmental niche, as their life stages face different abiotic constraints (e.g., larvae experiencing hot vs. cold conditions, see Santos et al. 2011b). One of the best-documented cases of true allochronic speciation sensu [START_REF] Taylor | The role of allochrony in speciation[END_REF] is the winter pine processionary moth (PPM), Thaumetopoea pityocampa Denis and Schiffermüller (Lepidoptera: Notodontidae). The PPM is a well-known defoliator occurring in all Mediterranean countries [START_REF] Kerdelhué | Quaternary history and contemporary patterns in a currently expanding species[END_REF], causing damages to Pinus species. In particular, the PPM has been reported almost everywhere in the Iberian Peninsula, as far as pines were present. The diversity of landscapes and climates over the whole region proved to be suitable for PPM, except the top of the mountains and the most arid locations in southeastern Spain [START_REF] Roques | Climate warming and past and present distribution of the processionary moths (Thaumetopoea spp.) in Europe, Asia Minor and North Africa[END_REF]. It has a typical one-year life cycle, with adults emerging and reproducing in the summer (end of July to mid-September in Portugal). Larvae hatch one month after egg laying and develop during fall and winter (October to end of February). The larvae then bury in the soil for pupation until adult emergence the following summer. This life cycle corresponds to the typical "winter phenology" observed throughout the PPM range. A unique population that does not follow the typical life cycle was discovered in 1997 in Portugal in a pure maritime pine Pinus pinaster Ait. forest. This peculiar and localized population is characterized by a "summer phenology", as adult emergence and reproduction take place in spring, from mid-April to early July, and larval development is completed in summer [START_REF] Pimentel | Establishment and expansion of a Thaumetopoea pityocampa (Den. and Schiff.) (Lep. Notodontidae) population with a shifted life cycle in a production pine forest, Central-Coastal Portugal[END_REF]). It will hereafter be referred to as the "summer population" (SP) as opposed to all other known "winter populations" (WPs). Previous monitoring data suggested that a temporal gap separated the reproduction period of the two co-occurring populations, the duration of which varied between years from 2 to 4 weeks, ensuring reproductive isolation [START_REF] Santos | Genetic isolation through time: allochronic differentiation of a phenologically atypical population of the pine processionary moth[END_REF](Santos et al. , 2011a)). Experimental rearing demonstrated that each type of phenology is highly heritable [START_REF] Branco | Experimental evidence for heritable reproductive time in 2 allochronic populations of pine processionary moth[END_REF]. The SP was first discovered in a ca. 900 ha zone in the south of the Mata Nacional de Leiria (MNL). It is noteworthy that the SP co-occurred at that time in the MNL with typical WP individuals, some trees being 1 3 attacked by both populations the same year [START_REF] Santos | Genetic isolation through time: allochronic differentiation of a phenologically atypical population of the pine processionary moth[END_REF]. Intensive field surveys in the following years allowed to evidence and follow its range expansion using nest occurrence data and pheromone trapping [START_REF] Godefroid | Climate constrains range expansion of an allochronic population of the pine processionary moth[END_REF]; its present-day distribution, much larger than the initial one in 1997, is restricted to the coastal area of Central Portugal, and spans over a 80 × 20 km area encompassing the whole MNL (Fig. 1).

The two populations remain syntopic, develop on the same hosts and both preferentially consume one-year old needles. No differences between the pheromones of the two populations were evidenced [START_REF] Kerdelhué | Genetic diversity and structure at different spatial scales in the processionary moths[END_REF]). Yet, the phenological shift induced major changes in the ecological constraints faced by the SP. It notably experiences much higher temperatures during early larval development, a stress that negatively affects WPs [START_REF] Roques | Climate warming and past and present distribution of the processionary moths (Thaumetopoea spp.) in Europe, Asia Minor and North Africa[END_REF]. The SP proved to be adapted to these new climatic conditions, as demonstrated by its higher larval thermal tolerances compared to local WPs (Santos et al. 2011b). Still, species distribution modeling suggested that climatic factors, in particular summer temperatures, impeded the SP expansion eastwards, and restricted the SP distribution along the coast [START_REF] Godefroid | Climate constrains range expansion of an allochronic population of the pine processionary moth[END_REF]). Yet, as known for the typical WP over its whole range and notably in Portugal [START_REF] Roques | Climate warming and past and present distribution of the processionary moths (Thaumetopoea spp.) in Europe, Asia Minor and North Africa[END_REF], SP population size can reach high densities, experiencing outbreak situations as recorded when it was discovered [START_REF] Pimentel | Establishment and expansion of a Thaumetopoea pityocampa (Den. and Schiff.) (Lep. Notodontidae) population with a shifted life cycle in a production pine forest, Central-Coastal Portugal[END_REF].

Previous results showed that the SP was genetically differentiated from the local WP based on nuclear markers, while the two sympatric populations shared a main mitochondrial haplotype [START_REF] Santos | Genetic isolation through time: allochronic differentiation of a phenologically atypical population of the pine processionary moth[END_REF](Santos et al. , 2011a;;[START_REF] Burban | Evidence for low-level hybridization between two allochronic populations of the pine processionary moth, Thaumetopoea pityocampa (Lepidoptera: Notodontidae)[END_REF]. Recent coalescent analyses based on Rad-seq data [START_REF] Leblois | Deciphering the demographic history of allochronic differentiation in the pine processionary moth Thaumetopoea pityocampa[END_REF]) established a convincing demographic historical scenario of the differentiation. From a large and stable ancestral population, SP and WP would have diverged during the Little Ice Age few centuries ago, both populations then experiencing a severe bottleneck. The SP was later subjected to a second bottleneck event ca. 70 years ago, before a last increase in population size lead to the present-day situation. Current gene flow was found to be still limited but higher than in the first step of the differentiation. This last point is consistent with some striking observations relative to the monitoring of adult emergence coupled with genetic assignments [START_REF] Burban | Evidence for low-level hybridization between two allochronic populations of the pine processionary moth, Thaumetopoea pityocampa (Lepidoptera: Notodontidae)[END_REF]):

(1) Rare males trapped during the WP adult emergence period proved to belong to the SP, and were called "LateSP" individuals. This result suggested that the time of adult activity could partially overlap between the two populations.

(2) Some individuals resulted from recent hybridization events between populations. In particular, a F1 hybrid was found at the southern edge of the SP range, trapped during the expected temporal gap between the SP and the WP emergence periods. However, a main limitation of the previous works mentioned above is that the SP was never sampled over its whole distribution range, but was rather studied in a few sites along the coast. In particular, very few males were genotyped from the range limits, and the easternmost area was never included.

The present study first aimed to analyze the intra-population genetic diversity and to look for spatial genetic patterns within the SP and the WP. In particular, as the SP expanded rapidly in the last decades, we determined if any genetic structuration could be identified, which would inform about the process of dispersal. Second, we used a wide sampling to determine if hybridization or introgression was actually more frequent at certain points in space and time, concerning notably range margins and intermediate dates of capture. Third, we wanted to test if LateSP (i.e., genetically SP individuals reproducing during the WP mating season) were found each year, and to determine if the categorization proposed by [START_REF] Burban | Evidence for low-level hybridization between two allochronic populations of the pine processionary moth, Thaumetopoea pityocampa (Lepidoptera: Notodontidae)[END_REF] that LateSP could actually be clustered in two groups based on their phenology was confirmed. Last, we wanted to use a spatially explicit sampling to determine the fine-scale pattern of co-occurrence of the two populations, and to determine if the presence of the SP could affect that of the WP. To do so, as it is possible to use pheromone traps to catch PPM males and get information about phenology (as adults are short-liver, male emergence is a good proxy for the date of reproduction), spatiality and genetic assignation, we performed monitoring along two transects within the whole SP range, completed by large scale sampling of the WP over the country.

Materials and methods

Sampling design

Pheromone traps were set up in 2012 and 2014, to take into account annual fluctuations regarding demography and phenology. They were installed along two transects within the documented distribution of the SP, hereafter called the Leiria region, following [START_REF] Godefroid | Climate constrains range expansion of an allochronic population of the pine processionary moth[END_REF]. Both transects were parallel, ca. 100 km long and positioned 20 km from each other. The first one corresponded to 9 sites chosen along the seashore (coastal transect), while the second comprised 7 sites located close to the eastern limit of the SP range (inner transect) (Fig. 1, Table S1). In each site, 2 funnel traps baited with synthetic PPM pheromone dispensers were hung on 2 different trees separated by ca. 10 meters. The pheromone lures were replaced every 6 weeks, following the methodology described in [START_REF] Santos | Genetic isolation through time: allochronic differentiation of a phenologically atypical population of the pine processionary moth[END_REF]Santos et al. ( , 2011a)). The traps installed in mid-April each year were checked fortnightly until the end of September to monitor the temporal dynamics of the male emergences during the entire period of adult flight of the two populations (Tables S1 andS2). Sampled individuals 1 3

were kept in 95% ethanol, and classified a priori as belonging either to the SP or the WP populations according to their date of capture: individuals caught before mid-July were a priori coded as SP and all other individuals as WP.

PPM individuals were also sampled from six localities over Portugal, outside of the Leiria region. In Fundão, Caparica, Grândola and Tavira, we used pheromone traps to catch males from mid-April to September as described above (Fig. 1, Table S1). In Varges, we sampled L1 larvae from 26 nests in Pinus pinaster, while the larvae from Viseu were those from [START_REF] Burban | Evidence for low-level hybridization between two allochronic populations of the pine processionary moth, Thaumetopoea pityocampa (Lepidoptera: Notodontidae)[END_REF]. All sampled individuals were immediately stored in 95% ethanol and kept at -18 °C until genetic analysis.

Microsatellite genotyping

A fraction of the individuals sampled in each site was genotyped using 20 microsatellite markers. We used 24-30 individuals sampled in each of the six localities over Portugal. Concerning the two transects in the Leiria region, we used a sub-sample of males from each site and each date of collect over the whole trapping season both in 2012 and 2014 (see Table S2 for details). However, a bad preservation of some individuals, notably due to predation in the traps, precluded further genetic analyses in a few cases (in particular, all samples from Leirosa in 2012, and the WP samples from Praia Vieira in 2012).

DNA was extracted from thoracic muscles (males) or from the whole body (larvae) using the DNeasy Blood and Tissue Kit (Qiagen), following the manufacturer's protocol. DNA was eluted in 100 µl of AE buffer, and its concentration and quality were assessed using a NanoDrop spectrophotometer. Samples were genotyped using twenty microsatellite loci, namely MS-Thpit01, MS-Thpit03, MS-Thpit04, MS-Thpit05, MS-Thpit06, MS-Thpit07, Ppit09,Ppit20,Ppit30,Ppit31,Ppit33,Ppit34,Ppit35,Ppit38,Ppit39,Ppit43,Ppit46,Ppit47,Ppit48. Detailed protocols for amplification and genotyping are described in [START_REF] Rousselet | Characterization of five microsatellite loci in the pine processionary moth Thaumetopoea pityocampa (Lepidoptera Notodontidae Thaumetopoeinae)[END_REF], A 'Hara et al. (2012) and [START_REF] Sauné | Isolation, characterization and PCR multiplexing of 17 microsatellite loci in the pine processionary moth Thaumetopoea pityocampa (Lepidoptera, Notodontidae)[END_REF]. PCR amplifications were performed using the Multiplex PCR Kit (Qiagen) and products were run and detected on an ABI 3500XL Genetic Analyzer (Applied Biosystems) in the GenSeq technological platform (Montpellier). Allele sizes were scored using the size standard GeneScan-500 LIZ (Applied Biosystems) and the GeneMapper 5.1 software (Applied Biosystems).

Genetic analyses

Genotyping data were organized in two datasets, corresponding to two spatial scales. The whole dataset corresponded to the six Portuguese populations together with the males sampled along the two transects in 2012 and 2014. In that case, we grouped the individuals caught in the Leiria region according to their phenology and year of capture, and considered 4 groups, namely Leiria WP 2012, WP 2014, SP 2012and SP 2014. The Leiria dataset focused on the two gradients and the individuals trapped in each site in 2012 and 2014 were considered separately. We first used the whole dataset to test Hardy-Weinberg equilibrium, linkage disequilibrium and to estimate the proportion of null alleles for each locus and each site. Deviations from Hardy to Weinberg equilibrium were tested with 1000 permutations and 100,000 steps in the Markov chain, while linkage-disequilibrium was tested for all pairs of loci with 10,000 permutations, both using Arlequin 3.11 [START_REF] Excoffier | Arlequin (version 3.0): an integrated software package for population genetics data analysis[END_REF]. Expected and observed heterozygosities were estimated using the same software.

3

The proportion of null alleles was estimated for each locus and population using the FreeNA package [START_REF] Chapuis | Microsatellite null alleles and estimation of population differentiation[END_REF]. All file format conversions were done using Convert [START_REF] Glaubitz | CONVERT: A user-friendly program to reformat diploid genotypic data for commonly used population genetic software packages[END_REF]).

Population genetic structure

We performed Principal Component Analyses (PCA) to explore the population genetic structure from the microsatellite datasets, using the R package Adegenet 1.4-2 [START_REF] Jombart | Comparative population genetic study of two oligophagous insects associated with the same hosts[END_REF]. Individuals were assigned to genetic clusters using the Bayesian clustering method implemented in Structure 2.3.4 under the model with admixture [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF]. We set K, the number of clusters, varying from 1 to 10 when analyzing the whole dataset, and from 1 to the respective number of localities sampled for each transect (i.e. 9 for the coastal transect and 7 for the inner transect). We used a burn-in of 500,000 steps followed by 50,000 iterations of the Markov Chain. To assess the consistency of results, we performed 10 independent runs for each value of K and compared the obtained individual Q-values that were subsequently analysed with Clumpack [START_REF] Kopelman | Clumpak: A program for identifying clustering modes and packaging population structure inferences across K[END_REF]. The number of clusters K that best explains the data was chosen by examining the curve of Log P(X|K) and using the ΔK method [START_REF] Evanno | Detecting the number of clusters of individuals using the software structure: a simulation study[END_REF] implemented in Structure har- vester [START_REF] Earl | Structure harvester: a website and program for visualizing structure output and implementing the Evanno method[END_REF]. Population structure was analysed through pairwise Fst using the excluding null alleles correction implemented in FreeNA to correct for the positive bias induced by the presence of null alleles [START_REF] Chapuis | Microsatellite null alleles and estimation of population differentiation[END_REF]. Bootstrapping 1000 times over loci provided the 95% confidence intervals.

Isolation by distance (IBD) results in an increase in genetic differentiation with geographic distance due to spatially limited dispersal between parents and offspring [START_REF] Wright | Isolation by distance[END_REF][START_REF] Rousset | Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance[END_REF]. IBD was characterized using GenePop v. 4.7.0 [START_REF] Rousset | Genepop'007: a complete reimplementation of the Genepop software for Windows and Linux[END_REF]) utilizing the regression of the pairwise genetic differentiation estimator a r calculated between individuals and the logarithm of geographic distances [START_REF] Rousset | Genetic differentiation between[END_REF]. The confidence intervals of the regression slopes and Mantel tests with 10,000 permutations assessed the significance of the IBD pattern. This analysis was first conducted on the WP samples of the whole dataset, and then separately for SP and WP using the Leiria dataset.

Hybrid detection

We took advantage of the procedure used in [START_REF] Burban | Evidence for low-level hybridization between two allochronic populations of the pine processionary moth, Thaumetopoea pityocampa (Lepidoptera: Notodontidae)[END_REF] to detect potential hybrids along the two transects. Briefly, from Structure analysis of the transect datasets, we retrieved for each locus the allele frequencies corresponding to the SP and WP clusters to generate parental simulated datasets: 1000 parental simulated SP individuals (pSPsim) and 1000 simulated WP individuals (pWPsim). We then used these simulated parents to generate 1000 individuals for each hybrid category (first and second generations of hybrids, namely F1sim and F2sim, as well as backcross SP and WP, namely bcSPsim and bcWPsim) with Hybridlab [START_REF] Nielsen | HYBRIDLAB (version 1.0): a program for generating simulated hybrids from population samples[END_REF]). The complete simulated dataset was analyzed with Introgress [START_REF] Gompert | Introgress: a software package for mapping components of isolation in hybrids[END_REF] to determine the range of hybrid index of each category, by estimating the individual proportion of ancestry corresponding to each parental pool. Then, NewHybrid was used to estimate the posterior probabilities of each simulated individual to be assigned mostly in one of the mentioned categories [START_REF] Anderson | A model-based method for identifying species hybrids using multilocus genetic data[END_REF]. This procedure further allows calculating the performance to detect each category [START_REF] Vähä | Efficiency of model-based Bayesian methods for detecting hybrid individuals under different hybridization scenarios and with different numbers of loci[END_REF].

3

The field dataset was then analyzed jointly with the simulated dataset with NewHybrid, to avoid biases linked to unbalanced category sizes. To test the consistency of the results, ten replicates of each analysis were performed during 100,000 sweeps once the burn-in step was stabilized, using Jeffreys prior for allele frequencies and mixing proportions. The field dataset was also analyzed with Introgress, using parental simulated data as parental reference. The field individuals for which the h-index fell outside the range of parental simulated data were then considered as hybrid, while the comparison with the type and rate of mis-assignments obtained for the simulated dataset was taken into account for NewHybrid, as in [START_REF] Burban | Evidence for low-level hybridization between two allochronic populations of the pine processionary moth, Thaumetopoea pityocampa (Lepidoptera: Notodontidae)[END_REF].

Co-occurrence patterns of SP and WP

To explore the fine scale patterns of co-occurrence between the two populations, and investigate the impact of the SP density on WP, we examined the relative number of SP and WP trapped in each site. The number of SP and WP individuals caught per site and per year were plotted against each other, and their relationship was tested through a Pearson correlation after logarithmic transformation of the data, considering only site/year for which both populations were present. Note that some LateSP individuals identified after analysis of genetic data (see Results) were grouped with the phenological SP individuals for this analysis. As we did not genotype all trapped individuals, a slight proportion of the phenological WP individuals that were not genotyped could actually correspond to unidentified LateSP. Given the low number of LateSP identified after genetic analyses, the corresponding bias remains very limited.

Results

Population genetic analyses

Considering the whole dataset, only few cases of departure from Hardy to Weinberg equilibrium were observed in more than one site (namely Ppit30 and Ppit38 in Leiria WP and SP for both years and Ppit46 in Varges and Viseu). The percentage of null alleles was over 5% only for Ppit33 and Ppit30 in 6 and 5 sites, respectively. No pair of loci was consistently in linkage disequilibrium for a majority of sampling sites. Observed (and expected) 1).

The three first axes of the PCA performed for the whole dataset explained 5.6%, 3.3% and 2% of the total variance, respectively (Fig. 2). Axis 1 mainly separated the SP from all other populations. Eight WP individuals clearly clustered within the SP group. Hereafter, they will be referred to as LateSP individuals, sensu [START_REF] Burban | Evidence for low-level hybridization between two allochronic populations of the pine processionary moth, Thaumetopoea pityocampa (Lepidoptera: Notodontidae)[END_REF]. All winter populations were distributed along axis 2 following an approximate east-west order. Axis 3 corresponded roughly to a north-south gradient.

The results obtained from Structure were consistent with those obtained from the PCA. The LogP (X|K) curve increased until K = 6 and then reached a plateau, but the variance increased for K ≥ 4, and ΔK showed a clear peak for K = 2 (Fig. S1). At K = 2, one cluster corresponded to all WP populations while the other corresponded to all SP individuals including the eight LateSP also identified on the PCA (Fig. 3, Table S3). For K = 3, Leiria WP and Caparica were grouped in the same cluster, while Tavira and Varges individuals were grouped in the 3rd cluster with high individual Q-values, following a geographic East-West pattern. Individuals from the other sites were assigned to both of these WP clusters (Fig. 3). For K = 4, 5 and 6, major runs (respectively 9/10, 6/10, 7/10) defined the additional clusters as corresponding to one population, respectively Tavira, Grândola and Varges. When performed separately within SP and within Leiria WP, both the PCA and Structure analyses did not reveal any particular pattern (not shown).

Note that the eight LateSP individuals found in the transects were removed from the WP samples in all population-based genetic analyses. As in previous studies, pairwise Fst indicated that Leiria SP was highly differentiated (0.21 < Fst < 0.28) from all WPs (Table 2). No inter-annual differentiation was found, neither in the SP nor in the Leiria WP. Pairwise differentiation between Portuguese WPs varied from 0.01 between the proximate sites of Fundão and Viseu to 0.21 between the more distant sites of Caparica and Varges. In the two transects, the differentiation between sites within the SP or within the WP was consistently low, reaching 0.03 at most (Table S4).

A clear and significant pattern of IBD was found for WP individuals all over Portugal (Fig. S2a) but with two different patterns depending on the geographical scale: at a shorter scale, e.g. individuals separated by less than ca. 100 km, the IBD slope was clearly weaker (0.031) than at a larger scale (0.15). Analyses on the transects, corresponding to a small geographical scale, also showed a significant IBD pattern for WP individuals, although with a slightly weaker slope (0.014) (Fig. S2b). Contrarily, no IBD pattern was evidenced when considering SP individuals (Fig. S2c). Mantel tests were all highly significant except for SP individuals.

Hybrid detection

Analyses of the simulated dataset with Introgress revealed that h-index ranged from 0 to 0.243 for pSPsim and from 0.642 to 1 for pWPsim. The ranges of h-index obtained for simulated hybrid categories greatly overlapped with each other and therefore could not be used to discriminate between F1, F2 and backcrosses; consequently, individuals could be identified only as parental versus hybrid in a broad sense. Introgress identified 7 field individuals as hybrids (Fig. S3, Table S5). These hybrid individuals corresponded to 1 WP trapped in São Martinho Porto in 2014, showing a h-index value of 0.640 and 6 SP (trapped in Nazaré and Monte de Bois in 2012, São Martinho Porto in 2012 and 2014, and São Gregório in 2014) with h-index values ranging from 0.244 to 0.303.

Table S6 summarizes the performances in hybrid detection with NewHybrid. Both efficiency and accuracy were high for parental categories while they were lower for hybrids. A small proportion of parental individuals could however be mis-assigned in their respective backcross categories, and backcross individuals to their respective parental category (Table S7). More, as few simulated parents were mis-assigned as backcrosses and exhibited a posterior probability of 1 to be in hybrid categories, no threshold could be defined to ascertain any backcross identity. Therefore, it would have been possible to ascertain the hybrid origin only for individual assigned as F1, but no such field individual was identified.

However, while the percentage of field individuals assigned as bcWP with NewHybrid was identical to the percentage of pWPsim mis-assigned as bcWP (12%), the proportion of field individuals assigned as bcSP (13.2%) was much higher than the proportion of pSPsim mis-assigned in this category (5.8%) (Table S7). This suggests the occurrence of bcSP or individuals with a more complex hybrid ancestry among the SP, notably in the vicinity of São Martinho Porto (Fig. 4), which is consistent with Introgress results.

Co-occurrence patterns of SP and WP

As expected, no male was caught during the SP flight period (April to mid July) in any of the 6 sites outside the Leiria region. In contrast, individuals were collected during the expected flight periods of the SP and the WP in both transects, and a temporal gap without any capture was found in each locality where both populations were collected the same year (Table S2).

The numbers of males from each population (identified solely from their date of emergence) varied strongly between sites and dates (Fig S4). In particular, in 2012 and 2014, the number of trapped SP males was higher than that of WP males in the coastal transect ( 4501 3 SP vs. 21 WP in 2012; 157 SP vs. 59 WP in 2014 over all sites); WP males were mostly trapped in the northernmost and southernmost sites, at the limits of the SP range, while almost none occurred in any other coastal site. The inner transect revealed an opposite pattern (43 SP vs. 62 WP in 2012 and55 SP vs. 126 WP in 2014). More SP than WP males were however trapped in its central part.

Among individuals genetically analyzed, a large majority (98%) of SP males emerged before mid-July, while 8 individuals collected after this date corresponded to LateSP sensu Burban et al. 2016 (i.e., phenotypically WP individuals collected after July 15th but assigned to the SP genetic cluster from PCA and Structure analyses) (Table S2). In the inner transect, the end of the SP flight period was even earlier, no male being caught after the end of June.

Regarding the eight LateSP individuals, we observed that they were all found in the coastal transect, mainly in localities where WP was rare (Nazaré) or absent (Pedrógão, Praia Vieira, T171) (Fig. 4). Three of them trapped in the end of July could be considered as representing the extreme end of the usual period of SP flight. Among them, one individual was surprisingly trapped in Cabeça Gorda, a southern locality outside the range of the SP; it could therefore correspond to a long distant migrant. The five remaining LateSP individuals were caught much later, at the end of the summer, corresponding to the end of the WP period of flight. They were thus phenologically shifted compared to the other 3, and probably corresponded to a different situation. We thus suggest restricting the term "LateSP" to the first group of individuals, whose emergence dates overlap with the beginning of the WP flight period, while the second group will now be referred to as "timelagged SP". Taking into account the site/year for which both populations were present, logtransformed data showed a significant negative correlation (Pearson r = -0.6899; p value = 0.003) reflecting a negative exponential relationship between SP and WP densities (Fig. 5).

Discussion

Phylogeographic context and genetic structure

Previous works revealed that PPM exhibits a high mitochondrial structuration across Portugal, in relation to the low mobility of the females [START_REF] Burban | Evidence for low-level hybridization between two allochronic populations of the pine processionary moth, Thaumetopoea pityocampa (Lepidoptera: Notodontidae)[END_REF]). The present study, conducted at a larger spatial scale as compared to [START_REF] Santos | Genetic isolation through time: allochronic differentiation of a phenologically atypical population of the pine processionary moth[END_REF]Santos et al. ( , 2011a)), evidenced also a nuclear genetic structure for the WPs sampled over Portugal. Both PCA and Structure analyses highlighted successively an east-west then a north-south pattern, while a general pattern of IBD was found for the first time over Portugal. IBD was previously evidenced in France over similar geographic distances when physical geographical barriers (namely the Massif Central) were taken into account [START_REF] Jombart | Comparative population genetic study of two oligophagous insects associated with the same hosts[END_REF]. Groups of sampled sites with very low pairwise differentiation nonetheless showed that gene flow is favoured in some regions (e.g., between Varges, Viseu and Fundão, the north-eastern sites; or between Caparica and Leiria WP, along the coast in the central part of Portugal). On the other hand, higher Fst values were obtained between the Leiria WP-Caparica group and the two most distant sites of Varges in the North and Tavira in the South of Portugal. It suggests that, additionally to the effect of distances, barriers to dispersal impedes PPM migration between some locations. Such a pattern of strong versus weak genetic structure in particular regions could be due to some landscape, land use and land cover features creating ecological discontinuities and limiting gene flow, while others favour connectivity of moth populations. Moreover, open habitats and agricultural lands, where isolated pine trees were proved to help stepping stone moth movements in France [START_REF] Rossi | Trees outside forest in agricultural landscapes: spatial distribution and impact on habitat connectivity for forest organisms[END_REF], can also bridge otherwise distant PPM populations.

In accordance with PCA and Structure results, and consistently with previous works based on lower numbers of sampling sites and genetic markers [START_REF] Santos | Genetic isolation through time: allochronic differentiation of a phenologically atypical population of the pine processionary moth[END_REF](Santos et al. , 2011a;;[START_REF] Burban | Evidence for low-level hybridization between two allochronic populations of the pine processionary moth, Thaumetopoea pityocampa (Lepidoptera: Notodontidae)[END_REF], the differentiation between the SP and any of the sampled WPs were very high. They were even actually of the same magnitude than between the most distant WP populations sampled in Portugal. However, the differentiation between the two phenological populations is due to their recent bottlenecks and reproductive isolation while the differentiation between WPs is mostly due to geographic distances and landscape features.

At the spatial scale of the SP distribution, the SP and the WP showed striking differences in their genetic characteristics. First, the SP proved to be genetically less diverse than the WP at the same scale, which is likely due to the bottlenecks evidenced by [START_REF] Leblois | Deciphering the demographic history of allochronic differentiation in the pine processionary moth Thaumetopoea pityocampa[END_REF] both at the origin of the SP and in the last decades. Contrarily, WP did not suffer such drastic bottleneck events and certainly beneficiated from gene flow from surrounding populations. Second, a significant IBD was evidenced for the WP but not for the SP across the transects. SP was highly homogeneous over its whole range, without significant differentiation neither between years nor between sites, even though it expanded out of the Leiria forest into more fragmented urbanized areas. The recent spatial expansion can actually have prevented the IBD pattern to emerge (but see [START_REF] Leblois | Influence of spatial and temporal heterogeneities of the estimation of demographic parameters in a continuous population using individual microsatellite data[END_REF], showing that an IBD pattern can rapidly emerge after a spatial expansion). Overall, these results suggest that the geographical expansion of the SP observed from 1997 to the present day occurred via gradual diffusive dispersion, without stepping-stone migration, and did not leave any particular signature in its genetic structure. The absence of IBD pattern for the SP may also be due to its higher population density and/or to an increase in dispersal, which necessitates further investigations.

Both allochrony and some spatial segregation limit gene flow between the SP and the WP

Extended SP flight period albeit associated to rarity of hybrid individuals

The systematic sampling we used here included sites with more diversified environmental conditions compared to previous studies, which could modulate phenology [START_REF] Démolin | Comportement des adultes de Thaumetopoea pityocampa Schiff. Dispersion spatiale, importance écologique[END_REF][START_REF] Robinet | Looking beyond the large scale effects of global change: local phenologies can result in critical heterogeneity in the pine processionary moth[END_REF]. It allowed a fine temporal monitoring over the whole SP range to trap and analyse the rare individuals flying at intermediate dates, i.e. in the second half of July. Our study showed that males flying in mid or late July were pure SP individuals emerging later than expected, rather than hybrids. Interestingly, monitoring and genetic assignments confirmed the presence of a temporal gap in the inner transect only, while the phenology of the SP in the coastal transect was partially overlapping with that of the WP, even though a large majority of SP individuals were indeed trapped before mid-July. If it seems striking to observe that LateSP individuals were only found in the coastal transect, this could be explained by the very high density of the SP locally, if the number of LateSP remains proportional to the number of SP. Additionally, we cannot exclude that climatic conditions along the coast favour a higher variability of SP phenology. This extended duration of the SP flight period, causing some overlap with the early WP male emergence, was expected to increase the probability of hybridization between the two populations but our results do not 1 3 support this hypothesis. We advocate below that fine-scale spatial segregation and densitydependent processes could further limit gene flow between the allochronic populations.

We still identified some few potential backcross individuals, mostly found in São Martinho Porto, the only locality where one F1 and one F2 were previously reported [START_REF] Burban | Evidence for low-level hybridization between two allochronic populations of the pine processionary moth, Thaumetopoea pityocampa (Lepidoptera: Notodontidae)[END_REF]. We here confirmed that the phenology of backcross individuals falls within the emergence of the closest parental species, which probably favours recurrent backcrossing over time.

It is interesting to note that we also identified five "time-lagged SP", i.e. genetically SP individuals emerging long after the end of the SP flight period. As hypothesized already in [START_REF] Burban | Evidence for low-level hybridization between two allochronic populations of the pine processionary moth, Thaumetopoea pityocampa (Lepidoptera: Notodontidae)[END_REF], they could result from an abnormal process of diapause termination as observed for rare individuals in experimental rearing [START_REF] Branco | Experimental evidence for heritable reproductive time in 2 allochronic populations of pine processionary moth[END_REF]. Even though they occur in low numbers, they were regularly found in previous studies (Santos et al. 2011a;[START_REF] Burban | Evidence for low-level hybridization between two allochronic populations of the pine processionary moth, Thaumetopoea pityocampa (Lepidoptera: Notodontidae)[END_REF]). More information concerning these peculiar individuals is needed, especially regarding their fitness and the heritability of their phenology. Yet, they could be implicated in genetic interferences with WP individuals still present or migrating in the core area of the SP [START_REF] Gröning | Reproductive interference between animal species[END_REF].

The WP seems to decrease in the main range of the SP Indeed, the successive studies focusing on the SP after its discovery in 1997 showed some gradual changes in the respective ranges of the SP and the local Leiria WP. In 1997, the SP was found in a very limited area and co-occurred with WP individuals in a single plot, both populations sometimes attacking the very same trees [START_REF] Pimentel | Establishment and expansion of a Thaumetopoea pityocampa (Den. and Schiff.) (Lep. Notodontidae) population with a shifted life cycle in a production pine forest, Central-Coastal Portugal[END_REF][START_REF] Santos | Genetic isolation through time: allochronic differentiation of a phenologically atypical population of the pine processionary moth[END_REF]. Later, the 2007-2010 field surveys, restricted to the coastal area, indicated that the SP had extended and reached high population levels in sites where the WP was henceforth absent (extreme south of the MNL, Nazaré). The two populations were at that time nonetheless found in sympatry at the northern and southern limits of the SP coastal range, respectively in São Martinho Porto and in the north of the MNL [START_REF] Burban | Evidence for low-level hybridization between two allochronic populations of the pine processionary moth, Thaumetopoea pityocampa (Lepidoptera: Notodontidae)[END_REF]. The 2012-2014 dataset presented here evidenced that the lower densities (or even the absence) of SP in the northernmost and southernmost sites of the transects correspond to the current limits of its ongoing latitudinal expansion. Interestingly, the SP densities along the inner transect do not reflect the densities found in the coastal transect at the same latitudes. Number of males trapped along the inner transect is actually in accordance with the SDM models predicting that the SP suitable area is wider in the south than in the north of the transect: the northernmost sites of the inner transect are actually located in areas with low predicted climate suitability for the SP (Fig. 1b; [START_REF] Godefroid | Climate constrains range expansion of an allochronic population of the pine processionary moth[END_REF]). Yet, WP is now absent or scarce in sites where it used to be found in sympatry with the SP few years ago in the coastal transect. Noticeably in the inner transect, the highest density of SP recorded in Monte de Bois also coincides with a lower density of the WP.

Even if the WP tends to strongly decrease along the coast, it is not strictly absent from all coastal sites where the SP occurs. In particular, even if no WP males were recorded in Nazaré from 2007 to 2010 and in 2014, one was trapped in 2012 in the same site. Similarly, WP males were found in São Martinho Porto in 2014, even though none was trapped in 2012. It could be hypothesized that WP could re-colonize areas from which it disappeared in the previous years. On the other hand, such records could correspond to rare long-distance dispersal events from surrounding stands as PPM males have good flying capacities and can be efficiently attracted by the pheromone lures [START_REF] Salvato | Do sexual pheromone traps provide biased information of the local gene pool in the pine processionary moth?[END_REF]. If adults monitoring were useful to precisely link phenology and genotype, recording SP 1 3 and WP nests (which would demonstrate the occurrence of effective reproduction) along the two transects will be necessary to infer in detail the fine-scale co-occurrence patterns of the two populations. Overall, the spatial dynamics of each population in relation with the SP expansion thus tends to shift from sympatry to parapatry. If gene flow was (and is still) mainly limited by allochrony, it now seems to be also by the current range and local densities of the WP and SP. Actually, the absence of hybrids in the regions where LateSP individuals occur and the SP reproductive period could overlap with that of the WP is certainly due to the absence of WP locally.

Competition could be the driver of the observed parapatry

The restricted climatic niche of the SP could be a key factor explaining its current distribution limited to a coastal fringe in Central Portugal [START_REF] Godefroid | Climate constrains range expansion of an allochronic population of the pine processionary moth[END_REF]. On the other hand, the WPs are expected to occur all over Portugal, as the climate is favourable to the PPM in the whole country [START_REF] Roques | Climate warming and past and present distribution of the processionary moths (Thaumetopoea spp.) in Europe, Asia Minor and North Africa[END_REF]. The apparent decrease of WP in the SP range could therefore result from competitive exclusion that might be especially intense during SP outbreaks as observed in the late nineties [START_REF] Pimentel | Establishment and expansion of a Thaumetopoea pityocampa (Den. and Schiff.) (Lep. Notodontidae) population with a shifted life cycle in a production pine forest, Central-Coastal Portugal[END_REF][START_REF] Santos | Genetic isolation through time: allochronic differentiation of a phenologically atypical population of the pine processionary moth[END_REF]. When two species have overlapping climatic niches, competition along an environmental gradient can lead to parapatry at the physiological limit of the dominant species [START_REF] Miller | Pattern and process in competition[END_REF][START_REF] Bull | Ecology of parapatric distributions[END_REF]. The realized niche could therefore be limited by interspecific competition, which is recognized to be a major interaction between phytophagous insects [START_REF] Denno | Interspecific interactions in phytophagous insects: competition reexamined and resurrected[END_REF], and more generally could explain the limited co-occurrence of closely related species [START_REF] Weinstein | The role of environment, dispersal and competition in explaining reduced co-occurrence among related species[END_REF].

Competition between SP and WP individuals could be direct, through interaction of larvae on the same host, or indirect through modification of host quality due to previous attacks. Concerning direct interaction, WP females are known to select hosts through visual [START_REF] Démolin | Comportement des adultes de Thaumetopoea pityocampa Schiff. Dispersion spatiale, importance écologique[END_REF]) and olfactory [START_REF] Paiva | Pine volatiles mediate host selection for oviposition by Thaumetopoea pityocampa (Lep., Notodontidae)[END_REF]) cues, often resulting in an aggregated distribution of PPM nests on some trees, preferentially along stand edges [START_REF] Samalens | Does landscape composition alter the spatiotemporal distribution of the pine processionary moth in a pine plantation forest?[END_REF]. Co-occurrence of several nests on a single tree can even be advantageous for WP individuals in regions with cold winter temperatures, because larvae can merge in a single tent, which enhances their resistance against adverse cold conditions [START_REF] Battisti | Natural history of the processionary moths (Thaumetopoea spp.): new insights in relation to climate change[END_REF]. However, this might not be the case in Portugal, where mild winters offer good climatic conditions for larval development. On the other hand, high nest density at the tree level tends to increase overall mortality of young larvae [START_REF] Pimentel | Spatial patterns at host and forest stand scale and population regulation of the pine processionary moth Thaumetopoea pityocampa[END_REF]). More, SP larvae are still in the very beginning of their development when WP females lay their eggs, but most SP larvae have reached L4 and L5 instars when WP individuals hatch. As both populations consume preferentially one-year-old needles, the competition would clearly be in favour of the SP if SP and WP nests co-exist on the same host.

Previous defoliations affect in various ways the fitness of phytophagous insects [START_REF] Nykänen | Damage-induced changes in woody plants and their effects on insect herbivore performance: a meta-analysis[END_REF]. Considering PPM and black pines, [START_REF] Battisti | Host-plant relationships and population dynamics of the pine processionary caterpillar Thaumetopoea pityocampa (Denis and Schiffermüller)[END_REF] and [START_REF] Tamburini | Effects of climate and density-dependent factors on population dynamics of the pine processionary moth in the Southern Alps[END_REF] suggested that a reduction of host suitability after defoliation was responsible for a negative density dependant feedback through lower larval survival and fecundity. Contrarily, [START_REF] Hodar | No evidence of induced defense after defoliation in three pine species against an expanding pest, the pine processionary moth[END_REF] and [START_REF] Lombardero | Effects of defoliation and site quality on growth and defenses of Pinus pinaster and P. radiate[END_REF] did not detect any impact of previous defoliation on larval performance through induced defences in various pine species. However, these last results were obtained from experimental defoliations occurring during the WP larval development, at the end of summer after the needle growth was achieved [START_REF] Bosc | Etude expérimentale du fonctionnement hydrique et carboné des organes aériens du Pin maritime (Pinus pinaster Ait[END_REF]. Conversely, induced accumulation of phenolics was recorded on Pinus sylvestris submitted to earlier artificial defoliation, which caused lower sawfly 1 3 performances the following year [START_REF] Roitto | Induced accumulation of phenolics and sawfly performance in Scots pine in response to previous defoliation[END_REF]. Therefore, a systemic negative impact of previous SP defoliation on pine suitability for WP individuals could occur and remains to be specifically investigated experimentally.

A new and original scenario of evolution from sympatry to parapatry during speciation

Parapatry represents a compelling species distribution, questioning the location and maintenance of the suture zone 1991). When sister taxa are concerned, different scenarios account for such a distribution [START_REF] Mayr | Systematics and the origin of species[END_REF][START_REF] Mayr | Animal Species and Evolution[END_REF][START_REF] Endler | Geographic variation, speciation, and clines[END_REF][START_REF] Hewitt | The subdivision of species by hybrid zones[END_REF][START_REF] Harrison | Hybrid zones: windows on evolutionary process[END_REF]). The more largely documented situation involves a first step of allopatric distribution, allowing divergence due to genetic drift and/or differential selection pressure, followed by a secondary contact [START_REF] Case | The community context of species' borders: ecological and evolutionary perspectives[END_REF][START_REF] Hewitt | Quaternary phylogeography: the roots of hybrid zones[END_REF][START_REF] Pettengill | Phylogeography of speciation: allopatric divergence and secondary contact between outcrossing and selfing Clarkia[END_REF]. It can then promote reinforcement of reproductive isolation through different types of intrinsic mechanisms [START_REF] Hoskin | Speciation via species interactions: the divergence of mating traits within species[END_REF][START_REF] Matute | Intrinsic reproductive isolation between two sister species of Drosophila[END_REF][START_REF] Canestrelli | The tangled evolutionary legacies of range expansion and hybridization[END_REF].

In a second scenario referring to parapatric speciation, the question of the origin and maintenance of a contact zone are more intricated. It supposes primary intergradation between neighbouring populations, linked to a patchy environment or to ecological clines [START_REF] Gavrilets | Patterns of parapatric speciation[END_REF][START_REF] Durrett | Spatial models for hybrid zones[END_REF]. Ultimately, parapatric distribution can arise by segregation of sympatric ecotypes, as hypothesized for salamander [START_REF] Fisher-Reid | Parapatric divergence of sympatric morphs in a salamander: incipient speciation on Long Island?[END_REF], and along fine spatial scale as for rockfish [START_REF] Ingram | Speciation along a depth gradient in a marine adaptive radiation[END_REF] or cyclids [START_REF] Seehausen | Speciation through sensory drive in cichlid fish[END_REF].

Convincing cases of parapatric distribution following sympatric differentiation concern polyploid speciation in plants [START_REF] Stuessy | Diploid and polyploid cytotype distribution in Melampodium cinereum and M. leucanthum (Asteraceae, Heliantheae)[END_REF][START_REF] Baack | Cytotype segregation on regional and microgeographic scales in snow buttercups (Ranunculus adoneus: Ranunculaceae)[END_REF]) for which differential ecological niche and reproductive interference through hybridization are accountable for the spatial segregation of different cytotypes [START_REF] Burton | Fitness differences among diploids, tetraploids and their triploid progeny in Chamerion angiistifolium (Onagraceae): mechanisms of inviability and implication for polyploid evolution[END_REF]. Differential ecological tolerance of polyploids, allowing colonization of new habitats is a major driver shaping their geographic distribution [START_REF] Weiss-Schneeweiss | Evolutionary consequences, constraints and potential of polyploidy in Plants[END_REF].

The SP/WP history could illustrate another scenario explaining a shift from sympatry to parapatry. Allochrony, acting as an automatic magic trait [START_REF] Servedio | Magic traits in speciation: 'magic' but not rare?[END_REF], would have immediately disrupted gene flow between the incipient populations that experienced a strong concomitant genetic drift [START_REF] Leblois | Deciphering the demographic history of allochronic differentiation in the pine processionary moth Thaumetopoea pityocampa[END_REF]. While these authors estimated a continuous presence of the ancestral WP population in the Iberian Peninsula during the last 900,000 years, the initial divergence likely arose in sympatry, even if the precise range of both SP and WP since this first step of differentiation remains unknown. It could have changed over time as [START_REF] Leblois | Deciphering the demographic history of allochronic differentiation in the pine processionary moth Thaumetopoea pityocampa[END_REF] highlighted a complex demographic history. Therefore, in the recent past, rapid range expansion of the SP and possible competitive exclusion of local WP might have led to the parapatric pattern observed nowadays.

Based on the knowledge gained so far, including the new information brought by the current study, we could represent the differentiation scenario of the SP and the local WP on a Dieckman's cube (Dieckman et al. 2004) (Fig. 6.) We hypothesized that (1) a strong shift in phenology was the primary cause of mating disruption [START_REF] Santos | Genetic isolation through time: allochronic differentiation of a phenologically atypical population of the pine processionary moth[END_REF](Santos et al. , 2011a)). As a result, the differentiating SP experienced radically different environmental conditions (2), leading to further adaptive differentiation between populations (e.g., proven increased tolerance of the larvae to high summer temperatures (Santos et al. 2011b); changes in reproductive traits and reproductive success [START_REF] Santos | Phenotypic divergence in reproductive traits of a moth population experiencing a phenological shift[END_REF]). The current step of this on-going differentiation scenario suggests that (3) the system is currently tending to move rapidly from sympatry to parapatry, which further impedes gene flow and limits 1 3 introgression between populations in spite of the occurrence of LateSP individuals. Yet, speciation is still not an achieved process. Future climate changes could notably further modulate the ecological and spatial differentiation, as well as phenology, the main trait implicated so far in reproductive isolation.

Fig. 6 Dieckman's cube illustrating the hypothesized scenario of SP differentiation during the speciation process regarding ecology, mating and space. 1-2: From the ancestral WP, a primary mutational step in the Leiria region was responsible for allochronic differentiation of the SP, the abrupt phenological shift acting as a prezygotic barrier and initiating the process on the "mating" axis. 2-3: The shift of phenology submitted the SP to new environmental conditions during its development, causing divergent selection, i.e. differentiation along the "ecology" axis. 3-4: Environmental characteristics delimit the expansion area of the SP, while competition excludes the WP from the core area of the SP, leading to an evolution from initial sympatry to parapatry on the "space" axis, reducing further the probability of introgression between the diverging populations

Fig. 1

 1 Fig.1Maps showing the sampling sites of T. pityocampa a over Portugal with the Leiria region highlighted in white b on the two transects studied in the Leiria region: black and blue dots correspond to the coastal and inner transects, respectively. The background shows the habitat suitability estimated by species distribution modelling[START_REF] Godefroid | Climate constrains range expansion of an allochronic population of the pine processionary moth[END_REF]. The values depict the output of a consensus model and represent the proportion of different models predicting the presence of the species

Fig. 2 3 Fig. 3

 233 Fig. 2 Graph of the principal component analysis (PCA) ran on microsatellite genotypes of all sampled individuals. a Axis 1-2, b axis 1-3

Fig. 4

 4 Fig. 4 Genetic assignments of 345 males trapped in the two transects in the Leiria region in 2012 and 2014. Only hybrids detected from both Introgress and NewHybrid analyses are highlighted. Colors refer to NewHybrid majority assignment and phenology (red: pSP, pink: Late-pSP, bleu: Time-lagged pSP, green: pWP, yellow: bcSP, brown bcWP)

Fig. 5

 5 Fig. 5 Relative numbers of SP and WP individuals trapped per site and year in the Leiria region, taking into account the genetic assignment for LateSP and time-lagged SP individuals. Grey dots correspond to sites for which one population was absent. Pearson correlation was estimated after Ln transformation of the data, considering only those for which both populations were present, i.e. excluding grey dots

Table 1

 1 heterozygosities ranged from 0.43 (0.45) in the SP both in 2012 and 2014 to 0.64 (0.68 and 0.70) in Viseu and Fundão (Table

	Number of genotyped individuals (N), number of	Populations	N	Na	Ho	He
	alleles (Na), observed (Ho) and expected (He) heterozygosities, after removing the LateSP	Varges Viseu	26 24	150 140	0.63 0.64	0.67 0.68
	individuals from the dataset	Fundão	30	155	0.64	0.70
		Caparica	29	91	0.47	0.51
		Grândola	30	137	0.61	0.63
		Tavira	30	139	0.57	0.64
		Leiria WP12	45	128	0.60	0.58
		Leiria SP12	100	103	0.43	0.45
		Leiria WP14	80	142	0.54	0.58
		Leiria SP14	112	106	0.43	0.45

Table 2

 2 Pairwise Fst between populations, after removing LateSP individuals from the dataset

		Var	Vis	Fun	Cap	Grâ	Tav	WP12	SP 12	WP 14
	Varges	-
	Viseu	0.04	-
	Fundão	0.05	0.01	-
	Caparica	0.21	0.10	0.10	-
	Grândola	0.11	0.06	0.07	0.08	-
	Tavira	0.10	0.09	0.10	0.18	0.10	-
	Leiria WP12	0.16	0.07	0.06	0.04	0.08	0.16	-
	Leiria SP12	0.26	0.24	0.25	0.27	0.22	0.28	0.24	-
	Leiria WP14	0.18	0.07	0.07	0.03	0.08	0.16	0.01	0.21	-
	Leiria SP14	0.26	0.24	0.24	0.26	0.22	0.27	0.23	0.01	0.21
	Significant values are presented in bold

Institut de Biologie Computationnelle, University of Montpellier, 34095 Montpellier, France

Acknowledgements

We thank Liliana Vilas Boas for her help with fieldwork and Helena Santos for cal advice during DNA extractions. We are grateful to Flora Abella who participated in DNA extractions. We acknowledge J. Rousselet (URZF, INRA Orléans, France) for providing the larvae from Varges. We are grateful to two anonymous reviewers for their helpful suggestions. Data used in this work were partly produced through the GenSeq molecular genetic analysis technical facilities of the Labex CeMEB, and ANR "Investissements d'Avenir" program (ANR-10-LABX-04-01). This study was partly financed by Fundação para a Ciência e Tecnologia, FCT-MCES, Portugal, (Project Pest-OE/AGR/UI0239/2011) and the Agence Nationale de la Recherche, France (Project GENO-SPACE ANR-16-CE02-0008). It was also supported by Susana Rocha Ph.D. scholarship with the reference SFRH/BD/90188/2012(FCT-MCES).

Compliance with ethical standards

Conflict of interest

The authors declare no conflict of interest.

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.