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. Reference Avise, Arnold, Ball, Bermingham, Lamb and Neigel1987; Avise Reference Avise2000). 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 Reference Funk and Omland2003; Ballard and Whitlock Reference Ballard and Whitlock2004; Rubinoff et al. Reference Rubinoff, Cameron and Will2006). 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 Reference Roe and Sperling2007; Galtier et al. Reference Galtier, Nabholz, Glemin and Hurst2009; Dupuis et al. Reference Dupuis, Roe and Sperling2012; Toews and Brelsford Reference Toews and Brelsford2012; Teske et al. Reference Teske, Golla, Sandoval-Castillo, Emami-Khoyi and van der Lingen2018; Campbell et al. Reference Campbell, MacDonald, Gage, Gage and Sperling2022). 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 Reference Hey and Nielsen2004; Carling and Brumfield Reference Carling and Brumfield2007; Gompert et al. Reference Gompert, Lucas, Buerkle, Forister, Fordyce and Nice2014; Garrick et al. Reference Garrick, Bonatelli, Hyseni, Morales, Pelletier and Perez2015; Satler et al. Reference Satler, Carstens, Garrick, Espíndola and Mardulyn2021). Recent technological advances provide access to single nucleotide polymorphisms (SNPs) spread throughout the genome (Baird et al. Reference Baird, Etter, Atwood, Currey, Shiver and Lewis2008; Peterson et al. Reference Peterson, Weber, Kay, Fisher and Hoekstra2012; McCormack et al. Reference McCormack, Hird, Zellmer, Carstens and Brumfield2013; Garrick et al. Reference Garrick, Bonatelli, Hyseni, Morales, Pelletier and Perez2015) and have been used successfully to resolve spatial genetic variation in many widespread, nonmodel organisms (e.g., Bagley et al. Reference Bagley, Sousa, Niemiller and Linnen2017; Lumley et al. Reference Lumley, Pouliot, Laroche, Boyle, Brunet and Levesque2020; MacDonald et al. Reference MacDonald, Dupuis, Davis, Acorn, Nielsen and Sperling2020; Cairns et al. Reference Cairns, Cicchino, Stewart, Austin and Lougheed2021; Polic et al. Reference Polic, Yildirim, Lee, Franzen, Mutanen, Vila and Forsman2022).
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. Reference Cooke, MacQuarrie and Lorenzetti2011; Schowalter Reference Schowalter2017). 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 Reference Fitzgerald1995, 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 Reference Lait and Hebert2018). Lait and Hebert (Reference Lait and Hebert2018) 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 Reference Lait and Hebert2018). 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. Reference Nicol, Arnason, Helson and Abou-Zaid1997; Parry and Goyer Reference Parry and Goyer2004). Recent work by MacDonald et al. (Reference MacDonald, Snape, Roe and Sperling2022) 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. (Reference MacDonald, Snape, Roe and Sperling2022) 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 Reference Lait and Hebert2018). 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 Reference Gray and Ostaff2012). 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. Reference Parry, Spence and Volney1998; Despland and Noseworthy Reference Despland and Noseworthy2006; McClure and Despland Reference McClure and Despland2010). Siblings form cohesive groups, foraging and resting together within the host plant canopy (Despland and Hamzeh Reference Despland and Hamzeh2004). Larvae remain on their natal host until their later instars and then may disperse to seek additional food sources (Batzer et al. Reference Batzer, Martin, Mattson and Miller1995). 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 Reference Struble1970), whereas females have a short obligatory flight before oviposition (Miller Reference Miller2006). The dispersal abilities of M. disstria are limited, with male moths capable of flying up to 3 km (Evenden et al. Reference Evenden, Whitehouse and Jones2015), but wind-assisted dispersal events over large distances (of more than 100 km) are possible (Brown Reference Brown1965).
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. (Reference MacDonald, Dupuis, Davis, Acorn, Nielsen and Sperling2020), generally following Peterson et al. (Reference Peterson, Weber, Kay, Fisher and Hoekstra2012) and Poland and Rife (Reference Poland and Rife2012). 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. (Reference MacDonald, Snape, Roe and Sperling2022). Briefly, we used Stacks 2.0 (Rochette et al. Reference Rochette, Rivera-Colon and Catchen2019) 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. Reference Martin, Dasmahapatra, Nadeau, Salazar, Walters and Simpson2013). 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 Reference Li and Durbin2009; Li Reference Li2013). 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. Reference Danecek, Auton, Abecasis, Albers, Banks and DePristo2011). 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. Reference O’Connell, Mulder, Maldonado, Currie and Ferraro2019). We used the R package SNPRelate, version 1.26.0, (Zheng et al. Reference Zheng, Levine, Shen, Gogarten, Laurie and Weir2012) 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. Reference Manichaikul, Mychaleckyj, Rich, Daly, Sale and Chen2010). 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. Reference Pritchard, Stephens and Donnelly2000). Principal component analyses were conducted in R using the adegenet package, version 2.1.2 (Jombart Reference Jombart2008), and visualised using ggplot, version 2 3.2.1 (Wickham Reference Wickham2016; 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. Reference Kopelman, Mayzel, Jakobsson, Rosenberg and Mayrose2015; http://clumpak.tau.ac.il/index.html). We identified the optimal value of K using the Pritchard (Pritchard et al. Reference Pritchard, Stephens and Donnelly2000) and Evanno methods (Evanno et al. Reference Evanno, Regnaut and Goudet2005). 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. (Reference Hebert, Cywinska, Ball and deWaard2003). 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 Reference Lait and Hebert2018) 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. Reference Katoh, Rozewicki and Yamada2019; https://mafft.cbrc.jp/alignment/server/; default settings), and trimmed sequences to 658 bp using Mesquite, version 3.6 (Maddison and Maddison Reference Maddison and Maddison2019; 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. Reference Templeton, Routman and Phillips1995; Clement et al. Reference Clement, Posada and Crandall2000) in PopART, version 1.7 (Bandelt et al. Reference Bandelt, Forstr and Röhl1999; Leigh and Bryant Reference Leigh and Bryant2015).
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 (Reference Lait and Hebert2018; 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 (Reference Lait and Hebert2018; 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. Reference Lumley, Pouliot, Laroche, Boyle, Brunet and Levesque2020). 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 Reference Lait and Hebert2018). 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 Reference Hewitt2004). 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 Reference Weir and Schluter2004; Jaramillo-Correa et al. Reference Jaramillo-Correa, Beaulieu, Khasa and Bousquet2009). 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. Reference Garrick, Hyseni, Arantes, Zachos, Zee, Oliver and Lozier2021). 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. Reference Lumley, Pouliot, Laroche, Boyle, Brunet and Levesque2020), Lynx canadensis (Carnivora: Felidae) (Stenseth et al. Reference Stenseth, Chan, Tong, Boonstra, Boutin and Krebs1999; Rueness et al. Reference Rueness, Stenseth, O’Donoghue, Boutin, Ellegren and Jakobsen2003), Ursus americanus (Carnivora: Ursidae) (Puckett et al. Reference Puckett, Etter, Johnson and Eggert2015; but see Bradburd et al. Reference Bradburd, Coop and Ralph2018), Campanula americana (Campanulaceae) (Prior et al. Reference Prior, Layman, Koski, Galloway and Busch2020), and Eptesicus fuscus (Chiroptera: Vespertilionidae) (Yi and Latch Reference Yi and Latch2022). Many boreal trees also show similar biogeographic patterns, including Picea glauca (Pinaceae) (de Lafontaine et al. Reference de Lafontaine, Turgeon and Carine2010), Abies balsamifera (Pinaceae) (Cinget et al. Reference Cinget, de Lafontaine, Gerardi and Bousquet2015), and two Betula species (Thomson et al Reference Thomson, Dick, Dayanandan and Carine2015). 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. Reference Smith, Tank, Godsoe, Levenick, Strand, Esque and Pellmyr2011; Satler and Carstens Reference Satler and Carstens2017). 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. Reference Lumley, Pouliot, Laroche, Boyle, Brunet and Levesque2020). A similar spatial structure was observed in P. glauca and A. balsamifera (de Lafontaine et al. Reference de Lafontaine, Turgeon and Carine2010; Cinget et al. Reference Cinget, de Lafontaine, Gerardi and Bousquet2015), 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. Reference Lumley, Pouliot, Laroche, Boyle, Brunet and Levesque2020). Recent analyses of P. tremuloides, an important host for M. disstria, showed a midcontinent division between northeastern and northwestern populations (Goessen et al. Reference Goessen, Isabel, Wehenkel, Pavy, Tischenko and Touchette2022), similar to M. disstria. Goessen et al. (Reference Goessen, Isabel, Wehenkel, Pavy, Tischenko and Touchette2022) 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. Reference Jaramillo-Correa, Beaulieu, Khasa and Bousquet2009; Ding et al. Reference Ding, Schreiber, Roberts, Hamann and Brouard2017). It is conceivable that M. disstria followed a similar evolutionary trajectory, dispersing from a single refugium (see also Lait and Hebert Reference Lait and Hebert2018) 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. Reference Bradburd, Coop and Ralph2018). 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 Reference Wright1943). This phenomenon is common in nature (Schwartz and McKelvey Reference Schwartz and McKelvey2008; Meirmans Reference Meirmans2012; Perez et al. Reference Perez, Franco, Bombonato, Bonatelli, Khan and Romeiro-Brito2018), more so if dispersal is limited (Wright Reference Wright1943; Slatkin Reference Slatkin1985). Many analytical approaches used to identify population structure can artificially delimit clusters in continuous genetic data (Perez et al. Reference Perez, Franco, Bombonato, Bonatelli, Khan and Romeiro-Brito2018), particularly with patchy or uneven sampling (Meirmans Reference Meirmans2019). This analytical conundrum has been termed the “cluster versus cline dilemma” (Guillot et al. Reference Guillot, Leblois, Coulon and Frantz2009). Recent work on Populus balsamifera clearly demonstrates this issue. Keller et al. (Reference Keller, Olson, Silim, Schroeder and Tiffin2010) 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. (Reference Meirmans, Godbout, Lamothe, Thompson and Isabel2017) 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 Reference Meirmans2015, Reference Meirmans2019). 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. Reference Goessen, Isabel, Wehenkel, Pavy, Tischenko and Touchette2022). 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. Reference MacDonald, Snape, Roe and Sperling2022). 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. Reference Bagley, Sousa, Niemiller and Linnen2017). Bagley et al. (Reference Bagley, Sousa, Niemiller and Linnen2017) showed that distinct geographic clusters varied in their patterns of host-associated divergence. Given that MacDonald et al. (Reference MacDonald, Snape, Roe and Sperling2022) 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 Reference Parry and Goyer2004), 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 Reference Avise2000). However, there is mounting evidence that mitochondrial variation is not ideal for inferring population-level variation (Ballard and Whitlock Reference Ballard and Whitlock2004; Toews and Brelsford Reference Toews and Brelsford2012; Teske et al. Reference Teske, Golla, Sandoval-Castillo, Emami-Khoyi and van der Lingen2018, Morón-López et al. Reference Morón-López, Vergara, Sato, Gajardo and Ueki2022) and that a shift towards unlinked genome-wide nuclear markers is occurring (Garrick et al. Reference Garrick, Bonatelli, Hyseni, Morales, Pelletier and Perez2015). 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 Reference Lait and Hebert2018); 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 Reference Toews and Brelsford2012) and has been described in a range of taxa (e.g., Wielstra and Arntzen Reference Wielstra and Arntzen2020; Cairns et al. Reference Cairns, Cicchino, Stewart, Austin and Lougheed2021; Morón-López et al. Reference Morón-López, Vergara, Sato, Gajardo and Ueki2022; Yi and Latch Reference Yi and Latch2022). Mitochondrial genomes, due to differences in inheritance, recombination, and effective population size, exhibit distinct responses to selection, demographic changes, and drift (Toews and Brelsford Reference Toews and Brelsford2012; Bonnet et al. Reference Bonnet, Leblois, Rousset and Crochet2017). 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 Reference Linnen and Farrell2007) and selective sweeps (Gompert et al. Reference Gompert, Forister, Fordyce and Nice2008; Wendt et al. Reference Wendt, Kulanek, Varga, Rakosy and Schmitt2022). 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 Reference Gillespie2000) or in species that have experienced rapid range expansions (Petit and Excoffier Reference Petit and Excoffier2009). Similar discordant results were detected in Speyeria aglaja (Lepidoptera: Nymphalidae) (Polic et al. Reference Polic, Yildirim, Lee, Franzen, Mutanen, Vila and Forsman2022), where distinct geographic lineages were detected with genome-wide SNPs but not with CO1 haplotypes. Polic et al. (Reference Polic, Yildirim, Lee, Franzen, Mutanen, Vila and Forsman2022) 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. Reference Bazin, Glemin and Galtier2006; Pavlova et al. Reference Pavlova, Amos, Joseph, Loynes, Austin and Keogh2013), sex-biased dispersal (Folt et al. Reference Folt, Bauder, Spear, Stevenson, Hoffman and Oaks2019), or cytoplasmic incompatibilities (Hurst and Jiggins Reference Hurst and Jiggins2005), 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. Reference James, Cooke, Brunet, Lumley, Sperling and Fortin2015) 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 Reference Lait and Hebert2018), 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 Reference Parry and Goyer2004) and show latitudinal clines in reproductive traits (Parry et al. Reference Parry, Goyer and Lenhard2001). Given that genomic differentiation was linked to both host association and variation in environmental conditions (MacDonald et al. Reference MacDonald, Snape, Roe and Sperling2022), 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.
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.4039/tce.2023.13.
Acknowledgements
The authors gratefully acknowledge all the collaborators and field crews who assisted with specimen collection for this project: Emma Despland, Vanessa Chaimbrone, Lia Fricano, Joshua Jarry, Mike Francis, Ariel Ilic, Kristin Hicks, Chris McVeety, Tyler Nelson, Christi Jaeger, Anne-Sophie Caron, Tyler Wist, and Reshma Jose. They also thank Sophie Dang, Erin Campbell, Tyler Nelson, Alice (Yuehong) Liu, Eric Lemieux, Kevin Ong, and Meng Zhang for assistance with laboratory procedures, bioinformatics, and larval rearing. Funding was provided by Natural Sciences and Engineering Council Discovery Grant RGPIN-2018-04920 to F.A.H.S, with support from Natural Resources Canada through A.D.R, and a NSERC Alexander Graham Bell Scholarship - Doctoral (CGS-D) and UCLA La Kretz Center for California Conservation Postdoctoral Fellowship to Z.G.M.
Open access
