Interactions between Climate and Stand Conditions Predict Pine Mortality during a Bark Beetle Outbreak

In temperate coniferous forests, biotic disturbances such as bark beetle outbreaks can result in widespread tree mortality. The characteristics of individual trees and stands, such as tree diameter and stand density, often influence the probability of tree mortality during a bark beetle outbreak. However, it is unclear if these relationships are mediated by climate. To test this, we assembled tree mortality data for over 3800 ponderosa pine trees from Forest Inventory and Analysis (FIA) plots measured before and after a mountain pine beetle outbreak in the Black Hills, South Dakota, USA. Logistic models were used to determine which tree, stand, and climate characteristics were associated with the probability of mortality. Interactions were tested between significant climate variables and significant tree/stand variables. Our analysis revealed that mortality rates were lower in trees with higher live crown ratios. Mortality rates rose in response to increasing tree diameter, stand basal area (both from ponderosa pine and non-ponderosa pine), and elevation. Below 1500 m, the mortality rate was ~1%, while above 1700 m, the rate increased to ~30%. However, the association between elevation and mortality risk was buffered by precipitation, such that relatively moist high-elevation stands experienced less mortality than relatively dry high-elevation stands. Tree diameter, crown ratio, and stand density affected tree mortality independent of precipitation. This study demonstrates that while stand characteristics affect tree susceptibility to bark beetles, these relationships may be mediated by climate. Thus, both site and stand level characteristics should be considered when implementing management treatments to reduce bark beetle susceptibility.


Introduction
Disturbances such as tree-killing bark beetles are a driver of tree mortality in North American forests and a concern for forest managers. Across North America, bark beetle outbreaks have resulted in widespread tree mortality in multiple forest types [1,2]. Understanding risk factors for individual trees and stands is important to determine where outbreaks are likely to occur and which trees may be more susceptible. However, the manner in which stand dynamics may interact with the site climate to influence mortality risk in a bark beetle outbreak is poorly understood.
Tree mortality as a result of bark beetle outbreaks is often correlated with site and forest structural characteristics, as well as individual tree characteristics within a stand [3]. Higher elevations, where extreme winter temperatures may limit the survival of bark beetle larvae, have previously been associated with lower tree mortality [4][5][6]. Dense stands [6][7][8][9][10][11] and low species diversity [12][13][14][15][16] may contribute to higher stand susceptibility due to a greater density of suitable host trees and a higher abundance of trees experiencing drought stress from overcrowding [17,18]. Similarly, trees with thinner crowns [19,20] and trees occupying subordinate positions in the canopy [7,21,22] are often more susceptible A subset of data was taken to include ponderosa pines measured within five South Dakota counties covering the Black Hills region (Fall River, Pennington, Custer, Lawrence, and Meade). Data were further subset to include only trees that were alive at initial measurement and were remeasured at least once. Trees were excluded from the analysis  FIA collects a wide range of variables at each fixed-radius plot, including tree diameter, height, canopy position (i.e., dominant, co-dominant, or overtopped), crown ratio (portion of the bole covered by live crown), and others. Trees are tagged to allow for long-term tracking of individuals. Stand-level variables such as live basal area and site index are derived from field measurements. Remeasured trees that died since the last inventory are assigned a cause of death, such as insect, fire, or silvicultural. For our purposes, we assumed that trees assigned a mortality code of "insect" were killed by bark beetles, as bark beetles are the only insect assemblage known to kill Black Hills ponderosa pine in significant numbers [34]. Additionally, all inventory plots occurred in or near known infestations, and this was verified against aerial detection survey data ( Figure 1) [35]. In the FIA plot design in South Dakota, all trees >12.7 cm in diameter at breast height (DBH) are measured on a 7.3 m-radius subplot, while all trees >2.5 cm in DBH are measured on a 2.1 m-radius microplot nested within the subplot. Four subplots/microplots are measured within each plot [32]. More information on FIA field protocol and data is available online at https://www.fia.fs.fed.us/library/, accessed on 13 April 2020. These inventory plots, measured at various points throughout the outbreak, provide an excellent resource to examine the ecological, geographical, and climatological factors influencing tree mortality during a bark beetle outbreak.
A subset of data was taken to include ponderosa pines measured within five South Dakota counties covering the Black Hills region (Fall River, Pennington, Custer, Lawrence, and Meade). Data were further subset to include only trees that were alive at initial measurement and were remeasured at least once. Trees were excluded from the analysis if they were harvested, killed for silvicultural reasons, or had a noninsect mortality agent (fire, weather, animal, etc.). Although most of the damage during our observation period was likely attributable to the mountain pine beetle (Dendroctonus ponderosae Hopkins), these events are often accompanied by outbreaks of secondary and related bark beetle species that co-attack susceptible hosts, and outbreaks of Ips bark beetles were also observed during the study period [36,37]. Therefore, when drawing general conclusions about the environmental factors affecting forest mortality during insect outbreaks, it may be more useful to focus on the combined mortality driven by species assemblages (i.e., both Dendroctonus and Ips spp.) rather than focusing on the effects of individual beetle species.
The 1981-2010 climate averages for each plot location were estimated from the PRISM database [38] at a pixel resolution of 800 m using bilinear interpolation. As the exact coordinates of FIA plots are kept confidential to preserve sampling integrity, we obtained modelled climate averages for each plot location from the Forest Service, which blurred the climate values of a random 10% subset of the data to within 20% of the actual value in order to preserve plot location confidentiality. Climate variables extracted from the PRISM database included average annual temperature and average annual precipitation. While both temperature and precipitation are correlated with elevation in the Black Hills, precipitation is less correlated than temperature, as precipitation is influenced by a general N-S gradient, in addition to an elevational gradient [38]. It is important to note that this study examines long-term regional climate averages, and the impact of drought (negative deviations from typical precipitation patterns) was not examined.

Analysis
Analysis was conducted in R version 3.5.2 [39]. We addressed each of our objectives with mixed-effects logistic regression using the "glmer" function from the "lme4" package [40]. A positive event was defined as a tree that was dead at remeasurement with a recorded cause of death of "insect." A negative event was defined as a tree that survived the entire study period. Each tree appears only once in the dataset-if a tree appeared in multiple sampling windows, these sampling windows were collapsed into one period for analysis. As we were primarily interested in the long-term effects of predictor variables on mortality risk, predictors were assigned based on the conditions that prevailed at the first time the tree was measured-i.e., if a tree was re-measured multiple times, the variables collected at the beginning of the study period were the variables used in the model. We retained individual trees as the unit of observation, which allowed us to parse out the effects of tree-and stand-level variables on mortality hazard. However, this approach necessitated a correction for the nonindependence of trees observed within the same plot. To address this, each model included a random effect of plot. Additionally, each model included a second covariate of total observation time (in years) to adjust for varying lengths of observation between plots. This correction was necessary to account for the higher likelihood of observing mortality during longer observation periods. Models initially included a second random effect of initial measurement year to account for potentially different levels of beetle pressure experienced by each time cohort, but the term did not significantly improve the model fit (Akaike information criterion (AIC) < 2) in any of the models and was consequently dropped. To reduce the impact of collinearity, all predictor variables were mean-centered and scaled by their standard deviation. Within each model, individual variable significance was assessed with a Wald test.
Our first objective-to identify the stand-and tree-level factors contributing to insectrelated mortality-was addressed by testing associations between 11 variables on the probability of insect-related mortality. Stand-level variables included slope ( • ), aspect ( • from NE), elevation (m), stand age (years), host (P. ponderosa) basal area (m 2 /hectare), non host basal area (m 2 /hectare), stand productivity (50-year site index), and heat load index. Heat load index is a composite of aspect, slope, and latitude that approximates the amount of solar radiation received by a particular site. Heat index was calculated using Equation where L = latitude, S = slope, and A = aspect folded about a southwest by northwest axis (all in radians) [41]. Trees located on plots with zero slope were assigned a neutral (90 • from NE) aspect. Tree-level variables included diameter at breast height (cm), crown ratio (% of bole containing live crown), and crown class-a number-based rating describing dominance in the canopy that ranges from 1 (completely open-growing) to 5 (overtopped; suppressed). The initial logistic model included all 11 predictor variables. A final model was obtained through backward stepwise model selection. Nonsignificant variables were eliminated from the saturated model one by one, beginning with the highest p-value (as determined by a Wald test) until only significant (p < 0.05) parameters remained. This selection procedure was necessary because several of our predictor variables were intrinsically correlated with each other, and we wanted to choose the "best" variable among these. While correlated, these variables may influence tree mortality through different mechanistic pathways. For instance, the slope is correlated with heat load index, as heat load index is itself calculated from slope, aspect, and latitude. Trees in stands with higher heat load indices may experience greater heat stress due to their topographic position, and trees on steeper slopes may similarly experience higher heat loads when the aspect is southerly. However, steep slopes also cause water to run off faster, causing trees to be more susceptible to drought stress and mortality [42][43][44]. Thus, while they may be correlated, the mechanisms by which slope and heat load index may influence mortality are subtly different. Consequently, using stepwise regression to identify the best variable (from a pool of correlated variables) can help identify mechanistic pathways in the system. One potential drawback to this approach is that the inclusion of correlated predictor variables in a model may artificially inflate p-values and generate false positives. To account for this, variance inflation factors (VIFs) were calculated for each predictor variable to identify collinearity in the model and verify that predictors were not contributing to artificially inflated p-values.
To address our second objective and determine how climate influenced mortality risk, we ran two additional models. The first tested the effect of average annual precipitation (mm/yr) and average annual temperature ( • C) on mortality risk in isolation. This allowed us to determine if climate, independent of stand-and tree-level factors, was predictive of mortality risk. In the second model, precipitation and temperature were analyzed in conjunction with the tree-and stand-level variables deemed significant in objective one. This allowed us to determine whether climate data were predictive of mortality risk when significant stand-and tree-level factors were accounted for.
To address our third objective and determine whether the relationships between stand/tree properties and mortality risk were mediated by climate, we included an interaction term between each significant stand/tree variable and each significant climate variable. Interaction terms were limited to variables that were proven to be significant in the additive model in order to limit exposure to type I errors. We applied a backward selection procedure to the interaction model identical to that used in Objective 1, eliminating nonsignificant interaction terms one at a time until only significant terms remained. Results reported for each variable include the p-value as determined by a Wald test and the change in log odds resulting from a one-standard-deviation increase in the predictor.
For each objective, statistics are reported for the final model containing significant predictors only. Variance inflation factors were calculated for each parameter from the final models within each objective.

Effect of Stand and Tree Variables on Mortality
First, we sought to determine which tree-and stand-level variables were associated with insect-related mortality risk, independent of climate. Five variables were significantly related to mortality (Figures 2 and 3, Table A2). As these variables were all analyzed in the same model, results can be interpreted as the impact of a given variable controlled for the other variables in the model. At the level of the individual tree, mortality was positively related to diameter (p < 0.001, ∆log odds = 0.599) and negatively related to crown ratio (p = 0.002, ∆log odds = −0.345). Accounting for other factors, the mortality risk for an individual tree roughly doubled for every 11 cm increase in tree diameter and halved for every 33% increase in crown ratio. At the stand level, tree mortality was positively associated with elevation (p < 0.001, ∆log odds = 2.082), host basal area (p < 0.001, ∆log odds = 1.306), and nonhost basal area (p = 0.039 ∆log odds = 0.446). Accounting for other factors, the mortality risk for an individual tree roughly doubled for every 83 m increase in elevation, every 5.3 m 2 /hectare increase in the host basal area of the stand, and every 4.6 m 2 /hectare increase in the nonhost basal area of the stand, though the difference between host vs. nonhost basal area was not statistically significant (p > 0.05).

Effect of Climate on Mortality
Next, we analyzed the impact of climatic variables on mortality, both in isolation and in conjunction with the significant stand/tree variables obtained in Objective 1. In the model containing climatic variables only, average annual temperature was inversely related to insect-caused mortality risk (p < 0.001, ∆log odds = −1.670, Table A3), while average annual precipitation was not correlated with mortality risk (p = 0.683). However, when treeand stand-level variables were included in the model, there was a marginally significant (p = 0.063, ∆log odds = −0.797, Table A4, Figures 2 and 3) inverse relationship between average annual precipitation and mortality risk, while average annual temperature was not correlated (p = 0.271). The significance or directionality of tree/stand variables found to be significant in Objective 1 was not influenced by the addition of climatic variables. every 33% increase in crown ratio. At the stand level, tree mortality was positively associated with elevation (p < 0.001, Δlog odds = 2.082), host basal area (p < 0.001, Δlog odds = 1.306), and nonhost basal area (p = 0.039 Δlog odds = 0.446). Accounting for other factors, the mortality risk for an individual tree roughly doubled for every 83 m increase in elevation, every 5.3 m 2 /hectare increase in the host basal area of the stand, and every 4.6 m 2 /hectare increase in the nonhost basal area of the stand, though the difference between host vs. nonhost basal area was not statistically significant (p > 0.05).

Effect of Climate on Mortality
Next, we analyzed the impact of climatic variables on mortality, both in isolation and in conjunction with the significant stand/tree variables obtained in Objective 1. In the model containing climatic variables only, average annual temperature was inversely re-

Interactions between Climate and Stand Condition
Finally, interactions between average annual precipitation and tree/stand variables were examined to determine if the relationship between tree/stand variables and mortality risk varied across a precipitation gradient. This analysis revealed a significant precipitation × elevation interaction (p = 0.013, ∆log odds = −1.112, Table 1, Figure 4). Elevation remained directly related to higher mortality risk, but the slope of the trendline became shallower with increasing precipitation ( Figure 5). When the precipitation × elevation interaction term was included in the model, the main effect of precipitation was no longer significant (p = 0.291, Table 1). The final model from Objective 3-which included diameter, crown ratio, elevation, host basal area, nonhost basal area, average annual precipitation, and an elevation × precipitation interaction term-was the best fit out of all the models tested ( Table 2). The marginal R 2 for the final model was 0.55. Calculated VIFs for each predictor were <3 (Table A6), indicating that p-values were not being artificially inflated by multicollinearity. Table 1. Output from the final model containing all significant stand-and tree-level variables and significant interactions with climate variables. Prior to modelling, variables were scaled by subtracting the mean and dividing by the standard deviation, so the coefficient can be interpreted as the change in log odds of mortality resulting from a one-standard-deviation increase in the predictor.  Table 1). The final model from Objective 3-which included diameter, crown ratio, elevation, host basal area, nonhost basal area, average annual precipitation, and an elevation × precipitation interaction term-was the best fit out of all the models tested ( Table 2). The marginal R 2 for the final model was 0.55. Calculated VIFs for each predictor were <3 (Table A6), indicating that p-values were not being artificially inflated by multicollinearity. Figure 4. Proportion mortality at high and low elevation, split by precipitation. Higher elevations were associated with higher mortality (p < 0.001), but this effect was reduced in stands receiving a higher average annual precipitation.  . Proportion mortality at high and low elevation, split by precipitation. Higher elevations were associated with higher mortality (p < 0.001), but this effect was reduced in stands receiving a higher average annual precipitation. In relatively dry stands, the risk of mortality rises quickly with elevation, but in relatively moist stands, the elevation-associated increase in mortality risk is less pronounced (precip × elevation, p = 0.0013).

<170>1700
Elevation (m) In relatively dry stands, the risk of mortality rises quickly with elevation, but in relatively moist stands, the elevation-associated increase in mortality risk is less pronounced (precip × elevation, p = 0.0013).

Discussion
Several factors contributed to insect-related tree mortality during the recent bark beetle outbreak in the Black Hills, South Dakota. Although we did find a significant interaction between elevation and precipitation, the effect of stand density, tree size, and crown ratio on mortality was not mediated by the relative moistness of the site. This indicates that while climate may influence some of the factors that predispose trees to beetle-caused mortality, many of those factors are climate-independent.
At the level of the individual tree, larger diameters and smaller crowns were associated with higher mortality rates (Figures 2 and 3). The association between larger trees and higher mortality risk is consistent with past studies [7,9,23,45,46] that have examined mountain pine beetle-related mortality in ponderosa pine. Although rarely examined in the context of beetle-related mortality, crown ratio is an oft-used proxy for tree vigor [19,20,24] and larger crowns have been associated with lower rates of general tree mortality [47], so the association with lower beetle-related mortality is not surprising. Despite past work associating canopy position with mortality risk [7], we did not find a significant link, possibly because canopy position was strongly correlated with diameter (Pearson's coefficient = 0.57), and diameter is critical in quantifying mountain pine beetle susceptibility.
At the stand level, higher elevations and denser stand conditions were associated with higher mortality rates (Figures 2 and 3). The positive relationship observed between live basal area and mortality hazard during a bark beetle outbreak is well-established for bark beetles generally [6,8,10,48] and D. ponderosae specifically [7,9,11,49]. Interestingly, our model revealed that it did not matter whether the live basal area was contributed by ponderosa pine or a nonhost species, such as aspen or white spruce-the addition of live basal from either group contributed equally to higher mortality risk. This indicates that individuals from pure stands were just as likely to experience mortality as individuals from relatively heterogenous, species-diverse stands of a similar density.
Prior research on the impact of nonhost trees on bark beetle outbreaks is mixed. Several studies indicate that host trees located in stands containing higher proportions of nonhost material suffer lower rates of mortality [50,51], an assumption that is maintained in the stand hazard rating developed for D. ponderosae-P. ponderosa in the Black Hills for the Forest Service's National Insect and Disease Risk Map [15]. However, other work from ponderosa pine forests in Colorado indicates that the presence of nonhost material in the stand can enhance mortality risk for host trees [7], which is consistent with our results. While nonhost species cannot amplify pest populations, they still compete with host individuals for moisture, which may explain why mortality risk from our study increased with nonhost basal area. However, conifer forests in the Black Hills are relatively homogeneous, so effects of species diversity on bark beetle-driven tree mortality may not be fully captured in this study.
While we identified a significant main effect of basal area, there was no evidence of a significant interaction with precipitation. This indicates that high-basal area stands are prone to beetle-related mortality across the full gradient of precipitation in the Black Hills, and that high-basal area stands are no more susceptible in dry environments than moist environments within the study region. Past research [52] reasoned that the live basal area could be increased in productive sites without an associated increase in mortality risk from beetles, but our data do not match this prediction, though productivity is influenced by more than just precipitation. Our results may support other hypotheses that dense stands promote beetle-caused mortality through moisture-independent pathways, such as by creating cooler microenvironments conducive to larval beetle development [8] and/or by suppressing winds that could interfere with aerial pheromones [31]. Alternatively, the lack of a precipitation x basal area interaction could be due to the relatively high intensity of the outbreak, which may have led to stands experiencing little drought stress being targeted by bark beetles at the same intensity as drought-stressed stands. The gradient of precipitation in the Black Hills may also not have been broad enough to detect an interaction, and analyses of other regions with different precipitation regimes may produce different results.
The strong relationship observed between higher elevation and increased mortality is one that, to our knowledge, has not been previously documented for mountain pine beetle outbreaks in ponderosa pine. Below 1500 m, the observed mortality rate during the study period was~1%, but above 1700 m, the mortality rate rose to over 30% (Figure 2). While elevation has been shown to influence bark beetle-related mortality in other areas, the effect has typically been in the opposite direction. Past work has documented reduced bark beetle-associated mortality at higher elevations for D. ponderosae in lodgepole pine [4,5] and Ips in ponderosa pine [6], but those studies were conducted at higher elevations than ours. Lower rates of beetle-caused tree mortality at high elevation are thought to be caused by extreme low winter temperatures, which can kill brood developing below the bark, or through increased winter precipitation inputs that improve tree access to water during summer months [4][5][6]. This is unlikely to happen in the Black Hills, as the highest elevations in the Black Hills do not experience low winter temperatures needed to kill larval beetles [5,38]. Extreme highs-and not extreme lows-offer one alternative explanation for the relative dearth of insect-related mortality observed at low elevation in the Black Hills, as the mortality of larval bark beetles has been documented at high temperatures [53]. Other nonclimatic factors, such as soil quality or relative topographic position (i.e., whether a stand is located at the top of a ridge vs. the bottom of a draw) could also provide mechanistic explanations for the higher mortality rates observed at high elevation.
Although observed mortality rates were greatest at high elevation, this effect was buffered by a higher average annual precipitation. As elevation increased, mortality risk increased faster in dry sites than in moist sites (Figures 4 and 5). This finding demonstrates that climate may influence the relationship between stand properties and mortality risk during a bark beetle outbreak. Stand hazard ratings have been developed for a multitude of pest-host-region combinations [12,13,15,50,54], but few models incorporate local climatic variation into their ratings. Our results suggest that such variations may be important in understanding tree susceptibility.
Despite these potential implications, it remains that most of the significant stand-and tree-level variables were not influenced by climate. This indicates that while the inclusion of climate variables could potentially improve current stand hazard ratings, the characteristics of an individual tree or stand are likely to remain the most important determinants for mortality risk. However, the relative importance of stand vs. climate variables (and their interactions) is likely to vary according to geographic scale. The Black Hills represents a relatively small land area, and our plots were separated by no more than 140 km ( Figure 1). If the study area were broadened to include a greater range of climatic variation, climate variables could play a more important role in modulating mortality risk.
Results from this study could help determine which stands are at the highest risk of mortality during a future bark beetle outbreak. While high-value trees, such as those located in campgrounds or tourist areas, can be protected using insecticides or pheromones [55,56], these treatments are typically expensive, and managers must prioritize which areas to treat. Our results show that within a stand, large trees with sparse crowns were most susceptible. Meanwhile, across the landscape, stands that were dry, dense, and at high elevation were most susceptible to insect-related mortality. Thus, treatments to reduce tree density may be most effective in relatively dry stands above 1700 m. Additionally, outbreak mitigation may be most effective in areas with no recent outbreaks, as many susceptible mature trees in post-outbreak stands have presumably been selected out of the population [57].
The scope of this study is limited by several factors. First, our observation window covers only a portion of the latest mountain pine beetle outbreak in the Black Hills. FIA plots in the Black Hills are currently transitioning from a five-to a seven-year resampling period. This means that some plots measured in 2013-2015 will not be remeasured again until 2020-2022, and are not represented in our dataset. The outbreak subsided around 2016, so additional mortality is likely to be observed as plots continue to be remeasured. Second, this is an observational study with many uncontrolled or unobserved variables potentially influencing results. While we report correlations between observed variables, it is difficult to draw causal linkages. Third, our study relies on modelled climate averages for each plot [38], which may have varying levels of accuracy. Additionally, climate averages for each plot were provided by the Forest Service, which fuzzed the value of a subset of plots to within 20% of the actual value in order to preserve plot location confidentiality. Due to this fuzzing, the strength of correlations between climate variables and mortality risk could be underestimated.
Our study highlights the importance of examining both tree-and stand-level variables that could contribute to bark beetle-related mortality. The results from our analysis can help managers identify potentially high-risk stands if another bark beetle outbreak were to occur in the Black Hills. Additionally, our results demonstrate that correlations between stand characteristics and mortality hazard can be strongly mediated by climatic variables. The role that average annual precipitation played in our analysis underscores the potential for bark beetle-driven patterns of tree mortality to be strongly affected by alterations in regional climate, and it will be important to further determine whether these effects were driven by summer or winter precipitation inputs. As precipitation patterns across the western U.S. shift due to climate change, longstanding correlations between stand characteristics and bark beetle susceptibility may change as well. Data Availability Statement: Publicly available datasets were analyzed in this study. This data can be found here: https://apps.fs.usda.gov/fia/datamart/ (accessed on 13 April 2020).