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

Christian Burban, Susana Rocha, Raphaël Leblois, Jean-Pierre Rossi, Laure

Saune, Manuela Branco, Carole Kerdelhue

► To cite this version:

Christian Burban, Susana Rocha, Raphaël Leblois, Jean-Pierre Rossi, Laure Saune, et al. From sympatry to parapatry: a rapid change in the spatial context of incipient allochronic speciation. Evolutionary Ecology, 2020, 34 (1), pp.101-121. 10.1007/s10682-019-10021-4. hal-02624010

HAL Id: hal-02624010

https://hal.inrae.fr/hal-02624010

Submitted on 2 Feb 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# From sympatry to parapatry: a rapid change in the spatial context of incipient allochronic speciation

Christian Burban¹ ID · Susana Rocha² ID · Raphaël Leblois³,⁴ ID · Jean-Pierre Rossi³ ID · Laure Sauné³ · Manuela Branco² ID · Carole Kerdelhué³ ID

Received: 17 April 2019 / Accepted: 21 November 2019 / Published online: 27 November 2019
© Springer Nature Switzerland AG 2019

## Abstract

Speciation is nowadays recognized as a dynamic process in which the respective roles of forces driving ecological differentiation and reproductive isolation can change through time and space. Incipient speciation events are particularly useful to follow such processes that are not tractable when considering well-differentiated taxa. A promising case study was discovered in the pine processionary moth, *Thaumetopoea pityocampa*, a Mediterranean defoliator of *Pinus* species, for which allochrony acting as an automatic magic trait was recognized as the major driver of an incipient speciation process. In Portugal, a unique population with a shifted phenology, known as the summer population (SP), co-occurs with a population following the typical life cycle, known as the winter population (WP). We monitored male activity of both populations in the Leiria region, i.e. over the whole SP distribution range using a systematic sampling along two transects, and studied Portuguese WPs at a larger geographical scale to explore their genetic diversity and spatial pattern of differentiation. Results showed that the WPs were genetically more diverse than the SP, with a strong pattern of isolation by distance both at large and small spatial scales, while the SP was very homogeneous over its whole range, without signature of its recent spatial expansion. Contrary to our expectations, no F1 hybrids were identified, even though we found an extended flight period of the SP, overlapping with the beginning of the WP reproductive period. Interestingly, the SP was found to be mostly limited to the sea shore where the WP is now scarce or absent, which could suggest competitive exclusion. Once clearly occurring in a sympatric context, the allochronic differentiation tends to develop nowadays in parapatry.

**Keywords** Allochronic speciation · Phenology · Sympatry · Parapatry · Hybridization · *Thaumetopoea pityocampa*

Christian Burban and Susana Rocha should be considered joint first author.

**Electronic supplementary material** The online version of this article (https://doi.org/10.1007/s10682-019-10021-4) contains supplementary material, which is available to authorized users.

<figure><img src="image_2.png" /></figure>

Christian Burban
christian.burban@inra.fr

Extended author information available on the last page of the article
Introduction

Automatic magic traits refer to traits subject to divergent selection, which at the same time
contribute to non-random mating (Servedio et al. 2011; Nosil 2012). Classical examples
of such traits include phenology when assortative mating occurs via temporal isolation or
habitat choice when mating occurs in the preferred habitat (Nosil 2012). 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 Gavrilets 2013). 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. Butlin et al. 2008). Yet, studying how
ecological and spatial patterns driving the differentiation interact all along the speciation
process remains highly relevant (Dieckman et al. 2004; Rundle and Nosil 2005; Kisel and
Barraclough 2010; Aguilée et al. 2011). Some key information are however hardly acces-
sible when considering past speciation events because only the final outcomes, poten-
tially modulated by biogeographic changes, are available (Gavrilets et al. 2000; Bolnik
and Fitzpatrick 2007). 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 repro-
ductive time and is recognized as a salient magic trait (Taylor and Friesen 2017). In case of
seasonal shift, the diverging populations can experience a significant shift in their environ-
mental niche, as their life stages face different abiotic constraints (e.g., larvae experienc-
ing hot vs. cold conditions, see Santos et al. 2011b). One of the best-documented cases
of true allochronic speciation sensu Taylor and Friesen (2017) is the winter pine proces-
sionary moth (PPM), *Thaumetopoea pityocampa* Denis and Schiffermüller (Lepidoptera:
Notodontidae). The PPM is a well-known defoliator occurring in all Mediterranean coun-
tries (Kerdelhué et al. 2009), 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 (Roques
et al. 2015). 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 Por-
tugal in a pure maritime pine *Pinus pinaster* Ait. forest. This peculiar and localized pop-
ulation 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 (Pimentel et al. 2006). 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 repro-
ductive isolation (Santos et al. 2007, 2011a). Experimental rearing demonstrated that each
type of phenology is highly heritable (Branco et al. 2017). 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
attacked by both populations the same year (Santos et al. 2007). Intensive field surveys in the following years allowed to evidence and follow its range expansion using nest occurrence data and pheromone trapping (Godefroid et al. 2016); 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 (Kerdelhué et al. 2015). 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 (Roques et al. 2015). 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 (Godefroid et al. 2016). Yet, as known for the typical WP over its whole range and notably in Portugal (Roques et al. 2015), SP population size can reach high densities, experiencing outbreak situations as recorded when it was discovered (Pimentel et al. 2006).

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 (Santos et al. 2007, 2011a; Burban et al. 2016). Recent coalescent analyses based on Rad-seq data (Leblois et al. 2018) 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

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Maps 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 (Godefroid et al. 2016). The values depict the output of a consensus model and represent the proportion of different models predicting the presence of the species</figcaption></figure>
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 (Burban et al. 2016):
(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 men-
tioned 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 Burban et al. (2016) 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 Godefroid et al. (2016). 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 Santos et al. (2007, 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 and S2). Sampled individuals
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 Burban et al. (2016). 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, MS-Thpit10, Ppit09, Ppit20, Ppit30, Ppit31, Ppit33, Ppit34, Ppit35, Ppit38, Ppit39, Ppit43, Ppit46, Ppit47, Ppit48. Detailed protocols for amplification and genotyping are described in Rousselet et al. (2004), A'Hara et al. (2012) and Sauné et al. (2015). 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 2012 and 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 (Excoffier et al. 2005). Expected and observed heterozygosities were estimated using the same software.
The proportion of null alleles was estimated for each locus and population using the *FREENA* package (Chapuis and Estoup 2007). All file format conversions were done using *CONVERT* (Glaubitz 2004).

## 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 (Jombart 2008). Individuals were assigned to genetic clusters using the Bayesian clustering method implemented in *STRUCTURE* 2.3.4 under the model with admixture (Pritchard et al. 2000). 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* (Kopelman et al. 2015). 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 (Evanno et al. 2005) implemented in *STRUCTURE* HARVESTER (Earl and von Holdt 2012). 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 (Chapuis and Estoup 2007). 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 (Wright 1943; Rousset 1997). IBD was characterized using *GENEPOP* v. 4.7.0 (Rousset 2008) utilizing the regression of the pairwise genetic differentiation estimator $a_t$ calculated between individuals and the logarithm of geographic distances (Rousset 2000). 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 Burban et al. (2016) 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* (Nielsen et al. 2006). The complete simulated dataset was analyzed with *INTROGRESS* (Gompert and Buerkle 2010) 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 (Anderson and Thompson 2002). This procedure further allows calculating the performance to detect each category (Vähä and Primmer 2006).
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 Burban et al. (2016).

## 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)

<table><caption>**Table 1** Number of genotyped individuals (*N*), number of alleles (*Na*), observed (*Ho*) and expected (*He*) heterozygosities, after removing the LateSP individuals from the dataset</caption><thead><tr><th>Populations</th><th><em>N</em></th><th><em>Na</em></th><th><em>Ho</em></th><th><em>He</em></th></tr></thead><tbody><tr><td>Varges</td><td>26</td><td>150</td><td>0.63</td><td>0.67</td></tr><tr><td>Viseu</td><td>24</td><td>140</td><td>0.64</td><td>0.68</td></tr><tr><td>Fundão</td><td>30</td><td>155</td><td>0.64</td><td>0.70</td></tr><tr><td>Caparica</td><td>29</td><td>91</td><td>0.47</td><td>0.51</td></tr><tr><td>Grândola</td><td>30</td><td>137</td><td>0.61</td><td>0.63</td></tr><tr><td>Tavira</td><td>30</td><td>139</td><td>0.57</td><td>0.64</td></tr><tr><td>Leiria WP12</td><td>45</td><td>128</td><td>0.60</td><td>0.58</td></tr><tr><td>Leiria SP12</td><td>100</td><td>103</td><td>0.43</td><td>0.45</td></tr><tr><td>Leiria WP14</td><td>80</td><td>142</td><td>0.54</td><td>0.58</td></tr><tr><td>Leiria SP14</td><td>112</td><td>106</td><td>0.43</td><td>0.45</td></tr></tbody></table>
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 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 Burban et al. (2016). 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 *F*st indicated that Leiria SP was highly differentiated (0.21<*F*st<0.28) from all WPs (Table 2). No inter-annual differentiation was found, neither in the SP nor in the Leiria WP. Pairwise

<figure><img src="image_4.png" /><figcaption>**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</figcaption></figure>

<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Graphic representation of the individual genetic assignments given by the individual Q-values inferred by STRUCTURE for the 10 populations studied. Individual assignments for K=2–6. Only major runs are shown for K≥4</figcaption></figure>

<table><caption>**Table 2** Pairwise *F*st between populations, after removing LateSP individuals from the dataset</caption><thead><tr><th></th><th>Var</th><th>Vis</th><th>Fun</th><th>Cap</th><th>Grâ</th><th>Tav</th><th>WP12</th><th>SP 12</th><th>WP 14</th></tr></thead><tbody><tr><td>Varges</td><td>–</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Viseu</td><td>0.04</td><td>–</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Fundão</td><td>0.05</td><td>0.01</td><td>–</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Caparica</td><td><strong>0.21</strong></td><td><strong>0.10</strong></td><td><strong>0.10</strong></td><td>–</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Grândola</td><td><strong>0.11</strong></td><td>0.06</td><td>0.07</td><td><strong>0.08</strong></td><td>–</td><td></td><td></td><td></td><td></td></tr><tr><td>Tavira</td><td><strong>0.10</strong></td><td><strong>0.09</strong></td><td><strong>0.10</strong></td><td><strong>0.18</strong></td><td><strong>0.10</strong></td><td>–</td><td></td><td></td><td></td></tr><tr><td>Leiria WP12</td><td><strong>0.16</strong></td><td>0.07</td><td>0.06</td><td>0.04</td><td><strong>0.08</strong></td><td><strong>0.16</strong></td><td>–</td><td></td><td></td></tr><tr><td>Leiria SP12</td><td><strong>0.26</strong></td><td><strong>0.24</strong></td><td><strong>0.25</strong></td><td><strong>0.27</strong></td><td><strong>0.22</strong></td><td><strong>0.28</strong></td><td><strong>0.24</strong></td><td>–</td><td></td></tr><tr><td>Leiria WP14</td><td><strong>0.18</strong></td><td>0.07</td><td>0.07</td><td>0.03</td><td><strong>0.08</strong></td><td><strong>0.16</strong></td><td>0.01</td><td><strong>0.21</strong></td><td>–</td></tr><tr><td>Leiria SP14</td><td><strong>0.26</strong></td><td><strong>0.24</strong></td><td><strong>0.24</strong></td><td><strong>0.26</strong></td><td><strong>0.22</strong></td><td><strong>0.27</strong></td><td><strong>0.23</strong></td><td>0.01</td><td><strong>0.21</strong></td></tr></tbody></table>

Significant values are presented in bold
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 (450
<figure><img src="image_7.png" /><figcaption>**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)</figcaption></figure>

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 and 55 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 “time-lagged SP”.
Taking into account the site/year for which both populations were present, log-transformed 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 structure across Portugal, in relation to the low mobility of the females (Burban et al. 2016). The present study, conducted at a larger spatial scale as compared to Santos et al. (2007, 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 (Kerdelhué et al. 2006). 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

<figure><img src="image_8.png" /><figcaption>**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</figcaption></figure>

<figure><img src="image_9.png" /></figure>
were proved to help stepping stone moth movements in France (Rossi et al. 2016), 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 (Santos et al. 2007, 2011a; Burban et al. 2016), 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 divergence, 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 Leblois et al. (2018) 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 Leblois et al. 2004, 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 (Démolin 1969; Robinet et al. 2015). 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
support this hypothesis. We advocate below that fine-scale spatial segregation and density-
dependent processes could further limit gene flow between the allochronic populations.

We still identified some few potential backcross individuals, mostly found in São Mar-
tinho Porto, the only locality where one F1 and one F2 were previously reported (Burban
et al. 2016). We here confirmed that the phenology of backcross individuals falls within the
emergence curve of the closest parental species, which probably favours recurrent back-
crossing 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
Burban et al. (2016), they could result from an abnormal process of diapause termination
as observed for rare individuals in experimental rearing (Branco et al. 2017). Even though
they occur in low numbers, they were regularly found in previous studies (Santos et al.
2011a; Burban et al. 2016). 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 (Gröning and Hochkirch 2008).

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 (Pimentel et al. 2006; Santos
et al. 2007). 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 hence-
forth 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 (Burban et al. 2016). 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 cur-
rent 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; Godefroid et al. 2016). 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. Simi-
larly, 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 disap-
peared 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 (Salvato et al. 2005).
If adults monitoring were useful to precisely link phenology and genotype, recording SP
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 affected 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 (Godefroid et al. 2016). 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 (Roques et al. 2015). 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 (Pimentel et al. 2006; Santos et al. 2007). When two species have overlapping climatic niches, competition along an environmental gradient can lead to parapatry at the physiological limit of the dominant species (Miller 1967; Bull 1991). The realized niche could therefore be limited by interspecific competition, which is recognized to be a major interaction between phytophagous insects (Denno et al. 1995), and more generally could explain the limited co-occurrence of closely related species (Weinstein et al. 2017).

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 (Démolin 1969) and olfactory (Paiva et al. 2011) cues, often resulting in an aggregated distribution of PPM nests on some trees, preferentially along stand edges (Samalens and Rossi 2011). 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 (Battisti et al. 2015). 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 (Pimentel et al. 2017). 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 (Nykänen and Koricheva 2004). Considering PPM and black pines, Battisti (1988) and Tamburini et al. (2013) suggested that a reduction of host suitability after defoliation was responsible for a negative density dependant feedback through lower larval survival and fecundity. Contrarily, Hodar et al. (2015) and Lombardero et al. (2016) 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 (Bosc 1999). Conversely, induced accumulation of phenolics was recorded on *Pinus sylvestris* submitted to earlier artificial defoliation, which caused lower sawfly
performances the following year (Roitto et al. 2009). 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 main-
tenance of the suture zone (Bull 1991). When sister taxa are concerned, different scenarios
account for such a distribution (Mayr 1942, 1963; Endler 1977; Hewitt 1989; Harrison
1990). 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 (Case et al. 2005; Hewitt 2011; Pettengill and Moeller 2012). It can
then promote reinforcement of reproductive isolation through different types of intrinsic
mechanisms (Moskin and Higgie 2010; Matute and Coyne 2010; Canestrelli et al. 2016).

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 (Gavrilets et al. 2000; Durrett et al. 2000). Ultimately, parapatric distribution can arise by segregation of sympatric ecotypes, as hypothesized for salamander (Fisher-Reid et al. 2013), and along fine spatial scale as for rockfish (Ingram 2011) or cyclids (Seehausen et al. 2008).

Convincing cases of parapatric distribution following sympatric differentiation concern polyploid speciation in plants (Stuessy et al. 2004; Baack 2004) for which differential ecological niche and reproductive interference through hybridization are accountable for the spatial segregation of different cytotypes (Burton and Husband 2000). Differential ecological tolerance of polyploids, allowing colonization of new habitats is a major driver shaping their geographic distribution (Weiss-Schneeweiss et al. 2013).

The SP/WP history could illustrate another scenario explaining a shift from sympatry to parapatry. Allochrony, acting as an automatic magic trait (Servedio et al. 2011), would have immediately disrupted gene flow between the incipient populations that experienced a strong concomitant genetic drift (Leblois et al. 2018). 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 Leblois et al. (2018) 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 (Santos et al. 2007, 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 (Santos et al. 2013)). 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
<figure><img src="image_10.png" /></figure>

**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

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.

**Acknowledgements** We thank Liliana Vilas Boas for her help with fieldwork and Helena Santos for technical 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.

**References**

A'Hara SW, Amouroux P, Argo EE et al (2012) Permanent genetic resources added to molecular ecology resources database 1 August 2011–30 September 2011. Mol Ecol Resour 12:185–189

Aguilée R, Lambert A, Claessen D (2011) Ecological speciation in dynamic landscapes. J Evol Biol 24:2663–2677

Anderson C, Thompson EA (2002) A model-based method for identifying species hybrids using multilocus genetic data. Genetics 160:1217–1229

Baack E (2004) Cytotype segregation on regional and microgeographic scales in snow buttercups (*Ranunculus adoneus*: Ranunculaceae). Am J Bot 91(11):1783–1788
Battisti A (1988) Host-plant relationships and population dynamics of the pine processionary caterpillar *Thaumetopoea pityocampa* (Denis and Schiffermüller). J Appl Entomol 105:393–402

Battisti A, Avci M, Avtzis DN et al (2015) Natural history of the processionary moths (*Thaumetopoea* spp.): new insights in relation to climate change. In: Roques A (ed) Processionary moths and climate change: an update. Springer, Dordrecht, pp 15–80

Bolnik DI, Fitzpatrick BM (2007) Sympatric speciation: models and empirical evidence. Annu Rev Ecol Evol Syst 38:459–487

Bosc A (1999) Etude expérimentale du fonctionnement hydrique et carboné des organes aériens du Pin maritime (*Pinus pinaster* Ait.). Ph.D. Thesis, Université de Bordeaux 2, Bordeaux, France

Branco M, Paiva M-R, Santos H, Burban C, Kerdelhué C (2017) Experimental evidence for heritable reproductive time in 2 allochronic populations of pine processionary moth. Insect Sci 24:325–335

Bull CM (1991) Ecology of parapatric distributions. Annu Rev Ecol Syst 22:19–36

Burban C, Gautier M, Leblois R, Landes J, Santos H, Paiva M-R, Branco M, Kerdelhué C (2016) Evidence for low-level hybridization between two allochronic populations of the pine processionary moth, *Thaumetopoea pityocampa* (Lepidoptera: Notodontidae). Biol J Linn Soc 199:311–328

Burton TL, Husband BC (2000) Fitness differences among diploids, tetraploids and their triploid progeny in *Chamerion angiistifolium* (Onagraceae): mechanisms of inviability and implication for polyploid evolution. Evolution 54(4):1182–1191

Butlin R, Galindo J, Grahame JW (2008) Sympatric, parapatric or allopatric: the most important way to classify speciation? Philos Trans R Soc B Biol Sci 363:2997–3007

Canestrelli D, Porreta D, Lowe WH, Bisconti R, Carere C, Nascetti G (2016) The tangled evolutionary legacies of range expansion and hybridization. Trends Ecol Evol 31(9):677–688

Case TJ, Holt RD, McPeek MA, Keitt TH (2005) The community context of species' borders: ecological and evolutionary perspectives. Oikos 108:28–40

Chapuis MP, Estoup A (2007) Microsatellite null alleles and estimation of population differentiation. Mol Biol Evol 24(3):621–631

Démolin G (1969) Comportement des adultes de *Thaumetopoea pityocampa* Schiff. Dispersion spatiale, importance écologique. Ann Sci For 26(1):81–102

Denno RF, McClure MS, Ott JR (1995) Interspecific interactions in phytophagous insects: competition reexamined and resurrected. Annu Rev Entomol 40:297–331

Dieckman U, Doebeli M, Metz JAJ, Tautz D (eds) (2004) Adaptive speciation. Cambridge University Press, Cambridge

Durrett R, Buttel L, Harrison R (2000) Spatial models for hybrid zones. Heredity 84:9–19

Earl DA, von Holdt BM (2012) Structure harvester: a website and program for visualizing structure output and implementing the Evanno method. Conserv Genet Resour 4:359–361

Endler JA (1977) Geographic variation, speciation, and clines. Princeton University Press, Princeton

Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software structure: a simulation study. Mol Ecol 14:2611–2620

Excoffier L, Laval G, Schneider S (2005) Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evol Bioinf Online 1:47–50

Fisher-Reid MC, Engstrom TN, Kuczinsky CA, Stephens PR, Wiens JJ (2013) Parapatric divergence of sympatric morphs in a salamander: incipient speciation on Long Island? Mol Ecol 22:4681–4694

Gavrilets S, Li H, Vose MD (2000) Patterns of parapatric speciation. Evolution 54(4):1126–1134

Glaubitz JC (2004) CONVERT: A user-friendly program to reformat diploid genotypic data for commonly used population genetic software packages. Mol Ecol Notes 4:309–310

Godefroid M, Rocha S, Santos H, Paiva M-R, Burban C, Kerdelhué C, Branco M, Rasplus J-Y, Rossi J-P (2016) Climate constrains range expansion of an allochronic population of the pine processionary moth. Divers Distrib 22(12):1288–1300

Gompert Z, Buerkle CA (2010) Introgress: a software package for mapping components of isolation in hybrids. Mol Ecol Resour 10:378–384

Gröning J, Hochkirch A (2008) Reproductive interference between animal species. Q Rev Biol 33(3):257–282

Harrison RG (1990) Hybrid zones: windows on evolutionary process. Oxford Surv Evol Biol 7:69–128

Hewitt GM (1989) The subdivision of species by hybrid zones. In: Otte D, Endler JA (eds) Speciation and its consequences. Sunderland, pp 85–110

Hewitt GM (2011) Quaternary phylogeography: the roots of hybrid zones. Genetica 139:617–638

Hodar JA, Torres-Muros L, Zamora R, Antonio J, Perez-Luque AJ, Senhadji K (2015) No evidence of induced defense after defoliation in three pine species against an expanding pest, the pine processionary moth. For Ecol Manag 356:166–172
Hoskin CJ, Higgie M (2010) Speciation via species interactions: the divergence of mating traits within species. Ecol Lett 13(4):409–420

Ingram T (2011) Speciation along a depth gradient in a marine adaptive radiation. Proc R Soc Lond B 278:613–618

Jombart T (2008) Adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics 24:1403–1405

Kerdelhué C, Magnoux E, Lieutier F, Roques A, Rousselet J (2006) Comparative population genetic study of two oligophagous insects associated with the same hosts. Heredity 97(1):38–45

Kerdelhué C, Zane L, Simonato M, Salvato P, Rousselet J, Roques A, Battisti A (2009) Quaternary history and contemporary patterns in a currently expanding species. BMC Evol Biol 9:220

Kerdelhué C, Battisti A, Burban C et al (2015) Genetic diversity and structure at different spatial scales in the processionary moths. In: Roques A (ed) Processionary moths and climate change: an update. Springer, Dordrecht, pp 163–226

Kisel Y, Barraclough TG (2010) Speciation has a spatial scale that depends on levels of gene flow. Am Nat 175(3):316–334

Kopelman NM, Mayzel J, Jakobsson M, Rosenberg NA, Mayrose I (2015) Clumpak: A program for identifying clustering modes and packaging population structure inferences across K. Mol Ecol Resour 15(5):1179–1191

Leblois R, Rousset F, Estoup A (2004) Influence of spatial and temporal heterogeneities of the estimation of demographic parameters in a continuous population using individual microsatellite data. Genetics 166:1081–1092

Leblois R, Gautier M, Rohfritsch A, Foucaud J, Burban C, Galan M, Loiseau A, Sauné L, Branco M, Gharbi K, Vitalis R, Kerdelhué C (2018) Deciphering the demographic history of allochronic differentiation in the pine processionary moth *Thaumetopoea pityocampa*. Mol Ecol 27:264–278

Lombardero MJ, Ayres MP, Bonello P, Cipollini D, Herms DA (2016) Effects of defoliation and site quality on growth and defenses of *Pinus pinaster* and *P. radiate*. Forest Ecol Manag 382:39–50

Matute DR, Coyne JA (2010) Intrinsic reproductive isolation between two sister species of *Drosophila*. Evolution 64(4):903–920

Mayr E (1942) Systematics and the origin of species. Columbia University Press, New York

Mayr E (1963) Animal Species and Evolution. Belknap Press, Cambridge

Miller RS (1967) Pattern and process in competition. Adv Ecol Res 4:1–74

Nielsen EE, Bach L, Kotlicki P (2006) HYBRIDLAB (version 1.0): a program for generating simulated hybrids from population samples. Mol Ecol Notes 6:971–973

Nosil P (2012) Ecological speciation. Oxford series in ecology and evolution. Oxford University Press, Oxford

Nykänen H, Koricheva J (2004) Damage-induced changes in woody plants and their effects on insect herbivore performance: a meta-analysis. Oikos 104:247–268

Paiva M-R, Mateus E, Santos MH, Branco MR (2011) Pine volatiles mediate host selection for oviposition by *Thaumetopoea pityocampa* (Lep., Notodontidae). J Appl Ent 135:195–203

Pettengill JB, Moeller DA (2012) Phylogeography of speciation: allopatric divergence and secondary contact between outcrossing and selfing *Clarkia*. Mol Ecol 21(18):4578–4592

Pimentel C, Calvão T, Santos M, Ferreira C, Neves M, Nilsson JA (2006) 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. For Ecol Manag 233:108–115

Pimentel CS, Ferreira C, Santos M, Calvão T (2017) Spatial patterns at host and forest stand scale and population regulation of the pine processionary moth *Thaumetopoea pityocampa*. Agric For Ento-mol 19:200–209

Pritchard JK, Stephens M, Donnely P (2000) Inference of population structure using multilocus genotype data. Genetics 155:945–959

Robinet C, Laparie M, Rousselet J (2015) Looking beyond the large scale effects of global change: local phenologies can result in critical heterogeneity in the pine processionary moth. Front Physiol 6:334

Roitto M, Rautio P, Markkola A, Julkunen-Tiitto R, Varama M, Saravesi K, Tuomi J (2009) Induced accumulation of phenolics and sawfly performance in Scots pine in response to previous defoliation. Tree Physiol 29:207–216

Roques A, Rousselet J, Avci M et al (2015) Climate warming and past and present distribution of the processionary moths (*Thaumetopoea spp.*) in Europe, Asia Minor and North Africa. In: Roques A (ed) Processionary moths and climate change: an update. Springer, Dordrecht, pp 81–161

Rossi J-P, Garcia J, Roques A, Rousselet J (2016) Trees outside forest in agricultural landscapes: spatial distribution and impact on habitat connectivity for forest organisms. Landsc Ecol 31:243–254
Rousselet J, Magnoux E, Kerdelhué C (2004) Characterization of five microsatellite loci in the pine pro- cessionary moth *Thaumetopoea pityocampa* (Lepidoptera Notodontidae Thaumetopoeinae). Mol Ecol Notes 4:213–214

Rousset F (1997) Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance. Genetics 145(4):1219–1228

Rousset F (2000) Genetic differentiation between individuals. J Evol Biol 13:58–62

Rousset F (2008) *Genepop'007: a complete reimplementation of the Genepop software for Windows and Linux*. Mol Ecol Resour 8:103–106

Rundle HD, Nosil P (2005) Ecological speciation. Ecol Lett 8:336–352

Salvato P, Simonato M, Zan L, Patarnello T, Masutti L, Battisti A (2005) Do sexual pheromone traps pro- vide biased information of the local gene pool in the pine processionary moth? Agric For Entomol 7(2):127–132

Samalens J-C, Rossi J-P (2011) Does landscape composition alter the spatiotemporal distribution of the pine processionary moth in a pine plantation forest? Popul Ecol 53:287–296

Santos H, Rousselet J, Magnoux E, Paiva M-R, Branco M, Kerdelhué C (2007) Genetic isolation through time: allochronic differentiation of a phenologically atypical population of the pine processionary moth. Proc R Soc Lond B 274:935–941

Santos H, Burban C, Rousselet J, Rossi J-P, Branco M, Kerdelhué C (2011a) Incipient allochronic speciation in the pine processionary moth *Thaumetopoea pityocampa* (Lepidoptera, Notodontidae). J Evol Biol 24:146–158

Santos H, Paiva M-R, Tavares C, Kerdelhué C, Branco M (2011b) Temperature niche shift observed in a Lepidoptera population under allochronic divergence. J Evol Biol 24:1897–1905

Santos HM, Paiva M-R, Rocha S, Kerdelhué C, Branco M (2013) Phenotypic divergence in reproductive traits of a moth population experiencing a phenological shift. Ecol Evol 3(15):5098–5108

Sauné L, Abella F, Kerdelhué C (2015) Isolation, characterization and PCR multiplexing of 17 microsatel- lite loci in the pine processionary moth *Thaumetopoea pityocampa* (Lepidoptera, Notodontidae). Con- serv Genet Resour 7:755–757

Seehausen O, Terai Y, Magalhaes IS et al (2008) Speciation through sensory drive in cichlid fish. Nature 455(2):620–626

Servedio MR, Van Doorn GS, Kopp M, Frame AM, Nosil P (2011) Magic traits in speciation: ‘magic’ but not rare? Trends Ecol Evol 26:389–397

Stuessy TF, Weiss-Schneeweiss H, Keil D (2004) Diploid and polyploid cytotype distribution in *Melampo- dium cinereum* and *M. leucanthum* (Asteraceae, Heliantheae). Am J Bot 91(6):889–898

Tamburini G, Marini L, Hellrigl K, Salvadori C, Battisti A (2013) Effects of climate and density-dependent factors on population dynamics of the pine processionary moth in the Southern Alps. Clim Change 121:701–712

Taylor RS, Friesen VL (2017) The role of allochrony in speciation. Mol Ecol 26:3330–3342

Thibert-Plante X, Gavrilets S (2013) Evolution of mate choice and the so-called magic traits in ecological speciation. Ecol Lett 16:1004–1013

Vähä J-P, Primmer CR (2006) Efficiency of model-based Bayesian methods for detecting hybrid individuals under different hybridization scenarios and with different numbers of loci. Mol Ecol 15:63–72

Weinstein BG, Graham CH, Parra J (2017) The role of environment, dispersal and competition in explaining reduced co-occurrence among related species. PLoS ONE 12(11):e0185493

Weiss-Schneeweiss H, Emadzade K, Jang T-S, Schneeweiss GM (2013) Evolutionary consequences, con- straints and potential of polyploidy in Plants. Cytogenet Genome Res 140:137–150

Wright S (1943) Isolation by distance. Genetics 28:114–138

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

Christian Burban¹ ID · Susana Rocha² ID · Raphaël Leblois³,⁴ ID · Jean-Pierre Rossi³ ID · Laure Sauné³ · Manuela Branco² ID · Carole Kerdelhué³ ID

1 BIOGECO, INRA, University of Bordeaux, 33610 Cestas, France

2 CEF, Instituto Superior de Agronomia (ISA), Universidade de Lisboa (ULisboa), 1349-017 Lisbon, Portugal

3 CBGP, INRA, CIRAD, IRD, Montpellier SupAgro, University of Montpellier, 34988 Montferriez-sur-Lez, France

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