Introduction
The mountain pine beetle, Dendroctonus ponderosae Hopkins (Coleoptera: Curculionidae), is a tree-killing bark beetle that has become one of the most significant destructive forest pests in western North America (Natural Resources Canada 2017). Historically, mountain pine beetle populations in Canada have been confined largely to British Columbia and southwestern Alberta within the range of their primary host, lodgepole pine (Pinaceae) (Safranyik et al. Reference Safranyik, Carroll, Régnière, Langor, Riel and Shore2010). However, in the early 2000s, populations expanded and established in northwestern Alberta (Natural Resources Canada 2017), aided by wind updrafts (Jackson et al. Reference Jackson, Straussfogel, Lindgren, Mitchell and Murphy2008; Janes et al. Reference Janes, Li, Keeling, Yuen, Boone and Cooke2014). The species’ eastward spread has continued, resulting in more than $Cdn 500 million being spent on mitigating damage in Alberta alone (Whittaker Reference Whittaker2018). Northeastern Alberta is of particular concern because this is where lodgepole pine transitions to jack pine, a naïve host species that provides a potential entry point into the boreal pine forests of eastern Canada and the northern United States of America (Cullingham et al. Reference Cullingham, Cooke, Dang, Cooke and Coltman2011; Bleiker et al. Reference Bleiker, O’Brien, Smith and Carroll2014).
Various genetic markers have been used to map mountain pine beetle population structure and range expansions in British Columbia (Mock et al. Reference Mock, Bentz, O’Neill, Chong, Orwin and Pfrender2007) and Alberta (Langor and Spence Reference Langor and Spence1991; Cullingham et al. Reference Cullingham, Roe, Sperling and Coltman2012; Samarasekera et al. Reference Samarasekera, Bartell, Lindgren, Cooke, Davis and James2012; Janes et al. Reference Janes, Li, Keeling, Yuen, Boone and Cooke2014, Reference Janes, Roe, Rice, Gorrell, Coltman, Langor and Sperling2016, Reference Janes, Worth, Batista and Sperling2018; Batista et al. Reference Batista, Janes, Boone, Murray and Sperling2016). These studies have shown that mountain pine beetle spread into Alberta was composed of two incursions that are genetically distinct (Samarasekera et al. Reference Samarasekera, Bartell, Lindgren, Cooke, Davis and James2012; Janes et al. Reference Janes, Li, Keeling, Yuen, Boone and Cooke2014; Batista et al. Reference Batista, Janes, Boone, Murray and Sperling2016) and a third incursion that is genetically intermediate (Trevoy et al. Reference Trevoy, Janes and Sperling2018). The southern resident population has ranged as far north as Banff and Canmore (Hopping and Mathers Reference Hopping and Mathers1945), and a more recently established northern population extends eastwards past Grande Prairie (Janes et al. Reference Janes, Li, Keeling, Yuen, Boone and Cooke2014). The intermediate population recently entered Alberta by way of Jasper National Park, and its genetic similarity to lab-produced hybrids between the northern and southern populations suggests that it resulted from admixture between them (Trevoy et al. Reference Trevoy, Janes and Sperling2018).
Epidemic populations have been identified in Hinton, west–central Alberta, and Slave Lake, indicating continued northeastward range expansion in Alberta. Mountain pine beetle has now also been detected at low densities on the eastern edge of Alberta near the border with Saskatchewan, Canada in novel jack pine habitat (MacCormick Reference MacCormick2020). It is uncertain if beetles detected near Saskatchewan came from low-density, locally established populations from previous long-range dispersal events or directly from larger outbreak populations to the west. These recently detected populations have not yet been included in demographic studies on mountain pine beetle spread. Determining their sources would allow optimisation of pheromone monitoring (Miller et al. Reference Miller, Lindgren and Borden2005), guide forest management practices to limit localised population increases (Safranyik and Carroll Reference Safranyik and Carroll2006), and better inform assessments of outbreak risk (Bleiker Reference Bleiker2019). In the present study, we use genome-wide single nucleotide polymorphisms to characterise population structure and identify dispersal and spread of mountain pine beetle in northern Alberta, particularly in newly detected northeastern populations.
Methods
Beetle sampling, DNA extraction, and sequencing
A total of 306 specimens were collected at 44 different sites from nine locations between 2014 and 2018, and specimens were assigned to populations based on proximity to other collection localities (Supplemental materials, Table S1). Beetles were preferentially taken from different trees or egg galleries to reduce the influence of genetic structuring due to the inclusion of siblings from single families. Once collected, beetles were preserved in 90–95% ethanol and placed in a freezer at −20 °C until DNA extraction.
DNA extraction used the Qiagen DNEasy Blood & Tissue kit (Qiagen, Germantown, Maryland, United States of America) with an optional RNAse A treatment. Library preparation of samples collected in 2014–2015 followed the two-enzyme genotyping-by-sequencing protocol of Poland et al. (Reference Poland, Brown, Sorrells and Jannink2012), and the library preparation of samples collected in 2016–2018 followed double digest restriction-site associated DNA sequencing protocol of Peterson et al. (Reference Peterson, Weber, Kay, Fisher and Hoekstra2012), according to the general procedure outlined in Campbell et al. (Reference Campbell, Davis, Dupuis, Muirhead and Sperling2017). Single-end sequencing of the 2014–2015 samples was performed on an Illumina HiSeq 2000 (Illumina, Inc., San Diego, California, United States of America) at the Institut de biologie intégrative et des systèmes (IBIS) at Université Laval, Quebec City, Québec, Canada, whereas single-end sequencing of the 2016–2018 samples was performed on an Illumina NextSeq 500 (Illumina, Inc.) at the University of Alberta Molecular Biology Service Unit, Edmonton, Alberta, Canada.
Data filtering and single nucleotide polymorphism identification
Initial data demultiplexing was conducted in Stacks, version 2.0b (Rochette et al. Reference Rochette, Rivera-Colón and Catchen2019) on the Graham cluster of Compute Canada (Toronto, Ontario, Canada). The PstI restriction sites on the 5′ end of each DNA sequence that was read were trimmed using cutadapt, version 1.9.1 (Martin Reference Martin2011). Samples sequenced on a NextSeq 500 had a final read length of 67 base pairs after trimming. Following Campbell et al. (Reference Campbell, Davis, Dupuis, Muirhead and Sperling2017), longer reads produced on a HiSeq 2000 (Illumina, Inc.) were additionally truncated on the 3′ end to match the length of the NextSeq-generated data. The trimmed sequence reads were then aligned to the female mountain pine beetle draft genome (Keeling et al. Reference Keeling, Yuen, Liao, Docking, Chan and Taylor2013) using the alignment algorithm Burrows–Wheeler Aligner Maximal Exact Match (BWA-MEM), version 0.7.17 (Li and Durbin Reference Li and Durbin2009). Alignment quality was checked using SAMtools, version 1.9 (Li et al. Reference Li, Handsaker, Wysoker, Fennell, Ruan and Homer2009).
Single nucleotide polymorphisms were identified using the ref_map.pl script in Stacks (Catchen et al. Reference Catchen, Amores, Hohenlohe, Cresko and Postlethwait2011; Rochette et al. Reference Rochette, Rivera-Colón and Catchen2019). We specified a minimum minor allele frequency of 0.01 and adhered to the r80 principle recommended by Paris et al. (Reference Paris, Stevens and Catchen2017), which requires a locus to be present in at least 80% of the individuals in a population in order to be retained in the final data set. At each locus, only a single, random single nucleotide polymorphism was retained. Further filtering of single nucleotide polymorphism data and samples was conducted in VCFtools, version 0.1.14 (Danecek et al. Reference Danecek, Auton, Abecasis, Albers, Banks and DePristo2011). Specifically, loci were excluded if they had a genotype quality score below 30 or if a single nucleotide polymorphism locus had more than 5% missing data globally. Initial exploratory principal component analysis recovered clustering along the first principal component that appeared to be driven by sex (consistent with Trevoy et al. Reference Trevoy, Janes, Muirhead and Sperling2019), which largely confounded population-level clustering (Supplemental materials, Fig. S2). To remove this signal, we filtered for linkage disequilibrium, following Abdellaoui et al. (Reference Abdellaoui, Hottenga, de Knijff, Nivard, Xiao and Scheet2013). Linkage disequilibrium was calculated for all pairwise single nucleotide polymorphism combinations using the dartR package (Gruber and Georges Reference Gruber and Georges2019) in R (R Core Team 2018). For any groups of single nucleotide polymorphisms that were in linkage disequilibrium (R 2 > 0.5), we retained only one random single nucleotide polymorphism from each to reduce the impact on population-clustering analyses. Finally, we did not distinguish between putatively adaptive and neutral loci for population-clustering analyses because a combination of both has been shown to more accurately reflect population structure in mountain pine beetle (Batista et al. Reference Batista, Janes, Boone, Murray and Sperling2016).
Principal component analysis groups specimens based on the extent of covariance in their associated data, and often individuals that are closely related can confound population structure due to their high rate of covariance (Alhusain and Hafez Reference Alhusain and Hafez2018). Although our sampling procedure was designed to reduce the incidence of multiple specimens from a single family, data visualisation using principal component analysis indicated tight clustering patterns for two groups, consisting of a total of seven individuals, that were consistent with family-level structure (Price et al. Reference Price, Zaitlen, Reich and Patterson2010; Supplemental materials, Fig. S1). Three clustered individuals came from a single tree (sibling group 1; Supplemental materials, Fig. S1), and four clustered individuals came from two close trees (sibling group 2; Supplemental materials, Fig. S1). No patterns in missing data or poor locus coverage were identified that could explain this structuring. Because female mountain pine beetles are capable of producing multiple egg galleries from a single mating (Reid Reference Reid1958), these clustered individuals were assumed to be siblings. To reduce the influence of such family structure on population-level analyses, all but one individual in each set of closely related specimens were removed from the data set.
Population structure analysis
To explore and identify potential signs of population structure across northern Alberta, we used several approaches. A principal component analysis on all samples that passed filtering was performed in the R package adegenet (Jombart Reference Jombart2008; Jombart and Ahmed Reference Jombart and Ahmed2011) and plotted with ggplot2 (Wickham Reference Wickham2016). This was used to detect single nucleotide polymorphism covariance between specimens and to identify population structure among the specimens without a priori population assignment. We also used this set of samples to calculate observed heterozygosity, expected heterozygosity, and pairwise population differentiation due to genetic structure (F ST) in the R package dartR, and the R package Genepop (Raymond and Rousset Reference Raymond and Rousset1995; Rousset Reference Rousset2008) was used to calculate population differentiation based on genetic diversity using an exact G-test of genic differentiation with a Bonferroni correction based on the 36 pairwise tests.
To assess fine-scale population structure across the northernmost locations, we created a northern data set by removing all southern- and intermediate-cluster population samples, including any putative dispersers. We performed a principal component analysis on this northern data set, then used discriminant analysis of principal components (Jombart et al. Reference Jombart, Devillard and Balloux2010) in adegenet to detect potential weak signals of population structure using a priori population assignments based on collection locality and verified by the prior principal component analysis. We then performed isolation-by-distance analysis using the R packages Geosphere (Hijmans Reference Hijmans2019) and SNA (Butts Reference Butts2016) and calculated a Mantel test in adegenet.
The program Structure, version 2.3.4 (Pritchard et al. Reference Pritchard, Stephens and Donnelly2000; Falush et al. Reference Falush, Stephens and Pritchard2003, Reference Falush, Stephens and Pritchard2007; Hubisz et al. Reference Hubisz, Falush, Stephens and Pritchard2009), was used to infer genetic population structure among all populations. We subsampled the full data set to 125 individuals that were arbitrarily selected to represent potential sample-site variation (10 to 16 individuals from each location) in order to reduce bias due to differences in sampled population sizes (Puechmaille Reference Puechmaille2016). In this data set, we attempted to include as many samples as possible for each population while maintaining a roughly equal sample size among the many sampling locations. We used the admixture model without specifying a priori population or location information, and we tested values of K from 1 to 10 with one million Markov chain Monte Carlo generations and a burn-in period of 250 000. We replicated the analysis for each value of K a total of 10 times and then used CLUMPAK (Kopelman et al. Reference Kopelman, Mayzel, Jakobsson, Rosenberg and Mayrose2015) to average the 10 independent replicates for each value of K. The optimal value of K was assessed using both the difference between the likelihoods of subsequent K values (ΔK) (Evanno et al. Reference Evanno, Regnaut and Goudet2005) and the mean estimated natural logarithm of the probability of each K value (Ln(PrK)) (Pritchard et al. Reference Pritchard, Stephens and Donnelly2000). We then ran a principal component analysis using this subset of data to ensure our Structure results were comparable to the principal component analysis of the full data set described above.
Results
After filtering, 299 individuals and 2872 genomic single nucleotide polymorphisms were retained for data analysis (Supplemental materials, Table S1). Average individual read depth was 58.3 (minimum depth = 12.6, maximum depth = 217.3; Supplemental materials, Table S1). Pairwise F ST between populations ranged from approximately 0 to 0.13 (Table 1), with Canmore found to be the most distinct population. Tests for differences in pairwise population genetic differentiation also identified Canmore as having allelic distributions that significantly differed from those of all other tested populations and found some significant differences between intermediate and northern locations (Table 1). Observed heterozygosity was similar among populations, ranging from 0.14 to 0.19, and differences between expected and observed heterozygosity were minimal (Table 1). Some of the most recently detected populations, such as Battle Lake and Lac La Biche, had higher heterozygosity than some older, established populations, such as Grande Prairie and Edson, but no clear trend was evident (Table 1).
Table 1. Pairwise population differentiation P-values (above diagonal), F ST (below diagonal), and observed (Ho) and expected (He) population heterozygosity.
Statistics were calculated using n individuals for each population group. Significantly differentiated population pairs (in bold) had genic differentiation values of P < 0.001 after Bonferroni correction.
In the principal component analysis for the full data set of 299 individuals (Fig. 1A), beetle sex had strong signal that overwhelmed population structure (Supplemental materials, Fig. S2). All but three specimens could be sexed using methods outlined in Trevoy et al. (Reference Trevoy, Janes, Muirhead and Sperling2019). For the three ambiguous specimens, we removed all loci containing missing data from sex-related scaffolds, but they remained intermediate between the male and female clusters (Supplemental materials, Fig. S2). Despite this, after filtering for linkage disequilibrium, these specimens clustered as expected, based on sampling location. To focus on geographic population structure, sex-based signal was successfully removed from the data set by filtering for linkage disequilibrium.
Fig. 1. Principal component analyses results. A, Results for all sampling locations and individuals (n = 299) and B, results for northern locations (Edson, Grande Prairie, Lac La Biche, Slave Lake, and Whitecourt), with all intermediate-cluster individuals removed (n = 201). In A, ellipses were added to emphasise the south, intermediate, and north clusters. These genetic population clusters were divided at PC1 = −5.5 and PC1 = −0.06 for assignment of individuals in Supplemental materials, Table S1.
After filtering, the first principal component (PC1), which often relates to geographic location in genetic assessments of population structure (see Abdellaoui et al. Reference Abdellaoui, Hottenga, de Knijff, Nivard, Xiao and Scheet2013), explained 3.6% of the variation in the single nucleotide polymorphism data and distinguished northern and southern populations as well as a distinct intermediate group between the two; for this study, we used PC1 = −5.5 and −0.06 to delimit the intermediate cluster. The second principal component explained 1.0% of the remaining variance and related to variation within, rather than between, populations.
No population structure was detected in either principal component analysis or discriminant analysis of principal components when we focussed on a northern data set that included 206 specimens, after removing individuals from the intermediate and southern clusters (Fig. 1B; Supplemental materials, Fig. S3). Isolation-by-distance analysis of the northern locations likewise did not indicate any patterns of population differentiation based on geographic distances (Mantel test: P = 0.75).
The ΔK values from the Structure analysis indicated an optimal K value of 2, and Ln(PrK) indicated an optimal K of 2 or 3 (Supplemental materials, Fig. S4). A K value of 2 distinguished northern and southern population clusters and identified a genetically and geographically intermediate population of specimens from Jasper, Battle Lake, and Hinton (Fig. 2); K values of three and higher neither further resolved this intermediate population nor indicated any other meaningful geographic clustering. Thus, K = 2 was considered to be optimal (Supplemental materials, Fig. S4). Structure analysis also indicated some individuals that did not cluster according to the locations in which they were collected; one individual that was collected near Edson and one that was collected near Slave Lake grouped with the intermediate population. Principle component analysis using this subset of 125 samples also produced similar divisions to the Structure analysis and closely resembled the patterns shown in Fig. 1A, indicating that the results of the Structure analysis were consistent with the population structure identified through principal component analysis of the full data set.
Fig. 2. Structure analysis for equivalently sized population subsets, totalling 125 individuals. K = 2 was the best supported number of populations. Based on individual assignments, the Edson and Slave Lake populations each include a single admixed individual with a roughly equal assignment to both the northern and southern genetic clusters.
Based on the principal component analysis of the full data set (n = 299), we identified several additional individuals that were collected in the northern populations but grouped with the intermediate population cluster, bringing the totals to three beetles from Slave Lake, five from Whitecourt, and one from Edson (Supplemental materials, Table S1 and Fig. S5). By chance, none of these Whitecourt individuals were sampled in the subset Structure analysis of 125 individuals; however, a larger Structure analysis of all 299 specimens indicated that these individuals were admixed (Supplemental materials, Fig. S4) and was largely consistent with the results of the principal component analysis presented in Fig. 1A. We also identified 11 specimens collected in Hinton and one collected in Jasper that grouped with the northern cluster. Of these putative dispersers, 81% were female and 19% were male (Supplemental materials, Table S1); the sex ratio of these geographically mismatched individuals is consistent with known sex biases in mountain pine beetle (Blomquist et al. Reference Blomquist, Figueroa-Teran, Aw, Song, Gorzalski and Abbott2010; James et al. Reference James, Janes, Roe and Cooke2016; Fig. 3). No putative dispersers from the intermediate population were found among the newly established Lac La Biche specimens.
Fig. 3. Population expansion dynamics inferred from genetic evidence based on 299 mountain pine beetle specimens collected from sites across Alberta, Canada. Pie charts for each sampling location are coloured by proportion of individuals assigned to the north (blue), intermediate (brown), and south (grey) population clusters, as determined using principal component analysis (Fig. 1A).
Discussion
Our study used population structure to investigate the dispersal dynamics of mountain pine beetle in central Alberta, identifying primary source populations as well as fine-scale dispersal from other regions. We found the same distinct north–south beetle genetic population structure that was previously described for Alberta (Janes et al. Reference Janes, Li, Keeling, Yuen, Boone and Cooke2014; Batista et al. Reference Batista, Janes, Boone, Murray and Sperling2016; Trevoy et al. Reference Trevoy, Janes and Sperling2018; Figs. 1 and 2). When we focussed on individuals from northern locations, however, we found no evidence of population structure (Fig. 1B). This finding is consistent with a high rate of dispersal and gene flow between these regions. The long-range dispersal capability of mountain pine beetle (Jackson et al. Reference Jackson, Straussfogel, Lindgren, Mitchell and Murphy2008) likely promotes genetic admixture, particularly from west to east, and may have inhibited the further development of genetically distinct populations in this region.
The extensive dispersal of mountain pine beetle also applies to mixing of population clusters. Although we saw no admixture involving the southern cluster, mingling was evident among the north and intermediate populations. We identified several putative dispersers with affinities to Jasper and Hinton in sites across the northern locations, including as far to the northeast as Slave Lake (Fig. 3; Supplemental materials, Table S1). This illustrates the considerable dispersal capability of these populations and their potential contribution to range expansion in the north. The majority of putative dispersers that we identified were female, which could potentially increase establishment success rate because female beetles can be mated before they disperse (Bleiker et al. Reference Bleiker, Heron, Braithwaite and Smith2013) and are the pioneering sex in D. ponderosae (Blomquist et al. Reference Blomquist, Figueroa-Teran, Aw, Song, Gorzalski and Abbott2010). Specimens collected in Hinton and one specimen at Jasper showed the reverse pattern, with some beetles having affinities to the northern cluster, indicating that several individuals likely came from a northern source (Fig. 3). These results highlight the importance of considering multiple source populations when determining infestation origins in newly established regions.
Admixture and high dispersal rates from the north and intermediate populations towards the northeast may have important genetic ramifications. In particular, dispersal from the genetically intermediate Jasper and Hinton region is likely to increase genetic diversity in the currently expanding eastern range of mountain pine beetle, which would otherwise be solely composed of the genetically homogenous northern population. Although this introduces the possibility of potential benefits from increased genetic diversity in the leading-edge populations, high rates of dispersal and gene flow may also prevent local adaptation (Lenormand Reference Lenormand2002). Regardless, genetic differences between regions allow us to monitor new populations and to identify the magnitude and composition of the spread of mountain pine beetle towards the northeast from different genetic sources.
Currently, the leading edge of the mountain pine beetle range is found in jack pine forests near Lac La Biche (MacCormick Reference MacCormick2020), and this population is genetically indistinguishable from other northern populations as far west as Grande Prairie. This uniformity can be explained by high rates of dispersal and regular waves of reinforcing conspecifics from other populations. In another forest pest, the spruce budworm (Choristoneura fumiferana Clemens) (Lepidoptera: Tortricidae), reinforcing waves of dispersers from neighbouring populations may allow endemic populations to become epidemic (Larroque et al. Reference Larroque, Johns, Canape, Morin and James2020). Although their life histories differ, mountain pine beetle and spruce budworm have similar overall outbreak dynamics, with large increases in population size causing epidemic conditions (Royama Reference Royama1984; Larroque et al. Reference Larroque, Johns, Canape, Morin and James2020). In jack pine forests, mountain pine beetle can face intense competition from heterospecifics, such as wood-boring beetles (Klutsch et al. Reference Klutsch, Najar, Cale and Erbilgin2016). This may serve to mediate local population growth, but populations near Slave Lake likely make contributions to the Lac La Biche infestation. Continued recruitment in this region risks establishing large, potentially epidemic populations in the eastern jack pine forests.
Mountain pine beetle dispersal has been difficult to predict, but our study has helped to identify large-scale dispersal patterns in central Alberta. We demonstrate that the current Hinton infestation was derived from two sources, these likely being Jasper and, to a lesser degree, Grande Prairie. Although we found no evidence of any dispersal from the southern population cluster, there was a general trend of population expansion from the Jasper region towards the northeast, with occasional dispersal as far northeast as Slave Lake. To date, control efforts in Alberta have been largely successful at slowing the spread of mountain pine beetle, but the beetle nonetheless remains capable of extensive and rapid range expansion. Long-range dispersal may limit population structure within the northernmost regions, but it also may contribute to epidemic population development. Continued focus on reducing current outbreaking populations may be necessary to minimise long-range dispersal and reinforcement of populations at the beetle’s leading edge in Alberta.
Acknowledgements
This research was supported by funding awarded to Felix Sperling from the Natural Science and Engineering Research Council of Canada (grant no. NET GP 434810-12) to the TRIA Network, with contributions from Alberta Agriculture and Forestry, fRI Research, Manitoba Conservation and Water Stewardship, Natural Resources Canada – Canadian Forest Service, Northwest Territories Environment and Natural Resources, Ontario Ministry of Natural Resources and Forestry, Saskatchewan Ministry of Environment, West Fraser, and Weyerhaeuser. An NSERC Discovery Grant to F. Sperling (RGPIN-2018-04920) supported V. Shegelski during preparation of this manuscript. This research was also funded by the Government of Alberta. The authors also thank Julian Dupuis, Stephane Bordeleau, Jackson Lai, Dylan Sjolie, Phil Batista, Sebastian Lackey, Devin Letourneau, Melodie Kunegel-Lion, John Haley, Andrew Sperling, William Sperling, Janet Sperling, Stephen Trevoy, Neil Thompson, Jennifer McCormmack, and Fraser McKee for their help with sample collection, logistics, and advice on statistical analyses. This research was enabled in part by support provided by Westgrid (www.westgrid.ca) and Compute Canada (www.computecanada.ca).
Data availability
Raw sequence data in the form of fastq files are available in the National Center for Biotechnology Information Sequence Read Archive (NCBI SRA) under the accession number PRJNA656769 (https://www.ncbi.nlm.nih.gov/sra/PRJNA656769).
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.4039/tce.2020.83.
