Postglacial recolonization shaped the genetic diversity of the winter moth (Operophtera brumata) in Europe

Abstract Changes in climate conditions, particularly during the Quaternary climatic oscillations, have long been recognized to be important for shaping patterns of species diversity. For species residing in the western Palearctic, two commonly observed genetic patterns resulting from these cycles are as follows: (1) that the numbers and distributions of genetic lineages correspond with the use of geographically distinct glacial refugia and (2) that southern populations are generally more diverse than northern populations (the “southern richness, northern purity” paradigm). To determine whether these patterns hold true for the widespread pest species the winter moth (Operophtera brumata), we genotyped 699 individual winter moths collected from 15 Eurasian countries with 24 polymorphic microsatellite loci. We find strong evidence for the presence of two major genetic clusters that diverged ~18 to ~22 ka, with evidence that secondary contact (i.e., hybridization) resumed ~ 5 ka along a well‐established hybrid zone in Central Europe. This pattern supports the hypothesis that contemporary populations descend from populations that resided in distinct glacial refugia. However, unlike many previous studies of postglacial recolonization, we found no evidence for the “southern richness, northern purity” paradigm. We also find evidence for ongoing gene flow between populations in adjacent Eurasian countries, suggesting that long‐distance dispersal plays an important part in shaping winter moth genetic diversity. In addition, we find that this gene flow is predominantly in a west‐to‐east direction, suggesting that recently debated reports of cyclical outbreaks of winter moth spreading from east to west across Europe are not the result of dispersal.

. For contemporary populations, the effects of these periods of isolation are often observed by the presence of distinct genetic lineages that correspond with the numbers and locations of the glacial refugia utilized by a particular species (Taberlet, Fumagalli, Wust-Saucy, & Cosson, 1998).
In addition to promoting patterns of diversification, the Quaternary climatic oscillations have also had profound impacts on genetic variation within populations. One common finding for postglacial recolonization studies is that populations in geographic regions that acted as glacial refugia during the glacial maximums are more genetically diverse than populations in regions that have been recolonized (i.e., the "southern richness, northern purity" paradigm; Hewitt, 2000;Canestrelli, Cimmaruta, Costantini, & Nascetti, 2006;Babik et al., 2009;Rousselet et al., 2010;Wielstra et al., 2013;Tison et al., 2014;Mezzasalma et al., 2015;Vitales, Garcia-Fernandez, Garnatje, Pellicer, & Valles, 2016). As populations of species recolonized temperate regions, they likely extended their distributions in a stepping-stone manner whereby multiple genetic bottleneck events were associated with continued northward range expansions (e.g., Tison et al., 2014).
However, during these cycles genetic lineages that were reproductively isolated during glacial maximums might also have subsequently come into secondary contact with one another during periods of recolonization and the signatures of these secondary contact events have been shown by the presence of detectible hybrid zones (Hewitt, 2000;Schmitt, 2007;Schmitt & Müller, 2007). Hybridization-here used in the context of admixture between distinct populations (Harrison & Larson, 2014)-has long been known to play an important role in shaping patterns of host-associations (e.g., Feder et al., 2003), and has the potential to greatly increase measures of genetic diversity (Verhoeven, Macel, Wolfe, & Biere, 2011) and even reverse the process of speciation (Seehausen, Takimoto, Roy, & Jokela, 2007). How hybridization might affect the "southern richness, northern purity" paradigm, however, is unclear. In part this may be because the strength of the signal of this paradigm may vary depending on the recolonization paths taken by a species, and the geographic locations of secondary contact zones (see Hewitt, 2000 andSchmitt, 2007).
Species of Lepidoptera have been crucial for examining both of the above biogeographic patterns. This is could be due to the fact that many species of Lepidoptera have the ability to disperse long distances and can quickly fill open niche space at the species-level, while dispersal distances for individuals are generally quite small, allowing for the preservation of biogeographic structure (Schmitt, 2007). These short dispersal distances may also be important for the promotion of local adaptation (Tison et al., 2014), and recent studies of Lepidoptera have also highlighted the importance of micro-refugia (von Reumont, Struwe, Schwarzer, & Misof, 2012), postcolonization hybridization (Schmitt & Müller, 2007), and colonization from refugia outside of Europe (Habel, Lens, Rodder, & Schmitt, 2011) for patterns of Eurasian biogeography. However, one species of Lepidoptera for which elucidating historical biogeographic patterns has been particularly elusive is the winter moth, Operophtera brumata (Lepidoptera: Geometridae). This species has a broad range of woody host plants including, but not limited to, oak (Quercus), maple (Acer), and birch (Betula) trees in Europe (Wint, 1983), and has long been studied as a model organism for studies of local adaptation (Tikkanen & Lyytikainen-Saarenmaa, 2002;Tikkanen, Woodcock, Watt, & Lock, 2006;Van Dongen, Matthysen, & Dhondt, 1996), and population ecology (Hassell, 1968;Macphee, Newton, & McRae, 1988;Varley & Gradwell, 1960, 1968, and has been at the center of an ongoing debate in regards to whether cyclical outbreaks of geometrid moths (including winter moth) move across western Eurasian from east to west approximately every 10 years (Tenow et al., 2013;but see Jepsen, Vindstad, Barraquand, Ims, &Yoccoz, 2016 andTenow, 2016). In addition, this species has also been the subject of much genetic research, including population structure (Leggett et al., 2011;Van Dongen, Backeljau, Matthysen, & Dhondt, 1998), hybridization (Elkinton, Liebhold, Boettner, & Sremac, 2014;Elkinton et al., 2010;Havill et al., 2017), and a draft genome for this species was recently published (Derks et al., 2015). Yet, the only continent-scale study of winter moth phylogeography (Gwiazdowski, Elkinton, Dewaard, & Sremac, 2013) found little evidence for geographically distinct genetic lineages when the mitochondrial locus cytochrome oxidase I (COI) was analyzed-although it did find support for a division between northern (i.e., Norway, Scotland, and Sweden) and southern European (i.e., all other sampled locations) populations.
This result was unexpected, given ample evidence of local adaptation by winter moth populations (Tikkanen & Lyytikainen-Saarenmaa, 2002;Tikkanen et al., 2006;Van Dongen et al., 1996), and the fact that females are flightless and that males are considered poor dispersers (Van Dongen et al., 1996).
Therefore, we were interested if a recently developed set of 24 highly variable microsatellite loci  could be used to further elucidate patterns of genetic diversity for populations of winter moth. Specifically we were interested the following: (1) whether individual winter moths could be assigned to distinct genetic clusters using Bayesian clustering analyses and genetic distance methods and whether these genetic groupings corresponded with geographically distinct glacial refugia, (2) whether divergence times between contemporary populations in regions that likely acted as glacial refugia corresponded with isolation during the last glacial maximum (LGM), (3) whether winter moth populations display evidence of decreasing genetic diversity with increases in latitude, and (4) whether contemporary gene flow could be detected and whether the direction of gene flow can inform the debate about cyclical outbreaks of geometrid moths.
For this study, male moths were collected at 44 locations in Eurasia (Table S1, Figure S1) using sex pheromone-baited traps deployed in tree canopies that attract winter moth and several closely related congeners (Elkinton, Lance, Boettner, Khrimian, & Leva, 2011;Elkinton et al., 2010). Moths were collected from the traps and then placed in glassine envelopes (Uline Corporation, USA), before storage at −80°C.

| Winter moth postglacial recolonization
After genotyping, individuals from which ≥20 microsatellite loci were successfully amplified were included in subsequent analyses.
Genotype data for each included individual is provided as a single-line formatted file titled (Appendix S2).

| Bayesian genetic clustering
Previous postglacial recolonization studies have found that the Quaternary climatic cycles have played an important role in shaping genetic diversity. This diversity is evident by the presence of two to three distinct genetic lineages in Eurasia (Hewitt, 2000;, with well-established hybrid zones where these lineages have come into secondary contact with each other (see Hewitt, 2000;Figure 3). Preliminary analyses of our dataset using RetEucetuED v.2.3.2 (Falush, Stephens, & Pritchard, 2003;Pritchard, Stephens, & Donnelly, 2000) indicated that individual winter moths could be assigned to an optimal number of two distinct genetic clusters (Appendix S1). Therefore, we estimated the probability of assignment (Z) that an individual could be classified as belonging to one of these two clusters or to hybrid categories (F1, F2, or backcrosses), using the software program ADwHybEiNR v.1.1.b3 (Anderson, 2008;Anderson & Thompson, 2002). Four independent runs, each of 1 million generations, discarding the first 100,000 burn-in generations, were analyzed using random starting values, and uniform priors for the estimates of θ (i.e., allele frequencies) and π (i.e., mixing proportions). Results were then averaged across runs, and visualized using EcMap v.10.3.1 (Esri, Redlands, CA, USA).

| Population genetic statistics and genetic distances
For all localities from which ≥10 individuals were successfully genotyped, we calculated standard population genetic summary statistics, including: the number of individuals genotyped (n), the average number of alleles per locus (N a ), the effective number of alleles (Eff_N a ), the average observed heterozygosity across loci (H o ), the average inbreeding coefficient (G IS ), and the presence of deviation from Hardy-Weinberg Equilibrium (HWA) using GDAoNivD v.2.0b27 (Meirmans & Van Tienderen, 2004). In addition, we tested each locus-pair within populations for evidence of Linkage-Disequilibrium (LD) using GDADPoP (Raymond & Rousset, 1995;Rousset, 2008), and we estimated the frequency of null alleles present for each locus and population using FEDDNA (Chapuis & Estoup, 2007). The degree of genetic differentiation (F ST ) among populations was calculated (accounting for the presence of null alleles) in FEDDNA, and Jost's estimator of actual differentiation (D est ) was calculated using RmoGN v.1.2.5 (Crawford, 2010). To examine whether genetic differentiation was correlated with geographic distances, the presence of isolation-by-distance was estimated by performing linear-regression of F ST / (1-F ST ), calculated in GDADPoP, versus the linear distance between populations (km) in R v. 3.1.3 (R Core Team 2015). Finally, to examine historical relationships among populations a "NeighborNet" network was reconstructed with the above null-allele corrected F ST estimates with the program RPaietRetEDD v.4.14.2 (Huson & Bryant, 2006).

| Historical demography
To determine whether divergence timing for contemporary Eurasian populations of winter moth correspond with genetic isolation during the LGM, population history parameters were estimated using Approximate Bayesian Computation (ABC) with the software program NiyABC v.2.1.0 (Cornuet et al., 2008(Cornuet et al., , 2014. For these analyses we included samples from Spain, Serbia, and Georgia as representative of three potential glacial refugia (Iberaian, Balkan, and Caucasus, respectively). Due to limited sampling from Italy (Italian refugium), we did not include samples from this region in our analysis. In addition, we included samples from Germany to represent the product of postglacial recolonization. We evaluated eight different recolonization scenarios ( Figure S2). Three scenarios included simple bifurcating trees that differed in whether the German population was derived from the Spanish, Serbian, or Georgian populations, and five scenarios that included the German population being derived from admixture. All scenarios included multiple population size parameters to allow for changes in population sizes following splitting and/or merging events.
Prior information for all parameters, including minimum and maximum values, as well as distributions, is presented in Table S2. Default mutation model parameters were used, except for the following: we set the minimum mean mutation rate to 1 × 10 −5 , and maximum values for the Mean and Individual locus coefficient P's to 1.0. Visualization of preliminary results based on principal component analyses (PCAs) within NiyABC indicated that these changes in the mutation rate and the Mean and Individual locus coefficient P's improved the shape of the cloud of simulated datasets. In addition, visualization of preliminary results based on PCAs indicated that the removal of four loci (02339, 00925, 02191, and 12042) also improved the shape of the cloud of simulated datasets. These four loci each had allelic ranges ≥120 repeat units (i.e., the number of tandem repeats). Given that lepidopteran genomes have high rates of duplication events (Edger et al., 2015;You et al., 2013), and that the development of microsatellite loci for these taxa have been difficult due to high frequencies of null alleles and associations with transposable elements (see citations in Sinama et al., 2011), it is possible that there might be an association between these four loci and transposable elements and/or duplication events. We will continue to monitor the appropriateness of using these four loci. We then generated a reference table by simulating 8 million datasets (1 million datasets per scenario). Pre-evaluations of the simulated datasets were conducted using PCA, and scenarios were compared using both the "Direct" (i.e., the proportion of simulated datasets from each scenario "closest" to the sample dataset) and "Logistic Regression" (i.e., a logistic regression analysis of the probability of deviations between summary statistics calculated for the simulated datasets and the sample dataset) tests in NiyABC.

| Genetic diversity and latitude
To examine whether measures of genetic diversity (N a , Eff_N a , H o , and G IS ) were correlated with latitude, we constructed generalized additive models (GAM) using the R package "mgcv" (Wood, 2011). Results were then visualized using the "plot.gam" function.

| Contemporary rates of gene flow
The rates of contemporary gene flow among populations in Eurasia were examined by estimating the proportion of migrant individuals in each population using the software program byDRRR v.4.0 (Wilson & Rannala, 2003). To reduce the number of pair-wise comparisons, individuals were grouped by country. Gene flow rates were then estimated for all 205 possible comparisons with four independent analyses, each with a random starting seed, runtimes of 10 million generations, burn-in periods of 1 million generations, and by setting all mixing parameters to 0.8. Results were then summarized across analyses using etEcDE v.1.6.0 (Rambaut & Drummond, 2007), and all migration rates whose 95% confidence intervals (mean ± [1.96 × SE]) did not include 0 are reported as significant. Average rates of gene flow between eastern (Austria, Czech Republic, Georgia, Poland, Serbia, Slovakia) and western (England, France, Germany, Italy, Norway, Scotland, Spain, Sweden, Switzerland) European countries were then calculated. The difference between these rates was used as the basis for simulations to determine if one or the other populations would eventually become ubiquitous and if so in how many generations. Simulations were carried out using the software program PoPG v. 4.03 (© University of Washington and Joseph Felsenstein), under a simple two allele model with equal starting frequencies and the observed difference in migration rates. Ten independent simulations were conducted; each with random starting seeds, and the number of generations to fixation of one of the alleles averaged across runs was recorded.

| Bayesian genetic clustering
After filtering the dataset to include only those individuals from which ≥20 microsatellite loci were genotyped, analyses included 669 individual moths. The proportional assignment of individuals by ADwHybEiNR is presented in Table S1 and summarized by country in Figure 1. These results showed that 145 individuals could be assigned to Cluster 1 (123 individuals with high probability of assignment Z ≥ 0.8, and 22 individuals with moderate probabilities of assignment 0.5 ≤ Z < 0.8).
Four hundred and thirty-four individuals were assigned to Cluster 2 Cluster 2 (386 with high probability of assignment and 47 with moderate probability of assignment). In addition, 78 individuals were identified as hybrids (13 with high probability of assignment and 65 with moderate probabilities of assignment). All hybrids were assigned to the F2 hybrid class. Thirteen individuals were classified as "unassigned" due to low probabilities of assignment (Z < 0.5 to any category). In general, assignments showed a clear distinction between eastern and western Europe, with moths from Georgia, Poland, and Serbia being primarily assigned to Cluster 1 (100%, 97%, and 88%, respectively), and individuals from Spain, Scotland, and England being primarily assigned to Cluster 2 (100%, 100%, and 96%, respectively).
In addition, while the overall number of hybrids was low, the greatest proportions of hybrid individuals were found in Austria, Sweden, and the Czech Republic (39%, 33%, and 32%, respectively), roughly corresponding with a previously identified hybrid zone in central Europe (see Figure 3 in Hewitt, 2000).

| Population genetic statistics and genetic distances
Population genetic statistics for 26 locations from which ≥10 individual moths were sampled showed that the average number of alleles per locus was N a = 8.63, the average effective number of alleles per locus was Eff_N a = 4.34 the average observed heterozygosity was H o = 0.58, and the average inbreeding coefficient was G IS = 0.17 (Table 1). The average presence of null alleles across loci was 6.7%.
All locations showed significant (p ≤ .05) deviations from HWE, while only 13 of the 276 pair-wise comparisons showed significant (p ≤ .05) evidence of LD between loci after using the Benjamini-Hochberg false-discovery method (Benjamini & Hochberg, 1995).
Among populations, the average values for pair-wise differentiation for F ST (corrected for the presence of null alleles) was 0.060 and for D est was 0.091 (Table S3). Genetic differentiation also showed strong patterns of isolation-by-distance (p = 2.2 × 10 −16 , r 2 = .5235, df = 323) among sampled populations. Network analysis indicated that most populations could be divided into two large subnetworks (Figure 2).
One of these subnetworks included populations from Austria, Czech Republic, Poland, Serbia, and Slovakia (labeled "Eastern Europe"), while the other subnetwork included populations from England, France, Germany, Italy, and Switzerland (labeled "Western Europe"). In addition, populations of moths from Scotland, Spain, Norway, and the Republic of Georgia were each placed at the ends of long branches.

| Historical demography
Comparisons of the posterior probabilities of the scenarios included in the NiyABC analysis indicated that Scenarios 4 and 7 were the best fit for the data ( Figure S3). Under the Direct approach simulations based on Scenario 7 were consistently the most supported, whereas under the Logistic Regression approach support moves from Scenario 4 to Scenario 7 with increasing distance from the sample data ( Figure S3).
These two scenarios differ in regards to whether the Serbian population is sister to the Georgian or the Spanish populations, yet both scenarios indicate that the German populations are derived from recent admixture between the populations that recolonized Europe from the Iberian and Balkan refugia. The branching pattern observed in Scenario 7 also matches the results observed from the Bayesian clustering and distance analyses (i.e., Georgia being most closely related to eastern European populations). Results for all estimated parameters for these two scenarios are presented in Table 2, and the mean values F I G U R E 1 Proportional assignment of individuals to hybrid categories, and contemporary gene flow rates among Eurasian countries. Pie charts, representing the proportion of individuals in each country assigned to different hybrid categories calculated with ADwHybEiNR (Z ≥ 0.5 to Cluster 1 = Black; Z ≥ 0.5 to Cluster 2 = White; Z ≥ 0.5 to F2 Hybrid = Dark Gray; Z < 0.5 to any category = Light Gray), were generated in EcMap v.10.3.1 and visualized using the Europe Albers Equal Area Conic projection. The location of each chart is placed approximately in the center of each country. The label for the Czech Republic has been abbreviated "C.R." Gene flow rates are presented as the proportion of individuals that are likely migrants from a different Eurasian country as calculated by byDRRR. The thickness of the curved lines represents the proportion of migrants (see legend, top right), and the arrowhead indicates the direction of gene flow. The distribution of glacial refugia during the LGM (diagonal shading) is drawn as per Figure 3 in Schmitt (2007) (2007) Fig. 3 Glacial Refugia: for divergence times and population sizes are presented graphically in Figure 3. In both scenarios, the German and Serbian effective population sizes are significantly larger (non-overlapping 95% CIs) than the Georgian and Spanish populations. The mean estimates for population sizes were all larger than the estimates for the ancestral population 2 (NA5), although again these were only significant for the German and Serbian populations, and there were no significant differences between the ancestral population 1 (NA3) and any contemporary popu-

| Genetic diversity and latitude
Results from the GAM analyses found no significant relationships between measures of genetic diversity (N a , GIS, and H o ) and latitude.

| Winter moth postglacial recolonization
We find that populations of winter moth can be broadly assigned to two geographically and genetically distinct populations (Figure 1), which diverged ~ 20 ka (Figure 3). While this finding differs slightly from previous evidence that suggested northern winter moth populations might be distinct from southern populations (Gwiazdowski et al., 2013), the observed results are consistent with many previous studies that have highlighted the importance of the Quaternary climate cycles in shaping patterns of genetic diversity and geographic distributions for other lepidopteran species in Europe and western Asia (e.g., Hewitt, 2000;Schmitt, 2007;Schmitt, Rober, & Seitz, 2005;Schmitt & Seitz, 2001a,b;Wu et al., 2015).
For winter moth, our findings suggest that populations became reproductively isolated in glacial refugia located both in the Iberian and divergence is similar to that which has been observed for other species of Lepidoptera (e.g., Schmitt & Seitz, 2001a; F ST = 0.060), but less than that which has been observed for species with recognized subspecies (e.g., Schmitt & Seitz, 2001b; F ST = 0.149). Beginning ~20 ka the glaciers began to retreat (Clark et al., 2009), and during the subsequent reforestation of the European continent, these two lineages expanded northward, coming into secondary contact ~ 5 ka along a well-documented hybrid zone in Central Europe (see Figure 3 in Hewitt, 2000).
Following the northward recolonization of the European continent, populations in Spain and Georgia-regions that acted as glacial refugia during the last glacial cycle-likely became reproductively isolated in these regions as evident by the long branches connecting these populations to other European populations (Figure 2). For the Spanish population, this result might highlight the importance of physical barriers, such as the Pyrenees, in promoting reproductive isolation, as seen for the pine processionary moth, Thaumetopoea pityocampa (Rousselet et al., 2010). Similarly, the distance of the branch leading to the Georgian population might also highlight the importance of geographic isolation, a pattern similar to that found for white oaks (Dumolin-Lapègue, Demesure, Fineschi, LeCorre, & Petit, 1997), one of the primary hosts for winter moth. The Caucasus region has long been recognized as having played an important role as a refugium during glacial cycles (e.g., Babik et al., 2005;Beck, Schmuths, & Schaal, 2007;Dubey, Zaitsev, Cosson, Abdukadier, & Vogel, 2006;Hampe, Arroyo, Jordano, & Petit, 2003;Jaarola & Searle, 2002), however the area connecting the Caucasus region to continental Europe has gone through significant changes in forest cover since the last glacial maximum (Adams, 1997), and it is possible that these changes eventually prevented gene flow between winter moth populations in eastern Europe and the Caucasus region.
Genetic distance analyses also revealed that populations from two countries, Norway and Scotland, that were recolonized after the retreat of the glaciers showed additional evidence of genetic differentiation ( Figure 2). In northern Norway, winter moth recently became established (~ 125 ya), and its introduction has caused extensive defoliations due to frequent outbreaks (Vinstad et al., 2013). populations is also interesting given that populations in Scotland have become a serious pest of Sitka spruce (Hunter, Watt, & Docherty, 1991;Stoakley, 1985), suggesting that this long branch might be the result of micro-evolutionary processes as a result of this significant host shift.

| Genetic diversity and latitude
One common finding for postglacial recolonization studies is that populations in geographic regions that acted as glacial refugia T A B L E 2 Parameter estimates (mean and 95% Confidence Intervals [CI]) from NiyABC for Scenarios 4 and 7 based on 1 million simulated datasets each. These estimates include; current effective population sizes in Georgia (N1), Serbia (N2), Spain (N3), and Germany (N4), as well as effective population sizes for ancestral population 1 (NA3) and ancestral population 2 (NA5). In addition, estimates are shown for the proportion of admixture (ra), time of divergence (t1, t2, and t3); as well as the average mutation rates (μ-mic), the average coefficient of P (p-mic), and the average SNI rates (sni-mic) for the microsatellite loci during the LGM are more genetically diverse than populations in recolonized regions (Babik et al., 2009;Canestrelli et al., 2006;Hewitt, 2000;Mezzasalma et al., 2015;Rousselet et al., 2010;Tison et al., 2014;Vitales et al., 2016;Wielstra et al., 2013). Here we find some evidence that measures of genetic diversity might be greater in middle European latitudes ( Figure S4), a region that is experiencing secondary contact (i.e., hybridization; Figure 1).
Given that hybridization has the potential to provide populations with additional genetic diversity, and that rates of hybridization are much higher than previously suspected (Allendorf, Leary, Spruell, & Wenburg, 2001;Harrison & Larson, 2014) particularly among insect taxa (Schwenk, Brede, & Streit, 2008), secondary contact following postglacial recolonization may play an important, and underreported, role in promoting the genetic diversity of populations in western Eurasia.

| Contemporary rates of gene flow
Another interesting finding from our study involves the level and direction of gene flow among contemporary European populations.
Recent studies have explored whether outbreaks of winter moth occur on a cyclical basis, suggesting a cycle length of approximately 10 years during which winter moths move across Europe predominantly from east to west (Tenow et al., 2013; but see Jepsen et al., 2016 andTenow, 2016). These studies have been important contributions to our understanding of the incidence and causes of cyclic population dynamics and of damaging insect outbreaks, however the factors responsible for these outbreaks remain unclear. Dispersal may be an important driver of outbreaks of herbivorous insects (Liebhold, Koenig, & Bjornstad, 2004;Peltonen, Liebhold, Bjornstad, & Williams, 2002), although it is but one of numerous biotic and F I G U R E 3 The top two scenarios for diversification and recolonization identified by simulations in NiyABC. Populations in each scenario were chosen to represent different glacial refugia (e.g., Spain = Iberian Refugium, Serbia = Balkan Refugium, Georgia = Caucasus Refugium), as well as a population that would have established after the northward retreat of the glaciers following the LGM (e.g., Germany). Both of these supported scenarios include an admixture event between the Serbian and Spanish populations giving rise to the German population. Mean estimates for merging (t1) and divergence (t2 and t3) events are shown along the y-axis along with 95% confidence intervals. The thicknesses of branches are drawn proportional to the mean estimated population size. Estimates for all parameters, including 95% confidence intervals, are presented in Table 2. All eight tested scenarios is show in Figure S2 Spain Germany Serbia Georgia Spain Germany Serbia Georgia Scenario 7 Scenario 4 1-ra ra ra 1-ra abiotic factors that can be responsible for the observed patterns (Hunter & Elkinton, 2000;Price & Hunter, 2015;Riolo, Rohani, & Hunter, 2015). Previous studies have highlighted the ability of winter moth populations to undergo long-distance dispersal (Leggett et al., 2011), likely through a process known as ballooning where larvae secrete a silky thread and are subsequently wind-dispersed (Elkinton et al., 2015). Therefore, long-distance dispersal by larvae may provide a mechanism that drives these outbreak waves. If this were the case, then we might expect levels of gene flow to be primarily in the east to west direction (the proposed direction of the outbreak waves), however; our results suggest that long-distance dispersal occurs primarily in the opposite direction (from west to east). In addition, if dispersal were responsible for these cyclical outbreaks, such frequent dispersal events would have long since obliterated the genetic distinction we discovered between eastern and western populations resulting in the predominance of the eastern genetic lineage. In contrast, our simulation results suggest that the western population will become predominant after another ~2,000 years of continued secondary contact, although populations of the eastern genetic lineage might likely still persist in the Caucasus region due to geographic isolation.