Introduction
The European winter moth, Operophtera brumata L. (Lepidoptera: Geometridae), is a non-native invasive pest in the Northeastern USA that is currently spreading and causing defoliation and tree mortality in the Northeastern USA (Elkinton et al., Reference Elkinton, Boettner, Sremac, Gwiazdowski, Hunkins, Callahan, Scheufele, Donahue, Porter, Khrimian, Whited and Campbell2010; Simmons et al., Reference Simmons, Lee, Ducey and Dodds2014). In its native range of Europe, winter moth undergoes 9–10-year cyclical outbreaks (Tenow et al., Reference Tenow, Nilssen, Bylund, Pettersson, Battisti, Bohn, Caroulle, Ciornei, Cs Ka, Delb, De Prins, Glavendekic, Gninenko, Hrasovec, Matosevic, Meshkova, Moraal, Netoiu, Pajares, Rubtsov, Tomescu, Utkina and Gurney2013). North America has experienced three major winter moth epidemics. The first one occurred in Nova Scotia, Canada, where it was detected in the 1930s and caused widespread damage to oak forests and apple orchards in the 1950s (Cuming, Reference Cuming1961). Two parasitoids from Europe were released as biological controls: a tachinid fly, Cyzenis albicans (Fallén), and an ichneumonid wasp, Agrypon flaveolatum (Gravenhorst). By 1961, winter moth populations collapsed, as parasitoid populations increased. To this day, winter moth remains at low density in Nova Scotia (Roland & Embree, Reference Roland and Embree1995). The second winter moth outbreak occurred during the 1970s in the Pacific Northwest of Canada and the USA, where it became a serious pest of apples, hazelnuts, cherries, and blueberries (Kimberling et al., Reference Kimberling, Miller and Penrose1986). This outbreak subsided after release of the same parasitoids used in the Nova Scotia outbreak (Roland & Embree, Reference Roland and Embree1995).
The most recent winter moth outbreak was first detected in Massachusetts in the late 1990s. It has since spread to adjacent states and is currently causing persistent and widespread defoliation of forest trees and crops such as apples, blueberries, and cranberries (Elkinton et al., Reference Elkinton, Boettner, Liebhold and Gwiazdowski2015). The expansion of its range in the Northeastern USA was unexpected, given successful control in other regions. The leading edge of spread is moving approximately 6 km year−1, and defoliation caused by winter moth is expanding (Elkinton et al., Reference Elkinton, Liebhold, Boettner and Sremac2014). It is not known whether this outbreak originated from an existing North American population, or if it represents a new introduction from Europe.
Winter moth is known to hybridize with a closely related North American species, O. bruceata (Hulst), the Bruce spanworm (Elkinton et al., Reference Elkinton, Boettner, Sremac, Gwiazdowski, Hunkins, Callahan, Scheufele, Donahue, Porter, Khrimian, Whited and Campbell2010). Bruce spanworm is common across the Northern USA and Canada, but it rarely experiences outbreaks (Rose & Lindquist, Reference Rose and Lindquist1982). Winter moth and Bruce spanworm are difficult to distinguish based on wing patterns, but the males have distinct genital morphology (Eidt et al., Reference Eidt, Embree and Smith1966). Pivnick et al. (Reference Pivnick, Barton, Millar and Underhill1988) observed males with intermediate genitalia at several locations in British Columbia, as did Troubridge & Fitzpatrick (Reference Troubridge and Fitzpatrick1993), who speculated that 0.4% of the moths they examined from British Columbia with intermediate morphological characters may be hybrids.
Elkinton et al. (Reference Elkinton, Boettner, Sremac, Gwiazdowski, Hunkins, Callahan, Scheufele, Donahue, Porter, Khrimian, Whited and Campbell2010) confirmed the presence of field collected hybrids in the Northeastern USA using DNA sequence data from one mitochondrial gene (cytochrome oxidase subunit I (COI)) and one nuclear gene (glucose-6-posphate dehydrogenase (G6PD)). Hybrids were heterozygous for G6PD alleles fixed in each species, or had a COI haplotype characteristic of one species, and a G6PD allele characteristic of the other. Since the discovery of field-collected hybrids using this method, yearly sampling using pheromone traps that attract both species and their hybrids found 32 confirmed hybrids out of 667 samples along an East–West transect across central Massachusetts (4.8%; Elkinton et al., Reference Elkinton, Liebhold, Boettner and Sremac2014). However, since this assay relies on only one nuclear gene, it cannot be determined whether these hybrids represent generations beyond F1, and if they do, whether there is asymmetry in the direction of backcrosses.
If multi-generational transfer of genes and adaptive traits between the introduced and native species is occurring (reviewed in Harrison & Larson, Reference Harrison and Larson2014), this could impact both the invasiveness of winter moth and the pest status of Bruce spanworm. Determining whether winter moth/Bruce spanworm F1 hybrids can reproduce will therefore help us decide whether more detailed studies of hybridization are warranted in this system.
In this study, we develop single nucleotide polymorphism (SNP) loci using Restriction-Associated DNA Sequencing (RAD-seq; Baird et al., Reference Baird, Etter, Atwood, Currey, Shiver, Lewis, Selker, Cresko and Johnson2008) and microsatellite loci from shotgun genomic reads. Microsatellite markers are historically difficult to develop for Lepidoptera because microsatellites and their flanking regions tend to be duplicated in the genome (Zhang, Reference Zhang2004). We resolve this by filtering thousands of microsatellites from genomic reads to identify likely single copy markers. Microsatellite loci were then evaluated for utility in intra- and inter-specific population genetic studies of winter moth and Bruce spanworm and used to determine whether hybrid generations beyond F1 occur in the field using a group of moths collected in the Northeastern USA that had been previously classified as hybrids using G6PD (Elkinton et al., Reference Elkinton, Liebhold, Boettner and Sremac2014). We also performed laboratory crosses within and between the species to compare egg viability for the pure species versus their hybrids.
Materials and methods
RAD-seq SNP development
A RAD-seq library was prepared using DNA extracted from the head and thorax of 48 male moth samples using the DNeasy Blood & Tissue Kit (Qiagen) following manufacturer instructions. Individually barcoded samples included 23 Bruce spanworms, 19 winter moths, four putative hybrids, and two Operophtera fagata Scharfenberg samples included as an outgroup (Table S1). Library preparation and sequencing were performed by Floragenex Inc. (Eugene, OR) using the methods of Baird et al. (Reference Baird, Etter, Atwood, Currey, Shiver, Lewis, Selker, Cresko and Johnson2008). Sequence reads for each sample were aligned with the software BWA (Li & Durbin, Reference Li and Durbin2009) to the recently published winter moth genome (Derks et al., Reference Derks, Smit, Salis, Schijlen, Bossers, Mateman, Pijl, de Ridder, Groenen, Visser and Megens2015) accessed from NCBI BioProject (Accession Number PRJNA263715). The resulting alignments were used to identify SNPs with the software Stacks (Catchen et al., Reference Catchen, Hohenlohe, Bassham, Amores and Cresko2013). Summary statistics were generated using a minimum read coverage of 4× and including only loci present in ≥75% of individuals. One random SNP was exported from each ‘stack’ for subsequent analysis. To evaluate whether individuals putatively identified as Bruce spanworm, winter moth, or hybrids could be identified with these SNPs, the program fastStructure (Raj et al., Reference Raj, Stephens and Pritchard2014) was used to calculate the probability of assignment (q) of individuals to k = 2 clusters.
To identify a subset of SNPs diagnostic for distinguishing winter moth and Bruce spanworm and their hybrids, the output from Stacks was filtered to exclude polymorphic sites that were within 50 bp of each other and sorted by allele frequencies to identify those that were fixed in Bruce spanworm and winter moth. For each candidate SNP, 300 bp of flanking sequence was exported and a single primer pair was designed for each locus using the Primer3 v. 2.3.4 plug-in (Untergasser et al., Reference Untergasser, Cutcutache, Koressaar, Ye, Faircloth, Remm and Rozen2012) in Geneious v. 8.1.6 (Kearse et al., Reference Kearse, Moir, Wilson, Stones-Havas, Cheung, Sturrock, Buxton, Cooper, Markowitz, Duran, Thierer, Ashton, Mentjies and Drummond2012), with the following specifications: primers at least 25 base pairs from the SNP, minimum product size of 100 bp, maximum product size of 300 bp, minimum T m of 52°C, and minimum 20% GC content.
Microsatellite development
Genomic reads for microsatellite discovery were generated using an Ion Torrent Personal Genome Machine Sequencer (ThermoFisher). DNA was extracted from the head and thorax of a single male winter moth collected in a pheromone trap (see Elkinton, et al., Reference Elkinton, Boettner, Sremac, Gwiazdowski, Hunkins, Callahan, Scheufele, Donahue, Porter, Khrimian, Whited and Campbell2010 for trapping methods) in Lugo, Spain in December 2008 using the DNeasy Blood & Tissue Kit. 100 ng of DNA was sheared to approximately 300–500 bp fragments using Ion Shear enzymatic fragmentation and a library was prepared using the Ion Xpress Plus Library Kit (Life Technologies). Size selection was performed using a Pippin-Prep Blue (Sage Science), targeting 475 bp fragments. The library was prepared using the Ion PGM Sequencing 400 kit (Life Technologies). Sequencing was performed using 850 flows on an Ion 318 sequencing chip. Base-calling was performed with TorrentSuite 2.2 (Life Technologies).
Microsatellite discovery and primer design were performed using QDD 3.1.2 (Meglecz et al., Reference Meglecz, Pech, Gilles, Dubut, Hingamp, Trilles, Grenier and Martin2014). The default parameters were used except that the minimum sequence length was set to 200 bp. Finally, sequences with successful primer design were compared with the RepBase database of known transposable elements (Jurka, Reference Jurka2000) to identify potentially duplicated genomic regions. The ‘best’ primer pair per target sequence was selected as described in Meglecz et al., (Reference Meglecz, Pech, Gilles, Dubut, Hingamp, Trilles, Grenier and Martin2014), based on primer alignment score, the distance from the microsatellite, and length of the PCR product. In addition, loci were selected that contained only pure microsatellites (not compound) with at least eight repeats, and fewer than ten reads in consensus sequences (to avoid duplicated loci).
The resulting primer pairs were tested for amplification success in four male winter moth samples collected in the Czech Republic, Spain, Scotland, and Italy, and four male Bruce spanworm samples from the USA collected in Washington, New Jersey, New York, and Pennsylvania. Reverse primers were modified with the addition of a 5′ pig-tail sequence (GTTT) to reduce stutter (Brownstein et al., Reference Brownstein, Carpten and Smith1996), and forward primers were modified with the addition of an M13 tail (TCCCAGTCACGACGT) to allow incorporation of a 6-FAM labeled M13 primer during PCR (Schuelke, Reference Schuelke2000). PCR reactions (10 µl) contained 1× PCR buffer, 1.0 µl dNTPs (10 mm each), 0.8 µl MgCl2 (25 mm), 0.1 µl BSA (10 mg ml−1), 0.025 µl of forward primer (10 mm), 0.25 µl of reverse primer (10 mm), 0.05 µl of fluorescently labeled M13 primer (100 mm), 0.10 µl Go Taq DNA polymerase (Promega), and 1.0 µl sample DNA template. PCRs were amplified using a touchdown program with an annealing temperature that started at 61°C and decreased 2°C for each of the first five cycles, then continued at 51°C for 30 cycles. PCR products were run on an ABI 3730 sequencer (Life Technologies) at the DNA Analysis Facility on Science Hill at Yale University with a Liz 500 internal size standard (Gel Company). Alleles were scored using the software Geneious.
Loci that had detectable peaks for at least three of four winter moth samples were further characterized using 24 male moths each from two winter moth populations collected with pheromone traps in Marcillac, France and Reinhardshagen, Germany. Primer pairs that also had readable peak patterns in Bruce spanworm were also characterized using 24 male moths each from two Bruce spanworm populations collected with pheromone traps in Harbor Creek, Pennsylvania and Oconto, Wisconsin. These populations were located approximately 600 and 1300 km, respectively, from the nearest known winter moth population (Elkinton et al., Reference Elkinton, Boettner, Sremac, Gwiazdowski, Hunkins, Callahan, Scheufele, Donahue, Porter, Khrimian, Whited and Campbell2010, Reference Elkinton, Liebhold, Boettner and Sremac2014).
ARLEQUIN 3.5 (Excoffier et al., Reference Excoffier, Laval and Schneider2005) was used to calculate allelic diversity, observed and expected heterozygosities, to test for Hardy–Weinberg equilibrium, to test for linkage disequilibrium, and to estimate differentiation between populations and species. Multiple comparisons were controlled using the method of Benjamini & Hochberg (Reference Benjamini and Hochberg1995) with a false discovery rate of 0.05. MICRO-CHECKER 2.2.3 (van Oosterhout et al., Reference van Oosterhout, Hutchinson, Wills and Shipley2004) was used to test for null alleles.
Finally, the samples from the four winter moth and Bruce spanworm populations were used to assess multiplexing combinations of 24 loci that amplified most consistently and were in Hardy–Weinberg equilibrium. The ten of these loci that amplified for both species were kept together for multiplexed combinations to allow for efficient genotyping of populations with potential hybrids.
Analysis of field-collected putative hybrids
To evaluate the ability to classify hybrids using the ten microsatellite loci that consistently co-amplified for winter moth and Bruce spanworm, hybrid genotypes were simulated with Hybridlab 1.0 (Nielsen et al., Reference Nielsen, Bach and Kotlicki2006). Parental genotypes from 43 known Bruce spanworms from Maine, New York, and Pennsylvania, and 43 known winter moths from Connecticut, Massachusetts, Rhode Island, New Hampshire, and New York were used to generate 100 genotypes each of pure winter moth, pure Bruce spanworm, F1 hybrids, F2 hybrids, winter moth backcrosses, and Bruce spanworm backcrosses. A data set that included the 86 known samples, plus the 600 simulated samples was analyzed with Structure 2.3.4 (Pritchard et al., Reference Pritchard, Stephens and Donnelly2000) and NewHybrids 1.1 (Anderson & Thompson, Reference Anderson and Thompson2002). Structure was run with k = 2 clusters, known individuals coded as learning samples for updating allele frequencies, and 20,000 burn-in iterations followed by 100,000 sample iterations. NewHybrids was run with 10,000 burn-in iterations followed by 100,000 sample iterations, with known individuals identified. The probability of assignment was estimated for four different genotype frequency categories: each pure species, F1 hybrids, F2 hybrids, and backcross to each species.
Fifty-four putatively hybrid male moths collected with pheromone traps in Massachusetts, Rhode Island, Connecticut, and New Hampshire that were previously classified using the G6PD gene (Elkinton et al., Reference Elkinton, Liebhold, Boettner and Sremac2014) were genotyped using the ten microsatellite loci. Genetic ancestry and hybrid class of the putatively hybrid moths were analyzed using Structure and NewHybrids with the same run conditions as above. The data set included the putative hybrids plus the 43 reference moths for each parental species.
For each putative hybrid, the 658 base pair ‘DNA Barcoding’ portion of the mitochondrial COI gene was amplified and sequenced using the primers LepF1 and LepR1 (Hebert et al., Reference Hebert, Penton, Burns, Janzen and Hallwachs2004). Sequences were compared with known samples of each species (Gwiazdowski et al., Reference Gwiazdowski, Elkinton, deWaard and Sremac2013) to determine the maternal lineage of each sample.
Laboratory crosses
Laboratory crosses were performed to produce pure winter moth, pure Bruce spanworm, and hybrid offspring with reciprocal female and male parents from each species. In 2014 and 2015, female moths were reared from field collected larvae, while males were either reared from larvae or collected as adults using pheromone traps. The larvae and adults were collected in Massachusetts on either side of the documented hybrid zone: winter moth in the east and Bruce spanworm in the west (Elkinton et al., Reference Elkinton, Boettner, Sremac, Gwiazdowski, Hunkins, Callahan, Scheufele, Donahue, Porter, Khrimian, Whited and Campbell2010). Larvae were reared in ventilated 5-gallon buckets with the tree foliage from which they were found. When they were in their last instar, sifted peat moss was added to the bottom of the buckets for pupation. The resulting pupae were sifted, sorted by sex, and stored at 10–12°C until adult emergence at the end of October to beginning of November.
Matings were set up in clear, 12-ounce cups lined with paper, in an incubator (Percival) set at 10–12°C with a 13.5 h dark to 10.5 h light schedule. One or two females were combined with 3–30 males, as available, as they emerged (median ratio of 2:5). Females laid eggs on the paper lining of the mating cups. The temperature in the incubator with the eggs was reduced to 9.5°C on 4 December, to 5°C on 12 December, and to 2°C on 18 January.
Egg condition was measured 4–5 weeks following mating. Eggs were scored according to color (orange or green) as an indication of viability (Peterson, Reference Peterson1962; du Merle & Brunet, Reference du Merle and Brunet1991). Orange eggs were considered to have been fertilized and were classified as viable while green eggs were considered unfertilized and were classified as non-viable. This proved true as green eggs inevitably shriveled and died. Egg viability data were analyzed using a Chi-square (χ2) test to determine differences between observed and expected proportions of viable eggs among the four cross types using RStudio 0.98.501 (RStudio Team, 2015).
Results
RAD-seq SNP development
De-multiplexed reads from 46 of the 48 moth samples (sequencing failed for one winter moth and one putative hybrid sample) were aligned to the reference winter moth genome. A total of 92,852,319 sequence reads (mean of 2,018,528 reads per sample with mean coverage of 258×; Table S1) were analyzed with Stacks, which identified 11,779 stacks with 104,033 polymorphic SNPs. Filtering these loci to include one SNP per stack, and those with at least 4X coverage in at least 75% of individuals resulted in 1216 polymorphic loci. Analysis of these loci with fastStructure correctly identified all pure species individuals and identified two of the three putative hybrids as having mixed ancestry (fig. 1). The third putative hybrid was identified as pure Bruce spanworm.

Fig. 1. Bar plot showing individual assignment to Bruce spanworm or winter moth for the known pure species and putative hybrids analyzed with 1216 SNPs using fastStructure (Raj et al., Reference Raj, Stephens and Pritchard2014). The height of each colored bar represents the proportion of an individual's genotype assigned to each species. The individual sample number is below each bar corresponding to Table S1.
Further filtering of the 1216 SNPs to those that were separated by >50 bp and fixed in each species identified 95 SNPs suitable for species-diagnostic marker development. A file with 300 bp of leading and trailing sequence centered on each SNP, and a list of the primer sequences with expected fragment size are available in the Supplementary Materials.
Microsatellite development
Genomic sequencing for microsatellite discovery resulted in 4,461,215 sequences after filtering for polyclonal and low-quality reads. Median sequence length was 373 bp. Raw sequence reads have been deposited in NCBI BioProject Accession Number PRJNA330716. A total of 3,551,925 sequences were at least 200 bp long, and 207,331 contained microsatellites. The sequence comparison and alignment step of analysis with QDD took 18.6 days to complete on a High Performance Computing Cluster at Yale University. This resulted in 14,491 singletons and 14,148 unique consensus sequences, for a total of 28,639 from which to design primer pairs. Primers were successfully designed for 17,141 of these. Stringent filtering resulted in 104 loci. Of these, 17 matched an arthropod transposable element (12 from Lepidoptera) and were discarded. An additional primer pair was also discarded because the reverse primer had a TTTCTTT sequence on the 5′ end that might have interfered with the pig-tail modification, leaving 86 remaining primer pairs to test. Forty-six of these had decipherable peak patterns, with one or two alleles in at least three of the four winter moth individuals tested, and were further characterized in two winter moth populations (tables 1 and 2). These loci included 22 dinucleotides, 20 trinucleotides, two tetranucleotides, and two pentanucleotides. Fourteen primer pairs also had decipherable peak patterns for at least three of the four Bruce spanworm individuals and were further characterized with two Bruce spanworm populations (table 2). Primer sequences are shown in Table S2.
Table 1. Characterization of 32 microsatellite loci tested with 24 individuals each from winter moth populations from Reinhardshagen, Germany and Marcillac, France.

Asterisks indicate significant deviation from Hardy–Weinberg equilibrium after correcting for false discovery rate.
Table 2. Characterization of 14 microsatellite loci tested with 24 individuals each from two winter moth populations from Reinhardshagen, Germany and Marcillac, France, and two Bruce spanworm populations from Harborcreek, Pennsylvania and Oconto, Wisconsin.

Asterisks indicate significant deviation from Hardy–Weinberg equilibrium using after correcting for false discovery rate.
The number of alleles per locus ranged from 4 to 30 in winter moth and 3 to 14 in Bruce spanworm (tables 1 and 2). Four loci deviated significantly from Hardy–Weinberg equilibrium in both of the winter moth populations, and ten loci deviated in one of the populations (tables 1 and 2), after controlling for false discovery rate. Two loci deviated from Hardy–Weinberg equilibrium in both of the Bruce spanworm populations, and two deviated in one of the populations (table 2). For winter moth, significant linkage was found for 19 out of 1035 pairwise comparisons of loci in the Marcillac winter moth population, and three in the Reinhardshagen population after controlling for false discovery rate. WM03296 and WM12042 was the only pair that was linked in both populations. No loci were linked in either Bruce spanworm population. Seven loci showed evidence for null alleles in both winter moth populations, five in Reinhardshagen only, and four in Marcillac only. One locus showed evidence of null alleles in both Bruce Spanworm populations, one in Harborcreek only, and one in Oconto only (Table S3).
The 24 loci that amplified most consistently for winter moth and were in Hardy–Weinberg equilibrium in at least one of the populations are indicated in Table S2. Using these 24 loci, genetic differentiation (F st) between the two winter moth populations was 0.0067. Using the ten loci that amplified consistently in both species, F st between the two winter moth populations was 0.0007, and between the two Bruce spanworm populations was 0.0130. Pairwise F st values between populations of different species were: 0.3824 (Marcillac versus Harborcreek), 0.4030 (Marcillac versus Oconto), 0.3917 (Reinhardshagen versus Harborcreek), and 0.4112 (Reinhardshagen versus Oconto). When populations were combined (N = 48), F st between species was 0.3952. Suitable PCR multiplex combinations for the 24 tested loci found are shown in Table S4.
Analysis of field-collected putative hybrids
Structure analysis of simulated pure parental and hybrid genotypes (fig. 2) resulted in a mean probability of assignment (q) to Bruce spanworm of 0.97 for simulated Bruce spanworm genotypes, and 0.03 for simulated winter moth winter moth genotypes. Mean q values for simulated F1, F2, Bruce spanworm backcross, and winter moth backcross genotypes were 0.74, 0.51, 0.50, and 0.25, respectively. NewHybrids analysis resulted in simulated genotypes being assigned with the highest probability to the correct hybrid class for 100% of simulated F1s (mean probability = 0.92), 62% of simulated F2s (mean probability = 0.64), for 91% of simulated Bruce spanworm backcrosses (mean probability = 0.84), and 95% of simulated winter moth backcrosses (mean probability = 0.84).

Fig. 2. Histograms of the probability of species assignment (q) using Structure (Pritchard et al., Reference Pritchard, Stephens and Donnelly2000) for simulated pure species and select hybrid classes using ten microsatellite loci.
Structure analysis of field-collected putative hybrids resulted in q values ranging from 0.009 to 0.991, indicating a mix of pure Bruce spanworm, hybrids, and pure winter moth (fig. 3; Table S5). Assignment of putative hybrids to pure species or hybrid classes using NewHybrids resulted in 6 moths assigned with highest probability as Bruce spanworm, three as winter moth, 24 as F1 hybrids, three as F2 hybrids, 18 as winter moth backcross, and none as Bruce spanworm backcrosses (Table S5).

Fig. 3. Bar plot showing individual assignment to Bruce spanworm or winter moth for the known pure species and putative hybrids analyzed with ten microsatellite loci using Structure (Pritchard et al., Reference Pritchard, Stephens and Donnelly2000). The height of each colored bar represents the proportion of an individual's genotype assigned to each species.
Of the 45 moths classified with highest probability to a hybrid class, 30 had COI mitochondrial haplotypes matching winter moth and 15 had haplotypes matching Bruce spanworm (Table S5). Sequences have been deposited in GenBank (Accessions KX591003–KX591056). Broken down by hybrid class, 18 out of 24 F1 hybrids (75%), one out of three F2 hybrids (33%), and 11 out of 18 (61%) winter moth backcrosses had winter moth COI haplotypes.
Laboratory crosses
There was no significant difference in the proportion of viable eggs between years, so the data were combined for analysis. Sixty-two females laid 2262 eggs in the pure winter moth crosses (mean = 42.9 eggs/female), 39 females laid 1869 eggs in the pure Bruce spanworm crosses (mean = 47.9 eggs/female), 358 females laid 14,023 eggs in the winter moth female × Bruce spanworm male crosses (mean = 39.2 eggs/female), and six females laid 344 eggs in the Bruce spanworm female × winter moth male crosses (mean = 57.3 eggs/female). The numbers of crosses that include Bruce spanworm females were limited by the difficulty in collecting Bruce spanworm larvae, which are rare. There was a strong effect of cross type on egg viability (χ2 = 2044.9, df = 3, P < 0.0001). Pure winter moth and pure Bruce spanworm crosses had 93.5 and 94.1% viable eggs, respectively, while crosses between winter moth females and Bruce spanworm males had 60.8% viable eggs and crosses between Bruce spanworm females and winter moth males had 22.1% viable eggs (fig. 4). The latter two proportions were significantly different (χ2 = 209.56, df = 1, P < 0.0001).

Fig. 4. Proportion of viable eggs by laboratory cross type between pure winter moths, pure Bruce spanworms, and their reciprocal female to male hybrids.
Discussion
The SNP and microsatellite loci that we developed will be useful for studying population dynamics of winter moth in Europe, the introduction history of winter moth in North America, and spatial patterns of hybridization between winter moth and Bruce spanworm. The suite of SNP loci, together with the available draft winter moth genome (Derks et al., Reference Derks, Smit, Salis, Schijlen, Bossers, Mateman, Pijl, de Ridder, Groenen, Visser and Megens2015), could be especially useful for tracking selective introgression of particular regions of the genome.
Isolation of microsatellites by stringent filtering of genomic DNA reads has successfully identified a suite of 24 loci useful for winter moth population genetics. Ten of these loci were also shown to be suitable for characterizing hybrids between winter moth and Bruce spanworm through simulation of hybrid genotypes. Additional suitable loci could undoubtedly be developed by using less stringent filtering parameters than we used in this study.
Our study confirms the presence of field collected F1 hybrid moths and is the first to report later hybrid generations (F2 and winter moth backcross individuals), which opens up the possibility that selective introgression between the species could be impacting invasion. For example, hybridizing with a ubiquitous native species could relax Allee effects (the difficulty of finding mates) that limit population growth and establishment (Elkinton et al., Reference Elkinton, Liebhold, Boettner and Sremac2014). In addition, since Bruce spanworm can persist in colder regions than winter moth (Roelofs et al., Reference Roelofs, Hill, Linn, Meinwald, Jain, Herbert and Smith1982; Elkinton et al., Reference Elkinton, Boettner, Sremac, Gwiazdowski, Hunkins, Callahan, Scheufele, Donahue, Porter, Khrimian, Whited and Campbell2010), the transfer of alleles that confer cold tolerance could allow winter moth to invade colder areas than if they were not hybridizing.
Finally, the microsatellite data indicate a strong asymmetry among winter moth/Bruce spanworm hybrids in the Northeastern USA. Forty percent of the hybrids were classified as winter moth backcrosses and none were Bruce spanworm backcrosses. In addition, 75% of the moths classified as F1 hybrids had winter moth COI haplotypes (Table S5) indicating that it is more likely for winter moth females to successfully reproduce with Bruce spanworm males than the reciprocal cross. This result is consistent with Smith & Ring (unpublished, cited in Underhill et al., Reference Underhill, Millar, Ring, Wong, Barton and Giblin1987), who found that winter moth females and Bruce spanworm males from British Columbia could successfully hybridize in cages to at least the F2 generation, but that the reciprocal cross between Bruce spanworm females and winter moth males did not occur. Our data did not indicate a similarly complete directional incompatibility. This discrepancy could be explained by the Bruce spanworms in British Columbia observed by Smith & Ring being a different subspecies (O. brumata occidentalis) than those found on the East Coast of North America (Gwiazdowski et al., Reference Gwiazdowski, Elkinton, deWaard and Sremac2013). The discrepancy could also be an artifact of the imperfect ability to assign individuals to a hybrid class. Among our simulated hybrid genotypes, 8% of the F2 hybrids were incorrectly classified as F1 hybrids (fig. 2), so some of the field collected moths that were classified as F1 hybrids could actually be F2. Laboratory crosses between winter moth females and Bruce spanworm males resulted in a higher proportion of viable eggs than the reciprocal cross (fig. 4), supporting this pattern, and implicating asymmetrical egg viability as an important source of mortality contributing to this pattern.
There are several possible causes for the observed asymmetrical, sex-biased hybridization between winter moth and Bruce spanworm. It could simply be a result of population demographics in the sampled region, where winter moths greatly outnumber Bruce spanworms during an outbreak (Elkinton et al., Reference Elkinton, Boettner, Liebhold and Gwiazdowski2015), making it less likely to sample Bruce spanworm backcrosses. Alternatively, there could be a sex-linked genetic incompatibility. Females are the heterogametic sex in Lepidoptera and there could be genes on the Bruce spanworm sex chromosome that are incompatible with genes on winter moth autosomes (Jiggins et al., Reference Jiggins, Linares, Naisbit, Salazar, Yang and Mallet2001). Finally, a symbiotic microorganism, such as Wolbachia, could be causing sex-biased directional incompatibility. Bruce spanworms could harbor a maternally transmitted symbiont that is incompatible with the winter moth genome or with a winter moth symbiont (Werren et al., Reference Werren, Baldo and Clark2008). Sequencing the genome of a winter moth collected in the Netherlands detected Wolbachia DNA (Derks et al., Reference Derks, Smit, Salis, Schijlen, Bossers, Mateman, Pijl, de Ridder, Groenen, Visser and Megens2015), but it is not yet known whether it is present in North American populations of winter moth or Bruce spanworm.
The observation of multi-generation, asymmetrical hybridization between non-native winter moths and native Bruce spanworms, highlights the need to further study the causes and consequences of this pattern. For example, samples from across the ranges of both species and through hybrid zones can be used to test whether hybridization could be contributing to the unexpected new outbreak of winter moth in the Northeastern USA.
Acknowledgements
The authors gratefully acknowledge funding from the USDA Forest Service, Forest Health Technology Enterprise Team, and The Norwegian Institute for Bioeconomy Research (NIBIO). We are grateful to Jason Boone and Jenna Gribbon (Floragenex) for RAD-seq support and Ben Evans (Yale University) who helped with preliminary data analysis. We thank Frank Krueger (Nordwestdeutsche Forstliche Versuchsanstalt), Sandy Gardosik (Pennsylvania Department of Agriculture), Linda Williams (Wisconsin Department of Natural Resources), Alain Roques (Institut National de la Recherche Agronomique), Marek Turcáni (Czech University of Life Sciences), Maria Lombardero (University of Santiago de Compostela), Adam Vanbergen (NERC Centre for Ecology & Hydrology), Chad Phillips (Washington State Department of Agriculture), Saul Vaiciunas (New Jersey Department of Agriculture), Kenneth Carnes (New York State Department of Agriculture & Markets), and Andrea Battisti (University of Padova) for collecting the samples used in this study. We thank Ida Fløystad, Siv Grethe Aarnes, and DeAdra Newman for technical assistance.
Supplementary material
The supplementary material for this article can be found at http://dx.doi.org/10.1017/S0007485316000857.