Whitebark pine mortality related to white pine blister rust, mountain pine beetle outbreak, and water availability

Whitebark pine (Pinus albicaulis) forests in the western United States have been adversely affected by an exotic pathogen (Cronartium ribicola, causal agent of white pine blister rust), insect outbreaks (Dendroctonus ponderosae, mountain pine beetle), and drought. We monitored individual trees from 2004 to 2013 and characterized stand-level biophysical conditions through a mountain pine beetle epidemic in the Greater Yellowstone Ecosystem. Specifically, we investigated associations between tree-level variables (duration and location of white pine blister rust infection, presence of mountain pine beetle, tree size, and potential interactions) with observations of individual whitebark pine tree mortality. Climate summaries indicated that cumulative growing degree days in years 2006–2008 likely contributed to a regionwide outbreak of mountain pine beetle prior to the observed peak in whitebark mortality in 2009. We show that larger whitebark pine trees were preferentially attacked and killed by mountain pine beetle and resulted in a regionwide shift to smaller size class trees. In addition, we found evidence that smaller size class trees with white pine blister rust infection experienced higher mortality than larger trees. This latter finding suggests that in the coming decades white pine blister rust may become the most probable cause of whitebark pine mortality. Our findings offered no evidence of an interactive effect of mountain pine beetle and white pine blister rust infection on whitebark pine mortality in the Greater Yellowstone Ecosystem. Interestingly, the probability of mortality was lower for larger trees attacked by mountain pine beetle in stands with higher evapotranspiration. Because evapotranspiration varies with climate and topoedaphic conditions across the region, we discuss the potential to use this improved understanding of biophysical influences on mortality to identify microrefugia that might contribute to successful whitebark pine conservation efforts. Using tree-level observations, the National Park Service-led Greater Yellowstone Interagency Whitebark Pine Long-term Monitoring Program provided important ecological insight on the size-dependent effects of white pine blister rust, mountain pine beetle, and water availability on whitebark pine mortality. This ongoing monitoring campaign will continue to offer observations that advance conservation in the Greater Yellowstone Ecosystem.


INTRODUCTION
Global forest area has declined by 40% over the last three centuries (Shvidenko et al. 2005), and continued decreases are projected over the next few decades (Pereira et al. 2010). Continued mortality of forests will contribute to further reductions in global carbon pools (Pan et al. 2011), diminishing biodiversity (Maron et al. 2015), and losses to watershed protection (Adams et al. 2012). While human conversion of forest to agricultural lands has historically reduced global forest area (Geist and Lambin 2002), increasing temperatures and drought stress have been implicated in the contemporary, large-scale die-offs of the world's forests (Allen et al. 2010). Although climate has always shaped forests, recent changes in temperature and precipitation patterns are interacting with primary forest disturbance agents (disease, insects, and drought) to further alter the distribution, abundance, and growth of forest tree species (Sturrock et al. 2011).
Increases in tree mortality have been well documented on the six forested continents (Allen et al. 2010, Carnicer et al. 2011, Das et al. 2013, ter Steege et al. 2015. Agents of forest mortality can indiscriminately remove all species within a stand (Westerling et al. 2011) or select against specific taxonomic groups (Anderegg et al. 2013). This latter phenomenon has occurred in several forested regions of western North America where numerous species of Pinus have experienced high rates of mortality in recent decades (Raffa et al. 2008). In the Greater Yellowstone Ecosystem (GYE), elevated levels of mortality occurring in whitebark pine forests (Pinus albicaulis) have been documented, and these shifts are emblematic of a forest species undergoing dynamic change (Logan et al. 2009).
Like many forest species that have undergone dramatic decreases in recent decades, there are multiple hypotheses about the singular and interactive causes contributing to whitebark pine decline. Cronartium ribicola, an exotic fungal pathogen that causes the disease white pine blister rust (blister rust), recent mountain pine beetle (Dendroctonus ponderosae) outbreaks, and decades of fire suppression have individually contributed to the observed decline of whitebark pine across its range (Arno 1986, Kendall and Arno 1990, Logan et al. 2010, Tomback and Achuff 2010, MacFarlane et al. 2013). However, it remains uncertain during an epidemic phase of mountain pine beetle whether the effects of these disturbance agents are additive or interactive-the presence of multiple agents increases the probability of mortality compared with an agent in isolation (Schwandt and Kegley 2004, Six and Adams 2007, Bockino and Tinker 2012, Dooley and Six 2015. Further complicating these relationships, drought and tree diameter may interact with mountain pine beetle and blister rust (van Mantgem et al. 2009, Logan et al. 2010, Allen et al. 2015. For example, mountain pine beetle preferentially select large-diameter pine trees (Cole andAmman 1969, Negron et al. 2008); the availability of larger whitebark pine trees promotes mountain pine beetle outbreaks and sustains epidemics under favorable climate conditions (Perkins and Roberts 2003, Safranyik and Carroll 2006, Raffa et al. 2013. In contrast, blister rust-related mortality is highest in small size classes (Tomback et al. 1995); however, the processes contributing to mortality are more nuanced. Blister rust mortality is dependent on infection location, where bole infections are often considered more immediately lethal (Koteen 2002, Newcomb 2003, and, although not formally assessed by others, on infection duration. Additionally, there are mounting concerns that changes in climate will further exacerbate the impacts of disturbance agents on whitebark pine (e.g., Powell 2004, Chang et al. 2014;Keane et al., in press) by amplifying host tree susceptibility to disease and insect outbreaks (sensu Sturrock et al. 2011). For instance, previous periods of warming in the 1930s, 1970s, and 1980s resulted in epidemic mountain pine beetle levels, with populations declining to endemic levels following the return to cooler temperatures (Logan et al. 2010). In addition, drought stress can impair the natural defense mechanisms of whitebark pine leaving them more vulnerable to successful mountain pine beetle attack (Raffa et al. 2013).
In the GYE, increases in temperatures and snowpack declines have been documented (Pederson et al. 2011, Sepulveda et al. 2015. These changes have occurred coincident with observations of widespread whitebark pine die-offs in the region, but to date, there have been few considerations of water availability that contributed to the mass mortality event (Buotte et al. 2016) and no formal consideration given to topoedaphic (terrain and soil) interactions. A water balance approach (Lutz et al. 2010) integrates temperature and precipitation with topoedaphic factors to generate biophysical metrics such as evapotranspiration and water deficit that may affect whitebark pine's ability to withstand disease agents (Millar et al. 2012). This approach also acknowledges that local conditions moderate climate (Dobrowski 2011) and contribute to spatiotemporal heterogeneity in the biophysical environment of whitebark pine stands. Importantly, biophysical estimates of water availability may prove useful in explaining patterns and mechanisms of observed whitebark pine decline.
As described, the ecological threats to whitebark pine are multiple, mechanistically complex, and play out over 7 million ha in the GYE. Given that management resources are limited, spatially explicit information on the environmental conditions related to tree mortality and survival is increasingly necessary for targeted management action. The National Park Service-led, long-term, ground-based monitoring of permanently tagged trees provides fundamental repeat observations to disentangle the effects of blister rust, mountain pine beetle, and water availability on whitebark pine mortality and offers new insight for managing protected areas. Here, we use long-term monitoring data from high-elevation whitebark pine stands across the GYE to assess the interplay between multiple disturbance agents and water availability on tree mortality during a mountain pine beetle epidemic. Specifically, we examine whether (1) larger trees with evidence of mountain pine beetle have greater probability of mortality, (2) increasing duration of blister rust infection in the canopy or bole results in increasing probability of mortality, and whether these relationships vary with tree diameter, (3) blister rust and mountain pine beetle interact to increase the probability of whitebark pine mortality, and (4) an interaction between mountain pine beetle attack and stand-level water availability influences mortality.

Study area
The study area is geographically defined as the Yellowstone Plateau volcanic fields and the 14 surrounding mountain ranges (Marston and Anderson 1991). Our study included whitebark pine stands located at elevations spanning 2400-3172 m in five national forests (NF), two national parks (NP), as well as state and private lands in portions of Wyoming, Montana, and Idaho ( Fig. 1). Whitebark pine is a keystone conifer and its forests are biologically significant components of high-elevation regions in the GYE and Northern Rocky Mountains. Blister rust was introduced to North America in 1910 and has resulted in severe declines in whitebark pine populations throughout its range (Kendall and Arno 1990, Keane and Arno 1993, Tomback and Achuff 2010. The blister rust life cycle is complex (McDonald and Hoff 2001), and infection remains ubiquitous and variable across the GYE; estimated overall infection rates vary from 20% to 30% for whitebark pine trees in this region (WPMWG 2014). Whitebark pine mortality and mountain pine beetle activity are also widespread and vary in severity across the GYE. Temperature strongly influences voltinism and adult survival as well as the synchrony and duration of mountain pine beetle emergence Schen-Langenheim 2007, Logan et al. 2010). Combined, these responses to warming temperatures at higher elevations contribute to elevated mountain pine beetle outbreaks in whitebark pine habitat.

Field sampling
Whitebark pine stands were defined as a contiguous area of forest with whitebark pine as the dominant or co-dominant component. Stands were identified using photo-interpretation (Dixon 1997) and stand composition cover maps (WPMWG 2011). Stands were selected randomly from a sample frame of approximately 10,770 mapped whitebark pine stands ≥2.0 ha (Dixon 1997, Landenburger et al. 2008; Fig. 1). We used a probabilistic, two-stage cluster design; whitebark pine stands were the primary sample units and 10 9 50 m transects (secondary sample units) were established within selected stands (Lohr 2010).
From 2004 to 2007, 176 permanent transects in 150 whitebark pine stands were established. All live, whitebark pine trees >1.4 m tall (4768 individual trees) located within the boundaries of transects were tagged and tree size (diameter at breast height in centimeters; dbh) was recorded.
Live trees were examined for blister rust infection, and the location of blister rust infections (canopy or bole) was recorded. A canopy canker was defined as any infection occurring within the canopy of a tree ≥5 cm distal from the main bole. A bole canker was an infection on the main trunk of the tree or within 5 cm of the bole along peripheral branches. Bole cankers are considered more immediately lethal and therefore more severe than canopy cankers (McDonald andHoff 2001, Schwandt andKegley 2004). At all subsequent site revisits, we documented tree status (live/dead), blister rust infection status and infection location for live trees, and mountain pine  beetle indicators (pitch tubes and frass in live trees and J-shaped galleries beneath the bark of dead trees). After 2007, stands were randomly assigned to one of four rotating panels. A panel is a grouping of 44 permanent transects that are revisited every two to four years. Over the course of the mountain pine beetle epidemic (2008)(2009)(2010)(2011)(2012)(2013), all transects were surveyed every two years to document mortality (Appendix S1: Table S1). Trees with all brown needles or devoid of needles were classified as dead in the year that this condition was observed (any year from 2008 to 2013).
Following initial blister rust infection, it takes approximately four years to visually detect infection (McDonald and Hoff 2001). Therefore, field crews evaluated trees for blister rust infection every four years (Appendix S1: Table S1). Tagged trees had two to five observations between 2004 and 2013.

Data analysis
We modeled the probability of whitebark pine mortality (i.e., our response variable was tree death) using a generalized linear mixed model with logit-link function. We included a random intercept for each transect (Zuur et al. 2009), which accounted for potential correlation among trees within a transect and for the cluster sampling design (Lohr 2010). We assumed no additional complexity for the correlation structure, which was supported by graphical displays that indicated the tree-level variable relationships with mortality did not vary spatially among transects. Graphical investigation of predictors also demonstrated that there were no outliers that would exert undue leverage in models.
To develop statistical models for understanding tree-level characteristics that influence probability of mortality, we simplified the repeated measurements of a tree into one record (n = 4768 unique trees). A tree was considered alive only if green needles were observed during our last survey of each transect (either 2012 or 2013). To characterize blister rust infection, we used information on the duration (number of years) and location (bole vs. canopy) as a proxy for severity, where infections in the bole were considered more severe than those in the canopy. Additionally, we considered infections of longer duration in the bole more imminently lethal than those of shorter duration. Duration of infection was calculated as the difference in years between the last and first time blister rust was recorded on a tree in a given location. For trees that died, blister rust infection location on trees was assumed unchanged since the previous survey. Both dbh and dbh-squared were centered (meaning the average dbh was subtracted from each tree's specific dbh value) for modeling to reduce variance inflation and the effects of collinearity. Predictors had variance inflation factors less than five, suggesting no evidence of multicollinearity (Marquardt 1970, Neter et al. 1996. To determine the best predictive statistical model for blister rust infection-related mortality, we fitted six a priori models: three with presence/absence predictors for blister rust and three that included duration of blister rust infection (Table 1). For each tree-level model, we included a variable that denoted whether a tree had evidence of mountain pine beetle or not (presence or absence), dbh, dbh-squared, and their twoand three-way interactions with each other and with the respective blister rust variables (Table 1). We selected a top model based on the strength of evidence from AIC (Burnham and Anderson 2002). Using this model as a starting point, we then used nested subset, backwards selection via drop-in-deviance tests (Cheng et al. 2010) to arrive at a more parsimonious model for ecological interpretation. Specifically, we assessed evidence of an interaction based on whether the probability of mortality was greater for trees with signs of both mountain pine beetle and blister rust infection compared to either disturbance agent alone, conditional on other included predictors.
We provide an assessment of the predictive performance of the top-ranked AIC model from the set of six a priori models compared to its more parsimonious counterpart in Appendix S2. Next, we incorporated climate predictors into our modeling process to evaluate whether they provided additional explanation of the observed variation in mortality and whether the climate effect was dependent on mountain pine beetle attack. Weather data for stands were obtained from single-pixel, daily time series of Daymet interpolated 1-km climate grids (Thornton et al. 2014). Biophysical parameters were calculated on a monthly time-step using a Thornthwaite-type water balance model to estimate water input to the soil, soil moisture storage, and loss of water via evapotranspiration. The monthly water balance equations were identical to Lutz et al. (2010) with these exceptions: Each year was calculated separately and December soil moisture and snowpack were carried over to the following January. The water balance model accounts for variation in stand-level slope and aspect that affect heat load and water-holding capacity of the soil. Soil water-holding capacity was obtained from soil components within soil map units containing each stand (USDA NRCS 2013).
We selected two climate variables, temperature and precipitation (denoted temp and ppt in Table 2) that coarsely represent water and energy for growth. We also included the interaction between temperature and precipitation to determine whether this simple variable was as informative as more complex water balance variables. We evaluated eight water balance variables that represent seasonal estimates of biophysical factors that may be more closely linked to tree condition than annual measures of temperature and precipitation (Stephenson 1998, Hu et al. 2010, Millar et al. 2012, Chang et al. 2014, Allen et al. 2015; Table 2). We used annual values summarized for each stand for all climate variables. Potential evapotranspiration (PET) is the amount of water that could be evapotranspired with available energy if water availability was unlimited. Actual evapotranspiration (AET) is the monthly upward loss of water from soil via transpiration and evaporation, which is limited by soil moisture availability. Water deficit, a positive valued measure of drought stress, was calculated as the difference between monthly PET and AET (Stephenson 1998, water deficit, Table 2). Water minus PET (W_PET) is snow melt plus rain that reaches the soil surface minus PET. We also include cumulative growing degree days (cumgdd, Table 2) as a measure of energy available for tree growth and difference between water reaching the soil surface and PET (W_PET) in case estimates of soil properties and their effect were not well characterized. We extended the parsimonious, tree-level statistical model for mortality by singly including each lag one climate variable in an effort to match the climate metrics with tree condition (stress) at the time of mountain pine beetle attack. A similar generalized linear mixed-model framework was used, as the binary response was still alive or dead, but we included a random intercept for each transect-year combination. We considered correlation among binary observations (lack of independence among trees within a transect for a given year) a nuisance and not of direct interest, so we specified an unstructured spatiotemporal random effect term in the models. This model structure assumes that there is an emergent climatic relationship with mortality in the GYE. A deeper investigation into whether there was among-transect variation in the climatic-mortality relationship was beyond the scope of the data available during this time period and was also complicated by the limited number of repeat observations available for each transect over time.
To compare non-nested models for whitebark pine mortality that included the lag one climate predictors, we used AIC (Burnham and Anderson 2001) to rank our models and also examined the adjusted R 2 (Nakagawa and Schielzeth 2013) as an indicator of goodness-of-fit. Adjusted R 2 can be considered as the proportion of variance explained by fixed factors alone (marginal adjusted R 2 ) or by both fixed and random factors (conditional adjusted R 2 ; Nakagawa and Schielzeth 2013). We used the function glmmML from the glmmML package (Brostr€ om 2013) available for estimating generalized linear mixed models in R (R Core Team 2014).

RESULTS
Of all the trees in the study through 2013, roughly one quarter had evidence of mountain pine beetle (1272 trees), while the remaining 3496 trees did not. Only 563 of 4768 trees had evidence of both mountain pine beetle attack and blister rust infection; 2129 trees lacked evidence of either stressor. Overall, trees with evidence of mountain pine beetle were larger in diameter than those without (Fig. 2). The majority of trees that had evidence of mountain pine beetle attack were classified as dead that same year or during the next revisit two years later. Even though this was an observational study, the observed size class distribution within a transect was similar to the among-transect distribution of tree sizes. Observed patterns in mortality, mountain pine beetle attack, and blister rust infection were suggestive of possible spatial processes for the spread across the region. However, with our dataset, we were unable to investigate spatial structure beyond incorporating a random intercept for each transect given the sampling design and available spatial predictors.

Tree-level predictors of mortality
The model used to describe whitebark pine mortality that included duration of canopy and bole infections had the greatest support (Table 3). Generally, there was more support for models that differentiated the severity of infection (bole vs. canopy infections) compared to models that included a variable for simply presence or absence of blister rust infection ( Table 3).
All of the models in the a priori set included mountain pine beetle, dbh (cm), dbh-squared, and their two-and three-way interactions (Table 3). We were able to reduce the complexity of the top-ranked AIC model substantially as there was no evidence for the more complicated interactions (results available in Appendix S2). We also found that the reduced model was essentially equivalent (in terms of predictive performance) to the top-ranked AIC model (results provided in Appendix S2). Using our reduced model (Appendix S2: Table S2), we found that for all size classes, trees without evidence of either mountain pine beetle or blister rust showed lower mortality levels (Fig. 3a). We found overwhelming evidence that smaller size class trees with longer duration of bole infection (≥5 years) had increased probability of mortality regardless of mountain pine beetle attack during the time frame of our study (Fig. 3c). There were only four trees with long durations of bole infections and evidence of mountain pine beetle accounting for the large uncertainty and potentially non-representative patterns (Fig. 3c). Interestingly, we found that for larger trees with evidence of mountain pine beetle attack, the probability of mortality was greater for those without bole infections compared to trees with longer duration of bole ❖ www.esajournals.org infections (Fig. 3, comparison across panels). This pattern suggests that there was no evidence for an interactive effect between primary disturbance agents that led to higher probability of mortality for larger trees.

Climate associations with probability of mortality
We found that adding AET to our final treelevel model was equivalent to adding PET (AET and PET within DAIC < 4 units, Table 4). For these top two, equivalent climate models, there Table 3. Rankings of six a priori models that provide alternative characterizations of tree-level blister rust infection impact on whitebark pine mortality. Notes: We explored location of infection (bole, canopy, or any) as a proxy for severity, where infections in the bole were considered more severe than those in the canopy. We compared models that specified duration of infection vs. using a variable that was simply presence or absence (P/A), where infections of longer duration are considered more severe than those of shorter duration. For each model, we included tree-level characteristics or variables for presence/absence of mountain pine beetle, dbh, dbh-squared, and their two-and three-way interactions with each other and with the respective blister rust variables. The best fitting model was selected based on the change in Akaike's information criterion (DAIC). The total coefficients column refers to the number of fixed parameter coefficients estimated for each model and is a measure of model complexity. was strong evidence for a cross-scale interaction (tagged tree to stand-level climate and water balance: MPB 9 AET P = 0.0001 from drop-indeviance test with 1 df; MPB 9 PET P < 0.0001 from drop-in-deviance test with 1 df). Given the similarity of AET and PET in models of mortality, hereafter we refer to these variables collectively as evapotranspiration. For larger trees attacked by mountain pine beetle, increased evapotranspiration was associated with a lower estimated probability of mortality, whereas trees not attacked by mountain pine beetle provided no evidence of a relationship between evapotranspiration and mortality (Fig. 4, patterns for AET were similar to PET displayed in Appendix S3: Fig. S1).
There was marginal, but inconclusive, support that areas with higher water deficit had an increased probability of mortality for whitebark pine attacked by mountain pine beetle (Appendix S3: Fig. S2). Transects with soils in the lowest quartile for water-holding capacity were consistently among those with the highest water deficit values for 2008-2013. There was no support that any of the remaining climate variables such as annual temperature or precipitation were informative (DAIC > 15 units, Table 4).
More of the variance in mortality observed during the mountain pine beetle epidemic was explained by incorporating the spatiotemporal random effect for each year 9 transect combination than by just a single climate variable (Table 4, comparison of marginal and conditional R 2 ). This suggests that in our analyses, there was a substantial amount of spatiotemporal variation that was left unexplained after accounting for climatic conditions (e.g., water deficit or evapotranspiration) and tree measures of blister rust and mountain pine beetle attack. It is possible that additional stand-level variables that are proxies for the spatial spread over time of the mountain pine beetle Fig. 3. Estimated probability of mortality given tree diameter at breast height (dbh) from the model selected by backward, subset elimination (nine coefficients; see Appendix S2: Table S2). Duration of blister rust bole infection was divided into three discrete groups: (a) no bole infection, 0 years; (b) short-duration bole infection, 1 ≤ years ≤ 4; and (c) long-duration bole infection, ≥5 years. Histograms correspond to counts of individual trees that were dead (top axis) or live (bottom axis) at last survey. Pink indicates that mountain pine beetle (MPB) were absent and green indicates that mountain pine beetle were present. Duration of canopy infection was assumed to be the mean for all trees (1.65 yr, SD = 2.8). Confidence interval is a pointwise, t-based, 90% interval. outbreak in the GYE would improve model fit (increase in marginal R 2 ).

DISCUSSION
The NPS-led Greater Yellowstone Interagency Whitebark Pine Long-term Monitoring Program chronicled the mountain pine beetle epidemic and blister rust-related mortality in whitebark pine from 2008 to 2013 within the GYE. Here, we demonstrate how ground-based observations combined with local climate and topoedaphic information can be used to clarify the interrelationships underpinning the dramatic decline of whitebark pine (>80% of canopy coverage; Logan et al. 2010, MacFarlane et al. 2013. We highlight the most prominent results related to our objectives, provide biological explanations for these findings, and discuss implications in the context of current and potential management actions.

Mountain pine beetle-mediated demographic shift
Across all stand-year combinations, the diameter of mountain pine beetle-attacked whitebark pine trees was consistently greater than that of un-attacked trees. As a consequence, the size class distribution of surviving whitebark pine within the GYE shifted to smaller-diameter trees in just six years (Appendix S4: Fig. S1). Others have documented mountain pine beetle selection preference for larger trees (Perkins and Roberts 2003, Raffa et al. 2013, Buotte et al. 2016. Larger trees also support greater numbers of emerging adult mountain pine beetles that are required to sustain an outbreak of epidemic levels (Cole and Amman 1969, Negron and Popp 2004, Negron et al. 2008. Following the observed mass mortality peak (in 2009), we documented a decline in the annual count of newly dead trees with evidence of mountain pine beetle starting in 2011 (WPMWG 2014). There are at least two possible explanations for the observed decline in whitebark pine mortality. First, there may have been a decrease in the density of large trees such that epidemic levels of mountain pine beetle could no longer be sustained. Another possible explanation is that the cold temperatures present in October 2009, which averaged six degrees colder than other Octobers in our time series, may have killed mountain pine beetles before they became cold-hardened (Dooley and Six 2015; Appendix S5: Fig. S1).

Size and blister rust mortality
While larger-diameter trees were more adversely impacted by mountain pine beetle in our study, we found that smaller-diameter trees (<20 cm dbh) had greater probability of mortality with evidence of blister rust infection. In addition, we determined that bole infections lasting ≥5 years were more imminently lethal in smaller size class trees. Smaller trees have fewer and shorter branches, and therefore, once a smaller tree is infected in the canopy, the distance an infection has to grow from a branch to the main bole is less than in larger-diameter trees. Another biological hypothesis for this pattern is that larger trees can resist infection by shedding branches or walling off infections (Tomback et al. 1995). In addition to losing future reproductive potential, smaller trees may die from blister rust Notes: Each model is a generalized linear mixed model with regression coefficients for 11 fixed effects: an intercept and regression coefficients for mountain pine beetle presence/ absence, duration of bole infection, duration of canopy infection, dbh-squared, and the specified climate variable; twoway interactions between mountain pine beetle and dbh, mountain pine beetle and dbh-squared, duration of bole infection and dbh, and mountain pine beetle and the climate variable. The row with the dash (Climate Variable column) indicates a model with no climate variable. Each model also included a random intercept for spatiotemporal variation for each transect-year combination. The marginal adjusted R 2 can be interpreted as the proportion of variance explained by fixed factors (climate variables), while the conditional R 2 can be interpreted as the proportion of variance explained by both fixed (climate variables) and random (transect) 9 year random intercepts factors. infection more rapidly than larger infected trees (Smith and Hoffman 2000, Koteen 2002, Newcomb 2003. With a shift toward smaller trees in the overall population, blister rust may once again pose a significant management challenge in the GYE.

Interactive effects of mountain pine beetle and blister rust on mortality
Other studies in the GYE have described an interactive association between mountain pine beetle and blister rust where blister rust-infected whitebark pine were preferentially selected by Fig. 4. Fitted estimates of whitebark pine mortality given continuous PET conditional on tree diameter and mountain pine beetle status (MPB, present/absent). Each graph provides conditional estimates of mortality by size class using the midpoint of each size class (dbh ≤ 2.5 cm = 2 cm, 2.5 < dbh ≤ 10 cm = 6 cm, 10 < dbh ≤ 30 cm = 18 cm, dbh > 30 cm = 40 cm). Estimates assume that duration of bole infection was the mean of all trees (0.8 yr, SD = 2.4 yr) and that duration of canopy infection was the mean of all trees (1.6 yr, SD = 2.8 yr). Histograms correspond to counts of individual trees (right margin) that were dead (top axis) or live (bottom axis). In (a) and (c), the green represents trees with dbh less than or equal to 2.5 cm and the pink represents the trees with dbh greater than 2.5 and less than or equal to 10 cm. In (b) and (d), the yellow represents the trees with dbh greater than 10 and less than or equal to 30 cm and the gray represents the trees with dbh greater than 30 cm. Confidence interval is a pointwise, t-based, 90% interval. mountain pine beetle for attack (Bockino andTinker 2012, Dooley andSix 2015). These studies used a narrower spatiotemporal sampling design. In contrast, our efforts in the GYE were conducted over the course of the mountain pine beetle outbreak (MacFarlane et al. 2013, Buotte et al. 2016) across a large, topographically, and climatologically heterogeneous region, where we found no support for an interaction on the probability of mortality from concurrent mountain pine beetle infestation and blister rust infection. We offer at least three explanations for the discrepancy between our findings and other whitebark pine studies in the GYE. First, our study sites were randomly selected over an extensive region where mountain pine beetle infestation and blister rust levels were both spatially and temporally diverse (WPMWG 2014). Other studies targeted sites for the prevalence of active mountain pine beetle infestations in stands with high levels of blister rust infection (Bockino andTinker 2012, Dooley andSix 2015). Given this fundamental difference in site selection, the presence or absence of an interactive effect may simply be an artifact of alternative study designs. A secondary explanation for the disparity is that the previous studies were much shorter in duration relative to ours and therefore were possibly conducted over a single phase of the mountain pine beetle outbreak. In contrast, our work chronicled whitebark pine mortality that resulted from changes in mountain pine beetle population dynamics (i.e., spatially and through endemic and epidemic levels). Lastly, above-average water deficit in some years preceding the epidemic in the GYE (Appendix S5: Fig. S2) may have contributed to whitebark pine stress that increased their susceptibility to even low or endemic levels of mountain pine beetle. This latter explanation is supported by work with lodgepole pine (Pinus contorta), ponderosa pine (Pinus ponderosa), and whitebark pine that show drought-stressed trees may have impaired defenses, ultimately increasing their susceptibly to mountain pine beetle during endemic levels (Boone et al. 2011, Raffa et al. 2013. Taken together, our decade-long monitoring of whitebark pine, mountain pine beetle, blister rust, and drought makes us uniquely suited to evaluate whether interactions described between mountain pine beetle and blister rust in whitebark pine are strongly dependent on mountain pine beetle population phases or a function of other tree stressors (Boone et al. 2011). Our inability to describe an interaction in the GYE demonstrates that the relationship between mountain pine beetle and blister rust is complex and that interactive effects must be contextualized through an understanding of stand characteristics, mountain pine population phase, and the severity of drought and other environmental pressures.
Mountain pine beetle-dependent climate associations with tree mortality Historically, mountain pine beetle were prevented from reaching epidemic levels at high elevations by cold temperatures (À18°to À40°C; Safranyik and Carroll 2006). Relative to mountain pine beetle developmental requirements, warm fall temperatures allow larva to synchronize with other development stages that may otherwise emerge one to three years post-oviposition Powell 2004, Bentz et al. 2015). Warm temperatures can also cause a shift from semivoltine (multiple years to complete life cycle) to univoltine development (single-year life cycle). These changes coupled with warmer winters and longer flight seasons determine when and where outbreaks occur if host resources are available (Logan and Bentz 1999, Hicke et al. 2006, Mitton and Ferrenberg 2012. Temperature limitations on mountain pine beetle development and population growth resulted in weakly co-evolved defenses against this native insect in whitebark pine (Raffa et al. 2013). Mountain pine beetles boring into trees sever resin ducts that may result in attacking mountain pine beetles being expelled (Safranyik and Carroll 2006). Secondary responses include additional resin response and production of chemicals toxic to mountain pine beetles (Raffa et al. 2013). If the collective effects of constitutive and inducible defenses are sufficient, a tree may survive beetle attack; otherwise, mountain pine beetles can exhaust and overcome tree defenses (Safranyik and Carroll 2006). At endemic levels, mountain pine beetles select trees weakened by drought, lightning, or other stressors (Raffa et al. 2013), but during epidemic phases, mountain pine beetle densities can overcome healthy trees. This positive feedback sustained by the availability of host resources and suitable temperature can result in exponential increases in tree mortality (Logan and Powell 2004).
Although there are several requirements that need to be met, our stand-level climate data show temperatures prior to the epidemic that are consistent with shifts to univoltinism and longer growing seasons (Mitton and Ferrenberg 2012). Specifically, years 2004Specifically, years -2007 were warmer than average and each of these years was successively warmer than the previous year. Immediately preceding the epidemic peak in 2009, mean cumulative calendar year growing degree days in 2006-2008 across our stands were greater than the threshold (833°C; Carroll et al. 2006) as the thermal requirement for a univoltine population (Appendix S5: Fig. S2). Although water deficit during the epidemic (for years 2006-2008; 36 mm) was greater than the 1980-2013 average (23 mm), by comparison it was much less than water deficit estimated for an historic drought year (1988; 67 mm). Suggestive, but inconclusive support for a relationship between mortality and water deficit in mountain pine beetle-attacked trees indicates that drought stress may have contributed to conditions that led to the observed widespread forest mortality, but more likely thermal release on mountain pine beetle development and population growth was the primary driver of the epidemic until host depletion and cold temperatures in October 2009 contributed to a disruption of the epidemic.
While acknowledging that temperature influences the timing and severity of mountain pine beetle epidemics, following attack, the probability of tree mortality was most strongly related to water availability. Across our high-elevation study sites, where temperature limits the process of evapotranspiration, we show that estimates of evapotranspiration (AET and PET) were similar (AET was approximately 80% of PET; Appendix S5: Fig. S2) and water deficit was low. We found evidence that the probability of mortality with evidence of mountain pine beetle is inversely related to evapotranspiration and directly related to water deficit. These associations suggest that defense mechanisms may be supported by greater water availability and hindered by drought stress. Other studies offer support for interactions between drought stress and insect-related mortality of trees (Breshears et al. 2005, van Mantgem and Stephenson 2007, Adams et al. 2012, Kaiser et al. 2012, Buotte et al. 2016. For example, Millar et al. (2012) found that whitebark pine mortality caused by mountain pine beetle in California was associated with a 2-yr lagged water deficit that may have predisposed trees to insect pests in the following years. A lack of co-evolution between whitebark pine and mountain pine beetle (Raffa et al. 2008) or the moderate water deficit conditions estimated during the epidemic may have contributed to the significant, but weak biological effect of water availability.
Climate change projections for this region (Westerling et al. 2011, Chang et al. 2014) suggest increasing water deficit (Alder and Hostetler 2013). We found that aspect and soil water-holding capacity were the two most important topoedaphic determinants of modeled water deficit within whitebark pine stands. Importantly, water deficits were greatest in soils with lowest water-holding capacity across all stand-year combinations, which suggests that conservation efforts should be considered in locations with aspect and soil properties that combine to produce the lowest water deficits. Another consideration is that north-facing aspects that reduce water deficit due to reduced heat load may be locations that predispose whitebark to blister rust due to higher humidity (Larson 2011). For these reasons, we suggest that low water deficit conditions found with finer-textured soils in less humid settings offer the best management option for planting to reduce susceptibility to either mountain pine beetle or blister rust. Ancillary conservation benefits of this approach would accrue because these locations with lower water deficit are also more resistant to drought and fire , Kane et al. 2015 and may be advantageous for germination and tree establishment (Gelderman et al. 2016).

Conservation implications for a long-lived species
Factors such as fire, competition with other tree species, resistance to blister rust, continued animal-assisted seed dispersal, and even human management actions (Keane et al., in press) not evaluated in our study will influence future whitebark pine densities in the GYE. However, our findings evaluate the principal contemporary threats to whitebark pine and their interactions. Through our ground-based, long-term monitoring program, we also show that these threats to whitebark pine are dynamic in space and time and, importantly, that these threats periodically surpass one another in their impact. Combined, our findings have significant implications for the conservation of whitebark pine in the GYE and highlight spatial and temporal aspects that previously were not well recognized in conservation strategies (GYCCWPS 2011).
Our study described a recent mountain pine beetle-induced shift in the size class distribution of surviving whitebark pine and presented evidence that blister rust may once again surpass mountain pine beetle as a more probable cause of mortality and damage to whitebark pine, particularly in the smaller size classes in the GYE. Given this finding, management actions should acknowledge the persistent threat of blister rust and its potential to affect recruitment of the more susceptible, smaller size class cohort. To promote regeneration of resistant whitebark pine trees, continued collection of cones from blister rust-resistant trees and the planting of blister rust-resistant stock should concentrate in areas projected to support whitebark pine in the future (Chang et al. 2014, Buermeyer et al. 2016. Although benefits of these efforts may take many human generations to be realized, active planting of blister rust-resistant trees combined with the dynamics of natural selection could enable future persistence in high-elevation areas of the GYE (Loehman et al. 2010, Mahalovich 2013Keane et al., in press).
With projected higher temperature in the GYE, large-scale outbreaks of mountain pine beetle should be anticipated (Buotte et al. 2016). Mitigating for this aggressive, native forest pest at broad scales is currently implausible using chemical means (Safranyik et al. 2010). However, our finding that mountain pine beetle infestations disproportionately affect large trees provides evidence that stand-level strategies should be prioritized. For example, strategies that promote and manage for a mosaic of tree sizes within whitebark pine stands and an increase in the heterogeneity of stand ages across the landscape may aid in disrupting the spread and severity of future mountain pine beetle outbreaks.
Although it was clear that the spatiotemporal pattern of the mountain pine beetle outbreak was the most important factor in observed mortality during this study, we found evidence that climatic factors also played a role. The models that included evapotranspiration which accounts for slope, aspect, and soil properties were marginally better (in terms of variance explained) than models using only simple measures of temperature and precipitation for understanding patterns of mortality. We suggest that managers can prioritize locations for conservation actions by considering patterns of soil, slope, and aspect that are likely better for establishment, provide a buffer against drought and fire, and may also reduce probability of mature tree mortality when attacked by mountain pine beetle.

ACKNOWLEDGMENTS
We would like to thank all of the field crew members who have worked on this monitoring effort over that past 12 years. We are especially grateful to Rob Daley for his expertise and assistance with data management, as well as graphical and GIS support. We would also like to thank the multiple collaborators who have helped to inform and sponsor this ongoing, long-term monitoring program including the National Park Service, US Forest Service, US Geological Survey, Bureau of Land Management, Montana State University, and American Forests. Participation by Siri Wilmoth and Kathryn Irvine was made possible through IA P12PG70586 and in-kind USGS funding. We thank Cynthia Hollimon for her efforts analyzing earlier whitebark datasets and Dr. Steve Cherry for his statistical consultation during various aspects of the monitoring effort. And finally, thank you to Dr. Kenneth Raffa for sharing his expertise and knowledge on the generalities of mountain pine beetle and whitebark pine health. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.