Neighborhood context mediates probability of host tree mortality in a severe bark beetle outbreak

. Understanding drivers of disturbances across scales is critical as environmental constraints change in a warming climate. Outbreaks of native bark beetles (Curculionidae: Scolytinae) are key natural disturbances that shape the structure and function of conifer forests across the northern hemisphere. While drivers of bark beetle outbreaks have been studied extensively at spatial scales ranging from stands to con-tinents, within-stand processes governing individual tree mortality in an outbreak are less well under-stood. Here, we use a spatially explicit long-term monitoring dataset of a lodgepole pine ( Pinus contorta var. latifolia ) forest ( > 9000 individually mapped trees in three 2-ha plots) impacted by a severe mountain pine beetle ( Dendroctonus ponderosae ) outbreak to explore interactions among ﬁ ne scale drivers of beetle-caused tree mortality. Using a Bayesian spatial modeling approach, we evaluated how tree scale and tree neighborhood scale characteristics interact with tree size to mediate host tree susceptibility to mountain pine beetle outbreak in the Southern Rocky Mountains (USA). We found evidence that both tree growth rate preceding the outbreak and neighborhood structure (within a 10 meter radius of the host tree) mediate the effect of tree size, and that the direction and magnitude of these mediating effects vary with tree size. Tree scale mortality probability increased with pre-outbreak growth rate for small to medium sized host trees ( ~ 10 – 25 cm diameter), but that same effect was not detected for large trees. Conversely, tree scale mortality probability increased with greater neighborhood density, with the most pronounced effects for medium to large sized host trees ( ~ 15 – 30 cm diameter). Within-stand topographic variability was not an important predictor of mortality probability; among stands, however, the stand in the driest topographic position experienced the greatest overall mortality. By explicitly considering how within-stand heterogeneity mediates individual tree scale susceptibility to bark beetle outbreak, our ﬁ ndings bridge an important gap in understanding multi-scale drivers of disturbance dynamics.


INTRODUCTION
Native bark beetles (Curculionidae: Scolytinae) are one of the primary disturbance agents affecting temperate forests worldwide (Kautz et al. 2017, Sommerfeld et al. 2018. In North America, historical (i.e., 19th century; Baker andVeblen 1990, Jarvis andKulakowski 2015) and contemporary (i.e., 1997Meddens et al. 2012) bark beetle outbreaks have impacted millions of hectares of forests, resulting in widespread tree mortality across multiple tree species. While such outbreaks are naturally occurring disturbances, climate warming is projected to alter outbreak dynamics (Bentz et al. 2010, Seidl et al. 2017, with important implications for forest structure and the ecosystem services that forests provide (Raffa et al. 2008, Thom and Seidl 2016, Seidl et al. 2016b, Negrón and Cain 2019. Outbreak patterns and processes have been studied extensively at broad scales (continental to regional; e.g., Meddens et al. 2012) and at the among-stand scale (aggregating all data to the stand; e.g., Simard et al. 2012). However, stand averages do not capture fine scale variability in stand conditions that may mediate tree scale susceptibility to bark beetle outbreak in important ways (Negrón et al. 2001). Interactions among factors that govern tree mortality at the tree scale and in the immediate environment of the tree are a key knowledge gap (Das et al. 2008), with clear implications for understanding how bark beetle outbreaks alter stand structure and thus dynamics of subsequent disturbances (e.g., Kashian et al. 2011).
Tree mortality often results from a complex combination of predisposing, inciting, and contributing factors (Manion 1981, Franklin et al. 1987. In the context of bark beetle outbreaks, long-term stress induced by less favorable edaphic or topographic position can predispose trees to mortality, and short-term stress induced by drought can incite mortality events through changes to both beetle demography and host tree physiology (McDowell et al. 2008, Raffa et al. 2008, Gaylord et al. 2013. Feedbacks across scales (e.g., host tree availability and susceptibility, beetle density, weather, and/or climate) can periodically trigger extensive beetle outbreaks (Raffa et al. 2008), which then act as the contributing factor resulting in tree death.
Tree characteristics (e.g., tree size and defensive capabilities) are key determinants of whether a host tree is attacked by bark beetles and if so, whether it survives. Large diameter host trees with thick phloem are preferred by beetles, as these trees increase the potential for beetles to reproduce and survive (Amman 1969). In non-outbreak conditions (beetle populations at low, endemic levels), typically only weakened or damaged host trees are successfully attacked and killed by beetles, while healthy trees can defend against an attack and survive. During bark beetle outbreaks, host tree defenses (e.g., resin defense traits; Ferrenberg et al. 2014) can be overwhelmed by high-density beetle populations capable of pheromone-mediated aggregation and mass attack (Raffa and Berryman 1983). With beetle populations unconstrained by tree defenses (Boone et al. 2011), host trees that are larger in size when outbreaks occur are preferentially mass attacked. Drought-induced stresses on tree physiology (e.g., hydraulic failure and carbon starvation) can both amplify and be amplified by bark beetle herbivory, with complex interactions among these factors contributing to host tree mortality when droughts and outbreaks coincide (McDowell et al. 2008, Gaylord et al. 2013. The preferential attack of large host trees by beetles, in combination with the greater susceptibility of large trees to hydraulic failure (Bennett et al. 2015), means that larger host trees experience the highest rates of mortality during bark beetle outbreak (e.g., Hopping and Beall 1948, Furniss 1965, Schmid and Hinds 1974, Amman and Cole 1983, Björklund and Lindgren 2009). Bark beetle herbivory may result in selective pressure for host trees to grow slowly and invest more in defenses at mature stages, counteracting selective pressure to grow quickly and invest less in defenses at juvenile stages and creating potential life-history trade-offs for host trees (de la Mata et al. 2017).
Spatial attributes of the neighborhood surrounding a host tree, in addition to characteristics of the tree itself, can also influence susceptibility to bark beetle attack (Bakaj et al. 2016). For example, while large diameter host trees are preferentially attacked during an outbreak, small diameter host trees may also be attacked when in close proximity to large diameter host trees (Preisler andMitchell 1993, Kashian et al. 2011). Further, prevalence and basal area of beetle-killed host trees may be greater in denser neighborhoods or clusters within a stand (Olsen et al. 1996, Negrón et al. 2001. Fine scale neighborhood effects are expected to be secondary to tree scale characteristics in determining host tree susceptibility to bark beetle outbreak (Bakaj et al. 2016). However, it remains unclear how tree neighborhood characteristics might mediate individual host tree susceptibility, and whether these effects vary with tree size. At watershed scales, outbreak severity can be driven by variation in topographic position (e.g., mesic or xeric v www.esajournals.org 2 August 2020 v Volume 11(8) v Article e03236 conditions; Vorster et al. 2017, Cooper et al. 2018. However, whether fine scale (e.g., submeter) variation in topography affects tree susceptibility has not been tested.
Here, we used a spatially explicit long-term monitoring dataset of a lodgepole pine (Pinus contorta var. latifolia) forest affected by a severe mountain pine beetle (MPB; Dendroctonus ponderosae) outbreak (in Fraser Experimental Forest, Colorado, USA) to characterize within-stand drivers of individual lodgepole pine mortality from bark beetles. MPB is the most widespread biotic agent of tree mortality in western North American pine forests, affecting tens of millions of hectares of forest during episodic outbreaks Carroll 2006, Graham et al. 2016). We studied three 2-ha plots in Fraser Experimental Forest in which every tree ≥5 cm in diameter (9357 trees) was censused, measured, and mapped leading up to (1989,2004) and following (2018) a severe MPB outbreak to address the following research question: How do tree scale (pre-outbreak growth) and tree neighborhood scale (neighborhood structure and topographic position within 10 meter radii) characteristics interact with tree size to mediate individual host tree susceptibility to bark beetle outbreak?
We expected that host tree diameter would be the primary driver of mortality probability, and that additional tree and tree neighborhood scale characteristics would interact with tree size to mediate mortality probability. At the tree scale, growth-defense trade-offs may result in a positive (de la Mata et al. 2017) or negligible (Boone et al. 2011) relationship between pre-outbreak growth and host tree mortality. At the tree neighborhood scale, we expected mortality probability would be greater in highdensity neighborhoods due to more favorable microclimate for MPB (e.g., temperature and wind speed; Bartos and Amman 1989) or preference of MPB for denser stands (Negrón 2019). We also expected mortality would be greater in tree neighborhoods dominated by lodgepole pine due to high beetle pressure resulting from the emergence of adult MPB from neighboring colonized trees. We expected that fine scale topographic position would further mediate mortality probability since topoedaphic outbreak refugia (i.e., locations where outbreak impacts are less severe relative to the surrounding area) exist among stands (Vorster et al. 2017, Cartwright 2018).

Study area
Fraser Experimental Forest is located within the USA Southern Rockies Ecoregion, approximately 80 kilometers northwest of Denver, Colorado, in the Arapaho National Forest (Appendix S1). The study stands are located in subalpine forest composed primarily of lodgepole pine and secondarily of subalpine fir (Abies lasiocarpa) and Engelmann spruce (Picea engelmannii), with elevation ranging from 2790 to 2970 m (Appendices S1 and S2). Stands are >300 years old with no known major disturbance since the stand-initiating fire

Sampling design
Three 2-ha plots separated by a minimum distance of 580 m were established in 1938 (Wilm and Dunford 1948), at which time a complete census of every living tree measuring at least 9 cm in diameter at breast height (DBH; measured 140 cm above the ground) was conducted. Subsequent surveys were conducted in 1940, 1946, 1960, 1989, 2004, and 2018. DBH was recorded for each live tree meeting the minimum diameter threshold, and each dead tree that was live at the time of the previous survey was examined for evidence of mortality agents. In 2004, the minimum diameter threshold was lowered to 5 cm. We conducted the first post-outbreak recensus of the plots in July and August 2018 and mapped the locations of all live and dead v www.esajournals.org surveyed trees (n = 9357). Stem-mapping was conducted using FieldMap (Hédl et al. 2009, www.fieldmap.cz), a mobile unit that combines a laser rangefinder (TruPulse 360, Laser Technology, USA) with geographic information system software. Standing trees were mapped in place, and uprooted trees were mapped at their previous rooted positions. Each tree was mapped with sub-meter accuracy to its horizontal coordinates (x, y) and its local elevation (z) relative to other trees.
While some mortality from MPB occurring prior to 2004 was noted in the 2004 survey, 90% of the MPB-caused tree mortality in the plots occurred after the 2004 survey. Therefore, we focused our analysis on the mortality that occurred between 2004 and 2018. We assumed that the survey conducted in 2004 provides a reasonable census of the potentially susceptible population of lodgepole pine (n = 5011), as well as the stand structure within each plot, at the start of the outbreak. The post-outbreak survey conducted in 2018 provides the individual tree scale mortality outcomes [killed by MPB or not killed by MPB (i.e., alive in 2018 or killed by another agent)] within each plot (Fig. 1). Dead trees were only considered killed by MPB if MPB adult and larval galleries in the vascular cambium could be positively identified, following established methods (Safranyik and Carroll 2006). MPB was considered the proximate cause of tree mortality in these cases, as their colonization and reproduction essentially girdles the cambium and outer xylem of host trees (Diskin et al. 2011, Pelz and Smith 2012, Simard et al. 2012).

Tree scale covariates
Tree size at the start of the outbreak was quantified using DBH recorded in the 2004 survey. Growth rate preceding the outbreak was quantified as mean annualized basal area increment (BAI) estimated from DBH measurements recorded in 1989 and 2004. While tree growth is expected to follow a sigmoidal curve both with increasing age and size, using BAI rather than radial growth minimizes the effect of increasing circumference (Weiner et al. 2001, Sullivan et al. 2016). We found a moderate positive correlation between DBH and BAI (Pearson r = 0.5) in the dataset. However, variance inflation factors for all covariates included our analysis were <2, suggesting there was not a high level of collinearity in the dataset (Zuur et al. 2010). BAI estimates are below zero in some cases (7% of trees; Table 1), which we attribute to measurement error from slight differences in positioning of diameter tapes among surveys. We made no adjustment to the BAI estimates to account for measurement error, as we assume these errors would have occurred in both directions and did not want to introduce bias to the data.

Neighborhood scale covariates
Within-stand neighborhood density was quantified using the sum of horizontal angles originating from each focal tree center and spanning the DBH of each neighbor tree. This neighborhood density index (NDI) is calculated within a predetermined neighborhood radius as follows: v www.esajournals.org where neighboring trees are indexed i = 1, . . ., n; DBH i is the DBH of neighboring tree i and distance i is the distance between the center of the focal tree and the center of neighboring tree i (Appendix S3). NDI increases with increasing large, numerous, and/or close proximity trees in the neighborhood of the focal tree. While NDI was originally proposed as a tree competition index to characterize tree vigor (Rouvinen and Kuuluvainen 1997), it also effectively serves as an index of relative stand density at the position inhabited by each tree. NDI was calculated using all neighboring trees, regardless of species. To quantify the relative density of host trees, we also calculated the proportion of total NDI contributed by lodgepole pine (pNDI).
Local topographic position was quantified using a topographic position index (TPI) that compares the elevation of each cell in a digital elevation model to the mean elevation of a specified neighborhood around that cell (Weiss 2001). TPI is calculated as follows: where z 0 is the elevation of the cell containing the focal tree, and z is the average elevation around the focal tree within a predetermined radius. TPI is positive (convex local topography) when the focal tree is rooted in a position higher than its neighbors and negative (concave local topography) when the focal tree is rooted in a position lower than its neighbors (Appendix S2). We calculated TPI using a 1-m resolution digital elevation model developed from the local elevation of each individual tree measured in the field.
Calculating neighborhood covariates requires specifying the radius of the neighborhood that is considered, which should ideally be informed by the ecological system. When characterizing tree competition effects, using a neighborhood radius approximately 3.5 times the average radius of tree crowns is appropriate (Contreras et al. 2011). We estimated that 10 m would be an appropriate neighborhood radius, as mature lodgepole pine crowns range from 5 to 8 m in diameter (Herman et al. 1996). To avoid boundary effects, trees located <10 m from the edge of each plot were excluded from analysis. To test the sensitivity of neighborhood size, 5 and 15 m neighborhoods were also considered and model results were compared.

Statistical analyses
To make inference regarding the potential influence of covariates on MPB-induced Notes: Dataset includes lodgepole pine ≥5 cm in diameter at breast height that were alive in 2004 (pre-outbreak) and re-surveyed in 2018 (post-outbreak). Trees without pre-outbreak growth data (i.e., trees that were not surveyed in 1989) and trees <10 m from the edge of each plot were excluded from the dataset.
v www.esajournals.org lodgepole pine mortality, we specified the following hierarchical spatial model. The binary response (1 = killed by MPB, 0 = not killed by MPB) for lodgepole pine tree i (i = 1, . . ., n) was modeled as a Bernoulli random variable Y i measured at location s i . The model takes the following form: where x ik (k = 1, . . ., K) are covariates, and u i is a spatial random effect modeled as a continuous Gaussian random field (GRF) measured at location s i , u i ðs i Þ ∼ GRFð0,∑Þ: We modeled the covariance matrix Σ using a Matérn covariance function parameterized by marginal standard deviation σ u and practical range r (the distance at which the correlation drops to approximately 0.1; Lindgren et al. 2011). We included tree scale covariates (DBH, BAI; Table 1), neighborhood scale covariates (NDI, pNDI, TPI; Table 1), and plot (B2, C3, D2) as fixed effects in the model. Because we expected that effects might vary with tree size, we also included interactions between DBH and all other tree and neighborhood scale covariates. To enable comparison of the magnitude of model coefficients, we standardized all covariates by subtracting their means and dividing by two times their standard deviations (Gelman 2008).
We used the integrated nested Laplace approximation (INLA) approach to implement this spatial model in a Bayesian framework. INLA is a deterministic algorithm developed to approximate posterior marginal distributions of model parameters (Rue et al. 2009) that offers substantial computational advantages when fitting spatial models (Blangiardo and Cameletti 2015). To estimate the GRF at irregularly spaced point locations, we used the stochastic partial differential equation (SPDE) approach, which approximates the continuous GRF as a discrete Gaussian Markov random field (Lindgren et al. 2011). The SPDE approach models the Gaussian Markov random field at the nodes of a triangulated grid (Appendix S4) and interpolates the Gaussian Markov random field throughout the study domain (Lindgren et al. 2011, Blangiardo andCameletti 2015).
Implementing the model in a Bayesian framework requires setting prior distributions for the model fixed effects and hyperparameters of the spatial random effect. We defined uninformative uniform and normal priors for the intercept and fixed effects, respectively. We set priors for the marginal standard deviation and practical range of the GRF using a penalized complexity prior approach, which avoids spatial overfitting by shrinking the marginal variance toward zero (Fuglstad et al. 2018). The priors were defined as follows: Pðσ u >0:1Þ ¼ 0:01 The prior for the practical range r was set based on the a priori hypothesis that r might be close to 30 m, the typical dispersal distance of MPB in lodgepole pine forests (Safranyik et al. 1992). Both the INLA and SPDE approaches were implemented in R (R Development Core Team 2019) using the R-INLA package (www.r-inla.org).
We evaluated model fit via the conditional predictive ordinate (CPO; Pettit 1990), which is a leave one out cross-validation score. CPO is calculated for each observation Y i by omitting Y i , fitting the model to all remaining observations Y −i , and then using the fitted model to predict the probability of Y i . The sum of the log of the CPO values is a useful summary of model fit, with larger sums indicating a better fitting model (Blangiardo and Cameletti 2015). Since our primary goal was to evaluate whether within-stand spatial context mediates MPB-induced lodgepole pine mortality, we evaluated the fit of our model against the fit of models (1) excluding neighborhood covariates, (2) excluding the random spatial effect, and (3) excluding both neighborhood covariates and the random spatial effect.
A covariate or interaction term x k was considered an important predictor in the model if the 95% credible interval for the coefficient β k did not include zero and a marginally important predictor in the model if the 90% credible interval for the coefficient β k did not contain zero. To illustrate how both the size and position of a tree with a stand might affect its probability of mortality due to a severe MPB outbreak, we used the model to predict and produce maps of mortality probability across each plot for a range of tree sizes.

RESULTS
Across all plots, 54% of lodgepole pine tree stems and 82% of lodgepole pine basal area live in 2004 were killed by MPB by 2018. At the individual tree scale, tree size (i.e., DBH) was the strongest predictor of MPB-induced lodgepole pine mortality (Fig. 2). Predicted mortality probability increased sharply for trees >15 cm DBH and the effect of DBH on mortality probability was approximately five times greater than other effects (Figs. 2, 3; Appendix S5). The mean DBH of lodgepole pine killed by MPB was 27.1 cm (SD = 7.3 cm), and the mean DBH of lodgepole pine not killed by MPB was 13.0 cm (SD = 5.9 cm). BAI was a marginally important predictor in the model ( Fig. 2; Appendix S5). Greater BAI increased predicted mortality probability for small to medium diameter host trees (10-25 cm DBH), while having a negligible effect for large diameter trees (Fig. 3).
At the tree neighborhood scale, the strength and direction of effects on mortality probability depended on the size of the focal tree. A positive effect of NDI and a positive DBH:NDI interaction term ( Fig. 2; Appendix S5) suggest that neighborhood density contributes to MPB-induced lodgepole pine mortality, but effects vary by tree size. Greater NDI increased predicted mortality probability for medium to large diameter trees (15-30 cm DBH) but had a negligible effect on mortality for small-diameter trees (Fig. 3). Host proportion (i.e., pNDI) was a marginally important predictor of mortality probability, and the DBH:pNDI interaction term was an important positive predictor ( Fig. 2; Appendix S5). Greater pNDI decreased predicted mortality probability for small-diameter trees (10-15 cm DBH) and increased predicted mortality probability for large diameter trees (25-30 cm DBH; Fig. 3). Neither TPI nor the DBH:TPI interaction term were important predictors ( Fig. 2; Appendix S5). At the plot scale, plot D2 had a higher mortality probability than plots B2 or C3 overall ( Fig. 2; Appendix S5).
The posterior median practical range of the spatial GRF was 36.7 m, with a posterior median marginal standard deviation of 0.46 (Appendix S5). Spatial patterns of predicted mortality probability were comparable across the three plots, with the highest variability evident for medium diameter trees (15-25 cm DBH) and probabilities approaching zero and one for small and large diameter trees, respectively (Fig. 4). Accounting for effects of neighborhood scale covariates and the spatial GRF resulted in substantially different predictions when compared to the model that considered only plot and tree scale covariates, with differences in posterior median predicted mortality probability as large as AE0.35 (Appendix S6).
The model including neighborhood scale covariates and the spatial GRF was the best fitting model that we evaluated (sum of log CPO = −1081). When building on a model that includes only plot and tree scale covariates (sum of log CPO = −1137), adding both neighborhood scale covariates and the spatial GRF increased model fit by a greater amount (+56) than adding either the spatial GRF (+40) or neighborhood scale covariates (+19) alone. Autocorrelation in the residuals of models without the spatial GRF indicated that including the spatial GRF was necessary to account for the spatial structure of the data. A sensitivity analysis of neighborhood sizes showed that the effect of increasing or decreasing the neighborhood radius by 5 m was negligible (Appendix S7).

DISCUSSION
Explicitly considering how fine scale patterns and processes can mediate disturbance dynamics is important as climate warms and disturbance activity changes (Das et al. 2008, Bentz et al. 2010). Our study demonstrates how within-stand heterogeneity can mediate individual tree scale susceptibility to MPB outbreak, helping to x-axis) and predicted probability of mortality (y-axis) across gradients of basal area increment (BAI; column 1), neighborhood density index (NDI; column 2), neighborhood host proportion (pNDI; column 3), and topographic position index (TPI; column 4). (b) Change in relationship between covariates (x-axis) and predicted probability of mortality (y-axis) across gradients of DBH. Solid lines represent posterior medians, dark shading represents 95% credible intervals for one location in plot C3, and light shading represents variability in 95% credible intervals for locations across plot C3 (i.e., variability modeled by the spatial GRF). Predictions consider each combination of covariates in isolation, holding all other covariates at zero (i.e., their average values). inform understanding of disturbance dynamics across scales. As expected, and consistent with previous research (Safranyik and Carroll 2006), we found tree size to be the primary predictor of host tree mortality probability in stands severely attacked by MPB. Since MPB require a minimum phloem thickness to construct their galleries, small-diameter trees are often not capable of supporting MPB survival and reproduction, and large diameter trees are preferred hosts (Amman and Cole 1983). Furthermore, large trees are more susceptible to hydraulic failure, which can be amplified by bark beetle attack (McDowell et al. 2008, Bennett et al. 2015. Our model predictions suggest that for sufficiently small or sufficiently large diameter host trees, tree size largely overrides any additional mediating effects, with mortality probability approaching zero and one at 10 and 30 cm DBH, respectively. For host trees of intermediate size, however, we found that pre-outbreak growth rate and neighborhood structure (i.e., density and host proportion) may mediate tree scale susceptibility to severe MPB outbreak.
Although growth rate preceding an outbreak (measured here as BAI) was only a marginally important predictor of mortality probability in our model, our findings suggest it interacts with tree size in key ways to mediate mortality probability. We found that predicted mortality probability increased with BAI for small to medium diameter trees (~10-25 cm DBH), but the change Fig. 4. Spatial predictions of mortality probability for a range of lodgepole pine DBH (10-30 cm). Predictions are posterior median predictions based on the specified DBH (columns) and the values of NDI, pNDI, and the spatial GRF for a hypothetical tree located within each grid cell in each plot (rows). BAI and TPI were held at zero (i.e., their average values).
in mortality probability with BAI appeared negligible for larger trees. Growth rate preceding an outbreak can increase mortality probability in lodgepole pine (Cooper et al. 2018), as BAI may be a strong indicator of phloem thickness for lodgepole pine trees (Shrimpton and Thomson 1985). It is possible that small-diameter lodgepole pine, which may otherwise not physically be able to support MPB survival and reproduction, are more likely to be killed in an outbreak if they exhibit faster pre-outbreak growth rates and thus potentially have thicker phloem. Conversely, fast growth rates may not affect mortality probability for large diameter trees, as large trees already possess sufficient phloem thickness. While our study does not account for tree age, it does support the idea that slower growth may increase tree survival at certain life stages (de la Mata et al. 2017). With more precise growth measurements (i.e., tree cores), it may be possible to find stronger evidence for the effect of growth rate while disentangling the effects of tree size and age.
At the neighborhood scale, we found that the probability of individual lodgepole pine mortality increased with neighborhood density, supporting our expectation. At the stand scale, dense stands can have microclimate characteristics (e.g., lower temperature and wind speeds) that are more favorable for, and therefore preferred by, MPB (Bartos andAmman 1989, Negrón 2019). The increase in predicted mortality probability with neighborhood density in our model suggests neighborhood density may mediate tree scale mortality rates during MPB outbreak. Within a stand, it is possible that variation in neighborhood density may affect MPB infestation dynamics by influencing the dispersal of MPB-produced pheromones. While the interaction with tree diameter suggests that the effect of neighborhood density is most pronounced for medium to large lodgepole pine (~15-30 cm DBH), being positioned in a lower density area may confer added tree scale resistance to MPBinduced mortality across tree sizes.
After accounting for overall neighborhood density, we found the proportion of neighborhood density contributed by host trees (i.e., lodgepole pine) decreased predicted mortality for small-diameter lodgepole pine (~10-15 cm DBH) and increased predicted mortality probability for large diameter lodgepole pine (~25-30 cm DBH). We expected to find a positive effect of host proportion due to increased beetle pressure in areas dominated by lodgepole pine. The negative effect of host proportion for smalldiameter trees suggests mortality probability for small trees is reduced if they are rooted in an area dominated by host trees, contrary to what has been found in previous studies (Preisler and Mitchell 1993). One possible explanation may be that increased intraspecific competition in these areas may lead to reduced growth rates and ultimately thinner phloem, as captured partially by the BAI covariate. However, the magnitude of the effect is modest when compared to the main effect of neighborhood density.
Our expectation that within-stand topographic position would affect tree scale mortality probability was not supported in this study. While topographic effects may not be important at the within-stand spatial scales (e.g., sub-meter) we analyzed, it is possible that local topographic position may mediate tree scale mortality probability in more topographically complex stands or at broader spatial scales (e.g., several hectares). For example, at broader scales (among stands), outbreak refugia in topographically concave areas may be linked to cool, moist conditions, which potentially relieve tree stress during drought conditions and reduce MPB infestation (Vorster et al. 2017, Cartwright 2018. Complex feedbacks make it difficult to disentangle hydraulic stress, which may vary with both topographic position and tree physiology, from beetle-induced tree mortality, particularly for larger trees (McDowell et al. 2008, Gaylord et al. 2013, Bennett et al. 2015. Our data do provide support for the effect of topography among stands, as mortality was greatest in the plot with the steepest slope and the most southwesterly aspect (plot D2, Appendix S1), characteristics that relate to faster runoff, less infiltration, warmer conditions, and drier soils (Cartwright 2018).
Given the stochastic nature of beetle dynamics and the switching that occurs from focal to neighboring trees during outbreak (Preisler andMitchell 1993, Safranyik andCarroll 2006), there is inherent clustering in MPB outbreak dynamics and resultant spatial patterns of tree mortality. The spatial random effect in our model accounts for deviations from the pattern of mortality and v www.esajournals.org survival predicted using fixed effects at the tree, neighborhood, and plot scale. These deviations may represent variation in the probability of mortality due to environmental factors or tree scale physiological factors (Ferrenberg et al. 2014, Zhao andErbilgin 2019) not considered in this study, or variation in the probability of mortality due to the beetle-generated semiochemical landscape.
The interactions among tree scale and tree neighborhood scale factors in our study have important implications for tree scale resistance to MPB outbreak. Because the biology of MPB places constraints on which trees are susceptible to attack, and these constraints are linked tightly to tree size (Safranyik and Carroll 2006), the strength and direction of potential mediating effects cannot be considered without accounting for tree size. We found that pre-outbreak growth rates and neighborhood structure had the greatest effect on predicted mortality probability for trees of intermediate size (e.g.,(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)Figs. 3,4). Previous research has suggested that lodgepole pine trees ≤25 cm in diameter are beetle sinks (i.e., more beetles attack than emerge; Safranyik et al. 1974). In this intermediate range, trees are large enough to physically support MPB survival and reproduction, yet not large enough that they would be considered the most desirable hosts. It is in these cases that factors such as pre-outbreak growth rate, neighborhood density, and neighborhood host proportion may have mediating effects.
Our findings that tree scale mortality probability can vary with pre-outbreak neighborhood density, and to a lesser degree, pre-outbreak growth rate, have potential implications for management treatments that seek to confer resistance to bark beetle outbreak. First, our study suggests that neighborhood structure and growth rate matter most for trees of intermediate size (~15-25 cm DBH), whereas mortality probability for smaller (<10 cm DBH) and larger (>30 cm DBH) trees is driven more exclusively by their size. Second, the direction of effects for neighborhood density and growth rate suggest a trade-off when considering thinning stands to foster resilience to outbreaks. Decreased neighborhood density in thinned stands may increase survival probability of residual trees, whereas higher growth rates may decrease survival probability of residual trees, especially for trees of small to intermediate size. While thinning has the potential to reduce the probability of outbreak development at the stand scale (Fettig et al. 2007), it remains unclear whether and how thinned stands confer resistance to MPB at the individual tree scale in the event MPB outbreak does occur. For example, in Oregon (USA), MPB infestation was delayed in thinned lodgepole pine stands compared to unthinned stands, yet colonization patterns were comparable across stands once infestation did occur (Preisler and Mitchell 1993). We emphasize that our study was conducted in plots without prior management treatments, and outcomes could differ if neighborhood density was intentionally manipulated. Further study is warranted to test the trade-offs between neighborhood density and tree growth, how they may affect mortality probability for residual trees in stands that have been thinned, and how long any potential effects may persist after treatment.
In this study, we modeled the probability that an individual lodgepole pine will be killed by MPB over the course of an outbreak, given that an outbreak has occurred. This probability differs from both (a) the probability that an individual lodgepole pine will be killed by MPB at endemic (i.e., background) population levels, where different processes and dynamics would be expected at both the tree and neighborhood scale, and (b) the probability of an outbreak developing within a stand, which is the focus of stand scale hazard and risk ratings for MPB (e.g., Shore et al. 2000, Negrón andPopp 2004). Tracking tree scale spatio-temporal dynamics over time periods that include both endemic and epidemic stages of MPB activity (e.g., Preisler andMitchell 1993, Logan et al. 1998), as well as linking tree scale data to a range of other spatial scales (e.g., stand to regional; Seidl et al. 2016a) are important additional dimensions to consider in future research.

CONCLUSION
Conifer forests have been shaped by periodic bark beetle outbreaks for millennia (Raffa et al. 2008, Jarvis andKulakowski 2015). As climate warms and broad scale environmental constraints change (Bentz et al. 2010), fine scale constraints may become increasingly important for v www.esajournals.org understanding both the causes and effects of outbreaks. Our findings suggest that while tree size is the primary factor driving host tree mortality during MPB outbreak, pre-outbreak growth and neighborhood structure can mediate mortality probability in meaningful ways. Further, the direction and magnitude of these mediating effects may vary with tree size. The interactions among tree scale and tree neighborhood scale factors in this study highlight the importance of natural within-stand heterogeneity in shaping tree scale susceptibility to disturbance.