Limited Genetic Structure of Gypsy Moth Populations Reflecting a Recent History in Europe

The gypsy moth, Lymantria dispar, a prominent polyphagous species native to Eurasia, causes severe impacts in deciduous forests during irregular periodical outbreaks. This study aimed to describe the genetic structure and diversity among European gypsy moth populations. Analysis of about 500 individuals using a partial region of the mitochondrial COI gene, L. dispar was characterized by low genetic diversity, limited population structure, and strong evidence that all extant haplogroups arose via a single Holocene population expansion event. Overall 60 haplotypes connected to a single parsimony network were detected and genetic diversity was highest for the coastal populations Croatia, Italy, and France, while lowest in continental populations. Phylogenetic reconstruction resulted in three groups that were geographically located in Central Europe, Dinaric Alps, and the Balkan Peninsula. In addition to recent events, the genetic structure reflects strong gene flow and the ability of gypsy moth to feed on about 400 deciduous and conifer species. Distinct genetic groups were detected in populations from Georgia. This remote population exhibited haplotypes intermediate to the European L. dispar dispar, Asian L. dispar asiatica, and L. dispar japonica clusters, highlighting this area as a possible hybridization zone of this species for future studies applying genomic approaches.


Introduction
The gypsy moth, Lymantria dispar L. (Lepidoptera, Erebidae), is a polyphagous pest feeding on more than 400 plant species with a preference for species within the genus Quercus [1]. Lymantria dispar is native to the Palearctic region [1,2] and became invasive in North America after its intentional introduction for hybridization experiments [3]. Currently, L. dispar is considered one of the most notorious forest pests worldwide [4], with outbreaks in its native range being common particularly in the oak forests of the Mediterranean region, while in Asia the broader spectrum of potential hosts renders gypsy moth capable of even more massive outbreaks [5]. Nonetheless, the potential of gypsy moth populations to erupt varies greatly [6] as in cases of recent invasion an outbreak can be sustained [7]. For all these reasons, gypsy moth, has been in the spotlight of various investigations ranging from population dynamics [8] and dispersal ecology [9,10] to the distribution of its genetic diversity at a global scale [11,12]. Three subspecies are recognized throughout the temperate part of the world due to their economic importance; European (L. dispar dispar) and two Asian (L. dispar asiatica in Russia and Asia and L. dispar japonica in Japan), and each of them is described separately based on morphological and behavioral traits [1]. Nonetheless, the status of the two Asian subspecies still remains disputable [13], as recent molecular investigations argue for an even more elaborate pattern among the far east populations [12]. However, female flight capability turned out to be the most important behavioral trait distinguishing Asian subspecies from European [10,14].
First genetic analysis based on DNA sequences resulted in developing an mtDNA marker that verified the close resemblance of North American and European population, that both were distinctly separated from Asian populations [11]. This pattern was confirmed in subsequent studies which showed that North American populations were related to French populations and additionally that Asian population were separated into a western and eastern genetic entity [14], where Japanese populations form a separate cluster differentiated from Western Asian populations [15].
Until recently, studies on European populations of gypsy moth were mostly integrated into broader investigations that aimed at a global overview of this species [14,15]. A study analyzing Croatian populations revealed a genetic split between coastal and continental populations caused by the Dinaric Alps [16]. In this study we sampled 38 locations from sixteen European countries, and additionally two Georgian populations to reveal insight into the phylogeography and demography of this species in Europe.

Materials and Methods
Samples (larvae) were collected from 36 localities in sixteen European countries and two additional sites from Georgia between 2010 and 2014 (Table 1), with specimens from the same locality/site constituting a geographic population. Larvae were taken from distant trees (not more than one larva per tree) so to avoid any bias induced by the use of mtDNA marker, and after collection they were immediately put into 99% ethanol and stored at −20 • C in laboratory. In total, 497 individuals were analyzed. One pair of abdominal legs was used for genomic DNA extraction from each specimen individually, using SIGMA Aldrich mammalian genomic DNA extraction kit (Sigma, St. Louis, MO, USA) and following manufacturer's instruction. DNA amplification was performed with primers UEA5 and UEA10 [17] in 20 µL reactions containing 2 mM of MgCl2, 100 mM of dNTPs (Fermentas, Vilnius, Lithuania), 0.5 mM of each primer 1U of Taq (Fermentas, Vilnius Lithuania). Thermocycling consisted of an initial denaturation step at 94 • C (2 min), which was followed by 33 cycles at 94 • C (30 s), 48 • C (60 s), and 68 • C (60 s), and a final extension step at 68 • C (10 min). All reactions were checked for amplification by gel electrophoresis, and confirmed products were purified using peqGOLD Cycle-Pure Kit (PeqLab, Erlangen, Germany) and directly sequenced with UEA10 at Eurofins MWG Operon (Mainz, Germany). Sequences were manually examined using Chromas Lite (Technelysium, Brisbane Australia) and aligned using ClustalW algorithm incorporated in BioEdit [18]. Upon alignment, ends of sequences were trimmed in order to obtain full overlap. Finally, this alignment was used for determination of haplotype sequences, using TCS 1.21 [19]. Singleton sequences were verified by an additional PCR and sequencing. Haplotype sequences were deposited in NCBI with accession numbers: KY628707 to KY628749.

Country
Abbrev. Host Lat. Long. N HT Hd ± SD π ± SD MNPD ± SD r (5) To separate the source of variation in three hierarchical levels (among groups of populations (F CT ), among populations within groups (F ST ), and within populations (F SC )), objective grouping and Analysis of Molecular Variance were performed using ARLEQUIN [19] based on the host of each population ( Table 1). The best-fit nucleotide substitution model of our data was selected among 88 different models (in this case GTR with Gamma distributed heterogeneity of rates) as implemented in jModeltest 2.1.7 [25], and this model was then used in MrBayes 3.2 [26]. Two separate chains were run in the Bayesian analysis, and the initial run lasted 2 × 10 6 generations, with a sampling frequency of 100 generations in order to define the trees to be discarded in the second run. As a plateau between standard deviations was reached after 5 × 10 5 generations, 500 trees were discarded from the estimation of the consensus tree. To better understand the assignment of haplotypes in clades and their geographic distribution, we also used sequences deposited in NCBI GenBank from previous investigations of L. dispar (HM013724 (Russia, Vostochnyy), HM013736 (Japan), HM013737 (Japan) [15]; KX436205 (Japan, Honshu), KX436223 (Japan, Nagano), KX436230 (Russia, Vladivostok) [27]; KY923059 (Russia, Primorski), KY923063 (Lithuania), KY923064 (Russia, Krasnoyarsk), KY923065 (Kazhakstan), GU994784 (Mongolia) [28]).

Results
In total, 497 sequences 676 bp long (spanning from 699 to 1375 bp of L. dispar COI sequence) were used for further analyses. Cladistic analysis resulted in 60 haplotypes ( Figure 1) distinguished by 25 singletons and 19 parsimony informative polymorphic sites, with the rarefaction curve of haplotypes reaching a plateau (data not shown) showing that the sampling effort was sufficient to describe intraspecific haplotype diversity [32]. Total haplotype diversity H T was 0.81 and intrapopulation diversity H S was 0.437 (data not shown). The lowest genetic diversity values were spotted in populations from the continental basin of the Balkan Peninsula between Dinaric Alps, Alps and Carpathian Alps and in northern Continental regions. On the other hand, the highest genetic diversity was observed in Coastal Croatia, Italy, and France. Populations from the Aegean coast and Black sea coast showed rather low diversity, while populations from Georgia showed intermediate values (Table 1).
The phylogenetic tree grouped European haplotypes in a similar pattern supporting the main clades with posterior probabilities higher than 80%, with the exception of green cluster that was the most common one in Central Europe. The inclusion of NCBI GenBank sequences from Asia, revealed the position of haplotypes LD52, LD53, and LD54 from Georgia diverged 1.46% from the European groups, while averagely diverged 0.86% from the cluster of Asian haplogroups L. dispar asiatica/L. japonica haplotypes (KX436223, KY923509, HM013724, GU994784, HM013737, HM013736, KX436205, KX436230) retrieved from NCBI ( Figure 3). Finally, while Siberian sequences integrated in the European L. dispar dispar cluster, Georgian haplotypes formed a distinct clade, separated from L. dispar japonica and L. dispar asiatica.   Table 1) and dimension of pies are proportional to sample size and haplotype abundance. Underlaying map represents biogeographic regions of Europe [33]. The figure was constructed using ArcMap 10.1 (ESRI, Redlands, CA, USA).   Table 1) and dimension of pies are proportional to sample size and haplotype abundance. Underlaying map represents biogeographic regions of Europe [33]. The figure was constructed using ArcMap 10.1 (ESRI, Redlands, CA, USA).

Figure 2.
Haplotype distribution of Lymantria dispar in Europe. Colored circles represent the clades according to parsimony network (Figure 1), acronyms indicate localities ( Table 1) and dimension of pies are proportional to sample size and haplotype abundance. Underlaying map represents biogeographic regions of Europe [33].

Discussion
Gypsy moth is a holarctic species, flexible in terms of temperature tolerance [34] and it has extremely wide spectrum of host plants [35]. This species successfully invaded a vast variety of ecosystems, causing disturbances on a wide array of tree and general plant species. Model predictions for L. dispar dispar support a range shift of up to 900 km northwards in Europe as a consequence of increase in average air temperature [36]. Therefore, it is important to understand the phylogeography of L. dispar dispar in Europe regarding management strategies [37].
Major [38][39][40] and minor refugia [27,41] during ice ages coupled with recent gene flow events [42] determined the current pattern of divergence. Though for some species it has been demonstrated that the most recent ice age had the most profound impact on intraspecific diversity [43], for several other species the origin of divergence goes further back before the last ice age [44,45], with Lymantria dispar being one of them [13]. Intraspecific divergence among haplotypes of European L. dispar dispar (0.5029%) is congruent with previous findings that this clade arose around 200 thousand years (kyr) [13]. In particular, the observed genetic pattern of L. dispar dispar in Europe seems to have been shaped by at least two refugia, one located near the Carpathian Mountains (green) spreading to Central

Discussion
Gypsy moth is a holarctic species, flexible in terms of temperature tolerance [34] and it has extremely wide spectrum of host plants [35]. This species successfully invaded a vast variety of ecosystems, causing disturbances on a wide array of tree and general plant species. Model predictions for L. dispar dispar support a range shift of up to 900 km northwards in Europe as a consequence of increase in average air temperature [36]. Therefore, it is important to understand the phylogeography of L. dispar dispar in Europe regarding management strategies [37].
Major [38][39][40] and minor refugia [27,41] during ice ages coupled with recent gene flow events [42] determined the current pattern of divergence. Though for some species it has been demonstrated that the most recent ice age had the most profound impact on intraspecific diversity [43], for several other species the origin of divergence goes further back before the last ice age [44,45], with Lymantria dispar being one of them [13]. Intraspecific divergence among haplotypes of European L. dispar dispar (0.5029%) is congruent with previous findings that this clade arose around 200 thousand years (kyr) [13].
In particular, the observed genetic pattern of L. dispar dispar in Europe seems to have been shaped by at least two refugia, one located near the Carpathian Mountains (green) spreading to Central Europe and the Pontic-Mediterranean region (blue) spreading to the Balkans and the north of Europe. Furthermore, it is likely that contemporary pattern of gypsy moth phylogeography in Europe was also influenced by anthropogenic migration. This has been already hypothesized in the recent study of the Asian gypsy moth based on COI barcode sequences [46] that revealed an unusually close relationship between geographically distant haplotypes which could not be explained merely by natural dispersal of gypsy moth.
Being a highly polyphagous species, L. dispar populations would have been expected to exhibit deep intraspecific divergence as different selection pressures are associated with different host plants, something that favors the emergence of barriers to gene flow [47,48]. Should these barriers be maintained for a sufficiently long period of time, they may lead to the formation of distinct host-associated lineages [49]. This concept was initially thought to be valid mostly for parthenogenetically reproducing species [50][51][52], yet recent studies have shown that it is also evident in sexually reproducing species [49,53]. However, the limited genetic structure revealed among gypsy moth populations in Europe coupled with the AMOVA outcome do not support this expectation. As with other polyphagous insect species such as Pityogenes chalcographus [54] or Heliconius melpomene [55,56] life-history traits might have evolved on one host plant upon secondary contact could have been transmitted to populations feeding on a different host species, and thus render the herbivore insect species capable to feed on different host species [54]. Furthermore, this effect is reinforced by the high gene flow of L. dispar as the combination of the migratory pathways (flying moths, wind-borne dispersal of neonate larvae and long distances anthropogenic movement [9]) could facilitate the redistribution and coalescence of populations [57].

Conclusions
Despite the relatively large sample, we were able to determine only limited inferences of phylogeographic structure among localities in Europe, with at least two refugial areas contributing to the current genetic structure of gypsy moth. The near panmictic status of haplotype LD4 also present in Georgia suggests that more complex forces have taken role in shaping of the contemporary phylogeographic pattern of gypsy moth, and one of likely explanations could be a man-aided dispersal. Such hypotheses have already been drawn on larger scale studies [46]. However, the separation of the Georgian populations from other European populations, but also from the Asian haplotypes L. dispar asiatica and L. dispar japonica, shows that there are likely many more genotypic entities than recently described [28]. In fact, as differences in biology and behavior are described between gypsy moths from Europe and Asia, i.e., flight capability of females [10,58], the region of the Caucasus should be studied more thoroughly in order to examine the taxonomic status of the genetic entities in that geographic area. Funding: This research was partly funded by Croatian Science Foundation project HRZZ 7616 DIFPEST, and Bilateral scientific and technical cooperation project funded by Croatian Ministry of Science and ÖAD. This study was also supported by the Austrian Science Fund FWF project I2604-B25.

Conflicts of Interest:
The authors declare no conflict of interest.