Improving Mountain Pine Beetle Survival Predictions Using Multi-Year Temperatures Across the Western USA

: Global climate change has led to an increase in large-scale bark beetle outbreaks in forests around the world, resulting in signiﬁcant impacts to forest ecosystems, timber economies, and forest-dependent communities. As such, prediction models that utilize temperature for estimating future bark beetle locations and consequential tree mortality are critical for informing forest management decision-making in an attempt to mitigate and adapt to pending and current outbreaks. This is especially true for physiological models that account for the e ﬀ ects of overwinter temperatures on bark beetle survival, as seasonal temperatures, speciﬁcally during winter months, exert the greatest impact on bark beetle mortality during various stages of life cycle development. Yet, how temperature observations are used to predict bark beetle survival can signiﬁcantly under- or over-estimate the role that temperature variability plays in annual tree mortality, especially under current climate change trajectories. This study evaluates how representations of winter temperature inﬂuence bark beetle survival estimates. Using the recent outbreak of mountain pine beetle ( Dendroctonus ponderosae Hopkins) across the western USA as a case study, single-year to decade-long winter temperature averages were used as inputs into a physiological beetle survival prediction model, the results of which were compared against beetle-induced tree mortality observations using temporal autoregressive models. Results show that using longer-term survival averages of seven to ten years signiﬁcantly increases the likelihood that temperature alone can predict general levels of beetle survival and hence beetle-induced mortality. These ﬁndings demonstrate the importance of considering the role of long-term temperature observations when forecasting bark beetle outbreaks, and that year-to-year temperature variability may be constrained in predicting beetle survival during outbreak periods.


Introduction
Climate change impacts have led to a rise in unprecedented native bark beetle epidemics in forests around the globe in recent decades. Warming temperatures and increasing prevalence of drought conditions have facilitated outbreaks of bark beetles across large regions, notably the coniferous forests of western North America [1][2][3][4], and across western and central Europe [5][6][7][8]. These forests are often critical sources of revenue for timber industries and forest-dependent communities, and as such large-scale epidemics put at risk economic and social systems, and drive significant changes to local and regional forest ecosystems [9]. Predictive tools are therefore critical in helping estimate the location and intensity of beetle-induced tree mortality from one year to the next in order to assist with mitigation measures and forest planning.
Of the many regions recently impacted by bark beetle outbreaks, the most widespread and enduring has been in the western United States and Canada, where nearly 50 million hectares of pine forest have been attacked by the mountain pine beetle (Dendroctonus ponderosae Hopkins) since the mid-1990s [10]. Mountain pine beetle (hereafter MPB) attack mostly lodgepole (Pinus contorta Douglas ex Loudon) and ponderosa pine (Pinus ponderosa Dougl. ex Laws.), the two most abundant pine species in these forests. The recent epidemic has impacted forests in eleven western USA states, resulting in millions of lost timber revenues [11], and substantial shifts in forest structure [12][13][14]. Recent studies have confirmed that increased aridity and annual winter temperatures, especially when trends occur for consecutive years, have allowed MPB to survive in areas beyond the species' historical geographic range and at higher altitude [15][16][17].
Like many Scolytidae species, MPB populations are constrained to endemic population levels by cold winter temperatures, which, along with other biotic as well as abiotic factors, result in larval mortality of approximately 98% [18]. However, warmer winter temperatures, especially in the critical early and late winter periods, have led to a decline in larval mortality and a subsequent increase in adult beetles that are available to attack trees during summer months [19,20]. MPB are not tolerant to body tissue freezing, and thus avoid freezing by voiding their digestive tracts in the fall to eliminate ice-nucleating material and by secreting antifreeze compounds in their cells during the winter in response to temperature changes [21]. In response to daily changes in phloem temperature, MPB larvae progressively adapt the concentration of these antifreeze compounds to adjust the supercooling point, above which ice crystal nucleation cannot occur. Since MPB's supercooling point tracks daily temperature changes, rapid temperature fluctuations such as sudden late fall and early spring cold snaps can push temperatures below the current supercooling point, resulting in ice crystallization and high mortality [21]. MPB larvae can potentially reduce their supercooling point to −40 • C prolonged winter temperatures below this value have been observed to result in complete MPB mortality [22].
Observed relationships between MPB physiology and winter temperatures have led to the development of an array of MPB prediction models in which seasonal temperature variability plays a central role in estimating MPB population levels and consequential tree mortality [23][24][25][26][27][28][29]. MPB frost resistance continually responds to winter temperature fluctuations. As such, utilizing daily temperature observations or simulations is appropriate for predicting daily survival levels. However, using such observations for a single-year, as most models do, can potentially mask the effects of temperature trends on beetle survival rates. For example, multiple winters of above-normal temperatures can facilitate MPB population growth from endemic to epidemic levels. At higher population densities, MPB can successfully mass-attack larger, healthier trees which provide more thermal protection than small trees [30]. Thus, early or late season cold snaps as well as enduring mid-season cold temperatures could have a dampened effect on mortality of larval MPB in large trees [18]. Therefore, longer term temperature averages would implicitly capture multi-year trends of reduced morality in successive warm winters, especially for populations in incipient-epidemic or epidemic phases. This is especially pertinent in the majority of the MPB species range in the USA where seasonal temperature variability is relatively high compared to surrounding areas. Figure 1 illustrates the mean minimum winter temperature from 1980 to 2016, temperature variability during this period (as shown by standard deviation of temperatures at each location over the same time period), and total tree mortality. It can be noted that most areas with highest tree mortality experience some of the coldest winter temperatures but also the highest year-to-year variability. This begs the question: Do multi-year temperature averages provide improved predictive capability for bark beetle-induced tree mortality than single-year temperatures? If the answer to this question is yes, many existing physiological bark beetle models, and insect models in general, would benefit from considering how such measures are influencing the prevalence of insects in forests over time. The objective of this study is to provide insight into how utilizing multiple years of temperature data can potentially improve the ability of models to predict overwinter beetle survival. To accomplish this objective, we compare single-year versus multi-year estimates of temperature-mediated MPB survival as predictors of MPB-induced tree mortality. Focusing efforts on the eleven western USA states, we utilize an existing physiological model that estimates tree mortality as a function of overwinter MPB survival that is calculated by inferring how daily winter temperatures impact MPB mortality [31]. We compare single-year and multi-year estimates of MPB survival from this model for predicting tree mortality using an autoregressive integrated moving average (ARIMA) modeling approach, the outputs from which allow us to determine the degree to which using longer time periods of MPB survival improve model predictive capabilities. The goal here is not to evaluate the Régnière and Bentz [31] per se, but rather to understand how the results of this model, and potentially others, change when utilizing multiple years of temperature data for estimating overwinter beetle survival.

Study Area and Data
The study area covers the native range of P. ponderosa and P. contorta within national forests of the eleven western states in the USA ( Figure 2). These species are dominant in low-to mid-elevation forests in this region, and are the main hosts for MPB. We acquired interpolated daily minimum and maximum temperature data covering the study area for the years 1980-2016 from the Daily Surface Weather And Climatological Summaries (DAYMET model version 3) [32]. The DAYMET data are distributed as gridded tiles in a Lambert Conformal Conic projection with 1 km × 1 km grid cells with a coverage area consisting of North America and a portion of western Greenland.
Tree mortality data for the years 1997-2010 were obtained from Meddens and Hicke [33], which contains information on beetle-killed trees from the western USA and British Columbia. The data consist of gridded layers for each MPB host species, transformed from spatial polygons that were originally derived from regularized photographic aerial survey data provided through the United States Department of Agriculture's Forest Service Pest Portal [34]. Raster cell values represent the estimated number of MPB-killed trees in that cell. For analysis, we created a composite raster layer from the layers for P. ponderosa and P. contorta var. latifolia in which cell values represented the sum of trees killed of both species.
To determine host species range, we used spatial polygon data for individual tree species derived from Little Tree Atlas [35], while spatial data for national forest boundaries were obtained from the United States Forest Service GeoData Clearinghouse [36]. The pine native range was calculated as the spatial union of the polygons for P. ponderosa and P. contorta. To define the final study area, we calculated the spatial intersection of the native pine range with the national forest boundaries ( Figure 2).
The study site spatial polygons were used as a mask for the gridded data set in order to determine which cells to use for analysis, resulting in a study area of approximately 460,000 square kilometers.

MPB Overwinter Survival Model
To estimate the potential MPB overwinter survival rates, we used a physiological process-based model of MPB freeze avoidance created by Régnière and Bentz [31]. The model uses daily temperature data to estimate changes in physiological state of MPB populations via the gain and loss of cold hardening, i.e., the capacity of MPB larvae to avoid ice crystal nucleation in temperatures below the normal freezing point of water. In the model, MPB larvae exist in three possible states. State 1 (no cold hardening) is a summer feeding state in which they are highly susceptible to freezing. In State 2, first phase of cold hardening, the MPB larvae have voided ice nucleating material from their gut. In State 3, winter high cold hardening, MPB larvae tissues contain antifreeze compounds. Each state is characterized by a symmetrical distribution of individual beetles' supercooling points (SCP) below which ice crystals form inside the larvae, causing mortality. The MPB population is characterized by the proportion of larvae in each of the states, from which an overall population SCP is calculated. The population-level SCP represents the temperature below which 50% of the population experience freeze-related mortality. Changes of the population proportions of each state are mediated by daily gains or losses of cold hardening in response to temperature changes. The model structure introduces a lag between changes in temperature and changes in population proportions. The time structure of the model allows it to respond to absolute low temperatures as well as to differentiate between slow and rapid cooling events. Unexpected cold snaps can result in high mortality if the temperature suddenly falls below the SCP, while gradual cooling allows the SCP to decrease, resulting in lower MPB mortality.
We implemented a spatially-explicit version of this model in the Java programming language, version 1.8 (Oracle, Redwood City, CA, USA). The performance of our model implementation was validated (Table 1) by comparing model outputs using daily temperature data from the weather stations to those found in Table 1 of Régnière and Bentz [31]. For all weather stations, our model survival predictions differed by no more than 9% from those provided in Table 2

Analysis
Yearly tree mortalities and modeled overwinter survival for all grid cells were averaged for the entire study area, resulting in a single value for each variable for each year. The tree mortality data reflects trees in the red needle stage, which occurs in the summer following the seasons when MPB-induced tree mortality takes place. For example, we defined red-stage pines observed in the summer of 2001 as having a 0-year time lag from the 2000-2001 winter and a 1-year time lag from winter 1999-2000. Next, we evaluated the relationship between tree mortality for a given year, t 0 , to MPB survival for year t 0 , the mean MPB survival for t 0 and t −1 , the mean MPB survival for t 0 , t −1 , and t −2 , and so forth up until the mean MPB survival for t 0 , t −1 , t −2 , . . . t −16 resulting in mean estimated survival or periods of 1-16 years.
Linear regression models were fit to the data using mean yearly MPB survival over all raster cells as a predictor for the logarithm of pine mortality. A true mechanistic description of MPB survival and tree mortality may be highly nonlinear and would include terms representing factors such as forest density and initial MPB population density, for which datasets with continental spatial extents do not exist. We have opted for a purely descriptive, phenomenological approach in order to use datasets with such a large spatial extent. For the purpose of meeting the assumptions implicit in linear Ordinary Least Squares (OLS) models, a logarithmic transformation of the response (tree mortality) was used to linearize the relationship and ameliorate heteroscedasticity.
Due to the possibility of temporal autocorrelation, i.e., observations in one year being correlated to observations in the preceding years, we used autoregressive integrated moving average (ARIMA) models to account for possible autocorrelation in the model residuals. Two candidate sets of models were inspected and used to select ARIMA model structures for analyses. The first were built using an autoregressive order of 1 year (AR1 correlation structure). Next, a set of models were built using a procedure that searches among different correlation structures to find an optimal model for a data set using the Akaike Information Criterion, the auto.arima() function in package forecast in R (R Core Team, Vienna, Austria). The second set included models with different correlation structures. A test of normality of residuals (shapiro.test()) revealed that a single AR1 model and 6 auto fitted models had significantly non-normal residuals (p < 0.05). We used a set of models with AR1 structures because we felt that the Akaike Information Criterion (AIC) was not the optimal model selection criterion for comparing models with different input. We also felt that using identical correlation structures in all models would facilitate model comparisons. All scripts for conducting analysis on model outputs are available here: https://github.com/michaelfrancenelson/Process_MPB_Overwinter_Spatial_Data.

Results
Separate models were fitted for mean MPB survival for periods of 1 to 16 years prior to the year that MPB-induced tree mortality was observed. Scatterplots of MPB survival and tree mortality, log-transformed to linearize and reduce heteroscedasticity, show an increasing convergence to a linear relationship as the length of the averaging period increases. Comparing the regression results ( Figure 3) revealed that single-year and short-term survival averages (<7 years) are not significant predictors of tree mortality, while models with longer-term mean survivals (≥7 years) are significant predictors of tree mortality. Shorter-term models all had low slope coefficients and insignificant p-values, and slope coefficients generally increased in the 1-year to 8-year models then gradually plateaued. The model goodness-of-fit measure, root mean square error (RMSE), was also higher for the short-term models and particularly low for models with 9-year and 16-year temperature means. The model coefficients ranged from approximately 0.52 to 0.57 with 8-year to 16-year temperature means ( Figure 4). Since the response was log-transformed, the coefficient is interpreted in terms of the percent increase in pine mortality associated with a 1% increase in MPB survival. For the longer-term models, the coefficients corresponded to approximately a 70% increase in pine mortality with 1% survival increase, 200% increase in mortality with 2% increase in survival, and 1500% increase in kills with a 5% increase in mean survival.  Figure 5). The longer-term average survival corresponds more strongly with high tree mortality areas, but there are some regions where low estimated beetle survival co-occurs with high tree mortality.

Discussion
A substantial body of research demonstrates that increasing temperatures have a direct impact on MPB survival and consequential tree mortality in the pine forests of western North America [1,3,9,16,19,20,37]. Models for predicting year-to-year changes in MPB population, both in terms of spatial extent and population levels, typically rely on single-year overwinter temperatures to estimate MPB survival and the potential impact on tree mortality. Yet, we find that multi-year MPB survival averages, calculated by way of daily temperature changes, provide improved predictive capability than single-year measures at our scale and resolution of analysis. This is particularly evident from the decrease in RMSE, a measure of model goodness-of-fit, averages of seven or more years are used in the analysis. Furthermore, when longer-term MPB survival means are used, coefficient values are notably higher, indicating that small temperature changes over multiple years can lead to a substantially increased level of tree mortality.
The results from this study raise questions as to why single-year MPB survival estimates, or estimates based on averages involving relatively few years of MPB survival, are poorer predictors of tree mortality despite overwhelming evidence of heightened tree mortality with increasing temperature. One potential explanation is that it is difficult to know the extent to which annual temperatures impact MPB populations. For one, the magnitude and abruptness of sudden temperature drops in early and late spring, when MPB larvae have lower levels of antifreeze compounds, required for significant mortality in MPB populations remain uncertain. Similarly, the duration of mid-winter cold snaps that are required to have a negative impact on MPB populations remains vague. Therefore, applying multiple years of survival predictions are more likely to be reliable for predicting tree mortality, because unusually cold winters or spring/fall cold snaps that occur several years in a row may have a much greater impact on the population-level beetle mortality than a single unfavorable year. Likewise, for areas with initially low MPB population densities, such as populations in the endemic phase, consecutive years with favorable temperatures may be necessary for large increases in population densities. This study, therefore, has found that longer term survival averages are more suitable for capturing such trends.
Another challenge with using single-year survival estimates is that it is difficult to associate specific temperatures with MPB over-winter survival across different beetle population phases. At low endemic populations, MPB are more susceptible to cold temperatures and rapid temperature changes as they typically reside in smaller trees with compromised health, such as those affected by disease or environmental constraints (e.g., localized drought conditions) [38,39]. These compromised trees have thinner cambium layers, providing a poor food source, and thinner bark, which provides only a weak buffer against rapid temperature changes. Thus, single-year MPB survival models may potentially overestimate tree mortality at low MPB populations because only compromised trees are at risk of being attacked. Alternatively, epidemic MPB populations attack larger diameter trees, which provide better shelter and food resources through thicker bark and cambium that improve MPB vigor while moving through various life stages in the tree over the winter. As such, when MPB populations are in an epidemic state, single-year survival estimates may underestimate survival because MPB can better withstand cold temperatures due to the protection offered by larger trees.
The issue of spatial scale is important to consider in studies covering a large spatial extent. In some areas, local-or regional-scale patterns may not reflect a continental-wide average. This can be observed in Figure 5 where some regions, for example in central Montana, with low 10-year mean estimated MPB survival overlapped with high tree mortality while the continent-wide average trend showed an increase in tree mortality with increased MPB survival. Thus, a clear trend at a continental spatial extent may be a poor tool for explaining or prediction patterns in individual forests. Such local-scale results may be more highly influenced by factors, such as local trends in MPB populations, not detected in averages at broader scales. Another local-scale factor not considered in our study is the presence of physical stand characteristics that affect the amount of thermal buffering in the phloem. Such differences may cause underestimation of potential larval survival in dense forest stands in which favorable conditions for MPB epidemics occur. The model of Regniere and Bentz [31] that we used to estimate MPB survival uses a simple adjustment to account for differences between air and cambium temperatures. For more detailed use at smaller spatial scales, this model could be expanded to include factors known to influence temperature discrepancies such as the physical characteristics of individual tree stands [40].
The findings from this study result from aggregating the relationship between MPB survival and tree mortality across all locations in the study area. It is likely that spatial variation exists in the accuracy of the prediction, as some locations may not demonstrate significance with multi-year averages, while others show significant relationships. It is therefore important not to assume that the findings are applicable across the MPB species range, and further work is needed to examine how multi-year estimate perform among latitudinal gradients, different tree species, elevation variation, and other variables important in the MPB-host tree relationship. Yet, this study has shown that it is critical to consider the role of multi-year temperature averages, particularly under processes of climate change, in determining the potential for MPB-induced tree mortality.
Unlike several other predictive models, we did not include MPB population pressure (i.e., the location of existing beetles) in our estimates of tree mortality. Population pressure is often used to represent the presence of MPB in the landscape, such that higher levels of estimated population would result in higher levels of MPB overwinter survival and subsequent tree mortality. Instead, we were most interested in how temperature-mediated MPB survival alone influences overall tree mortality predictions, and thus did not consider beetle population levels in our predictive modelling.
One challenge with predicting MPB survival and annual tree mortality is the presence of temporal autocorrelation, which is evident in the plots in Figure 4. For most plots, particularly those demonstrating a significant fit between the two variables (i.e., the 6-year to 10-year averages), the location of points appears dependent on the particular year for which the data was predicted and observed. That is, the first year of the study, 1997, represents low beetle survival and low tree mortality, while later years, particularly 2008 to 2010, demonstrate high values for both variables, and the remaining years occupy a position in the plot in near chronological order. The violation of the assumption of independence of observations poses a challenge for statistical analyses, which can be addressed with analyses that account for serial correlation in model residuals, such as the AR1 models.
Insect outbreak prediction models are an essential tool for assisting forest management with mitigating and adapting to insect epidemics, especially considering predictions of future climate change that will leave an increasing area of the world's forests susceptible to insect outbreaks. As regions enter into different climate trajectories including prolonged periods of above normal temperatures and drought conditions, it is important to evaluate how current models and model parameterization based on conventional knowledge of insect physiology and insect-host interactions perform under such emerging conditions.

Conclusions
This study evaluates how representations of winter temperature influence bark beetle survival estimates. The results from our study show that using longer-term survival averages of seven to ten years significantly increases the likelihood that temperature can predict general levels of beetle survival and hence beetle-induced mortality. These findings demonstrate the importance of considering the role of long-term temperature observations when forecasting bark beetle outbreaks, and that year-to-year temperature variability may be constrained in predicting beetle survival during outbreak periods.
Author Contributions: C.B. was the principal investigator on the project funding this work. His role on this manuscript was developing research context for this study, assisting M.F.N. with analyzing model outputs, and writing the manuscript. M.F.N.'s role on this manuscript was to retrieve spatial data, from multiple sources, implement the MPB overwintering model, and analyze model outputs.