Molecular Ecology (2015) 24, 1275–1291

doi: 10.1111/mec.13103

Genetic structure, admixture and invasion success in a Holarctic defoliator, the gypsy moth (Lymantria dispar, Lepidoptera: Erebidae) YUNKE WU,*† JOHN J. MOLONGOSKI,* DEBORAH F. WINOGRAD,* STEVEN M. B O G D A N O W I C Z , † A R T E M I S S . L O U Y A K I S , ‡ D A V I D R . L A N C E , * V I C T O R C . M A S T R O * and RICHARD G. HARRISON† *Otis CPHST Lab, Joint Base Cape Cod, United States Department of Agriculture, 1398 West Truck Road, Buzzards Bay, MA 02542, USA, †Department of Ecology and Evolutionary Biology, Cornell University, Corson Hall, Ithaca, NY 14853, USA, ‡Department of Microbiology and Cell Science, Space Life Sciences Laboratory, Kennedy Space Center, University of Florida, Gainesville, FL 32611, USA

Abstract Characterizing the current population structure of potentially invasive species provides a critical context for identifying source populations and for understanding why invasions are successful. Non-native populations inevitably lose genetic diversity during initial colonization events, but subsequent admixture among independently introduced lineages may increase both genetic variation and adaptive potential. Here we characterize the population structure of the gypsy moth (Lymantria dispar Linnaeus), one of the world’s most destructive forest pests. Native to Eurasia and recently introduced to North America, the current distribution of gypsy moth includes forests throughout the temperate region of the northern hemisphere. Analyses of microsatellite loci and mitochondrial DNA sequences for 1738 individuals identified four genetic clusters within L. dispar. Three of these clusters correspond to the three named subspecies; North American populations represent a distinct fourth cluster, presumably a consequence of the population bottleneck and allele frequency change that accompanied introduction. We find no evidence that admixture has been an important catalyst of the successful invasion and range expansion in North America. However, we do find evidence of ongoing hybridization between subspecies and increased genetic variation in gypsy moth populations from Eastern Asia, populations that now pose a threat of further human-mediated introductions. Finally, we show that current patterns of variation can be explained in terms of climate and habitat changes during the Pleistocene, a time when temperate forests expanded and contracted. Deeply diverged matrilines in Europe imply that gypsy moths have been there for a long time and are not recent arrivals from Asia. Keywords: global phylogeography, hybridization, mitochondrial DNA, Pleistocene glaciation

invasion

success,

microsatellite,

Received 23 September 2014; revision received 29 January 2015; accepted 2 February 2015

Introduction Patterns of genetic variation within and among conspecific populations are a product of both history and current ecology. Analyses of such patterns allow inferences Correspondence: Richard G. Harrison, Fax: (607) 255 8088; E-mail: [email protected] © 2015 John Wiley & Sons Ltd

about gene flow between populations and how gene flow has changed over time. With sufficiently intensive sampling, it becomes possible to identify regions of recent or ongoing admixture between differentiated populations. Amounts and patterns of variation likewise provide information about recent range expansions, including human-mediated introductions.

1276 Y . W U E T A L . For invasive species, defining current population structure provides the context for identifying source populations and tracing the probable paths and genetic consequences of invasion (Bock et al. 2015). Historical records can often tell us how, when and from where invasive organisms arrived in new areas, but for many organisms such records are not complete or do not exist. Characterizing allele and haplotype frequencies for molecular markers provides direct estimates of both the amount and spatial pattern of genetic variation in native and introduced populations. These data can be used for testing hypotheses about the loss of genetic variation during invasions, the role of multiple introductions in augmenting diversity, and the importance of admixture in catalysing colonization success (Dlugosch & Parker 2008; Rius & Darling 2014). Recent studies have suggested that high levels of genetic diversity may contribute to invasion success and that successful invasions often result from admixture in colonizing populations (Kolbe et al. 2008; Crawford & Whitney 2010; Krehenwinkel & Tautz 2013; Bock et al. 2015). However, the importance of admixture has also been questioned, given that bottlenecked, lowdiversity populations are also successful colonists. Thus, a review of molecular diversity in native and introduced populations of 80 species concluded that introduced populations often showed substantial losses in both allelic richness and heterozygosity compared to their probable source populations (Dlugosch & Parker 2008). Admixture of different source populations has often been observed in natural range expansions (Vila et al. 2005; Sakaguchi et al. 2011; Krehenwinkel & Tautz 2013); indeed, the range expansion of modern humans has been accompanied by admixture with both Neanderthals and Denisovans (Pr€ ufer et al. 2013), and some evidence suggests that alleles derived from interbreeding may have allowed human populations to adapt to more northern climates (Sankararaman et al. 2014). Rius & Darling (2014) have reviewed the evidence that admixture has a causal role in successful colonization by non-native species, concluding that it remains unclear the extent to which admixture is a driver rather than simply a correlate of that success. For example, they suggest that high propagule pressure will enhance the likelihood of invasion success and also increase the probability of admixture. A polyphagous pest of deciduous forests and trees in urban areas, the gypsy moth (Lymantria dispar) is one of the most destructive introduced insects in North America. Indigenous populations of L. dispar are found throughout the Japanese archipelago, across much of Asia, throughout Europe and in parts of North Africa. In 1868 or 1869, a handful of adult moths accidentally escaped from a house in Medford, Massachusetts,

where an amateur entomologist, Leopold Trouvelot, was rearing gypsy moths (Forbush & Fernald 1896). Although it has been suggested that Trouvelot was attempting to hybridize gypsy moths with native silk moths (to produce a more ‘robust’ silk moth), there is, in fact, no evidence to support this claim, nor any reason to believe that crosses between such distant relatives would actually be successful. Despite early attempts at eradication, gypsy moth populations became established in eastern Massachusetts and gradually expanded their range, which now extends from the Maritime provinces in Canada south to North Carolina and west to Wisconsin and Minnesota. Additional introductions may have occurred early in the 20th century (Headlee 1921), but the historical evidence is not compelling. In contrast, it is clear that gypsy moths from East Asia have arrived (but not established) in North America in recent decades; the sources of most of these moths have been ships from the Russian Far East (Animal and Plant Health Inspection Service [APHIS] 2014). The population structure of gypsy moth is of particular interest because of its wide geographic distribution, potential for invasion and great economic importance. Patterns of morphological variation in this insect were described in great detail by the evolutionary biologist Richard Goldschmidt. Goldschmidt (1934) suggested that L. dispar had its origins in Japan and that the moth then spread from Japan and/or Korea across Eurasia. Currently, three subspecies of L. dispar are recognized, primarily based on morphology: L. d. asiatica from Continental Asia, L. d. dispar from Europe and North America and L. d. japonica from the Japanese archipelago (Pogue & Schaefer 2007). This classification broadly corresponds to geographic region of origin, except for the grouping of European and North American populations, which is in accord with historical accounts of the introduction to North America. Asian gypsy moths differ in a number of ways from their European and North American relatives. First, unlike most European female moths, Asian females are capable of powered flight (but see Reineke & Zebitz 1998; Keena et al. 2008). This ability may accelerate the spread of populations; in the absence of female flight, dispersal occurs primarily through the ballooning of first instar larvae. Asian moths also have different host plant preferences (Baranchikov 1988), which may allow them to utilize resources largely unaffected by European moths. Previous studies that have characterized the population structure of gypsy moths have used small sample sizes and/or few molecular markers (Bogdanowicz et al. 2000; Keena et al. 2008; deWaard et al. 2010). Here, we investigate the global pattern of genetic variation in L. dispar, using allele frequencies for nine highly polymorphic microsatellite loci and ~2900 bp of mitochondrial DNA © 2015 John Wiley & Sons Ltd

G L O B A L P H Y L O G E O G R A P H Y O F G Y P S Y M O T H 1277 (mtDNA) sequences. We have access to a large collection of specimens obtained over 20 years (1992–present); these samples represent the entire Holarctic distribution of gypsy moth from Asia, Europe and North America. Our aim was to use this robust data set to characterize worldwide population structure, to examine evidence for admixture in natural populations of this species and to assess the genetic consequences of introduction(s) of this moth to North America. We also evaluate the extent to which Pleistocene glacial/interglacial cycles influenced lineage diversification in gypsy moths.

Materials and methods Sample collection and DNA extraction We included representatives of the three currently recognized subspecies, Lymantria dispar asiatica, L. d. dispar and L. d. japonica (Fig. 1; Table 1). Specimens collected from west of the Ural Mountains were considered to be L. d. dispar; specimens collected from east of the Ural Mountains were identified as L. d. asiatica (Pogue & Schaefer 2007). Sampling effort was most intense in Asia, which is believed to be the centre of origin for the species and where moths exhibit the most genetic variation and the most variation in morphological and life history traits (Goldschmidt 1934; Harrison et al. 1983; Bogdanowicz et al. 2000). The number of moths collected from each location depended on local population density and ranged from 3 to 117 (average of 35.5 moths per location). All specimens were stored dry at 20 or 70 °C. Genomic DNA was extracted using a Proteinase K protocol (Maniatis et al. 1982). A leg or antenna of each moth was placed in 500 lL extraction buffer (5% Proteinase K, 0.1% Tergitol, 19 TE buffer) and incubated overnight. The extraction was deactivated by heating at 75 °C for 30 min and then stored at 20 °C.

Microsatellite loci development and amplification Methods for developing microsatellite loci are described in Bogdanowicz et al. (1997). In addition to the three loci used by Bogdanowicz et al. (1997), we identified six additional loci from the L. dispar genome that were polymorphic among populations. Locus name, motifs and flanking primers for these six loci are described in Table 2. Briefly, genomic DNA was digested with the restriction enzyme Sau3AI. Fragments of 400–600 bp in size were isolated by agarose gel electrophoresis and electroelution. Genomic fragments were ligated to the plasmid pUC 18 that was previously digested with BamHI and dephosphorylated. Products of ligation were used to chemically transform DH5-a cells. Aliqu© 2015 John Wiley & Sons Ltd

ots of the transformation were grown on Luria–Bertani (LB) plates with ampicillin. Colonies were transferred to nylon membranes and probed with radiolabelled (32P) repeat oligos, representing 10 trimers and two dimers. Hybridizations occurred in TMAC buffer to negate variation in melting temperatures among probes. Plasmid DNA from positive colonies was isolated by a rapid lysis-by-boiling procedure and sequenced with M13 primers that flanked the inserts. Amplification of microsatellite loci was conducted in polymerase chain reactions (PCR, 10 lL total volume) containing 3.5 lL molecular grade water, 1 lL 109 PCR buffer without MgCl2, 1.4 lL MgCl2 (25 mM), 1.6 lL dNTP solution (1.25 mM), 0.2 lL primers (10 pmol/lL) and 0.1 lL of JumpStart Taq DNA Polymerase (2.5 units/lL). Two microlitres DNA was used as template. Forward primers were labelled with fluorescent WellRed dyes (Sigma Aldrich, St. Louis, MO, USA). For microsatellite loci 10F1, 49, 106 and 202, PCR cycling conditions consisted of an initial 1 min denaturation at 95 °C followed by 40 cycles of 95 °C for 15 s, 55 °C for 15 s and 72 °C for 30 s, and a final extension of 72 °C for 10 min. Cycling conditions for loci 198, 9C2, 238, 254 and 101 were the same as above except that the annealing temperature was 60 °C. After amplification, fragments were separated and analysed on a Beckman Coulter CEQ 8000 DNA Analyzer. Allele sizes were manually checked and scored based on trace data. Presence of null alleles was investigated in MICROCHECKER 2.2.3 (van Oosterhout et al. 2004). Because missing data can affect the estimate of null allele percentage, for this analysis, we included only the 17 locations that had 25 or more specimens with complete data (specimens with missing data were discarded).

Summary statistics of population differentiation Genetic diversity of L. dispar was assessed using the nine microsatellite loci. We used CONVERTER 1.31 (Glaubitz 2004) to format the microsatellite data and count the number of private alleles over all loci. Linkage disequilibrium and deviation from Hardy–Weinberg (HW) equilibrium were calculated with GENEPOP 3.4 (Raymond & Rousset 1995). Linkage disequilibrium was tested between each pair of loci for 17 locations that had 25 or more specimens with complete data (the same data set used for MICROCHECKER) with 1000 dememorizations in 1 million Markov chain approximations. Sequential Bonferroni correction was used to account for multiple pairwise comparisons. Deviation from HW equilibrium was assessed for each locus at each of the 49 locations using an exact test. We used FSTAT 2.9.3.2 (Goudet 1995) to estimate mean allele richness (RS) for each location after correcting for sample

1278 Y . W U E T A L .

16 6 7 5

14

2 1

3

17

20

8

10

35

19 28 26 25

15

27

13 4

18

29

11

31

23

21

30 22

24

9 12

38

40

44

41 42

37

49 45

43 48

36

39 34

32

33

46

47 Fig. 1 Sampling of Lymantria dispar throughout the Holarctic region (upper figure). Red dots denote sampling locations in North America, Asia and Europe (middle and lower figures). Locations names are given in Table 1. Two locations (33 and 34) that lack detailed locality information are placed in respective country where large forests exist.

size differences. Samples with fewer than seven moths or having missing data for an entire locus were omitted from the RS estimation. Mean expected heterozygosity (HE) and observed heterozygosity (HO) and genetic diversity (D) were estimated in ARLEQUIN 3.5 (Excoffier & Lischer 2010). We compared RS and D values (see Table 1) among groups using one-way analysis of variance (ANOVA) in SPSS 19. If the null hypothesis

of equal group variances was rejected, a nonparametric Welch’s test was used. Distribution of genetic variance between and within genetic groups was assessed through a hierarchical analysis of molecular variance (AMOVA) using ARLEQUIN 3.5. Genetic groups were initially classified according to subspecies designations and then modified based on cluster analysis (see below). Statistical significance of © 2015 John Wiley & Sons Ltd

G L O B A L P H Y L O G E O G R A P H Y O F G Y P S Y M O T H 1279 Table 1 Sampling locations and microsatellite diversity indices No.

Sampling localities

North America 1 Connecticut, USA† 2 New Jersey, USA† 3 West Virginia, USA 4 North Carolina, USA 5 Michigan, USA 6 Minnesota, USA† 7 Wisconsin, USA† Asia 8 Enshi, Hubei, China 9 Beibei, Chongqing, China 10 Liu’an, Anhui, China† 11 Funing, Jiangsu, China 12 Yichun, Jiangxi, China 13 Tengzhou, Shandong, China† 14 Beijing, China† 15 Changli, Hebei, China 16 A’ershan, Nei Mongol, China 17 Tianjin, China 18 Suihua, Heilongjiang, China 19 Liaoyuan, Jilin, China 20 Jianping, Liaoning, China† 21 Honshu, Japan† 22 Hiroshima, Japan 23 Hokkaido, Japan† 24 Chiba, Japan 25 Donghae-si, Gangwon-do, South Korea† 26 Joan-myeon, Gyeonggi-do, South Korea 27 Incheon, South Korea 28 Seoul, South Korea 29 Busan, South Korea 30 Pohang-si, Gyeongsangbuk-do, South Korea 31 Ulsan, South Korea 32 Almaty, Kazakhstan 33 Unknown, Kazakhstan 34 Unknown, Kyrgyzstan 35 Vladivostok, Primorsky Krai, Russia† 36 Novosibirsk, Russia 37 Tomsk, Russia 38 Yekaterinburg, Sverdlovsk, Russia Europe 39 Krasnodar, Krasnodar Krai, Russia 40 Indre-Et-Loire, Verneuil-Sur-Indre, France 41 Near Germany, France 42 Petit Landau, Haut Rhin, France 43 Pino, Corsica, France† 44 Langen, Hesse, Germany† 45 Sopron, Gy} or-Moson-Sopron, Hungary 46 Sardinia, Italy† 47 Alcala de los Gazules, Cadiz, Spain† 48 El Casta~ no, Andalucıa, Spain†  49 Banska Stiavnica, Bansk a Bystrica, Slovakia†

n

AN

AP

Rs

D

HE

HO

91 75 60 38 17 79 69

3.33 3.44 1.67 2.67 2.78 3.00 1.11

0 0 0 0 1 1 0

2.53 2.81 n.a. 2.27 2.50 1.89 n.a.

0.457 0.446 0.384 0.462 0.416 0.253 0.354

0.486 0.526 0.239 0.475 0.472 0.274 0.213

0.452* 0.495* 0.189* 0.384* 0.446 n.s. 0.210* 0.193*

4 3 26 5 6 40 47 26 15 73 55 19 105 80 51 30 31 70 24 34 91 9 10 4 17 11 28 57 9 10 10

3.56 2.78 6.67 3.89 4.78 8.00 8.44 7.56 6.67 9.67 10.33 7.11 11.44 10.56 9.44 7.00 7.11 11.22 9.11 10.67 12.00 6.00 5.33 3.78 5.11 3.56 4.67 11.33 3.89 4.00 3.78

0 0 0 0 0 1 1 0 0 1 2 0 2 0 1 0 1 0 0 3 3 0 0 0 1 0 0 3 0 0 0

n.a. n.a. 4.53 n.a. n.a. 4.74 4.54 4.71 4.43 4.42 5.08 5.02 4.49 4.76 4.65 4.65 4.41 5.41 5.63 5.55 5.43 5.18 4.60 n.a. 3.92 3.19 3.05 5.28 3.42 3.56 3.17

0.655 0.496 0.609 0.610 0.544 0.592 0.625 0.589 0.614 0.619 0.668 0.686 0.584 0.621 0.554 0.555 0.516 0.707 0.693 0.741 0.680 0.677 0.589 0.615 0.611 0.503 0.467 0.556 0.560 0.514 0.474

0.696 0.730 0.719 0.690 0.770 0.720 0.693 0.728 0.648 0.671 0.744 0.758 0.675 0.678 0.670 0.715 0.653 0.764 0.762 0.779 0.760 0.773 0.663 0.701 0.649 0.568 0.494 0.753 0.560 0.599 0.520

0.611 n.s. 0.389 n.s. 0.524* 0.498 n.s. 0.607 n.s. 0.570* 0.532* 0.560* 0.445* 0.519* 0.574* 0.605* 0.518* 0.511* 0.501* 0.507* 0.518* 0.576* 0.670 n.s. 0.632* 0.586* 0.642 n.s. 0.467* 0.528 n.s. 0.527* 0.483* 0.425* 0.627* 0.469 n.s. 0.567 n.s. 0.404 n.s.

10 48 16 8 8 117 29 42 11 7 13

3.22 5.89 4.11 2.56 1.78 7.78 5.44 5.44 3.00 2.44 4.78

0 3 0 1 0 1 0 0 0 0 0

2.80 3.67 3.09 2.39 1.77 3.82 3.76 3.19 2.44 2.38 3.78

0.385 0.549 0.395 0.387 0.309 0.549 0.429 0.412 0.320 0.365 0.611

0.420 0.579 0.457 0.387 0.309 0.604 0.616 0.507 0.334 0.386 0.632

0.347 n.s. 0.431* 0.352* 0.389 n.s. 0.347 n.s. 0.480* 0.492* 0.431* 0.252 n.s. 0.370 n.s. 0.462*

Number of individual sampled (n); mean number of alleles per locus per locality (AN); number of private alleles over all the loci (AP); mean allelic richness corrected for sample size of 7 individuals (RS), localities with fewer than 7 individuals were omitted; Average gene diversity over loci (D); expected heterozygosity (HE) under Hardy–Weinberg (HW) equilibrium; observed heterozygosity (HO) and the statistical significance for departure from HW equilibrium (*significant after sequential Bonferroni correction for P < 0.05; n.s., nonsignificant). Localities also sampled for mitochondrial DNA sequences are denoted by †. © 2015 John Wiley & Sons Ltd

1280 Y . W U E T A L . Table 2 Microsatellite and mitochondrial DNA primers used in this study Locus name

Motif

Microsatellite 10F1 (AC) 49 (TGA) 198 (AC) 9C2 (TGTC) 238 (CA) 254 (CA) 101 (TGA) 106 (TGA) 202 (AC) Mitochondrial DNA ND2 ATP6/ATP8 ND6/cytb

Forward primer

Reverse primer

CGC ACA AAG CTC TCA GAT GA See Bogdanowicz et al. (1997) See Bogdanowicz et al. (1997) GTT TAT GGT TTG TAA TTA TTT AAC ATT CA ACT GTT CGT TTA TTC AAT AGT GTT GG TAC TGT TTG AAG TCG GTT TTG C AAT TTA CCC TTC CGT TAT GTA GAC AGG CTC GAT GCC AGT AGT GG See Bogdanowicz et al. (1997)

CGT TAC CGC GTG TCT ACA TT

TGG ATG TTG AAT TGG GTT AGA GCG GAA CTA ACC ACA GAT TT AAC CCC CTT TCT ATA GGA TT

covariance associated with each hierarchical level was calculated with 1000 permutations.

Bayesian cluster analysis Population structure of L. dispar was inferred from microsatellite data using a Bayesian clustering method implemented in STRUCTURE 2.3.4 (Pritchard et al. 2000). We used the admixture model to allow individuals to have mixed ancestry from multiple clusters in L. dispar. Sampling locations were used as prior information to assist clustering. Allele frequencies were assumed to be correlated among populations because we expected some interpopulation gene flow. A series of Markov chain Monte Carlo (MCMC) simulations were run from K = 2 to K = 20 with 20 replicates for each value of K. The chain was run for 500 000 generations after a burnin of 100 000 generations. The most likely value of K was calculated following Evanno et al. (2005) in STRUCTURE HARVESTER (Earl & vonHoldt 2012). After choosing the optimal value of K, the associated 20 MCMC replicates were aligned under the greedy algorithm of CLUMP 1.1.2 (Jakobsson & Rosenberg 2007). A mean matrix was computed as the mean of 20 individual matrices after columns were aligned according to the permutation with the greatest H-value (average pairwise similarity). This mean matrix was plotted in DISTRUCT 1.1 (Rosenberg 2004). We followed Weisrock et al. (2010) in discarding individual membership coefficients with posterior probabilities <0.05 and assigned corresponding values proportionally to other membership coefficients so that coefficients for each individual still summed to one.

Microsatellite network We constructed a network to visualize relationships among populations. Microsatellite data were converted

AAA ATC AGT TCA GTG GTT TAG ACG ATA TCC CTT AGT CGC CTT TTA CG GAT GAC TAG GGT ATT CAA TAC GCA ACA TAT TCG AAC AGT TGT TTC ATA A ACA AAG CCA ATC GGA TAG AAC A

ATT GCA AAT TTT AAG GAG TAT TT TGA TTG GAT AAC CGC AAC TG TGA TCC AGT TTG ATG AAG AA

from number of nucleotides to number of repeats in 1.31. Nei’s (1978) unbiased standard genetic distance was calculated across loci among sample locations in MICROSAT 1.5b (Minch et al. 1996) with the stepwise-mutation model. We performed 1000 bootstrap replicates by resampling loci with replacement. The mean distance matrix was analysed with the NeighborNet algorithm in SPLITSTREE4 4.13.1 (Huson & Bryant 2006), because Neighbor-Net tends to construct networks that are better resolved than those from other algorithms (Bryant & Moulton 2004). The network was displayed using the EQUALANGLE algorithm. CONVERTER

Isolation by distance For populations of L. dispar that occupy separate continents or islands, water barriers limit gene flow. But for populations that are continuously distributed across the Eurasian continent, it is unclear where and how differentiation occurs. We tested for isolation by distance (IBD) among continental Eurasian populations using microsatellite data. We used a Mantel test to determine whether there is a positive association between genetic distance and geographic distance. Isolation by distance predicts that populations that are geographically close are more genetically similar than those that are at greater geographic distances (Wright 1943). Geographic distances were calculated by GEOGRAPHIC DISTANCE MATRIX GENERATOR 1.2.3 (Ersts 2014). Pairwise FST values (Weir & Cockerham 1984) calculated by GENEPOP were used for the genetic distance. The Mantel test was performed with 1000 random permutations to obtain statistical inferences at a = 5% using the program ARLEQUIN. Because tests of IBD can be confounded by the presence of population structure (Meirmans 2012), we performed an additional partial Mantel test, which incorporates the Bayesian cluster assignment as a © 2015 John Wiley & Sons Ltd

G L O B A L P H Y L O G E O G R A P H Y O F G Y P S Y M O T H 1281 covariate to the genetic and geographic distance matrices. Statistical inferences at a = 5% were achieved through 1000 random permutations.

MtDNA genealogy and divergence dates To evaluate the role of Pleistocene climate cycles in shaping the genetic structure of L. dispar, we used mitochondrial DNA to date intraspecific divergence events. We amplified and sequenced five mitochondrial protein-coding genes (ND2, ND6, ATP6, ATP8 and cytb) to reconstruct the matrilineal genealogy. These five genes exhibit greater variability within Lepidoptera than the commonly used COI gene and thus can improve phylogenetic resolution (Cameron & Whiting 2008). DNA sequence data were obtained from a subset of samples that included 68 moths from 18 locations (Table 1). Lymantria albescens was chosen as the out-group based on its close relationship with L. dispar (deWaard et al. 2010). We designed PCR primers (Table 2) to amply three mitochondrial fragments that covered the five genes. Cycling conditions consisted of initial denaturation at 94 °C for 3 min, followed by 30 cycles of denaturation at 94 °C for 30 s, annealing at 55 °C for 30 s, extension at 60 °C for 3 min and a final extension at 60 °C for 10 min. Negative controls were performed to monitor contamination. Amplified PCR products were examined on 3% agarose gels and subsequently purified by ExoSAP-IT (Affymetrix, Cleveland, OH, USA) and then sequenced on an ABI 3730XL. Sequences were aligned by the MUSCLE algorithm in GENEIOUS 5.6.6 with default parameters. The alignment was manually adjusted when insertions–deletions (indels) were not aligned properly. Because of the absence of recombination in most mitochondrial genomes, we concatenated the five mitochondrial genes. We calculated pairwise uncorrected p-distance among sample locations and the mean p-distance among genetic groups in MEGA 5.20 (Tamura et al. 2011). Portions of alignment where indels are present were eliminated from the calculation. To estimate divergence times within L. dispar, we first tested the adherence of the phylogeny to a molecular clock using Tajima’s one-degree-of-freedom method (Tajima 1993) implemented in MEGA. Lymantria albescens was chosen as the out-group for the taxon triplet. The test approximately followed the chi-square distribution with one degree of freedom (d.f.). If the clock hypothesis was not rejected, we used the widely cited lepidopteran mtDNA substitution rate of 1.1–1.2 9 10 8 per site per year (Brower 1994; Maroja et al. 2007) to estimate dates of cladogenetic events under a strict clock model in BEAST 1.7.5 (Drummond & Rambaut 2007). Data were bi-partitioned into nucleotides and indels, which were analysed as standard binary data. We did © 2015 John Wiley & Sons Ltd

not further partition nucleotides by genes or codon positions due to low variability of mtDNA within L. dispar. A substitution model for nucleotides was selected based on the Akaike Information Criterion obtained from JMODELTEST 2.1.4 (Darriba et al. 2012). A simple model without site heterogeneity was selected for the indel partition. A Bayesian MCMC chain was run for 50 million generations sampled every 1000 generations. Chain convergence was assessed by effective sample size (ESS) values of the parameters in TRACER 1.5 (Rambaut & Drummond 2007) and posterior distributions in AWTY (Wilgenbusch et al. 2004). The first 25% of trees were discarded as burn-in.

Results Microsatellite data analysis A total of 1738 moths were genotyped for nine microsatellite loci (Table 1). Several loci showed evidence of null alleles in the MICROCHECKER analysis, but estimated percentages of null alleles were generally low (<5%; Appendix S1, Supporting information). Deviations from HW equilibrium (HO < HE) were observed at seven loci, but usually only at a few locations (Table 1; Appendix S2, Supporting information). The exception was locus 9C2, for which a deficiency of heterozygotes was observed at most sites. Comparison of deviations from HW equilibrium with results from the MICROCHECKER analysis suggested that most deviations are due to the presence of null alleles. Given the low percentage of potential null alleles, we retained the full data set for subsequent analyses. In tests of linkage equilibrium, only 15 of 612 pairwise comparisons showed evidence of disequilibrium after a sequential Bonferroni correction (Appendix S3, Supporting information). For the microsatellite loci, Lymantria dispar was characterized by high numbers of alleles ranging from 11 (locus 49) to 30 (loci 10F1 and 238), with a mean of 19.8 alleles per locus (Table 2). Private alleles were recovered in all three continents and subspecies. Although gypsy moths arrived in North America only very recently, populations in Minnesota and Michigan harbour two alleles at locus 254 that were not recovered at any other locations. The average number of alleles per locus per location was overall high in populations from Asia, intermediate in European populations and relatively low in North American populations. After correction for differences in sample sizes, we still observed a decline of allelic diversity from East Asia to Europe, with a further decrease in North America. Across the entire global distribution, moths from northeastern China, the Russian Far East and the Korean Peninsula have the highest allelic richness (all show RS > 5;

21

H

Fig. 2 Estimated population structure of Lymantria dispar based on STRUCTURE analysis of microsatellite data. Four genetic clusters are recognized. A total of 1738 moths are represented by thin vertical lines with their membership coefficients in four clusters colour coded. Thick black vertical lines delimit geographic areas; from left to right: North America, Europe, Middle Asia, Mainland East Asia and Japanese archipelago. Location numbers and subspecies assignments are shown below.

on

s

Ja

hu ,

pa

n

S So H uth ok K ka or e i d o, a Ja pa

28

23 29 30 31

japonica

22

24

1282 Y . W U E T A L .

n

ea

K

ut h

N

So

25

26

27

or

ia

E

Ru

a

in

Ch

N

Ch

in

a

40

44

47

48

46

13 34 33 32 10 45 49 8 1112 9 38 37 36 39

14

15 16

17

asiatica

NE

20

18

19

35

ss

Sp a in Ita ly H un Sl ga o v ry C aki & a W K yr Rus g K yzs sia az ta ak n hs Ch tan in a C

G er m an Co

nn

ec tic ut ,U

2

se y,

Ne w

3 4 5

a,

ni

US

M in n e so ta ,U S W isc o n sin M ,U ic S h No igan ,U rt h Ca S ro W lin a, tV

7

S

US

Je r

irg i

6

dispar

43 1

y

Fr an ce

41

42

Table 1). Allelic richness (RS) and genetic diversity (D) values were compared among Asia, Europe and North America by ANOVA (null hypothesis of equal group variances was not rejected, both P > 0.05). For this analysis, three locations (36–38) from Central Russia were included with European populations, given results from Bayesian cluster analysis (see below). Both RS and D revealed that Asian populations had significantly higher values than European and North American populations (both P < 0.0005), which were not significantly different from each other (RS: P = 0.124; D: P = 0.483). However, the latter two groups show distinct allele frequencies at multiple loci. For example, at locus 198, the 151 bp allele was the most common allele in North American populations (frequency >60% in all populations except for one at 31%), while this allele was relatively rare (<9%) in European populations (Table 2; Appendix S4, Supporting information). When we define populations according to their assignment to one of the three morphological subspecies, an AMOVA showed that most of the genetic variance (80.73%, P < 0.001) is found within populations, with 6.60% (P < 0.001) of the variance among locations within subspecies and 12.67% (P < 0.001) of the total variance among subspecies. If locations are grouped into the four well-defined Bayesian clusters (see below), AMOVA gave similar results, with 14.80% (P < 0.001) of the variance among clusters and 81.40% (P < 0.001) of the variance within locations.

Bayesian cluster analysis Following Evanno et al.’s (2005) method, STRUCTURE analysis identified four well-differentiated genetic clusters (Fig. 2; Appendix S5, Supporting information). When the locus most affected by null alleles was discarded and the analysis rerun, the result was the same (Appendix S6, Supporting information). Three of the four clusters defined by STRUCTURE broadly correspond to the three morphological subspecies, with the fourth cluster representing North American populations. However, boundaries between putative subspecies are not discrete. Moths from northeastern China, the Russian Far East and the Korean Peninsula have been identified © 2015 John Wiley & Sons Ltd

es

US

G L O B A L P H Y L O G E O G R A P H Y O F G Y P S Y M O T H 1283 as L. d. asiatica, but based on microsatellite data, most of these moths are of mixed ancestry, with approximately equal individual membership coefficients from L. d. asiatica and L. d. japonica. In fact, the STRUCTURE analysis suggests a pattern of clinal variation along a latitudinal gradient in this region. Alleles of L. d. asiatica predominate in northeastern China, whereas alleles of L. d. japonica predominate in populations in the southern part of the Korean Peninsula (locations 29–31). These southern Korean populations were assigned to L. d. japonica for subsequent network analyses. Pure L. d. japonica is found on the island of Honshu of Japan, but moths from the northern island of Hokkaido have a small fraction of their genome from L. d. asiatica and thus resemble those from the southern part of the Korean Peninsula in genomic make-up. Pure L. d. asiatica was found across China. In Central Asia, populations also appear to exhibit mixed ancestry. Moths from Kazakhstan (locations 32 and 33) and nearby Kyrgyzstan (location 34) both have a large portion of their genomes from L. d. asiatica, but the former population has alleles from L. d. japonica and the latter has alleles from L. d. dispar. Cluster analysis also suggested that Central Russian populations collected from east of the Ural Mountains (locations 36–38) do not group with L. d. asiatica but instead are predominantly L. d. dispar, the subspecies that occurs throughout Europe. North American populations, although recently derived from western European moths, were assigned to a genetic cluster distinct from European populations (Fig. 2). European and North America moths remain distinct even when we choose different values of K (K = 2, K = 3, and K = 5; Appendix S7 Supporting information). Among North American populations, moths from two locations [Connecticut (1) and New Jersey (2)], which are geographically closest to the original introduction site in Massachusetts, displayed the highest proportion of ‘European ancestry’. Some moths from Michigan showed low levels (around 10%) of admixture from L. d. japonica.

Microsatellite network Based on Nei’s (1978) distance estimated from microsatellite data, SPLITSTREE produced a tree (Fig. 3) that was congruent with the genetic clusters recovered by STRUCTURE. The deepest divisions are between populations from Asia and those from Europe/North America. Samples from Central Asia are placed near the middle of the two major splits. Consistent with the cluster analysis, the network grouped the three locations from the southern end of the Korean Peninsula with Japanese populations. Clusters representing L. d. asiatica and © 2015 John Wiley & Sons Ltd

L. d. japonica were connected by samples from the eastern Asian hybrid zone. The three Central Russian locations collected from east of the Ural Mountains were indeed nested within European populations. Moths from Spain were characterized by large genetic distances from other parts of Europe. STRUCTURE revealed European populations and North American populations to be two distinct clusters, but the SPLITSTREE network suggests that North American samples derive from French populations, although the branches joining them are long. Moths from New Jersey and Connecticut are more similar to moths from French populations than moths from other North American locations, consistent with results from the STRUCTURE analysis and the documented introduction of L. d. dispar from Western Europe into Massachusetts in 1868 or 1869 (Forbush & Fernald 1896).

Isolation by distance Pairwise FST values estimated from microsatellite data (Appendix S8, Supporting information) were high between locations from different continents (mean between continent FST = 0.27) and substantially lower between populations within continent (Asia, FST = 0.09; Europe, FST = 0.15; North America, FST = 0.09). We detected a strong signal of IBD (correlation coefficient r = 0.6438, P < 0.001) when the Eurasian populations (excluding Japan) were analysed (Appendix S9, Supporting information). Even after the effect of population structure (i.e. differentiation between L. d. dispar and L. d. asiatica) was controlled by a partial Mantel test, a significant relationship between FST and geographic distance was still found among Eurasian populations (r = 0.3570, P < 0.001).

MtDNA genealogy and divergence dates For the five mitochondrial genes, we generated a total of 2901–2932 bp of DNA sequence for sampled individuals. One 6-bp indel was found in the coding sequence of ATP8. Another 5-bp indel occurs in a noncoding region between ND6 and cytb. In the same region, we discovered a microsatellite locus with dinucleotide repeat motif (TA). Allelic variation among individuals ranged from 9 to 22 repeats. This microsatellite locus also was present in the out-group species Lymantria albescens. Excluding indels (but including the microsatellite locus), 85 sites were variable within L. dispar, of which 53 were parsimony-informative. Intraspecific uncorrected p-distance averaged 0.47% across 18 sampling locations (Appendix S10, Supporting information). Average sequence divergences (mean p-distances) between populations of L. d. asiatica and L. d. japonica

1284 Y . W U E T A L .

7

Fig. 3 Neighbor-Net network based on Nei’s (1978) unbiased standard genetic distance calculated from microsatellite data. Nodes show 47 sampling locations, which are grouped into clusters defined by STRUCTURE analysis (colour code same as Fig. 2). The network clearly shows that North American populations are derived from European populations.

0.05

6

3 5

4 1

2 39 42

48 41 40 45

37 44 36 38

47

43 46 49

34 11

32

18 12 27 35 19 28 25 15

33

29 24 21 22

26

16 17 14 20

13

10

23

9 30 8 31

and between European and North American populations of L. d. dispar were about 0.2%. Mean p-distances were about 0.8% when comparing populations of either Asian subspecies to those of L. d. dispar. Phylogenetic analysis of mtDNA identified TrN93+ I as the optimal substitution model for the concatenated sequence data after indels were excluded. We coded indels as presence–absence binaries and analysed them with the sequence data in BEAST. Including or excluding indels in the analysis resulted in identical topologies (Fig. 4); nodal support values were slightly higher when indels were included. We recovered a well-supported dichotomy between the two Asian subspecies

and L. d. dispar. Within the Asian clade, the two subspecies are reciprocally monophyletic, but only L. d. japonica received 100% bootstrap support. Within the L. d. dispar clade, some locations have deeply divergent matrilines, and individuals from the same location did not always form monophyletic groups (e.g. Slovakia and Italy). Consistent with the microsatellite network, moths from North American locations were nested within clades from Western Europe. Despite the recent and presumably small introduction, the 16 North American moths sampled included four mtDNA haplotypes (mtDNA microsatellite region excluded). © 2015 John Wiley & Sons Ltd

G L O B A L P H Y L O G E O G R A P H Y O F G Y P S Y M O T H 1285 MIS 11 (MBE)

MIS 9

MIS 7

MIS 5

North America nested in this clade

dispar 198–215 kya

asiatica

Lymantria dispar 354–386 kya

113–123 kya

LR04 δ18Obenthic(‰)

japonica

Connecticut_10_46_03 NewJersey_07_2_01 Minnesota_12_77_01 Wisconsin_12_1_01 Wisconsin_07_54_02 Wisconsin_12_1_02 Minnesota_10_4_03 Minnesota_12_77_02 Wisconsin_07_54_01 Minnesota_10_4_01 Germany_1054_02 Germany_1054_03 Germany_1054_04 Connecticut10_45_11 Connecticut10_46_06 Italy_81_17 France_73_10 France_73_12 France_73_09 France_73_11 Italy_81_18 Germany_1054_01 NewJersey07_1_04 Connecticut10_49_01 NewJersey07_1_01 NewJersey07_1_02 Italy_81_19 Italy_81_16 Spain_1_02 Spain_1_03 Spain_1_01 Spain_2_01 Slovakia_13_03 Slovakia_13_02 Slovakia_13_01 Slovakia_13_04 China_74_02 Russia_339_06 China_20_02 Korea_67_19 Russia_339_01 Russia_339_02 China_74_03 Korea_67_20 Russia_339_05 China_20_14 China_20_16 Korea_67_26 Korea_67_27 China_22_10 China_20_01 China_29_12 China_22_12 China_22_11 China_74_04 China_29_11 China_29_13 China_29_14 China_22_13 China_74_01 Japan_315_01 Japan_321_02 Japan_321_01 Japan_344_11 Japan_321_04 Japan_344_12 Japan_344_13 Japan_344_14 Lymantria albescens

2.5

2.5

3.0

3.0

3.5

3.5

4.0

4.0

4.5

4.5 400

300

200

100

0 kya

Fig. 4 MtDNA genealogy and divergence dates of Lymantria dispar derived from five mitochondrial genes and assuming substitution rates of 1.1–1.2 9 10 8 per site per year. Lymantria albescens was used as the out-group. The three subspecies of L. dispar each represent a clade. Solid dots denote posterior probability greater than 99%. Shaded areas show the upper and lower estimates of divergence dates. In the background, dark grey columns indicate glacial periods and light grey columns indicate interglacial periods (Lisiecki & Raymo 2005). Bottom of the figure: d18O benthic record from the LR04 stack, modified from Tzedakis et al. (2009). MIS, marine oxygen-isotope stage; MBE, Mid-Brunhes Event.

Tajima’s one-degree-of-freedom method did not reject the hypothesis that mtDNA of L. dispar evolved under a molecular clock (v2 = 0.73, d.f. = 1, P = 0.39). Applying © 2015 John Wiley & Sons Ltd

the substitution rate of 1.1–1.2 9 10 8 per site per year, the basal dichotomy of mtDNA haplotypes in L. dispar was dated to 354–386 thousand years ago (Fig. 4), at a

1286 Y . W U E T A L . time of transition from a prolonged interglacial period [marine oxygen-isotope stage (MIS) 11] to the next glacial period, according to the LR04 Benthic Stack by Lisiecki & Raymo (2005). This divergence date closely followed the so-called Mid-Brunhes Event (c. 410– 430 thousand years ago), which marks the onset of a greater amplitude of glacial–interglacial variability after the Middle Pleistocene (Jansen et al. 1986; Tzedakis et al. 2009). Initial diversification within the European/North American clade (interglacial MIS 7) pre-dated that between the two Asian subspecies, which occurred during interglacial MIS 5 (Fig. 4).

Discussion Detailed analyses of patterns of genetic variation in Lymantria dispar confirm the presence of distinct subspecies but also reveal extensive natural hybridization and introgression at the boundaries between these forms. Environmental changes during the Pleistocene led to episodes of vicariance and secondary contact, resulting in genetically distinct geographic populations. Subsequent dispersal and gene flow have blurred the boundaries where secondary contact has occurred. Nonetheless, observed patterns of genetic variation allow identification of source populations for introductions. Gypsy moths introduced to North America are clearly derived from Western Europe; these populations exhibit reduced allelic diversity, and admixture between different subspecies has not played a role in the success of this invasive species.

Genetic structure of L. dispar populations and genetic consequences of the introduction into North America Recent morphological studies have recognized three subspecies of L. dispar, which represent populations from three major geographic areas where the moth is endemic (Continental Asia, Japan, Europe) and one area where the moth is introduced (North America; Pogue & Schaefer 2007). Two former subspecies of L. dispar, each restricted to a single island in the Japanese archipelago, have been determined to be distinct species (Lymantria umbrosa and Lymantria albescens; Pogue & Schaefer 2007), and these populations indeed exhibit substantial mtDNA sequence divergence from each other and from L. dispar (Bogdanowicz et al. 2000). Our molecular data, both mtDNA sequences and microsatellite allele frequencies, are generally consistent with the morphological analysis, but reveal that boundaries between subspecies are not discrete and that admixed populations are common in eastern and central Asia. Although AMOVA reveals that the majority of microsatellite variation is found within populations, both FST based on

microsatellite allele frequencies and mtDNA sequence divergence provide evidence of substantial genetic differences among subspecies/geographic regions. Of the four genetic clusters defined in STRUCTURE, three roughly correspond to L. d. asiatica, L. d. dispar from Europe and L. d. japonica. The fourth cluster includes all sample locations from North America. Although these moths were introduced to North America from Western Europe <150 years ago (Forbush & Fernald 1896), the founder event that occurred during the initial introduction appears to have substantially altered the allelic composition of North American populations, such that they do not cluster with European populations in STRUCTURE analyses regardless of the value of K. Measures of allelic richness and heterozygosity are lower in North America than in France or Germany, but the differences are not statistically significant (Table 1). Although allele frequencies in North America appear altered by the bottleneck at introduction, both the microsatellite network and the mtDNA gene tree clearly reveal the ancestry of North American gypsy moths to be in Western Europe. In both phylogenetic trees (Figs 3 and 4), North American populations are nested within a European cluster. It has long been thought that Trouvelot’s French origins make France the probable source of the North American introduction, but the exact source in Western Europe remains unknown (Forbush & Fernald 1896; Harrison et al. 1983). The microsatellite network is consistent with French ancestry, because two locations in Eastern France (41 and 42) are most similar to North American populations based on Nei’s distance. A Western Russian location (39) also shows a small Nei’s distance from North American populations, but it is unlikely to be the source population given what we know about the history of the release (Forbush & Fernald 1896). The mtDNA gene tree groups North American moths within the European clade, but again the pattern of ancestry is not entirely clear. Haplotypes from Connecticut and New Jersey are sister to haplotypes from Germany, but other haplotypes from these North American sites, as well as haplotypes from Wisconsin and Minnesota form a separate clade (although still within the larger European clade). In addition to the major allele frequency shifts that occurred during the founder event, we also find moths from Michigan and Minnesota that have alleles not found in other parts of the species’ distribution, despite extensive sampling. These may represent rare alleles in Europe or new mutations that have arisen after the introduction. Our data suggest that introduction and establishment of gypsy moths subsequent to the documented 1869 event have either not occurred or have had little impact on the genetic composition of gypsy moth populations in North America. Subsequent © 2015 John Wiley & Sons Ltd

G L O B A L P H Y L O G E O G R A P H Y O F G Y P S Y M O T H 1287 introductions from Europe might not be easily recognized, given the relative homogeneity of microsatellite allele frequencies across Western Europe. Regardless of whether secondary introductions have occurred, there is no evidence that admixture between genetically differentiated populations has played a role in the success of gypsy moths in North America (Rius & Darling 2014). An ongoing discussion about the genetics of introduced and invasive species focuses on the role of admixture, which occurs when multiple introductions of the ‘same’ species result in hybridization between two or more previously independent lineages. Admixture results in increases in standing variation (enabling adaptation) and can give rise to transgressive segregation (producing phenotypes outside of the ranges of both parents). It is evident that admixture is not required for successful invasion, but not yet clear whether admixture enhances the probability of success. Often ignored in the invasion literature is recognition that admixture is common in natural populations and that hybrid zones have been carefully studied for more than a century. Most hybrid zones are the result of secondary contact between previously isolated populations and therefore are useful models for admixture of introduced populations. Although hybridization has long been recognized by botanists as a major source of new variants (frequently, but not always involving polyploidy), the role of hybridization in producing novel and successful lineages in animals remains controversial (Barton 2001; Schumer et al. 2014). Our concern with admixture as a catalyst for the origin of new variants and successful invasion must be calibrated based on our understanding of admixture in natural populations. If admixture does increase colonization success (Kolbe et al. 2008; Crawford & Whitney 2010; Krehenwinkel & Tautz 2013; Bock et al. 2015), then natural hybridization should lead to source populations that have greater probability of invading other regions. For this reason, detailed analyses of patterns of genetic variation in natural populations of potentially invasive species are an essential context for assessing the probable consequences of an introduction. In gypsy moths, hybridization between the subspecies L. d. asiatica and L. d. japonica appears to be extensive, and hybrid zone populations have the highest genetic diversity across the entire global distribution. The proportion of ancestry assigned to the two subspecies varies in a northsouth latitudinal cline, with the proportion of the genome attributed to L. d. asiatica increasing with latitude. At the southern end of the Korean Peninsula, individuals have more than 70% of their genome descended from L. d. japonica and thus are, on average, genetically closer to L. d. japonica than to L. d. asiatica. © 2015 John Wiley & Sons Ltd

Regardless of subspecific assignment, East Asian populations are both highly variable and have been the probable source of moths found on ships entering North American ports and of moths trapped in the vicinity of these ports (APHIS 2014). In the face of these threats, it would seem that the protective measures implemented by North American countries to prevent new human-mediated introductions of L. dispar (APHIS 2014) have generally been successful. However, future introductions from nonEuropean source populations might substantially change the genetic composition of North American populations, given the high levels of genetic variation in Asian L. dispar and their substantial differences (both genotypic and phenotypic) from European moths. In such a scenario, North American gypsy moth populations might gain greater adaptive potential, through enhanced dispersal ability and broader host plant choice (Kolbe et al. 2008; Rius & Darling 2014).

Dispersal and gene flow In Asia, both gypsy moth males and females are capable of long-distance flight (Rozkhov & Vasilyeva 1982), whereas in most of Europe only males can fly. In all settings, there is wind-aided larval dispersal, potentially over long distances (Baranchikov 1988). Although potential dispersal distances are large, most caterpillars carried by wind do not move more than 100 m (Minott 1922; Mason & McManus 1981), and the movements of male and female moths usually do not exceed 3–5 km (Baranchikov 1988). Thus, isolation by distance across suitable habitat may be evident. Indeed, we see a strong signal of isolation by distance across the Eurasian continent. In contrast, the Ural Mountains, recognized as a phylogeographic break between sister groups for a number of taxa from Europe and Asia (Hewitt 1996; Brunhoff et al. 2003) and suggested as the barrier that separates L. d. asiatica and L. d. dispar (Baranchikov 1988; Pogue & Schaefer 2007), do not appear to be a major extrinsic barrier to gene flow for gypsy moths. Individuals collected from three Russian locations east of the mountain range are assigned to the European cluster in the STRUCTURE analysis and are grouped with European populations in the network. Here ecological factors [abundant availability of a desirable food source—oaks (Quercus spp.)] may facilitate movement of L. d. dispar across the mountains. Narrow ocean barriers can also substantially limit gene flow, which will resume if such barriers disappear. That changes in ocean barriers have influenced gypsy moth population structure is evident from the distribution of haplotypes and alleles in Asia. The Korea Strait (including the Tsushima Strait), which isolates Korea

1288 Y . W U E T A L . and Japan, was formed 1.5–1.7 million years ago (Kitamura et al. 2001), prior to the split date between L. d. asiatica and L. d. japonica inferred from mtDNA. If the two Asian subspecies were formed by allopatric separation, then the common ancestor of the two subspecies apparently crossed the strait prior to subspecies divergence. The split between L. d. asiatica and L. d. japonica is dated to early MIS 5, an interglacial period with high sea levels due to the melting of the Greenland and Antarctica ice sheets (Tzedakis et al. 2009). This divergence time is consistent with allopatric separation, because high sea levels submerged any possible land bridge formed in previous glacial periods and reduced faunal connections between Japan and the Korean Peninsula. However, alleles from L. d. japonica found at the southern end of the Korean Peninsula most likely trace back to moths that dispersed across the narrow Korea Strait from Japan. The Korea Strait is less than 130 m deep, and maximum sea level lowering in the Middle and Late Pleistocene was 130–140 m. Periods of low sea level would have either completely closed the strait or narrowed it significantly (MillienParra & Jaeger 1999), allowing for increased gene exchange. Studies of Japanese sawfly species reveal gene flow with Korean populations around the same time (Sakurai et al. 2009).

Diversification and Pleistocene glaciation cycles Patterns of ongoing gene flow and recent introductions certainly contribute to current population structure in L. dispar, but historical events and environmental changes have also had a major influence. MtDNA dating under a strict clock model suggests that L. dispar split into European and Asian lineages in the Middle Pleistocene. MIS 11 has been identified as one of the key interglacial periods over the past 800 000 years (EPICA Community Member 2004). It is an exceptionally warm and long interglacial. More importantly, MIS 11 is characterized by the Mid-Brunhes Event (c. 410–430 thousand years ago), after which the amplitude of glacial–interglacial variability increased significantly (Jansen et al. 1986; Tzedakis et al. 2009). Greater variability could drive larger fluctuations in forest extent, with forest expansion during interglacial periods and contraction during glacial periods. The population size and geographic range of L. dispar probably fluctuated with the extent of deciduous forest, because the most preferred host—Quercus spp.—is a major component of temperate forests in both Europe and Asia (Jin et al. 1987; Tong et al. 1988, 1992; Okuda et al. 2001; Tzedakis et al. 2006; Koutsodendris et al. 2010). Re treat to isolated refugia (e.g. in southern Europe; Hewitt 1996, 2000, 2001) would lead to divergence and could explain patterns that we see.

Evidence for multiple glacial refugia in Europe is found in our mtDNA data. The most deeply diverged (basal) matrilines are found in Slovakia, Spain and Italy, perhaps reflecting relatively early (100–200 thousand years ago) population subdivision in southern refugia. Because most females of the European gypsy moth (L. d. dispar) do not fly, and given that mtDNA is maternally inherited, we expect persistence of ancestral lineages in former glacial refugia. Southern Europe (e.g. Italy, Spain) served as refuges during periods of glacial advance (Hewitt 1996, 2000, 2001), and our findings are consistent with the presence of those isolated refugia. Furthermore, the deep divergence of mtDNA haplotypes suggests that L. d. dispar has had a long-term presence in Europe, contrary to the views of Goldschmidt (1934), who believed that gypsy moth had only recently colonized Europe. Both early molecular studies of gypsy moth populations (Harrison et al. 1983; Harrison & ODell 1989) and our current estimates of allelic richness reveal a pattern in which haplotype diversity or allelic richness is highest in East Asia, gradually decreases across Central Asia and Europe and is lowest in North America due to the recent introduction. Such a pattern has been explained by invoking an East Asian origin for gypsy moth and a recent arrival in Europe (Goldschmidt 1934; Harrison & ODell 1989; Bogdanowicz et al. 2000). However, given that observed mitochondrial diversification within the European clade pre-dates that within the Asian clade and that the relative age of each clade is not affected by our choice of model for nucleotide substitution rates, patterns of mtDNA variation suggest that reduced variation at microsatellite loci in Europe requires an alternative explanation. Reduced nuclear diversity may be a consequence of repeated bottlenecks rather than recent arrival.

Conclusions Gypsy moths in North America are clearly descended from source populations in Western Europe, and there is no evidence that enhanced variability due to admixture between subspecies was important in the success of the invasion. Indeed, North American populations exhibit reduced variability and differences in allele frequency compared with populations in Europe, consistent with a population bottleneck. However, if admixture does enhance the probability of invasion success, then East Asian populations constitute a significant threat. Diversification of Lymantria dispar lineages appears to have been driven by events during the Middle and Late Pleistocene, with diversification occurring independently and contemporaneously in Europe and Asia. Contrary to the assertions of Goldschmidt, © 2015 John Wiley & Sons Ltd

G L O B A L P H Y L O G E O G R A P H Y O F G Y P S Y M O T H 1289 European populations of gypsy moths do not represent recent arrivals. Furthermore, the three recognized subspecies, which may well have diverged in allopatry, have not remained distinct, with extensive hybridization between L. d. asiatica and L. d. japonica occurring in eastern Asia, and some gene exchange between L. d. asiatica and L. d. dispar in central Asia. Despite hybridization at the boundaries between subspecies, the subspecies have distinctive mtDNA haplotypes and microsatellite allele frequencies. Using molecular markers that allow identification of region of origin for moths arriving in North America will provide critical information for decisions about future management of this destructive pest.

Acknowledgements We thank Baode Wang, Heung-Sik Lee, Alamaz Orozumbekov, Yasutomo Higashiura, Yuri Baranchikov, Tamara Freyman, Vasily Ponomarev and Paul Schaefer for assistance in obtaining gypsy moth specimens from throughout Asia and Russia. We acknowledge Melody Keena, Milan Zubrik, Frank Kruger, Matthias Kolb, Horst Marohn and Mike South for their similar assistance in procuring specimens from Europe. We also thank innumerable additional foreign and domestic cooperators who have provided gypsy moths from throughout the world to the Otis CPHST laboratory over the past twenty years for use in scientific research. We thank Peggy Elder for technical assistance and Rodger A. Gwiazdowski for assistance with data analysis and discussion on the manuscript. Two anonymous reviewers and the Associate Editor (Rosemary Gillespie) provided comments that much improved the manuscript.

References Animal and Plant Health Inspection Service [APHIS] (2014) Asian gypsy moth. Plant protection and quarantine fact sheet. http://www.aphis.usda.gov/publications/plant_health/ content/printable_version/fs_phasiangm.pdf Baranchikov YN (1988) Ecological basis of the evolution of host relationships in Eurasian gypsy moth populations. In: Proceedings of the Lymantriidae: A Comparison of Features of New and Old World Tussock Moths (eds Wallner WE, McManus KA), pp. 319–338. USDA Forest Service Northeastern Forest Experiment Station, Broomall, PA. Barton NH (2001) The role of hybridization in evolution. Molecular Ecology, 10, 551–568. Bock DG, Caseys C, Cousens RD et al. (2015) What we still don’t know about invasion genetics. Molecular Ecology. doi:10.1111/mec.13032. Bogdanowicz SM, Mastro VC, Prasher DC, Harrison RG (1997) Microsatellite DNA variation among Asian and North American gypsy moths (Lepidoptera: Lymantriidae). Annals of the Entomological Society of America, 90, 768–775. Bogdanowicz SM, Schaefer PW, Harrison RG (2000) Mitochondrial DNA variation among worldwide populations of gypsy moths, Lymantria dispar. Molecular Phylogenetics and Evolution, 15, 487–495. © 2015 John Wiley & Sons Ltd

Brower AVZ (1994) Rapid morphological radiation and convergence among races of the butterfly Heliconius erato inferred from patterns of mitochondrial DNA evolution. Proceedings of the National Academy of Sciences, USA, 91, 6491–6495. Brunhoff C, Albreath KEG, Fedorov VB, Ook JAC, Jaarola A (2003) Holarctic phylogeography of the root vole (Microtus oeconomus): implications for late Quaternary biogeography of high latitudes. Molecular Ecology, 12, 957–968. Bryant D, Moulton V (2004) Neighbor-Net: an agglomerative method for the construction of phylogenetic networks. Molecular Biology and Evolution, 21, 255–265. Cameron SL, Whiting MF (2008) The complete mitochondrial genome of the tobacco hornworm, Manduca sexta (Insecta: Lepidoptera: Sphingidae), and an examination of mitochondrial gene variability within butterflies and moths. Gene, 408, 112–123. Crawford KM, Whitney KD (2010) Population genetic diversity influences colonization success. Molecular Ecology, 19, 1253– 1263. Darriba D, Taboada GL, Doallo R, Posada D (2012) jModelTest 2: more models, new heuristics and parallel computing. Nature Methods, 9, 772. Dlugosch KM, Parker IM (2008) Founding events in species invasions: genetic variation, adaptive evolution, and the role of multiple introductions. Molecular Ecology, 17, 431–449. Drummond AJ, Rambaut A (2007) BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evolutionary Biology, 7, 214. Earl D, vonHoldt BM (2012) STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conservation Genetics Resources, 4, 359–361. EPICA Community Member (2004) Eight glacial cycles from an Antarctic ice core. Nature, 429, 623–628. Ersts PJ (2014) Geographic Distance Matrix Generator (version 1.2.3). American Museum of Natural History, Center for Biodiversity and Conservation. Available from http://biodiversityinformatics.amnh.org/open_source/gdmg. Accessed August 19, 2014. Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software Structure: a simulation study. Molecular Ecology, 14, 2611–2620. Excoffier L, Lischer HE (2010) Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Molecular Ecology Resources, 10, 564–567. Forbush EH, Fernald CH (1896) The Gypsy Moth and the Brown Tail Moth. Wright and Potter Press, Boston, MA. Glaubitz JC (2004) CONVERT: A user-friendly program to reformat diploid genotypic data for commonly used population genetic software packages. Molecular Ecology Notes, 4, 309–310. Goldschmidt R (1934) Lymantria. Bibliographia. Genetica, 11, 1– 185. Goudet J (1995) FSTAT (Version 1.2): a computer program to calculate F-statistics. Journal of Heredity, 86, 485–486. Harrison RG, ODell TM (1989) Mitochondrial DNA as a tracer of gypsy moth origins. In: Proceedings of the Lymantriidae: A Comparison of Features of New and Old World Tussock Moths (eds Wallner WE, McManus KA), pp. 265–273. USDA Forest

1290 Y . W U E T A L . Service Northeastern Forest Experiment Station, Broomall, PA. Harrison RG, Wintermeyer SF, Odell TM (1983) Patterns of genetic variation within and among gypsy moth, Lymantria dispar (Lepidoptera: Lymantriidae), populations. Annals of the Entomological Society of America, 76, 652–656. Headlee TJ (1921) The present status of gypsy moth in New Jersey. Journal of Economic Entomology, 14, 172–177. Hewitt GM (1996) Some genetic consequences of ice ages, and their role in divergence and speciation. Biological Journal of the Linnean Society, 58, 247–276. Hewitt GM (2000) The genetic legacy of the Quaternary ice ages. Nature, 405, 907–913. Hewitt GM (2001) Speciation, hybrid zones and phylogeography: or seeing genes in space and time. Molecular Ecology, 10, 537–549. Huson DH, Bryant D (2006) Application of phylogenetic networks in evolutionary studies. Molecular Biology and Evolution, 23, 254–267. Jakobsson M, Rosenberg NA (2007) CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics, 23, 1801–1806. Jansen JHF, Kuijpers A, Troelstra SR (1986) A mid-Brunhes climatic event: long-term changes in global atmospheric and ocean circulation. Science, 232, 619–622. Jin Q, Wang P, Wang S (1987) The sporo-pollen assemblages and paleoclimate in the middle course area of the Huaihe River during Late Cenozoic. Marine Geology and Quaternary Geology, 7, 93–109. Keena MA, C^ ote MJ, Grinberg PS, Wallner WE (2008) World distribution of female flight and genetic variation in Lymantria dispar (Lepidoptera: Lymantriidae). Environmental Entomology, 37, 636–649. Kitamura A, Takano O, Takata H, Omote H (2001) Late Pliocene-early Pleistocene paleoceanographic evolution of the Sea of Japan. Palaeogeography, Palaeoclimatology, Palaeoecology, 172, 81–98. Kolbe JJ, Larson A, Losos JB, de Queiroz K (2008) Admixture determines genetic diversity and population differentiation in the biological invasion of a lizard species. Biology Letters, 4, 434–437. Koutsodendris A, M€ uller UC, Pross J, Brauer A, Kotthoff U, Lotter AF (2010) Vegetation dynamics and climate variability during the Holsteinian interglacial based on a pollen record from Dethlingen (northern Germany). Quaternary Science Reviews, 29, 3298–3307. Krehenwinkel H, Tautz D (2013) Northern range expansion of European populations of the wasp spider Argiope bruennichi is associated with global warming–correlated genetic admixture and population-specific temperature adaptations. Molecular Ecology, 22, 2232–2248. Lisiecki LE, Raymo ME (2005) A Pliocene-Pleistocene stack of 57 globally distributed benthic d18O records. Paleoceanography, 20, PA1003. Maniatis T, Fritsch EF, Sambrook J (1982) Molecular Cloning: A Laboratory Manual. Cold Spring Harbor Laboratory, Cold Spring Harbor, NY. Maroja LS, Bogdanowicz SM, Wallin KF, Raffa KF, Harrison RG (2007) Phylogeography of spruce beetles (Dendroctonus

rufipennis Kirby) (Curculionidae: Scolytinae) in North America. Molecular Ecology, 16, 2560–2573. Mason CJ, McManus ML (1981) The role of dispersal in the natural spread of the gypsy moth. In: Proceedings of the Second IUFRO Conference on Dispersal Of Forest Insects: Evaluation, Theory and Management Implications (eds Berryman AA, Safranyik L), pp. 94–115. Washington State University, Cooperative Extension, Pullman, WA. Meirmans PG (2012) The trouble with isolation by distance. Molecular Ecology, 21, 2839–2846. Millien-Parra V, Jaeger J (1999) Island biogeography of the Japanese terrestrial mammal assemblages: an example of a relict fauna. Journal of Biogeography, 26, 959–972. Minch E, Ruiz-Linares A, Goldstein D, Feldman M, CavalliSforza LL (1996) Microsat (version 1.5b): a computer program for calculating various statistics on microsatellite allele data. Available from http://genetics.stanford.edu/hpgl/projects/microsat/. Accessed August 19, 2014. Minott CW (1922) The Gypsy Moth on Cranberry Bogs. Kessinger Publishing, LLC, Whitefish, MT. Nei M (1978) Estimation of average heterozygosity and genetic distance from a small number of individuals. Genetics, 89, 583–590. Okuda M, Yasuda Y, Setoguchi T (2001) Middle to Late Pleistocene vegetation history and climatic changes at Lake Kopais, Southeast Greece. Boreas, 30, 73–82. van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P (2004) MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Molecular Ecology Notes, 4, 535–538. Pogue MG, Schaefer PW (2007) A Review of Selected Species of Lymantria H€ ubner [1819] (Lepidoptera: Noctuidae: Lymantriinae) from Subtropical and Temperate Regions of Asia Including the Description of Three New Species, Some Potentially Invasive to North America. Forest Health Technology Enterprise Team, Fort Collins, CO. Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics, 155, 945–959. Pr€ ufer K, Racimo F, Patterson N et al. (2013) The complete genome sequence of a Nenaderthal from the Altai Mountains. Nature, 505, 43–49. Rambaut A, Drummond AJ (2007) Tracer v1.4, Available from http://beast.bio.ed.ac.uk/Tracer Accessed August 19, 2014. Raymond M, Rousset F (1995) genepop (version 12): Population genetics software for exact tests and ecumenism. Journal of Heredity, 86, 248–249. Reineke A, Zebitz CPW (1998) Flight ability of gypsy moth females (Lymanfria dispar L.) (Lep., Lymantriidae): a behavioural feature characterizing moths from Asia? Journal of Applied Entomology, 122, 307–310. Rius M, Darling JA (2014) How important is intraspecific genetic admixture to the success of colonizing populations? Trends in Ecology and Evolution, 29, 233–242. Rosenberg NA (2004) Distruct: a program for the graphical display of population structure. Molecular Ecology Notes, 4, 137–138. Rozkhov AS, Vasilyeva TG (1982) Gypsy moth in central and eastern Siberia. Neparny Shelkoprayo v Srednai I Vostochnoi Sibiti, pp. 4–19. Russian Academy of Sciences, Siberian Branch, Novosibirsk, Russia.

© 2015 John Wiley & Sons Ltd

G L O B A L P H Y L O G E O G R A P H Y O F G Y P S Y M O T H 1291 Sakaguchi S, Takeuchi Y, Ymasaki M, Sakurai S, Isagi Y (2011) Lineage admixture during postglacial range expansion is responsible for the increased gene diversity of Kalopanax septemlobus in a recently colonised territory. Heredity, 107, 338–348. Sakurai A, Hamanishi Y, Maeto K, Naito C (2009) Molecular biogeography of two sibling species of the sawfly genus Macrophya (Hymenoptera: Tenthredinidae) in Japan. Zoological Science, 26, 325–329. Sankararaman S, Mallick S, Dannemann M et al. (2014) The genomic landscape of Neanderthal ancestry in present-day humans. Nature, 507, 354–357. Schumer M, Rosenthal GG, Andolfatto P (2014) How common is homoploid hybrid speciation? Evolution, 68, 1553–1560. Tajima F (1993) Simple methods for testing molecular clock hypothesis. Genetics, 135, 599–607. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S (2011) MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Molecular Biology and Evolution, 28, 2731–2739. Tong G, Yu S, Zhang J, Han S, Li Y (1988) On analysis of quantitative geology for Quaternary sporopollen in the Northern China Plain. Acta Bontanica Sinica, 30, 325–332. Tong G, Zhang J, Fan S (1992) Distribution of Quaternary palynoflora in China. Marine Geology and Quaternary Geology, 12, 45–56. Tzedakis PC, Hooghiemstra H, P€alike H (2006) The last 1.35 million years at Tenaghi Philippon: revised chronostratigraphy and longterm vegetation trends. Quaternary Science Reviews, 25, 3416–3430. Tzedakis PC, Raynaud D, McManus JF, Berger A, Brovkin V, Kiefer T (2009) Interglacial diversity. Nature Geoscience, 2, 751–755. orklund M (2005) The importance Vila M, Vidal-Romanı JR, Bj€ of time scale and multiple refugia: incipient speciation and admixture of lineages in the butterfly Erebia triaria (Nymphalidae). Molecular Phylogenetics and Evolution, 36, 249–260. deWaard JR, Mitchell A, Keena MA et al. (2010) Towards a global barcode library for Lymantria (Lepidoptera: Lymantriinae) tussock moths of biosecurity concern. PLoS One, 5, e14280. Weir BS, Cockerham CC (1984) Estimating F-statistics for the analysis of population structure. Evolution, 38, 1358–1370. Weisrock DW, Rasoloarison RM, Fiorentino I et al. (2010) Delimiting species without nuclear monophyly in Madagascar’s mouse lemurs. PLoS One, 5, e9883. Wilgenbusch JC, Warren DL, Swofford DL (2004) AWTY: A system for graphical exploration of MCMC convergence in Bayesian phylogenetic inference. http://ceb.csit.fsu.edu/ awty. Accessed August 19, 2014. Wright S (1943) Isolation by distance. Genetics, 28, 114–138.

Y.W., J.J.M., D.R.L., V.C.M. and R.G.H. planned the study. Y.W., J.J.M., D.F.W., S.M.B. and A.S.L. conducted laboratory work and data analyses. Y.W., J.J.M., S.M.B. and R.G.H. wrote the study. All authors contributed to revisions.

© 2015 John Wiley & Sons Ltd

Data accessibility DNA sequences are available in GenBank (KP705282– KP705488). The sequence alignment, mtDNA tree file, microsatellite genotype data, STRUCTURE input files and pairwise geographic distance matrix are available in DRYAD (doi:10.5061/dryad.v020j).

Supporting information Additional supporting information may be found in the online version of this article. Appendix S1 Estimated percentage of null alleles for microsatellite loci among 17 locations that had 25 or more specimens with complete data (specimens with missing data were discarded). Appendix S2 Statistical significance of deviation from Hardy– Weinberg (HW) equilibrium for microsatellite data calculated per locus per location. Appendix S3 Statistical test of linkage disequilibrium for microsatellite loci, using 17 locations that had 25 or more specimens with complete data (specimens with missing data were discarded). Appendix S4 Frequency of allele 151 bp at locus 198 between North American and European populations. Appendix S5 Evanno’s DK distribution (Evanno et al. 2005) calculated based on the complete microsatellite data set. Appendix S6 STRUCTURE analysis on microsatellite data with locus 9C2 excluded. Appendix S7 Genetic clustering of Lymantria dispar when different values of K are chosen in STRUCTURE analysis. Population names are given in Appendix S4. Appendix S8 Pairwise FST values of microsatellite data for 47 sampling locations. Appendix S9 Mantel test for correlation between FST and geographic distances for Eurasian populations (excluding Japan). The slope of the linear regression is (3.007e-05)  (9.447e-07), and the intercept is (3.625e-03)  (5.293e-03). Appendix S10 Uncorrected p-distance of mtDNA among three subspecies of Lymantria dispar.