Is Resistance to Mountain Pine Beetle Associated with Genetic Resistance to White Pine Blister Rust in Limber Pine?

Limber pine (Pinus flexilis James) co-evolved with the mountain pine beetle (Dendroctonus ponderosae Hopkins; MPB) and is now also challenged by the non-native pathogen Cronartium ribicola (J.C. Fisch.) that causes the lethal disease white pine blister rust (WPBR). Previous research suggests that trees infected with WPBR can be preferred hosts for MPB. Using resin duct traits associated with MPB resistance, we tested for a relationship between resistance to MPB and WPBR in limber pine, in the absence of either biological agent. These analyses will help evaluate if MPB historically may have contributed to natural selection for WPBR resistance in advance of WPBR invasion, and could help explain the unusually high frequency of the dominant Cr4 allele for complete resistance to WPBR in limber pine populations of the Southern Rocky Mountains. Resin duct production, density and relative duct area did not differ between healthy trees previously inferred to carry the dominant Cr4 allele and trees that lack it at 22 sites, though some duct traits varied with elevation. MPB resistance does not appear to have played an evolutionary role in contributing to the high frequency of Cr4 in naïve populations, however, MPB may affect the future evolution of resistance to WPBR in the pines where the two pests coincide and WPBR will affect forest recovery after MPB epidemics. MPB-WPBR interactions in a changing climate will affect the future trajectory of limber pine.


Introduction
The southern Rocky Mountains are being challenged by both native and invasive pests and pathogens with increasing intensity as temperatures continue to rise and drought increases [1]. As a native pest, the mountain pine beetle (Dendroctonus ponderosae Hopkins, MPB) has evolved alongside its hosts colonizing all native pines within its range except Jeffrey pine (Pinus jeffreyi Balf) [2,3]. With suitable host species to the north and south of the current MPB geographic range, there is speculation that climate (minimum temperature) has been a main determinant in preventing range expansion [4]. Within the MPB distribution in the western United States, the coldest night of the year has warmed approximately 4 • C since 1960 [5,6]. Since air temperatures have been increasingly more favorable for MPB overwinter survival, recent outbreaks are putting high elevation pines such as limber pine and whitebark pine at risk of population decline [7][8][9][10][11][12]. Limber pine is a preferred host of MPB for reproduction, development, and survival [13]. Not only does limber pine have a thicker phloem which is desirable in a suitable host, but limber pine phloem may also have more nutritional value than other species [14][15][16][17]. As a native species, MPB is an important driver of ecosystem dynamics and the recent outbreak was more severe than any documented over the past century [9,17,18].
Through the coevolution of pest and host, certain pine species colonized by MPB have developed resistance mechanisms and resin duct traits are proving to be of more influence on inhibiting herbivores than previously thought [19][20][21]. Resin can be both toxic to attacking beetles and can physically encumber them [22], and trees with higher resin duct numbers are considered to be more highly defended [21,23]. Similarly, constitutive resin flow is a first line of defense against attacks in pines, and the best predictor of resin flow is resin duct number, size and total duct area [24]. In limber pine, MPB resistant trees (those trees that survived recent mass attack) were found to have an average of 37% more resin ducts in their most recent 5 years of radial growth than susceptible trees (those trees that were killed by recent mass attack) [21].
The non-native pathogen, Cronartium ribicola (J.C. Fisch), which causes the lethal disease white pine blister rust (WPBR), was introduced to western North America over a century ago, and all nine North American five-needle white pine species are highly susceptible [25]. The spreading WPBR infection front is now in the Southern Rocky Mountains on limber pine (Pinus flexilis James) and Rocky Mountain bristlecone pine (Pinus aristata Engelm.) with some infections further south in New Mexico and Arizona on southwestern white pine (Pinus strobiformis Engelm.) [26][27][28][29]. Limber pine has been listed as endangered in Alberta Canada [30]; conservation strategies are in place for the Southern Rockies [31,32], and populations continue to decline due to this invasive fungus, MPB and a changing climate [33][34][35].
Introduced pests and pathogens can have devastating effects on native tree species. Some host species have genetic resistance to the invader while others have little and are now functionally extinct [36][37][38]. Although North American five-needle pines did not co-evolve with C. ribicola, genetic disease resistance has been identified in each of the nine pine host species [25]. Qualitative resistance to WPBR is conferred by a dominant R gene allele (aka major gene resistance or Cr alleles) and imparts immunity to the disease. Qualitative resistance has been reported in sugar pine (P. lambertiana Dougl.; Cr1), western white pine (P. monticola Dougl. ex D. Don; Cr2), southwestern white pine (P. strobiformis Engelm.; Cr3), and, most recently, limber pine (Cr4) [39][40][41][42][43]. Cr1, Cr2 and Cr4 are unique and occur on different chromosomes in their respective hosts [44]. The Cr4 allele has been inferred to be present in limber pine stands not yet infected with WPBR and at unexpectedly high frequencies [39]. In the Southern Rocky Mountains the frequency of the Cr4 allele in limber pine is estimated at 5.0% and shows greater range among sites (0% to 13.9%) [39] than current data on Cr2 in western white pine (<0.1%; range of 0% to 0.8%) [45] and Cr1 in sugar pine range-wide (2.2%; range of 0% to 8.9%) [46]. Without past exposure to and selection pressure by WPBR, the origins of the Cr alleles remain enigmatic [39,40,45,47,48].
There is evidence to suggest that MPB preferentially attack trees with WPBR and that the likelihood of MPB attack increases with the severity of WPBR infection [49,50]. It is possible that trees weakened by other agents are simply easier for MPB to colonize [51] yet not all studies have shown a linear relationship between MPB infestation and WPBR infection. Dooley and Six [52] demonstrated that whitebark pine (P. albicaulis Engelm.) trees with very high WPBR disease severity were less likely to be attacked by MPB, even when baited, while Schwandt and Kegley [53] observed that MPB preferences vary with MPB population levels. One hypothesis is that once disease severity is high, the monoterpene profile emitted by the tree is altered and less attractive to the beetle [52].
It is unknown whether or not MPB selects WPBR infected trees because they have been weakened by the disease, or if there are underlying resistance mechanisms that deter both pest and pathogen. Here we test the latter case. Prior off-site controlled C. ribicola inoculation studies have identified forest trees inferred to have or lack the Cr4 resistance allele to WPBR in stands that have not yet been invaded by WPBR [39,54]. Using resin duct traits associated with MPB resistance (e.g., [20,21,24]), we tested for a relationship between susceptibility to WPBR and MPB in limber pine in the absence of either biological agent. We sampled healthy trees previously inferred to carry the dominant Cr4 allele for complete resistance to WPBR and trees that lack the dominant allele. All trees had been treated with an anti-aggregate pheromone (verbenone) to deter MPB attack during the recent epidemic; none of the trees showed signs or symptoms of WPBR. If trees that lack the dominant Cr4 resistance allele have traits that make them more susceptible to mortality by MPB than those trees with the dominant Cr4 allele, MPB could have historically contributed to natural selection against WPBR susceptibility in advance of WPBR's invasion, and could help explain the unusually high frequency of Cr4 resistance in the southern Rocky Mountains. We also discuss how MPB may affect the future evolution of resistance to WPBR in the pines where the two pests coincide and how WPBR may affect post-MPB epidemic forest recovery processes.

Materials and Methods
In 2016, 93 trees were sampled at 22 sites in Southern Wyoming and Colorado with elevations ranging from approximately 2350 to 3350 m ( Table 1). As stated above, we sampled healthy trees formerly inferred to carry the dominant Cr4 allele for complete resistance to WPBR (R trees) and WPBR susceptible trees that lack the dominant allele (S trees). Mean stem diameter at 1.3 m (DBH) was 27.1 cm (2.0 cm standard error (SE)) and 23.6 cm (1.1 cm SE) for R and S trees, respectively. Progeny from each tree sampled had been previously tested for genetic resistance to WPBR as described by Schoettle et al. [39]. Two 5 cm long xylem cores were extracted from each tree-one from the west side and one from the east side-at DBH using a 4 mm diameter increment borer; metrics from the two cores were combined for each sample tree. At two sites a 5 mm diameter increment borer was used-these trees are omitted from non-standardized calculations resulting in different sampling values for these metrics. Cores were prepared using techniques described in Kane and Kolb [20]. Each core was mounted on a wooden block and sanded with a belt sander to create a flat cross-section for analysis.

Resin Duct Measurements
Using a combination of traits and methods described by Ferrenberg et al. [21] and Hood and Sala [24], extracted cores were analyzed using a Dino-lite pro digital microscope model AM-413ZT (AnMo Electronics Corporation, New Taipei, Taiwan, 2000). We counted the number of resin ducts (no. year −1 ) ( Figure 1) and measured resin duct size (mm 2 ), ring width (mm), and ring area (area of ring on the core; mm 2 ) for 20 years of growth on the core; previous research suggests that vertical resin ducts produced over 25 years continue to contribute to resin production [23]. We calculated total duct area (sum of duct size per annual ring; mm 2 yr −1 ), duct density (total number of ducts per annual ring divided by ring area; no. mm 2 yr −1 ), and relative duct area (total duct area divided by ring area × 100) for each tree. In addition to analyzing all 20 years of rings (Table S1), we analyzed a 10-year interval from 2001 to 2010 which captures dates of a significant drought both before and during the recent MPB epidemic. We also calculated the metrics for the entire 10-or 20-year period following Ferrenberg et al. [21] methods; the results yielded the same conclusions as the annual averages so we report the averages by year as recommended by Hood and Sala [24], unless otherwise stated.

Data Analysis
Resin duct metrics were compared between R and S trees using t-tests; to improve normality, metrics were log-transformed prior to analyses. To adjust for multiple comparisons, we used a Benjamin-Hotchburg adjusted p-value for significance. Regression analyses were used to examine relationships between growth and duct traits with elevation. All analyses were computed in the statistical software package R version 3.4.3 [55]. . We counted the number of resin ducts (no. year −1 ) ( Figure 1) and measured resin duct size (mm 2 ), ring width (mm), and ring area (area of ring on the core; mm 2 ) for 20 years of growth on the core; previous research suggests that vertical resin ducts produced over 25 years continue to contribute to resin production [23]. We calculated total duct area (sum of duct size per annual ring; mm 2 yr −1 ), duct density (total number of ducts per annual ring divided by ring area; no. mm 2 yr −1 ), and relative duct area (total duct area divided by ring area × 100) for each tree. In addition to analyzing all 20 years of rings (Table S1), we analyzed a 10year interval from 2001 to 2010 which captures dates of a significant drought both before and during the recent MPB epidemic. We also calculated the metrics for the entire 10-or 20-year period following Ferrenberg et al. [21] methods; the results yielded the same conclusions as the annual averages so we report the averages by year as recommended by Hood and Sala [24], unless otherwise stated.

Data Analysis
Resin duct metrics were compared between R and S trees using t-tests; to improve normality, metrics were log-transformed prior to analyses. To adjust for multiple comparisons, we used a Benjamin-Hotchburg adjusted p-value for significance. Regression analyses were used to examine relationships between growth and duct traits with elevation. All analyses were computed in the statistical software package R version 3.4.3 [55].

Results
There were no significant differences between limber pine trees previously inferred to carry the dominant Cr4 allele for complete resistance to WPBR (R trees) and those susceptible to WPBR (S trees) with respect to mean radial ring width, number of vertical resin ducts, total duct area, duct density or relative duct area averaged per year over the 2001-2010, 10-year interval in the sampled wood cores (Table 2, Figure 2).

Results
There were no significant differences between limber pine trees previously inferred to carry the dominant Cr4 allele for complete resistance to WPBR (R trees) and those susceptible to WPBR (S trees) with respect to mean radial ring width, number of vertical resin ducts, total duct area, duct density or relative duct area averaged per year over the 2001-2010, 10-year interval in the sampled wood cores (Table 2, Figure 2).   The values of each metric ranged approximately four-fold among the sample trees (see Figure  3). The total number of resin ducts produced during the 10-year period and the duct area was positively correlated with radial growth during the same period; the relationships were the same for R and S trees ( Figure 3A,B). Resin duct density and relative duct area significantly decreased with radial growth and again were the same for R and S trees ( Figure 3C,D).
Radial growth (p = 0.0003, n = 93), duct production (p < 0.0001, n = 83), and duct area (p = 0.0011, n = 83; data in table S2) decreased with increasing source elevation ( Figure 4A,B); the relationships were significant for R, S and all of the trees (Table S2). Duct density (p = 0.2072, n = 93) and relative duct area (p = 0.3599, n = 93) were not significantly correlated with source elevation (Figure 4C,D) when all sample trees were included. Relative duct density was positively correlated with source elevation for S trees (p = 0.0303, n = 44) and not for R trees (p = 0.1650, n = 49) (Table S2). The values of each metric ranged approximately four-fold among the sample trees (see Figure 3). The total number of resin ducts produced during the 10-year period and the duct area was positively correlated with radial growth during the same period; the relationships were the same for R and S trees ( Figure 3A,B). Resin duct density and relative duct area significantly decreased with radial growth and again were the same for R and S trees ( Figure 3C,D).
Radial growth (p = 0.0003, n = 93), duct production (p < 0.0001, n = 83), and duct area (p = 0.0011, n = 83; data in Table S2) decreased with increasing source elevation ( Figure 4A,B); the relationships were significant for R, S and all of the trees (Table S2). Duct density (p = 0.2072, n = 93) and relative duct area (p = 0.3599, n = 93) were not significantly correlated with source elevation (Figure 4C,D) when all sample trees were included. Relative duct density was positively correlated with source elevation for S trees (p = 0.0303, n = 44) and not for R trees (p = 0.1650, n = 49) (Table S2). were positively correlated with radial growth. Resin duct density (C) and relative duct area (D) significantly decreased with radial growth. All metrics were calculated using the extracted cores, not the entire tree stem. Black circles are trees inferred in previous studies to carry the dominant Cr4 allele for complete resistance to WPBR; red triangles are for trees that lack the dominant Cr4 allele. Regression lines are shown for those relationships that are significant (p < 0.05); regression relationships for each metric did not differ between R and S trees.  Resin duct production (A) and resin duct area (B) were positively correlated with radial growth. Resin duct density (C) and relative duct area (D) significantly decreased with radial growth. All metrics were calculated using the extracted cores, not the entire tree stem. Black circles are trees inferred in previous studies to carry the dominant Cr4 allele for complete resistance to WPBR; red triangles are for trees that lack the dominant Cr4 allele. Regression lines are shown for those relationships that are significant (p < 0.05); regression relationships for each metric did not differ between R and S trees.  were positively correlated with radial growth. Resin duct density (C) and relative duct area (D) significantly decreased with radial growth. All metrics were calculated using the extracted cores, not the entire tree stem. Black circles are trees inferred in previous studies to carry the dominant Cr4 allele for complete resistance to WPBR; red triangles are for trees that lack the dominant Cr4 allele. Regression lines are shown for those relationships that are significant (p < 0.05); regression relationships for each metric did not differ between R and S trees.  Radial growth (A) and resin duct production (B) significantly decreased with increasing source elevation. Duct density (C) was not significantly correlated with elevation and only relative duct area (D) of susceptible trees was significantly correlated with source elevation. All metrics were calculated using the extracted cores, not the entire tree stem. Symbols are as shown in Figure 3. Black circles are trees inferred in previous studies to carry the dominant Cr4 allele for complete resistance to WPBR; red triangles are for trees that lack the dominant Cr4 allele. Regression lines are shown for those relationships that are significant (p < 0.05; see Table S2).

Discussion
Genetic resistance in a species to a novel biotic stress for which it has not co-evolved may be neutral or carry a benefit to the host via a function that offers adaptive benefit for environmental factor(s), biotic or abiotic, under which it co-evolved [47]. Knowing the complete resistance status of trees in populations naïve to WPBR enables us to assess if adaptations to other stresses are associated with complete resistance and may provide insights into the evolutionary history of this resistance trait [47,48]. We hypothesized that MPB may be such a stress for limber pine, however, our data on resin duct traits suggest otherwise. Resin duct traits of mature limber pine trees that are inferred to carry the dominant Cr4 allele for complete resistance to WPBR did not constitutively differ from those of trees susceptible to WPBR. We conclude that the Cr4 allele does not predispose the trees to be more or less resistant to MPB (when predicted by resin duct traits) in the absence of WPBR. Consequently, these findings suggest historical selection imposed on limber pine by MPB, prior to the invasion of C. ribicola, which was neutral with respect to the Cr4 status of the tree and, therefore, is not likely to have contributed to the unusually high frequency of Cr4 in the Southern Rocky Mountains. Further study of resin chemistry and the induction response of duct traits and chemistry to MPB attack could provide other insights.
The density of vertical ducts observed here are similar to those reported by Ferrenberg et al. [21] for limber pine and lodgepole pine (P. contorta var. latifolia Dougl. ex Loudon). Likewise, the positive correlation between duct number and total duct area with radial growth is consistent with previously reported patterns in this species and ponderosa pine (P. ponderosa Lawson & C. Lawson) [21,24]. We found a significant negative correlation of duct density and relative duct area with radial growth in limber pine while Hood and Sala [24] found this relationship was less reliable in ponderosa pine.
Resin duct number decreased and relative duct area increased in the wood of trees growing at increasing elevation ( Figure 4) consistent with Wimmer and Grabner's [56] report of fewer resin ducts of Norway spruce from cooler habitats (Picea abies L.). Likewise, the increase in relative duct area with increasing elevation (and therefore lower air temperature; [57]), is similar to the linear relationship between relative duct area in needles and minimum air temperature in Scots pine (P. sylvestris L.) [58] yet different from the curvilinear response of needle duct traits to elevation of P. taiwanensis Hayata [59]. Limber pine invests more in resin duct defenses per unit radial wood growth with increasing elevation, however, resin flow is better predicted by resin duct number and area [24] which decreased in our study with increasing elevation. Furthermore, monoterpene concentration and diversity in resin decreased with increasing elevation in limber pine [60]. In total, these findings suggest that high elevation limber pine may be structurally and chemically less defended against MPB than those trees from lower elevations. This is consistent with the elevation gradient in the plant defense hypothesis [60] that posits greater pest pressure, and therefore evolution of greater defense, at lower elevations. However, these results are only partially consistent with the resource availability and plant defense hypothesis [61] that suggests an evolutionary response to low resource availability selects for slow growth rates and larger investments in defense. This discrepancy may be an outcome of the historical thermal limitation of MPB to lower elevations in limber pine ecosystems and, therefore, selection for defense to MPB has not been a significant force in the evolution of the high elevation populations. Both genetic and environmental factors contribute to constitutive and inducible vertical resin duct and chemical defenses [62,63]. The relative sensitivity and plasticity of the (1) defense response in limber pine and (2) MPB population dynamics to climate change will contribute to the outcome of the continued limber pine-MPB interaction in mountain-top forests.
Our data suggest that MPB has not historically selected directly for or against genetic resistance to WPBR, and therefore the current field observations of MPB preference for WPBR infected trees is likely because the trees have been weakened by the disease, and not because of underlying resistance mechanisms that deter both pest and pathogen. However, MPB preference for WPBR weakened trees, and therefore WPBR susceptible trees, can accelerate the selection against WPBR susceptibility thereby increasing the proportion of the remaining trees that may have resistance to WPBR [64]. This may occur in populations with trees with complete resistance but is less likely for trees with quantitative resistance. Complete (qualitative) resistance (i.e., Cr4) prevents WPBR from developing, however, there is still a carbon cost to the tree for expressing the immunity response upon infection by C. ribicola [48] which may be sufficiently high to affect MPB attack. Trees with quantitative resistance still develop the WPBR disease, albeit more slowly [65], and may therefore be preferentially attacked by MPB. It is unknown if trees that are susceptible to WPBR and trees with quantitative WPBR resistance that have a similar disease severity are equally attractive to MPB attack and mortality. Quantitative resistance mechanisms are polygenic and therefore are less likely to be overcome and rendered ineffective by mutation in the rust fungi (i.e., is more durable) than qualitative resistance. Increasing the frequency of quantitative resistance to WPBR in five-needle pine populations is a cornerstone of restoration and conservation strategies to sustain five-needle pine populations [32,66,67]. In the presence of WPBR, MPB may (1) accelerate selection against trees fully susceptible to and symptomatic with WPBR and for trees with complete resistance to WPBR, and (2) be neutral or select against infected trees with quantitative resistance to WPBR.
The five-needle pines are early seral species and the canopy openings generated in a post-MPB epidemic landscape can provide regeneration opportunities for the pines (e.g., [68,69]). However, if mature tree mortality is too high due to both pest and pathogen, the surviving population may be too few to support forest recovery [70,71]. Furthermore, WPBR kills trees of all ages and kills young seedlings especially rapidly [72], so even if sufficient seed sources persist, forest regeneration and development in the presence of WPBR is likely to be slowed due to seedling and sapling mortality by WPBR. If the frequency of heritable resistance is high enough in the surviving overstory population to generate viable seedling and sapling cohorts, the rapid mortality by WPBR in young age classes may accelerate selection for genetic resistance to WPBR [73]. If the surviving overstory population lacks heritable resistance, the population could be at risk for extirpation if supplemental planting of seedlings with genetic resistance is not implemented [74]. The regeneration dynamic of these species is already slow (e.g., References [75,76]) and a further delay in site occupation may change forest succession trajectories.

Conclusions
These findings suggest historical selection imposed on limber pine by MPB, prior to the invasion of C. ribicola, was neutral with respect to the Cr4 status of the tree and therefore is not likely to have contributed to the unusually high frequency of Cr4 resistance to WPBR in the Southern Rocky Mountains. Therefore, the current field observations of MPB preference for WPBR infected trees is likely because the trees have been weakened by the disease, and not because of underlying resistance mechanisms that deter both pest and pathogen. However, MPB and WPBR impacts on five-needle pine populations together are likely to be greater than additive. MPB may affect the future evolution of resistance to WPBR in the pines where the two pests coincide and WPBR will affect forest recovery after MPB epidemics. The interactive effects of MPB and WPBR in a changing climate put the health of high elevation western North American ecosystems at risk.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1999-4907/9/10/595/s1, Table S1: Mean metrics per year over a 20 year interval from 1996 to 2015 for R and S trees from the extracted limber pine cores, Table S2: Regression statistics for growth and resin duct traits over a 10 year interval from 2001 to 2010 for R and S trees from the extracted limber pine cores as a function of elevation.