Identification of winter moth (Operophtera brumata) refugia in North Africa and the Italian Peninsula during the last glacial maximum

Abstract Numerous studies have shown that the genetic diversity of species inhabiting temperate regions has been shaped by changes in their distributions during the Quaternary climatic oscillations. For some species, the genetic distinctness of isolated populations is maintained during secondary contact, while for others, admixture is frequently observed. For the winter moth (Operophtera brumata), an important defoliator of oak forests across Europe and northern Africa, we previously determined that contemporary populations correspond to genetic diversity obtained during the last glacial maximum (LGM) through the use of refugia in the Iberian and Aegean peninsulas, and to a lesser extent the Caucasus region. Missing from this sampling were populations from the Italian peninsula and from North Africa, both regions known to have played important roles as glacial refugia for other species. Therefore, we genotyped field‐collected winter moth individuals from southern Italy and northwestern Tunisia—the latter a region where severe oak forest defoliation by winter moth has recently been reported—using polymorphic microsatellite. We reconstructed the genetic relationships of these populations in comparison to moths previously sampled from the Iberian and Aegean peninsulas, the Caucasus region, and western Europe using genetic distance, Bayesian clustering, and approximate Bayesian computation (ABC) methods. Our results indicate that both the southern Italian and the Tunisian populations are genetically distinct from other sampled populations, and likely originated in their respective refugium during the LGM after diverging from a population that eventually settled in the Iberian refugium. These suggest that winter moth populations persisted in at least five Mediterranean LGM refugia. Finally, we comment that outbreaks by winter moth in northwestern Tunisia are not the result of a recent introduction of a nonnative species, but rather are most likely due to land use or environmental changes.

One species that frequently reaches outbreak densities across much of its native range of Europe and western Asia (Elkinton, Boettener, Liebhold, & Gwiazdowski, 2015) is the winter moth, Operophtera brumata (Lepidoptera: Geometridae). Throughout its distribution, winter moth is a serious pest of many broadleaf tree and shrub species, including oak (Quercus spp.), maple (Acer spp.), and birch (Betula spp.) trees (Wint, 1983). Winter moth is also present in North Africa, where it was first reported in the 1970s (Ferguson, 1978;Hausmann & Viidalepp, 2012). More recently it has been observed at outbreak densities in oak forests (Quercus spp.) in northwestern Tunisia beginning in 2009 (Mannai, Ezzine, Nouira, & Ben Jamâa, 2015), though the reasons for these outbreaks are currently unknown. In this region, mature oak forests extend mainly to the regions of Kroumirie and Mogods, from the Algerian-Tunisian border in the west, to Lake Ichkeul in the east (Debazac, 1959;DGF, 2005), and are home to five oak species: Quercus suber L., Quercus canariensis Willd, Quercus ilex L., Quercus coccifera L., as well as a Tunisian-Algerian endemic species Quercus afares Pomel (Boudy, 1950). In Tunisia, winter moth has mainly been observed on Q. canariensis in both mixed oak (Aïn Zena,36.435°N,8.515°E) and pure Q. canariensis oak (Mzara Forest,36.769°N,8.72°E) forests, with lower densities being found on Q. afares and Q. suber compared with Q. canariensis in Aïn Zena. These differences may be the result of differences in budburst timing of different oak species (Mannai, Ezzine, Hausmann, Nouira, & Ben Jamâa, 2017).
As a result of the movement of nursery stocks, invasive populations of winter moth have also been introduced to North America (Ferguson, 1978), with populations established in Nova Scotia in the 1920s (MacPhee, 1967), Oregon in the 1950s (Kimberling, Miller, & Penrose, 1986), British Columbia in the 1970s (Gillespie, Finlayson, Tonks, & Ross, 1978), and in the northeastern United States in the 1990s (Elkinton et al., 2015). Therefore, we were curious as to whether or not the recent reports of outbreaks in northwestern Tunisian oak forests were the result of a recent introduction (as has been the case in North America), or whether a native population has achieved outbreak levels for some other reason. Previously, Mannai et al. (2017) explored the native/nonnative status of outbreaking Tunisian winter moth populations by sequencing a fragment of the mitochondrial (mtDNA) locus cytochrome oxidase I (COI) from field-collected larvae in Tunisian oak forests, and compared these sequences to those published in NCBI GenBank and the Barcode of Life Data Systems (BOLD) database. This work confirmed the field identifications of winter moth in northwestern Tunisia.
Unfortunately, due to the overall lack of genetic diversity of this COI fragment, the native/nonnative status of these populations could not be determined. This low diversity in COI was similarly observed by Gwiazdowski, Elkinton, deWaard, and Sremac (2013) where winter moth samples collected from across Europe were predominantly assigned to two haplotypes that differed by only two substitutions (~0.5% divergence), thereby hampering the ability of COI data alone to determine relationships among populations.
Previously, we examined the use of glacial refugia during the last glacial maximum (LGM) by winter moth .
In that work, we genotyped field-collected individuals from across Europe with 24 highly polymorhic microsatellie loci. In contrast to mtDNA sequence data that showed limited biogeographic structure between winter moth populations, we found evidence that individuals could be assigned to either an eastern or a western European genetic lineage, with a clear hybrid zone in central Europe . Using approximate Bayesian computation (ABC) methods, we also hypothesized that the contemporary genetic diversity of European populations of this species was largely shaped by its retreat to refugia in the Iberian, Balkan, and Caucasus regions, and by hybridization between the Iberian and Balkan lineages postrecolonization.
Absent from the sampling in Andersen et al. (2017), however, were moths from the Italian peninsula and from North Africa, both regions known to have acted as glacial refugia for other species (see Hewitt, 2000;Schmitt, 2007). Therefore, in this study we obtain new winter moth samples from Tunisia, Italy, and additional specimens from Spain, to further explore the biogeographic history of winter moth and its use of potential glacial refugia during the LGM. LGM or whether they were recently introduced to the region, and if so from where.

| Data collection
Adult male winter moths were collected using pheromone-baited traps (Elkinton et al., 2010;Elkinton, Lance, Boettner, Khrimian, & Leva, 2011) deployed in Barcelona, Catalonia, Spain, and in Ortì, Calabria, Italy. After moths were collected from the traps, they were placed in glassine envelopes (Uline Corporation, USA), and stored at −80°C. In addition, we sampled larval winter moth collected from trees in Mzara Forest, Tunisia (pure forest of Q. canariensis, alt 653 m.). Once at the laboratory, larvae were placed in 96% ethanol, and then stored at −80°C prior to molecular analyses. See Table 1 for complete locality information. To this dataset, we added genotype scores from individuals presented in Andersen et al. (2017), including all 28 individuals from Tbilisi, Republic of Georgia, 30 randomly selected individuals from Germany, 30 randomly selected individuals from Serbia, and 20 randomly selected individuals from Spain (bringing the total for Spanish samples to 30).

| DNA extraction and microsatellite amplification
Whole genomic DNA was extracted using the Omega Bio-tek E.Z.N.A. ® Tissue DNA Kit (Omega Bio-tek) following the manufacturer's instructions. For adult males, prior to DNA extraction, the uncus and wings were removed and stored as morphological vouchers. Vouchers were deposited at the Yale Peabody Museum of Natural History, New Haven, Connecticut, USA. Specimen was then homogenized using 3/16″ stainless steel beads (GlenMills Inc.) using a FastPrep-24 Sample Homogenizer (MP Biomedicals). Twenty-four polymorphic microsatellite loci were amplified following protocols presented in Havill et al. (2017), and genotyping was performed at the DNA Analysis Facility on Science Hill at Yale University in comparison to the GeneScan 500 LIZ size standard (Thermo Fisher Scientific) using a 3730xl DNA Analyzer (Thermo Fisher Scientific).
Fragment lengths were scored using the microsatellite plugin in the software program Geneious v. R11 (https ://www.genei ous.com). To address potential variation in allele fragment lengths generated during different genotyping runs (see Hoffman & Amos, 2005), the same laboratory methods, personnel, sequencing instrument, and allele bins in Geneious were used for all runs, including those performed by Andersen et al. (2017). Individuals for which ≥20 microsatellite loci were successfully amplified were included in subsequent analyses.
This subset included up to 30 randomly selected individuals from TA B L E 1 Locality collection information including accession numbers, the number of samples collected (N), the collector, collection date (Date), and geographical coordinates (i.e., latitude and longitude) in decimal-degree format

| Bayesian genetic clustering
The assignment probabilities (Q) of each individual to distinct genetic clusters (K) were estimated using the Bayesian genetic clustering software structure v.2.3.2 (Falush, Stephens, & Pritchard, 2003;Pritchard, Stephens, & Donnelly, 2000) for values ranging from K = 1 to K = 8. For each value of K, ten independent analysis were run, each with random starting values, a burn-in period of 100,000 generations, and a total runtime of 1,000,000 generations. For each run, we utilized the admixture model, with correlated allele frequencies using default settings for both. The optimal number of clusters was then estimated using the Delta-K method of Evanno, Regnaut, and Goudet (2005) as implemented in the software package structureHarvester (Earl & vonHoldt, 2012), and values for Q were averaged across runs using Clumpak (Kopelman, Mayzel, Jakobsson, Rosenberg, & Mayrose, 2015).  (Chapuis & Estoup, 2007), the subsequent pairwise corrected F ST distance matrix was then used to calculate a "NeighborNet" network with the program SplitsTree v.4.14.2 (Huson & Bryant, 2006). The probability of population differentiation was then estimated using GenePop (Raymond & Rousset, 1995;Rousset, 2008).

TA B L E 2
Population genetic diversity including the number of individuals genotyped (n), the average number of alleles per locus (Na), the effective number of alleles (Eff_Na), the average observed heterozygosity across loci (H o ), the average expected heterozygosity within each population (H s ), the total expected heterozygosity (H t ), the average inbreeding coefficient (G IS ), the presence of deviation from Hardy-Weinberg Equilibrium (HWE) using GenoDive, and the average null allele frequency (Null) estimated using FreeNA

| Historical demography
To estimate the divergence time between the Tunisian population and Eurasian populations of winter moth using ABC as implemented in the software DiyABC v.2.1.0 (Cornuet et al., 2014(Cornuet et al., , 2008, we used a two-step approach to limit the number of potential scenarios (see For both sets of simulations, the following summary statistics were used: one sample; mean number of alleles, mean genic diversity, mean size variance: two sample; F ST , classification index, and (dμ) 2 distance. Distributions for parameter priors are presented in Table 3. The scenarios were then compared using and the "Logistic We then estimated error rates by implementing the "Evaluate confidence in scenario choice" analysis in DiyABC.

| Data collection, DNA extraction, and microsatellite amplification
After filtering individuals from which less than 20 of 24 polymorphic microsatellites were amplified, ten winter moth individuals

| Bayesian genetic clustering
The probability of assignment (Q) of individuals to major genetic clusters (K) for K = 1 through K = 8 after averaging across independent Structure runs is presented in Figure 1. Assignments for major and minor clusters, as well as the number of runs supporting each cluster, after summary in Clumpak are presented in Figure S3.
Results from structureHarvester indicated that the optimal number of genetic clusters (K) present in the dataset was K = 6 ( Figure   S4). At K = 6, all ten independent Structure runs were assigned to one major cluster using Clumpak, and these genetic clusters predominantly corresponded to the geographical origin of samples (Figures 1 and 2).

| Genetic diversity and differentiation
Population genetic statistics for each population are presented in  Lines widths are drawn relatively proportional to the mean estimated population size. At each node, the mean divergence/ merger time is shown along the y-axis (1 ka = 1,000 years). Ninetyfive percent confidence intervals for population sizes, divergence/ merger times, and for the amount of admixture (ra) are presented in Table 3 Spain Germany Serbia Georgia Values for all pairwise comparisons are presented in Table S1.
The NeighborNet analyses indicated that all sampled populations showed a high degree of genetic distinction, with all but the German population being placed at the end of a long branch (Figure 3).

| Historical demography
Results from simulations conducted during the first step in our DiyaBc analyses indicated that the Ortì, Italy population was most closely related to the Spanish population ( Figure S5 and Table S2).
Model-checking analyses for the posterior distributions for summary statistics based on this dataset formed a distinct cloud around the observed dataset ( Figure S6). Evaluations of the confidence in the scenario choice indicated a posterior predictive error rate of 0.064.
Results for simulations conducted as part of our second step of the analyses also indicated that the Tunisian population was most closely related to the Spanish one (Figure 4), and this result was strongly supported by the logistic regression test ( Figure S7 and Table S3).
Model-checking analyses for the posterior distributions for summary statistics based on this dataset formed a distinct cloud around the observed dataset ( Figure S8). Evaluations of the confidence in the scenario choice indicated a posterior predictive error rate of 0.287. Based on this scenario (Figure 4), the Tunisian population likely diverged from the Spanish population ~16,200 years ago (95% CI = 10,100-19,800 years ago; Table 3). All estimates for demographic parameters for the most likely scenarios were calculated based on the analysis of 1,000,000 simulated datasets and are presented in Figures S9 and S10, for the scenarios estimating the relationship of the southern Italian population and the Tunisian populations, respectively, and summarized for the Tunisian scenario in Table 4.

| D ISCUSS I ON
We previously demonstrated that winter moth persisted in at least three glacial refugia (Iberian, Aegean, and Caucasus) during the LGM . By analyzing newly sampled populations in North Africa and the Italian Peninsula, we provide evidence for two additional LGM refugia. Contemporary populations sampled in Tunisia and southern Italy differ from all other European populations (Figures 1 and   2) suggesting that the lineages that derived from these refugia did not recolonize Europe following the retreat of glaciers following the LGM Note: Mean, median, mode, and 95% confidence interval (CI) for each parameter are presented based on the summary of 1 million simulated datasets for the best supported scenario (Scenario 4).

TA B L E 4
Demographic parameter estimates from DiyABC for the best supported scenario presented in Figure  3 representing the relationship of the Mzara Forest, Tunisia population to other sampled localities as presented in Figure 3 ( Figure 4). For the southern Italian lineage, it is possible that the Alps presented a significant geographic barrier to their recolonization as has been seen in other systems (see Hewitt, 2000 Geometridae), and possibly winter moth to grow to outbreak densities (Mannai, 2017). In addition, climate change in northern Africa might have influenced the population dynamics of winter moth, as already demonstrated in northern Scandinavia (Jepsen, Hagen, Ims, & Yoccoz, 2008). Indeed, northern Tunisia saw dramatic changes in mean annual temperature during the 20th century, with an increase of 2.5-3°C on average, and future climate models predicting an additional increase of 3-4°C coupled with decreased rainfall in the 21st century (Hulme, Doherty, Ngara, New, & Lister, 2001). Environmental factors can play an imortant role in shaping plant-herbivore interactions (Kuglerová, Skuhrovec, & Münzbergová, 2019). It is therefore possible that changes in temperature and/or precipitation might promote winter moth populations either directly by reducing developmental times and/or indirectly by limiting the amount of secondary compounds produced by stressed oak trees. It is also possible that human-mediated disturbances in the region might be limiting the impact of native natural enemies (such as pathogens, parasitoids, and predators).
Further work in this region is needed to understand what factors are resulting in the outbreaks of winter moth in order to preserve these culturally and ecologically important native oak stands.

| CON CLUS IONS
Here, we find the use of two additional glacial refugia, North Africa and southern Italy, by winter moth during the LGM. Our results indicate that populations that utilized both refugia during the LGM diverged from a population that migrated into the Iberian Peninsula, and that subsequently each population became genetically distinct. Curiously, our findings continue to support the results from Andersen et al. (2017) that southern populations appear to be less genetically distinct than northern populations. As such, the biogeographical history of winter moth might represent an important case study for the role of secondary contact in promoting the genetic diversity of outbreaking pest species. Finally, we hope that this new genetic data, and the discovery of additional biogeographic structuring, can be used in conjunction with the recently published winter moth genome (Derks et al., 2015), to promote the use of winter moth as a model system for future research questions, including studies of local adaptation-a subject for which winter moth has previously been studied in the pregenomics era (Buse & Good, 1996;van Dongen, Backeljau, Matthysen, & Dhondt, 1997;Holliday, 1977; Mohamed Lahbib Ben Jamâa https://orcid. org/0000-0003-1608-2867 Adalgisa Caccone https://orcid.org/0000-0002-8949-9260 Joseph S. Elkinton https://orcid.org/0000-0001-9874-0426

O PE N R E S E A RCH BA D G E S
This article has earned an Open Data Badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at https ://osf.io/yfzh3/ ?view_ only=11a9e 5ca6c 2743b 38abf 5aa5f d2bec5f.

DATA AVA I L A B I L I T Y S TAT E M E N T
Genotype data for this study is provided as a tab-delimited file titled "WinterMothTunisiaStructure.txt" uploaded as supplemental material.