Mountain pine beetle host-range expansion threatens the boreal forest

1Increasing global temperatures have lead to elevational and/or latitudinal shifts in species ranges. This is especially true for insects because they are poikilothermic organisms and therefore quickly respond to changes in their thermal environment. For forest insect pests this is expected to result in negative economic and ecological impacts. Range shifts can expose novel habitats, where naïve hosts may not have evolved appropriate defences to ward off attack. Changes in climate suitability may go beyond the range of the host species, which could mean a dead end for the insect. However hybrid zones between forest species can provide a phenotypic stepping-stone which can help mediate range expansion into a new host.

2The mountain pine beetle (MPB; Dendroctonus ponderosae Hopkins) is a bark beetle indigenous to western North America that primarily feeds on lodgepole pine (Pinus contorta Dougl. ex Loud. var. latifolia), but also feeds on sugar pine (P. lambertiana Dougl.), western white pine (P. monticola Dougl. Ed. D. Don) and ponderosa pine (P. ponderosa, P. Laws. Ex C. Laws; ). Beetle populations typically infest damaged trees or trees with compromised defence capacity; however, given the right conditions they can erupt into large-scale outbreaks and cause significant losses of mature healthy stands. The most recent outbreak has affected over 14 million hectares of forest land in western Canada with considerable losses recorded in the USA. This is the largest outbreak that has been documented since record taking began approximately 125 years ago and brings negative economic impacts, inefficient nutrient cycling and carbon sequestration, and reduced biodiversity.

3Until recently, the range of MPB in Canada has been primarily restricted to the British Columbia interior due to physiological restrictions for MPB. Colder temperatures at higher elevations and increasing latitudes can affect their synchrony and over-winter survival, suppressing population growth at their range limit. However, a suite of studies examining future climate change scenarios have predicted a geographic range shift into previously marginal habitat. As predicted, the beetle has traversed the Rocky Mountains and has reached the eastern edge of the range of lodgepole pine in north-central Alberta where forest composition shifts to jack pine (P. banksiana, Lamb) through a hybrid zone. This hybrid zone could conceivably help mediate MPB host-range expansion into jack pine. Jack pine is a boreal species whose range extends from Alberta to Nova Scotia. There has been no record of MPB infection in natural hybrid or jack pine stands. However, there is considerable evidence to suggest hybrids and jack pine would be suitable hosts for MPB reproduction (reviewed in ). Given the close evolutionary relationship between lodgepole and jack pine, together with the ability of MPB to attack more distantly related pine and instances of MPB attack on jack pine in nursery settings, hybrids and jack pine are likely to be compatible hosts for MPB. Consequently, the arrival of MPB in north-central Alberta has raised major concerns regarding the potential for further expansion to the vast boreal forest.

4Evidence for host-range expansion by MPB to jack pine and hybrids has been difficult since it has not been possible to reliably distinguish lodgepole and jack pine from their hybrids. Seed and cone morphometry are the current criteria for distinguishing species; however, considerable morphological variation in the hybrid zone makes positive identification difficult. Previous efforts to identify hybrids using molecular approaches including chemical profiles, allozymes, organellar DNA, random amplified polymorphism and restriction fragment length polymorphisms have also been ineffective in resolving hybrids from parents in the hybrid zone. Microsatellites have been effectively utilized for a variety of taxa where hybrids have been difficult to distinguish from parentals using other characters. For that reason, our first objective was to develop a panel of microsatellite markers that amplify in both species and test their efficacy to distinguish jack pine, lodgepole pine and their hybrids. We verified our ability to accomplish this by analysing simulated datasets containing multiple levels of admixture generated using pure lodgepole pine (central and southern British Columbia) and jack pine (Saskatchewan and Ontario) as benchmarks. Due to the level of differentiation described between these species (FST = 0.108,; GST = 0.247, ) we expected accurate resolution of first and second generation hybrids, however, we anticipated diminishing power with advanced generations of backcrossing given these species are closely related. Our second objective was to analyse samples collected from MPB attacked and un-attacked trees from the region of MPB expansion including the leading edge in north-central Alberta where the ranges of lodgepole pine and jack pine overlap (Fig. 1). We classified the hybrid and species status of the trees, focusing on MPB-attacked individuals sampled at the leading edge of the expansion based on the information from our simulations and species assignment analysis. We hypothesized that some of the attacked trees in these regions would be jack pine because of the geographic location and the reported susceptibility of jack pine to MPB in nursery settings.

5Foliage samples were collected from lodgepole pine within 10 locations in British Columbia (n= 160), from 14 jack pine locations in Ontario and Minnesota (n = 70) and one location in Saskatchewan (n = 43), and from lodgepole pine, jack pine and putative hybrids within 12 locations in Alberta (n = 381) (Table 1, Fig. 1). In British Columbia and Alberta both MPB attacked and non-attacked trees were sampled. The majority of the foliage samples were collected during two time periods, from February 2007 to May 2007 and September 2007 to April 2008. Stands for sampling were selected on the basis of aerial and ground survey data of MPB attack provided by Alberta Sustainable Resource Development (ASRD). Aerial surveys were conducted in the late summer and early fall each year by ASRD to identify newly attacked trees by the presence of fading foliage, an indicator of MPB infestation. The majority of aerially identified trees were subsequently ground-truthed by ASRD. For each sampling period, waypoints from the most recently available aerial survey were used to locate candidate trees for sampling. Based on 2009 survey data, additional foliar samples from putative jack pine trees were collected in April and May 2010 in the easternmost extent of detected MPB attack by ASRD (n = 6, ‘Smith’–Table 1, Fig. 1) and by the Canadian Forest Service (n = 18, ‘Wildwood’–Table 1, Fig. 1).

6Prior to sampling, each tree was first confirmed for MPB attack and colonization by identification of diagnostic entrance holes in the bark, followed by bark removal around the area of these entrance holes to confirm the presence of MPB larval galleries, as opposed to galleries created by other species such as Ips pini Say, which co-occur in this region. For trees sampled in 2010, the existence of pupal chambers was specifically noted as proof that the host is suitable for completing development of all feeding larval stages. Further development of these pupae had likely not occurred because of the time of year that the samples had been taken. In all cases, foliage was collected from the crown using pole pruners or a shotgun. For samples collected in 2007 and 2008, MPB (typically as larvae) from attacked trees were sampled simultaneously for genetic analysis (Samerasekera et al., unpublished). All sampled trees were geo-referenced using Garmin GPS units (Garmin International, Olathe, KS, USA). Samples were stored in coolers until they could be brought to the laboratory, where they were processed and stored at –20 °C or –80 °C until DNA extraction could be performed.

7Genomic DNA was isolated from ground needle tissue using a CTAB (hexadecyl trimethyl ammonium bromide) protocol modified from, according to. Additional changes include 5 μL of RNase A (70 units/mg protein, Sigma) being added to each sample along with the CTAB buffer. As well, incubation time at 65 °C was increased to 2 h and all centrifugation steps were performed at 5800 g with the duration of the final two centrifugation steps being increased to 5 min. Pellets were resuspended in 125 μL Milli-Q water.

8We screened twenty-three microsatellite loci that were previously tested in lodgepole pine and found 10 to be polymorphic in a panel of eight individuals of each species (Appendix I). Additional loci were isolated from a microsatellite enriched (GTn/CTn) library constructed from a single lodgepole pine individual using the methods of. We selected 368 clones and obtained bi-directional sequences using T3 and T7 primers using BigDye v3.1 sequencing chemistry (Applied Biosystems, Carlsbad, CA, USA) and resolved on a 3730 DNA Analyzer (Applied Biosystems). Sequences were aligned in SeqMan (LaserGene; DNASTAR, Madison, WI, USA) resulting in 222 contigs. Twenty-five contigs were selected that contained long uninterrupted (11–33) repeats with sufficient flanking sequence for primer design. Amplification primers were designed using Primer3 with the default parameters except the optimal Tm was set to 56 °C and the maximum Tm difference between pairs of primers was restricted to 1 °C. Two loci were retained (GENBANK accession numbers: HQ404301– 2) that were polymorphic in both species and had clean amplification products (Appendix I).

9Genotyping was completed for all individuals at 12 microsatellite loci. These loci were amplified in two single, and five multiplex 15 μL reactions (A–G; Appendix I) containing: ∼400 ng DNA (∼200 ng DNA for the single reactions), 1X PCR Buffer, 160 μM each dNTP, 1% dimethyl sulfoxide (volume to volume), 1U Taq DNA polymerase (AB), and optimized MgCl2 and primer amounts (Appendix I). Amplifications were completed using an Eppendorf Mastercycler using the following cycling parameters: 94 °C for 5 min, 33 cycles of 94 °C for 30 s, 56 °C for 30 s, and 72 °C for 15 s, and a final extension at 72 °C for 30 min. These reactions were co-loaded into three injections (1–3; Appendix I) on an ABI 3730 DNA Analyzer and genotyped using GeneMapper software (Applied Biosystems) with allele sizes being determined relative to GeneScan-500LIZ (Applied Biosystems). Genotyping error rate was quantified by running duplicate genotypes for 46 samples.

10We assessed our microsatellite scoring for stutter errors, large-allele drop-out and null alleles using microchecker; any loci with scoring issues were removed from further analyses. Hardy–Weinberg equilibrium and linkage disequilibrium was assessed in genepop 4.0 (; web version, http://genepop.curtin.edu.au/) within locations (where n ≥ 20) for each species (hybrids excluded) as defined by analysis of outputs from Newhybrids and Structure (see Results section). Significance was assessed using Bonferroni corrected alpha values for multiple comparisons (α = 0.05, ). Standard measures of allelic diversity including number of alleles, effective number of alleles (defined as the number of alleles with equal frequency that would achieve the observed level of diversity, ) and observed and expected heterozygosities (HO and HE, respectively) were calculated in GenAlEx 6.0. Allelic richness and private allelic richness corrected for sample size differences were calculated in hp-rare 1.0. All measures were estimated for the entire data-set, for each species (hybrids excluded) once identified, and basic diversity measures were calculated for each location. In addition, we calculated the level of differentiation (FST) at several hierarchical levels: among locations in the entire dataset, between species and among locations within species, using GenAlEx.

11The detection of hybrid classes will depend on the degree of differentiation between the parental species and the loci used. Previous studies that have assessed hybrid zones and the resolving capacity of microsatellites have included only first (F1) and second generation hybrids (F2, F1 backcrosses) in their simulations. We also included a third generation of hybrids (F2 backcrosses and F1 double backcrosses) to assess our ability to resolve advanced introgression. We simulated five datasets using Hybridlab ver. 1.0. Hybridlab was developed to create artificial parental and hybrid genotypes to evaluate the power to correctly identify hybrids. As input, we used genetic profiles from 100 jack pine from Ontario and Saskatchewan and 100 lodgepole pine from British Columbia (we selected samples far removed from the hybrid zone) to represent the microsatellite allele frequency variation for each species. For each dataset we simulated profiles for 300 jack pine and 300 lodgepole pine; we chose 300 as we felt this would simulate our dataset closely. Using these 600 simulated genotypes we generated 100 F1 hybrids. With these hybrid profiles we were able to simulate F2 hybrids, F1 × jack pine (F1Pj), F1 × lodgepole pine (F1Pl), F2 × jack pine (F2Pj), F2 × lodgepole pine (F2Pl), F1-jackpine × jack pine (F1Pj-Pj), and F1-lodgepole pine × lodgepole pine (F1Pl-Pl) (Fig. 3). Each dataset for hybrid analysis was comprised of simulated genotypes for 300 jack pine and 300 lodgepole pine and 10 individuals from each hybrid class: F1, F2, F1Pj, F1Pl, F2Pj, F2Pl, F1Pj-Pj and F1Pl-Pl. We then analysed these datasets with two different Bayesian methods (Structure and Newhybrids) to establish a threshold (QT) for assigning parental and hybrid status. We used both Bayesian methods in our approach to maximize the accuracy of assigning individuals.

12In Newhybrids v. 1.0, we implemented the model with uninformative priors for both the allele frequency and admixture distributions. For each individual Newhybrids calculates the probability of belonging to the parental categories and each specific hybrid classes (e.g. F1, F2, F1 backcross). However, assignment to a specific hybrid class may be uncertain in which case the probability is often split among the hybrid classes. Therefore, to assign hybrids we used two different approaches; if the assignment to a hybrid class was of high probability (i.e. an assignment to the F1 hybrid class with over 0.9 probability) we took the assignment. However, if we found the probability was divided among the hybrid classes, we summed the probability over classes, and if the sum was ≥ 0.9 we considered this a hybrid with class undefined. To assign a QT for the parental categories we looked at the range and average assignment across all individuals. We ran each data-set through ten simulations with a burn-in of 50 000 and 500 000 Markov-chain Monte Carlo sweeps for data collection and we did not use prior population information as species designations were not included for the collected samples.

13In Structure 2.3.1, we set K = 2, and ran the admixture model with uncorrelated allele frequencies, inferring lambda for each population 10 times. The algorithm was run for 550 000 Markov-chain Monte Carlo sweeps with a burn-in of 50 000 and 500 000 for data collection. We used Clumpp to summarize the ten iterations for each of the five simulated data-sets then looked at the Q values for each cluster to determine the most appropriate QT. Again, we looked at the range and average of Q values for each class as we did with Newhybrids. We assessed the standard deviation among individuals across the 10 iterations for each of the five simulated datasets to determine the margin of error for both programs. We also compared the outputs from both programs for each simulated dataset to assess the level of agreement between the methods.

14We analysed our 678 genotyped samples in Newhybrids and Structure using the same parameters and summary methods as our simulations. Our choice of the most appropriate QT for assigning species class was based on the simulations (see Results section). The results from the two programs were compared and a final species class was assigned based on a combination of the two. If the assignment between the two programs did not agree we made a decision based on the Q values from the simulations (see Results section).

15Six-hundred and seventy-eight individuals were used for genotyping (Dryad entry doi:10.5061/dryad.8677), 154 of these represented individuals that had been successfully attacked by MPB. In the putative jack pine sampled in 2010, signs of successful reproduction were observed where pupal chambers were present (Fig. 2), which indicates the eggs hatched and the insects completed all larval stages.

16Individuals were genotyped at 12 microsatellite loci with less than 2% missing data (555 complete profiles, 101 missing one locus and 22 missing two loci). Our genotyping error rate was very low at 0.8%. Some loci were not in Hardy–Weinberg equilibrium at the location level (Appendix II) with a homozygote excess, however only the locus Pcon54 was consistently and extremely out of Hardy–Weinberg equilibrium, therefore we removed this locus from analyses. After Bonferroni correction we did not find any pairs of loci in linkage disequilibrium in lodgepole pine, and only two comparisons for jack pine. Genetic diversity was high among locations, average HO = 0.683 (Table 1). Genetic diversity measured for each species across loci by HO, HE, allelic richness and private allelic richness was higher in lodgepole pine (Appendix III). Across the entire dataset, differentiation among locations (FST = 0.125) and between species (FST = 0.133) was high. Within species however, differentiation among locations was very low (FSTlodgepole = 0.033, FSTjack = 0.016). However, there were some significant pair-wise comparisons; most notably for jack pine, where all locations were differentiated from Ontario/Minnesota.

17Across the five simulated datasets the assignment of individuals to their correct class was never less than 96% using either assignment method (Fig. 4, Table 2). The calculated proportion of ancestry for both programs were highly consistent, the standard deviation among individuals across the ten iterations for each simulated dataset were extremely low (Structure = 0.00014, Newhybrids = 0.00003). As well, the consistency in assignment across the programs was also very high with few discrepancies (Table 3). For both methods detection of 1st generation hybrids was 100%, for 2nd generation hybrids this decreased to 87% and 81% and for 3rd generation hybrids this decreased further to 68% and 63%, in Newhybrids and Structure, respectively (Table 2). Based on the results from Newhybrids the most accurate method to assign hybrids was to sum the estimated proportions across the hybrid categories. When we used a QT of 0.9 for the parental and hybrid classes in newhybrids, similar to other studies, we found a large number of individuals that did not assign to any category. However, over 98% of these individuals with QT < 0.90 were hybrids, therefore we used a QT ≥ 0.9 to assign pure species, and all other individuals were assigned hybrid status.

18Using Structure, we found that QT ≥ 0.9 gave the most accurate pure assignment for simulated data, but cases of advanced introgression were difficult to detect (Fig. 4). The average disagreement between Structure and Newhybrids across the five simulations was 1.4%. All of these discrepancies were either hybrids correctly assigned by Newhybrids, or pure individuals correctly assigned by Structure. Based on these discrepancies we developed a set of rules for assigning individuals when the methods had conflicting results: (i) if the probability of being a hybrid in Newhybrids is ≥ 0.9 and the probability of being a parental in Structure is > 0.9 but < 0.95 the individual was assigned hybrid status. (ii) If the probability is most likely to be a pure parental in Newhybrids but < 0.9 and in Structure is > 0.95, then the individual is assigned to the parental category.

19We were able to clearly delineate the two species and our power to detect hybrid individuals was 0.74 averaged across three generations of hybrids in the simulated data (Table 2). We therefore used a QT ≥ 0.9 in both Structure and Newhybrids for the assignment of pure lodgepole pine and jack pine, and considered the unassigned individuals as hybrids. The two assignment methods agreed for 668 of the 678 samples analysed (Table 3). The 10 disagreements were resolved using the rules that we developed based on the simulations (see above). The final breakdown of assignment for our sample data was 87 hybrids, 301 jack pine and 290 lodgepole pine (Figs 5 and 6). Ancestry of the hybrid trees was predominately lodgepole pine (Fig. 5). Of the trees sampled that were designated as attacked there were eight jack pine (five from the samples from Smith, and three from Wildwood), 127 lodgepole pine, and 19 hybrids. The eight trees assigned as jack pine had Q values > 0.99.

20Using simulated microsatellite genotypes and Bayesian cluster analyses, we were able to test the power of a set of microsatellites to distinguish lodgepole pine, jack pine and variable levels of hybrid ancestry. We found high agreement between the two Bayesian algorithms and low assignment error in the analysis of simulated data. Where discrepancies arose between the two methods, we generally found that Newhybrids had better power to detect hybrids while Structure performed better with pure individuals, highlighting the complementarily of the methods. While our microsatellite panel had excellent resolving power for parentals and recent hybrids, resolution declined with increasing generations of hybrid back-crossing. We expected decreased resolution based on the close evolutionary relationship between lodgepole pine and jack pine, however this may not be a major issue for our dataset. Given the colonization times for lodgepole pine (7000 YBP;;; ) and jack pine (6000 YBP;; ) in central Alberta, the long generation times for pine and low hybrid vigour, the potential for geographically spread, advanced introgression is not likely.

21Our ability to distinguish species classes enabled us to address the second objective of our study, namely to determine whether jack pine have been successfully attacked by MPB. We identified eight pure jack pine trees collected from MPB attacked stands at the edge of MPB range expansion in Alberta (Figs 1 and 2). In addition, we also identified 19 hybrid trees that have also been attacked confirming their susceptibility within the hybrid zone.

22The discovery that MPB has expanded its host repertoire to include jack pine has prompted us to consider whether MPB will be able to sustain eruptive populations locally and thus spread further eastward into the boreal forest. Host and beetle interaction is influenced by physiology, population dynamics and environment making this a complex system where many biotic and abiotic factors need to be considered concomitantly. This study demonstrates that both hybrids and pure jack pine are susceptible to MPB within the hybrid zone. It is currently unknown whether hybrids and jack pine have different susceptibilities to MPB attack relative to lodgepole pine. Differential susceptibility of jack pine and hybrids to MPB is plausible. For example these two species have different susceptibilities to western gall rust fungus. Also, lodgepole pine and MPB appear to share a long co-evolutionary history that has presumably allowed this tree species to adapt at some level to MPB. In fact, recently showed that naïve lodgepole pine stands had higher MPB reproductive output than stands that experienced epidemic outbreaks. This suggests that MPB may have higher reproductive success in hybrid and jack pine trees which would further sustain the epidemic and support its eastward expansion. There is a potential for introgression of lodgepole pine genes – including genes that condition defence – into the jack pine genome through historical hybridization. However, widespread hybridization and introgression of defence genes seems unlikely based on our findings, therefore any potential benefits of introgression would likely be limited to the hybrid zone. As well, those trees identified as hybrids in our study exhibited higher lodgepole than jack pine ancestry (Fig. 5) suggesting a higher percentage of lodgepole pine backcrosses and limited introgression of lodgepole genes into the jack pine range. Further, chemical defences produced by lodgepole pine differ in their composition from jack pine and hybrids. For example, found α-pinene, a chemical that may facilitate a successful mass attack by MPB, to be at considerably higher concentrations in jack pine than in lodgepole pine. and found that measures of beetle performance were similar between lodgepole pine and jack pine, suggesting the potential of jack pine to sustain populations in a manner similar to lodgepole pine. As well, there are MPB fungal associates that are part of the beetle invasion process, and it has been shown that hybrids and jack pine are susceptible. It should be noted that many of the published experiments were carried out using cut bolts and artificially infested rather than naturally infested live trees, limiting extrapolation to the natural boreal forest.

23While tree-level defences that are partly under genetic control will conceivably contribute to the probability of MPB eruptive population dynamics in jack pine, there also needs to be a sufficient density of available hosts. The severity of the recent MPB outbreak has been partly attributed to the present-day distribution of lodgepole pine, where the continuity of even-aged stands resulting from forest management practices has been ideal for maintaining beetle populations. In contrast, jack pine is not uniformly distributed in the boreal forest, typically occurring in a patchy distribution. and have analysed the potential for MPB spread in the boreal forest region and have suggested that Alberta populations are susceptible; however, jack pine occurrence may be too fragmented east of Alberta to sustain the type of epidemic that British Columbia has experienced.

24Climate has influenced the severity of the recent MPB outbreak, and will likely play an important role in determining whether this epidemic will be maintained in the boreal forest. Temperature has received the most attention and from this there are two considerations: (i) in the past, the summer climate in northern Alberta has not been suitable to sustain synchronous beetle populations which are necessary to maintain eruptive populations and (ii) beetles are not completely cold tolerant and can incur high mortality within the jack pine range from cold exposure. However, changes in global temperatures have improved the climatic suitability for both summer and winter seasons promoting conditions for eruptive behaviour.

25The discovery of successful MPB attack and evidence of completed larval development in jack pine is a critical first step in assessing future impacts of this destructive forest pest on the boreal forest and how climate change may affect the system. We have considered whether MPB populations can be maintained in the boreal forest and there are many factors that need to be met for continued MPB expansion and population growth. If jack pine can sustain endemic populations and thus maintain this host-range expansion it is critical that forest management incorporates these considerations in their future planning. MPB is not endemic to the boreal forest and therefore should be considered an invasive species and managed as such. Forest ecosystems in North America have already been challenged with numerous pest invasions that represent a considerable threat. When we factor in climate change, the vulnerability of ecosystems such as the boreal forest to disturbance is further increased putting an extremely important ecosystem in jeopardy.