Genome-wide markers show continental structuring and mitonuclear discordance in the forest tent caterpillar (Malacosoma disstria Hübner) (Lepidoptera: Lasiocampidae)

Abstract

The forest tent caterpillar, Malacosoma disstria Hübner (Lepidoptera: Lasiocampidae), is an irruptive forest pest found throughout North America. Widespread species such as M. disstria are exposed to historical and contemporary processes that are not uniform and can generate regionally distinct genomic variation. Previous analyses used a short mitochondrial fragment to infer broad-scale phylogeographic patterns in M. disstria, whereas nuclear markers have been previously applied only in a smaller geographic region. In this study, we quantified M. disstria population variation with genome-wide single nucleotide polymorphisms and cytochrome c oxidase from mitochondrial DNA. Using highly variable genome-wide markers, we resolved clear genomic differences among populations of M. disstria east of the Rocky Mountains that were not detected using mitochondrial variation alone. We also did not detect host-associated divergence in either our genomic or mitochondrial data. Our results highlight the utility of genome-wide markers to resolve intraspecific population structure within a widespread species and support the need for further biogeographic sampling of this forest insect pest.

Introduction

Insect populations are structured by historical and contemporary eco-evolutionary processes interacting at a range of spatial scales. Delimiting and describing spatial genetic structure are fundamental steps towards understanding and disentangling these processes. Mitochondrial DNA has been used extensively to resolve relationships within and among closely related species and was foundational to the field of phylogeography (Avise et al. 1987; Avise 2000). However, due to its reduced effective population size and unique pattern of maternal inheritance, this single genetic region can be shaped by neutral or selective forces that do not influence other (i.e., nuclear) genetic markers (Funk and Omland 2003; Ballard and Whitlock 2004; Rubinoff et al. 2006). Hence, the evolutionary history of the mitochondrial genome may not reflect the overall history of a species and may provide limited or contradictory information relative to the nuclear genome for inter- and intraspecific genomic variation, leading to observations of mitonuclear discordance (Roe and Sperling 2007; Galtier et al. 2009; Dupuis et al. 2012; Toews and Brelsford 2012; Teske et al. 2018; Campbell et al. 2022). Detecting discordance between the mitochondrial and nuclear genomes requires a complementary approach, with both mitochondrial and nuclear markers sampled from the same individuals. Intraspecific variation, especially if it is recently derived, often necessitates large numbers of independent and variable loci to accurately infer population and demographic histories and to resolve spatial genomic variation (Hey and Nielsen 2004; Carling and Brumfield 2007; Gompert et al. 2014; Garrick et al. 2015; Satler et al. 2021). Recent technological advances provide access to single nucleotide polymorphisms (SNPs) spread throughout the genome (Baird et al. 2008; Peterson et al. 2012; McCormack et al. 2013; Garrick et al. 2015) and have been used successfully to resolve spatial genetic variation in many widespread, nonmodel organisms (e.g., Bagley et al. 2017; Lumley et al. 2020; MacDonald et al. 2020; Cairns et al. 2021; Polic et al. 2022).

The forest tent caterpillar, Malacosoma disstria Hübner (Lepidoptera: Lasiocampidae), is widely distributed throughout North America, where it exhibits cyclical outbreaks and is an important defoliator of deciduous trees (Cooke et al. 2011; Schowalter 2017). The species feeds on a number of different host plants, including aspen (Populus spp.) (Salicaceae), maple (Acer spp., except A. rubrum) (Sapindaceae), oak (Quercus spp.) (Fagaceae), and birch (Betula spp.) (Betulaceae) in northern parts of its range, along with five additional host genera in the southern United States of America (Fitzgerald 1995, Table 5.1). Previous genetic surveys using mitochondrial DNA detected complex phylogeographic patterns within and among M. disstria populations, but underlying mechanisms were not identified (Lait and Hebert 2018). Lait and Hebert (2018) detected a distinct population west of the Rocky Mountains and a variable but unstructured population east of this barrier. Whereas the Rocky Mountains were hypothesised to be a significant barrier to gene flow, spatial genetic variation east of this barrier could not be explained by geographic factors alone (Lait and Hebert 2018). Furthermore, larval host plant was not explored as a potential factor structuring genetic diversity, despite M. disstria populations exhibiting significant performance differences on different host plants (e.g., Nicol et al. 1997; Parry and Goyer 2004). Recent work by MacDonald et al. (2022) used thousands of nuclear genomic markers to identify spatial, environmental, and host-associated divergence among regional populations of M. disstria in eastern Canada. This work explored the fine-scale genomic structuring among M. disstria collected on distinct host plants across a latitudinal gradient. MacDonald et al. (2022) showed that host association, as well as temperature and isolation by distance, collectively shaped genomic variation among eastern M. disstria populations. Because this work focused only on populations in eastern Canada, the authors were unable to assess whether wider spatial variation may exist among M. disstria populations, nor were they able link their work to the existing mitochondrial DNA data.

Here, we describe the population structure of M. disstria in northern boreal and mixed-hardwood forests of Canada east of the Rocky Mountains. We used reduced-representation sequencing to generate thousands of SNPs throughout the genome. Concurrently, we generated mitochondrial sequence data (cytochrome c oxidase 1 (CO1)) from the same individuals and combined these data with previously published data (Lait and Hebert 2018). Together, these nuclear and mitochondrial sequence data were used to: (1) describe the genomic population structure of M. disstria throughout the northern boreal and mixed-hardwood forests of Canada; (2) determine whether mitochondrial complexity is related to host association; (3) assess concordance between genome-wide SNPs and mitochondrial DNA, and (4) conclude whether M. disstria in northern boreal and mixed-hardwood forests represent a single variable population or distinct regional groups with discrete genetic clustering. Our genome-wide markers show that M. disstria exhibit regional population structure east of the Rocky Mountains, but this structure was not related to host association. We also demonstrated the importance of using genome-wide markers when describing intraspecific variation in this widespread species and highlighted the discordance observed between genomic and mitochondrial variation.

Methods

Study organism

Malacosoma disstria has a univoltine life cycle and overwinters as pharate larvae in egg masses laid on host plants during the previous summer (Gray and Ostaff 2012). Females deposit a single egg mass onto a selected host, and larvae emerge in the spring as a family group to feed on new leaves. Synchrony of larval emergence and leaf flush in their host plant is critical for survival (Parry et al. 1998; Despland and Noseworthy 2006; McClure and Despland 2010). Siblings form cohesive groups, foraging and resting together within the host plant canopy (Despland and Hamzeh 2004). Larvae remain on their natal host until their later instars and then may disperse to seek additional food sources (Batzer et al. 1995). Adult M. disstria are capital breeders (i.e., all or most adult energy is obtained as a larvae) and are relatively short lived. Dispersal by flight is male-biased (Struble 1970), whereas females have a short obligatory flight before oviposition (Miller 2006). The dispersal abilities of M. disstria are limited, with male moths capable of flying up to 3 km (Evenden et al. 2015), but wind-assisted dispersal events over large distances (of more than 100 km) are possible (Brown 1965).

Sample collection and DNA extraction

We sampled M. disstria individuals from Alberta (n = 18) and Saskatchewan (n = 6) in central Canada and from Ontario (n = 139) and Quebec (n = 10) in eastern Canada. Egg bands or early instar larvae were collected from the four principal host plants in the region: trembling aspen, Populus tremuloides Michaux (n = 131), sugar maple, Acer saccharum Marshall (n = 32), white birch, Betula papyrifera Marshall (n = 4), and northern red oak, Quercus rubra Linnaeus (n = 131). These hosts were not evenly distributed throughout the range, with only trembling aspen located in central Canada and boreal sites in eastern Canada (Table 1). We sampled egg masses and early instar larvae (at least the fourth instar) because host identity represents the ovipositional choice of females. Sampling of individuals across central and eastern Canada was completed in the summer of 2018, and egg bands were shipped to the Great Lakes Forestry Centre, Sault Ste. Marie, Ontario. Upon arrival, insects were allowed to emerge and were reared at 27° C, 55% relative humidity, and 16:8-hour light:dark photoperiod in a Z-labtech Flex environmental chamber (Z-Sciences Corp., Westmount, Quebec) in the Insect Production and Quarantine Facility. We maintained individual families in separate rearing containers and fed larvae locally collected foliage from their natal hosts. We sampled third to fifth instars from each family (in case of colony collapse due to disease) and reared the remaining larvae to adults when possible. We did in fact lose a number of families to disease, so the larvae were the only samples we had left for extraction from many locations. Adults were frozen at –20° C, and larval specimens were preserved in 100% ethanol before storage at –20° C. Voucher specimens are maintained in a frozen tissue repository at the Great Lakes Forestry Centre.

Table 1. Summary of Malacosoma disstria collection details. Sample size represents the number of individuals used in each set of analyses.

* Final sample size following quality control filtering.

ddRADseq, double-digest restriction site–associated DNA library preparation and sequencing; mt, mitochondrial; AB, Alberta; ON, Ontario; SK, Saskatchewan; QC, Quebec; n/a, location data not available on GenBank.

At each location, we tried to collect egg bands or larvae from 10 to 20 separate trees (Table 1). However, collecting intensity and shipping varied between collaborators, and our final sample size therefore varied among sites, with fewer sampling sites in central Canada than in eastern Canada. Also, some larvae hatched in transit and intermixed with unrelated larvae from other egg bands. This could have resulted in the accidental inclusion of siblings within our analyses, so we performed a kinship analysis on our genomic data to screen for related individuals (see the “Data filtering and SNP identification” section below).

We dissected thorax tissue from adults (n = 63) and the head capsule and upper thorax from larvae (n = 110; Table 1). For each larva, we removed the digestive tract to reduce plant or microbial contamination. We extracted genomic DNA from each specimen using a Qiagen DNeasy Blood and Tissue Extraction kit (QIAGEN, Hilden, Germany), following the manufacturer’s protocols, with the addition of bovine pancreatic ribonuclease A (RNaseA, 4 uL at 100 mg/m; Sigma-Aldrich Canada Co., Oakville, Ontario), followed by ethanol precipitation and resuspension in Millipore water. Final extracts were stored at –20° C. The quality and quantity of each sample were measured using a Nanodrop 1000 spectrophotometer (Thermo Fisher Scientific, Mississauga, Ontario) and Qubit fluorometer (Thermo Fisher Scientific), respectively.

Genome-wide SNPs and analysis

Double-digest restriction site sequencing

We submitted 173 specimens to the Molecular Biology Service Unit (University of Alberta, Edmonton, Alberta, Canada) for double-digest restriction site–associated DNA library preparation and sequencing. Libraries were prepared using 200 ng of genomic DNA, using PstI and MspI restriction enzymes as in MacDonald et al. (2020), generally following Peterson et al. (2012) and Poland and Rife (2012). We sequenced our pooled library of individuals using single-end, 75-bp sequencing on an Illumina NextSeq 500 (Illumina, San Diego, California, United States of America) platform in one high-output flowcell.

Data filtering and SNP identification

Bioinformatic processing of raw sequence reads followed MacDonald et al. (2022). Briefly, we used Stacks 2.0 (Rochette et al. 2019) to demultiplex Illumina single-end, 75-bp reads, removing those with quality scores of less than 20 within a sliding window equal to 15% of the read length. Illumina index sequences and the first 5-bp from the 5′ end of each read (corresponding to the PstI restriction site) were removed using Cutadapt 1.9.1 (Martin et al. 2013). The resulting 62-bp reads were then aligned to a M. disstria genome assembly (National Center for Biotechnology Information Accession PRJNA824522; 1090 scaffolds, each over 10 kb) using Burrows–Wheeler Aligner 0.7.17 (BWA-MEM; Li and Durbin 2009; Li 2013). Following alignment, we called SNPs among sequenced individuals using Stacks 2.0, stipulating a single population. Filtering of the resulting data was completed using VCFtools, version 0.1.14 (Danecek et al. 2011). We removed: (1) individuals with more than 25% missing data; (2) loci with a read depth less than five; (3) loci with minor allele frequencies of at least 0.05; (4) loci with percentages of more than 5% missing data; and (5) one locus for every pair of loci that were within at least 10 kb of each other.

We screened our data set for full-sibling relationships that could bias assessments of population structure (O’Connell et al. 2019). We used the R package SNPRelate, version 1.26.0, (Zheng et al. 2012) in the R statistical computing environment (R Core Team 2021) to estimate pairwise kinship coefficients among all sequenced individuals. For diploid organisms, the coefficient value expected for full siblings is 0.25 (Manichaikul et al. 2010). A natural break in coefficient values occurred at 0.22. Coefficient values above 0.22 were observed only for pairs of individuals collected from the same location and were therefore inferred to indicate full-sibling relationships. For each pair, we removed the individual with the greater percentage of missing data. Using this final set of individuals, we reverted to original binary alignment map files, recalled SNPs, and repeated filtering using the same parameters listed above.

Population structure

We quantified population structure among our M. disstria samples using two different approaches: principal component analysis and the software programme Structure (Pritchard et al. 2000). Principal component analyses were conducted in R using the adegenet package, version 2.1.2 (Jombart 2008), and visualised using ggplot, version 2 3.2.1 (Wickham 2016; https://ggplot2.tidyverse.org). For Structure analyses, 10 independent runs were completed for each value of K = 1:10, each stipulating the admixture model, correlated allele frequencies, a burn-in period of 100 000, and 1 000 000 Markov chain–Monte Carlo (MCMC) repetitions. Location priors (set to collection localities) were used to inform the MCMC algorithm without biasing the model. Runs were visualised in CLUMPAK (Kopelman et al. 2015; http://clumpak.tau.ac.il/index.html). We identified the optimal value of K using the Pritchard (Pritchard et al. 2000) and Evanno methods (Evanno et al. 2005). Structure admixture plots were visualised using ggplot2.

Mitochondrial phylogeography

We amplified a 658-bp CO1 fragment for 130 specimens using the LepF and LepR primers based on a modified protocol from Hebert et al. (2003). The majority of specimens were also used for double-digest restriction site–associated DNA library preparation and sequencing. However, due to stringent filtering for our downstream SNPs analyses, some individuals were removed, leaving specimens with only CO1 data (n = 8). Each polymerase chain reaction (PCR) contained 9.76 μL of dH2O, 2 μL of 10× buffer, 2 μL of MgCl2, 0.4 μL of deoxynucleoside triphosphate, 0.4 μL of LepF primer, 0.4 μL of LepR primer, 0.04 μL of TopTaq DNA polymerase (QIAGEN, Hilden, Germany), and 5 μL of template DNA. The PCR thermal cycle profile consisted of one cycle at 94° C for 2 minutes, 35 cycles of 94° C for 2 minutes, 45° C for 30 seconds, 72° C for 2 minutes, and a final cycle of 72° C for 5 minutes, followed by a 4° C hold. Afterward, all PCR samples were purified with an exonuclease-shrimp alkaline phosphatase protocol following the manufacturer’s instructions (Exo-SAP-IT, Thermo Fisher Scientific, Waltham, Massachusetts, United States of America). Mitochondrial fragments were submitted for Sanger sequencing to the Molecular Biology Service Unit (University of Alberta) on an Applied Biosystems 3730 (Mississauga, Ontario). All sequences were aligned and quality checked using Geneious Prime, version 2019.2.3 (https://www.geneious.com/). We combined our data with an additional 139 previously published sequences (Lait and Hebert 2018) downloaded from the Barcode of Life Data System (15 January 2020; http://www.boldsystems.org). These included additional material from British Columbia (western) and southern United States of America (southern). We aligned this combined data set using MAFFT online, version 7 (Katoh et al. 2019; https://mafft.cbrc.jp/alignment/server/; default settings), and trimmed sequences to 658 bp using Mesquite, version 3.6 (Maddison and Maddison 2019; https://www.mesquiteproject.org/), to match the previously published data. We constructed Templeton, Crandall, and Sing haplotype networks (also known as statistical parsimony; Templeton et al. 1995; Clement et al. 2000) in PopART, version 1.7 (Bandelt et al. 1999; Leigh and Bryant 2015).

Results

Genome-wide SNPs and analysis

Across all 173 individuals, a total of 186 588 563 75-bp reads were sequenced and passed Illumina quality filters. After filtering and removing adapters and restriction sites, 97.3% of the 62-bp reads were successfully aligned to the M. disstria reference genome. After calling and filtering a preliminary set of SNPs, we identified and removed nine individuals with more than 25% missing data and 42 individuals that shared a full-sibling relationship with at least one other individual. We then recalled SNPs and filtered the data based on read depth, minor allele frequency, missing data, and physical proximity. A total of 2828 SNPs with a mean read depth of 26.95 across 122 individuals comprised our final genomic data set.

We observed two distinct clusters of M. disstria individuals (n = 122) based on genome-wide SNPs (Fig. 1). Our principal component analysis suggested distinct clusters of central (n = 14) and eastern (n = 108) individuals that separated along principle component 1, which explained 3.55% of total genomic variation (Fig. 1A). Using Structure results, we inferred an optimal number of populations as K = 3 (Pritchard method) or K = 2 (Evanno method; Fig. 1B). Both K = 2 and K = 3 identified central samples as a distinct group. We did not, however, see separation of individuals based on host association, with larva from P. tremuloides appearing in multiple clusters (Fig. 1).

Figure 1. Population genetic structure of Malacosoma disstria using A, Principal component (PC) analysis and B, model-based clustering using Structure. In the principal component analysis, every point represents a single individual (n = 122) that has been coded for larval host plant. We compared K = 1–10 for our Structure analyses and identified K = 2* as optimal using the Evanno method, and K = 3 with the Pritchard method (see the “Results” section for details). In the admixture plot from Structure individuals are ordered by increasing longitude and we indicate the associated larval host in the bar plot above.

Mitochondrial phylogeography

We sequenced 658 bp of the CO1 gene, corresponding to the DNA barcode region, for 130 M. disstria specimens (National Center for Biotechnology Information Genbank Accession numbers MT791498-MT791627). We combined these with an additional 139 sequences previously published by Lait and Hebert (2018; Supplementary material, Table S1). We found complex geographic structuring within the combined data set. Western populations formed a distinct group (Fig. 2), but the mitochondrial haplotypes from central, eastern, and southern populations were variable. These formed an unstructured population that did not align with specific geographic regions (Fig. 2). Furthermore, we observed a number of geographically distant individuals that shared identical haplotypes (Fig. 2). We also did not observe any relationship between mitochondrial genetic variation and host association, with larvae on different hosts sharing identical haplotypes and with larvae on the same host having different haplotypes (Fig. 2).

Figure 2. Templeton, Crandall, and Sing haplotype network of 269 Malacosoma disstria mitochondrial cytochrome c oxidase subunit 1 (CO1) sequences from across Canada and the United States of America. Sequence data from this study (n = 130) were combined with data from Lait and Hebert (2018; n = 139). Each square represents a single individual, with colour representing a major geographic region and letters indicating larval host plant (when known). Individuals were clustered into distinct CO1 haplotypes, with black dots representing unsampled, potential haplotypes and small black lines representing additional mutational differences.

Discussion

Widespread species are shaped by geographically distinct ecological and evolutionary forces, often giving rise to regionally structured populations with unique genetic signatures. Mitochondrial variation is often insufficient for resolving complex genetic patterns, and genome-wide markers are needed to adequately quantify spatial variation. Using genome-wide SNPs, we show that M. disstria east of the Rocky Mountains in Canada are structured into two distinct genetic clusters, similar to other widespread forest species in Canada (e.g., Choristoneura fumiferana; Lumley et al. 2020). Although we observed geographically structured clusters in our genomic data, discontinuous sampling between clusters means that we cannot conclusively resolve whether they are connected by a cline of genomic divergence or are separated by a discrete barrier to gene flow. Regardless, our genomic results contradict mitochondrial inferences of homogeneous population structure east of the Rocky Mountains (Lait and Hebert 2018). The lack of concordance between mitochondrial variation and genome-wide markers is not uncommon and points to the importance of using thousands of unlinked markers to assess intraspecific variation in widespread species.

Population genomics of the forest tent caterpillar

Temperate forest species can have complex evolutionary histories arising from shifting distributions during the glacial and interglacial cycles of the Quaternary Period (Hewitt 2004). Although the boreal forest currently represents an extensive contiguous biome, glacial cycles fragmented these habitats and allowed populations to adapt and diverge in isolation and during subsequent postglacial spread (Weir and Schluter 2004; Jaramillo-Correa et al. 2009). Combinations of drift and adaptation in these regions have led to intra- and interspecific diversification in many forest species. Similar population structures among different forest species reflect the common demographic responses to historic bioclimatic change and ecological co-associations (Garrick et al. 2021). For example, midcontinent biogeographic breaks, similar to what we observed in M. disstria, have been detected in other forest animals, including Choristoneura fumiferana (Lepidoptera: Tortricidae) (Lumley et al. 2020), Lynx canadensis (Carnivora: Felidae) (Stenseth et al. 1999; Rueness et al. 2003), Ursus americanus (Carnivora: Ursidae) (Puckett et al. 2015; but see Bradburd et al. 2018), Campanula americana (Campanulaceae) (Prior et al. 2020), and Eptesicus fuscus (Chiroptera: Vespertilionidae) (Yi and Latch 2022). Many boreal trees also show similar biogeographic patterns, including Picea glauca (Pinaceae) (de Lafontaine et al. 2010), Abies balsamifera (Pinaceae) (Cinget et al. 2015), and two Betula species (Thomson et al 2015). Herbivorous insects, such as M. disstria, would have tracked the fragmentation and expansion of their host plants, leaving similar genetic signatures within their populations (Smith et al. 2011; Satler and Carstens 2017). Another lepidopteran forest pest, the eastern spruce budworm, Choristoneura fumiferana, has a midcontinental break that separated the species into central and eastern populations (Lumley et al. 2020). A similar spatial structure was observed in P. glauca and A. balsamifera (de Lafontaine et al. 2010; Cinget et al. 2015), the primary host plants for C. fumiferana. Each of these authors suggests that this biogeographic pattern is a result of migration from distinct glacial refugia, creating a common history of expansion and spread (Lumley et al. 2020). Recent analyses of P. tremuloides, an important host for M. disstria, showed a midcontinent division between northeastern and northwestern populations (Goessen et al. 2022), similar to M. disstria. Goessen et al. (2022) proposed that these two northern P. tremuloides populations arose following distinct colonisation routes out of a single common refugium, possibly south of the Great Lakes (Jaramillo-Correa et al. 2009; Ding et al. 2017). It is conceivable that M. disstria followed a similar evolutionary trajectory, dispersing from a single refugium (see also Lait and Hebert 2018) and then evolving into distinct regional populations.

Discriminating population structure in widespread species is not trivial. Although we observed a discrete break between central and eastern groups of M. disstria, these clusters might represent a genomic cline rather than discrete populations (Bradburd et al. 2018). Genomic clines may develop through drift and adaptation so that genetic variation becomes spatially autocorrelated; as samples get further apart, genetic differences between samples also increase (Wright 1943). This phenomenon is common in nature (Schwartz and McKelvey 2008; Meirmans 2012; Perez et al. 2018), more so if dispersal is limited (Wright 1943; Slatkin 1985). Many analytical approaches used to identify population structure can artificially delimit clusters in continuous genetic data (Perez et al. 2018), particularly with patchy or uneven sampling (Meirmans 2019). This analytical conundrum has been termed the “cluster versus cline dilemma” (Guillot et al. 2009). Recent work on Populus balsamifera clearly demonstrates this issue. Keller et al. (2010) described two discrete clusters of P. balsamifera in central and eastern Canada, similar to the pattern we observed in M. disstria. However, when Meirmans et al. (2017) increased sampling between these hypothesised groups, they resolved a genetic cline rather than two distinct clusters. Distinguishing between a genetic cline and discrete populations depends on sampling. Uneven sampling across a distribution can strongly influence the analytical delimitation of genomic variation (Meirmans 2015, 2019). Our sampling of M. disstria was patchy, particularly between our central and eastern clusters. Improved sampling of this region is required to resolve whether the observed genetic structure is in fact discrete or continuous. Interestingly, P. tremuloides, a host plant for M. disstria, does show discrete population structure across this part of the range (Goessen et al. 2022). We are interested to see, with increased sampling, whether M. disstria mirrors the evolutionary history of its host plant.

Our sampling reflects variation in the diversity of available hosts, with higher diversity occurring in eastern Canada (e.g., Acer and Quercus) than in central Canada (P. tremuloides). We did not observe any relationship between population structure and host association for either genome-wide SNPs or mitochondrial haplotypes. The majority of variation in our genomic data set was spatially structured (central versus eastern Canada), even among individuals sampled from a single host (P. tremuloides; Fig. 1). However, recent work on the eastern population of M. disstria identified signals of host-associated genomic divergence within M. disstria that was not apparent in patterns of population structuring (MacDonald et al. 2022). It is unclear whether M. disstria is currently undergoing divergence related to host association that will eventually be reflected in patterns of population clustering or whether gene flow among host-associated groups is simply too high to allow discrete host-associated structuring to precipitate. Given the polyphagous nature of this species, it would be worthwhile to examine the historical biogeography of M. disstria in relation to regional host association across its entire North American range using whole-genome sequence data. Another widespread forest insect pest, Neodiprion lecontei (Hymenoptera: Diprionidae), has distinct genetic structuring due to Pleistocene drift and host-associated divergence (Bagley et al. 2017). Bagley et al. (2017) showed that distinct geographic clusters varied in their patterns of host-associated divergence. Given that MacDonald et al. (2022) also detected host-associated divergence in the eastern cluster of M. disstria, it would be valuable to assess whether similar patterns exist in other populations at broader spatial scales, contrasting spatially structured and host-associated genomic divergences across the species’ entire range. Earlier reciprocal transplant experiments clearly showed that functional differences exist between southern and northern populations of M. disstria (Parry and Goyer 2004), suggesting that at least some host-associated divergence exists among populations. Therefore, it would be valuable to expand our genomic sampling to these regional populations to assess the impact of geography and host association on genomic divergence.

Mitonuclear discordance

Mitochondrial DNA has been widely used to explore species-level variation and biogeography for more than 20 years (Avise 2000). However, there is mounting evidence that mitochondrial variation is not ideal for inferring population-level variation (Ballard and Whitlock 2004; Toews and Brelsford 2012; Teske et al. 2018, Morón-López et al. 2022) and that a shift towards unlinked genome-wide nuclear markers is occurring (Garrick et al. 2015). We observed discordance between the population structure inferred from genome-wide SNPs and that shown in the mitochondrial CO1 data. We resolved clear differences between central and eastern populations (Fig. 1) using SNPs, but no such patterns were resolved among CO1 haplotypes. In fact, a number of CO1 haplotypes were shared between geographic regions (Fig. 2). Previous work described CO1 admixture in the eastern region (Lait and Hebert 2018); however, we found no evidence of this substructure in our genomic data. We also did not see correspondence between CO1 variation and host association.

Discordance between nuclear and mitochondrial genomes is a well-known phenomenon (reviewed in Toews and Brelsford 2012) and has been described in a range of taxa (e.g., Wielstra and Arntzen 2020; Cairns et al. 2021; Morón-López et al. 2022; Yi and Latch 2022). Mitochondrial genomes, due to differences in inheritance, recombination, and effective population size, exhibit distinct responses to selection, demographic changes, and drift (Toews and Brelsford 2012; Bonnet et al. 2017). Because the mitochondrial genome is maternally inherited as a single independent unit with little to no recombination, it is susceptible to asymmetrical introgression (Linnen and Farrell 2007) and selective sweeps (Gompert et al. 2008; Wendt et al. 2022). These evolutionary events may limit the power of mitochondrial DNA to resolve intraspecific variation or to create a distinct evolutionary history when compared with the nuclear genome. This may be particularly prevalent in species with large populations (Gillespie 2000) or in species that have experienced rapid range expansions (Petit and Excoffier 2009). Similar discordant results were detected in Speyeria aglaja (Lepidoptera: Nymphalidae) (Polic et al. 2022), where distinct geographic lineages were detected with genome-wide SNPs but not with CO1 haplotypes. Polic et al. (2022) suggested that gene flow or selective sweeps homogenised haplotype diversity and erased the spatial variation in the mitochondrial genome. We hypothesise that similar events may have occurred within M. disstria. We acknowledge, however, that other evolutionary processes, such as adaptive selection (Bazin et al. 2006; Pavlova et al. 2013), sex-biased dispersal (Folt et al. 2019), or cytoplasmic incompatibilities (Hurst and Jiggins 2005), may have also led to incongruent mitochondrial and nuclear variation.

Future directions

Malacosoma disstria is a widespread species found throughout Canada and the United States of America, and our sampling represents only a portion of its total distribution. The combination of postglacial history, host association, and contemporary dispersal associated with irruptive outbreaks (e.g., James et al. 2015) could lead to a complex population structure within this species. We know that distinct mitochondrial haplotypes persist in M. disstria populations west of the Rocky Mountains (Lait and Hebert 2018), and this structure was not captured in our study due to the spatial extent of our sampling. Previous work has quantified functional variation among other M. disstria populations. For example, southern populations of M. disstria have distinct regional host associations (Parry and Goyer 2004) and show latitudinal clines in reproductive traits (Parry et al. 2001). Given that genomic differentiation was linked to both host association and variation in environmental conditions (MacDonald et al. 2022), it is plausible that additional population structure remains undetected throughout the range of M. disstria. Increased sampling across the entire range of M. disstria is needed to fully describe the population variation in this species and would provide insight to the historical biogeography and evolutionary history of this important forest pest.