Bayesian Predictions of Bark Beetle Attack and Mortality of Three Conifer Species During Epidemic and Endemic Population Stages

Bark beetles naturally inhabit forests and can cause large-scale tree mortality when they reach epidemic population numbers. A recent epidemic (1990s–2010s), primarily driven by mountain pine beetles (Dendroctonus ponderosae), was a leading mortality agent in western United States forests. Predictive models of beetle populations and their impact on forests largely depend on host related parameters, such as stand age, basal area, and density. We hypothesized that bark beetle attack patterns are also dependent on inferred beetle population densities: large epidemic populations of beetles will preferentially attack large-diameter trees, and successfully kill them with overwhelming numbers. Conversely, small endemic beetle populations will opportunistically attack stressed and small trees. We tested this hypothesis using 12 years of repeated field observations of three dominant forest species (lodgepole pine Pinus contorta, Engelmann spruce Picea engelmannii, and subalpine fir Abies lasiocarpa) in subalpine forests of southeastern Wyoming paired with a Bayesian modeling approach. The models provide probabilistic predictions of beetle attack patterns that are free of assumptions required by frequentist models that are often violated in these data sets. Furthermore, we assessed seedling/sapling regeneration in response to overstory mortality and hypothesized that higher seedling/sapling establishment occurs in areas with highest overstory mortality because resources are freed from competing trees. Our results indicate that large-diameter trees were more likely to be attacked and killed by bark beetles than small-diameter trees during epidemic years for all species, but there was no shift toward preferentially attacking small-diameter trees in post-epidemic years. However, probabilities of bark beetle attack and mortality increased for small diameter lodgepole pine and Engelmann spruce trees in post-epidemic years compared to epidemic years. We also show an increase in overall understory growth (graminoids, forbs, and shrubs) and seedling/sapling establishment in response to beetle-caused overstory mortality, especially in lodgepole pine dominated stands. Our observations provide evidence of the trajectories of attack and mortality as well as early forest regrowth of three common tree species during the transition from epidemic to post-epidemic stages of bark beetle populations in the field.

Bark beetles naturally inhabit forests and can cause large-scale tree mortality when they reach epidemic population numbers. A recent epidemic (1990s-2010s), primarily driven by mountain pine beetles (Dendroctonus ponderosae), was a leading mortality agent in western United States forests. Predictive models of beetle populations and their impact on forests largely depend on host related parameters, such as stand age, basal area, and density. We hypothesized that bark beetle attack patterns are also dependent on inferred beetle population densities: large epidemic populations of beetles will preferentially attack large-diameter trees, and successfully kill them with overwhelming numbers. Conversely, small endemic beetle populations will opportunistically attack stressed and small trees. We tested this hypothesis using 12 years of repeated field observations of three dominant forest species (lodgepole pine Pinus contorta, Engelmann spruce Picea engelmannii, and subalpine fir Abies lasiocarpa) in subalpine forests of southeastern Wyoming paired with a Bayesian modeling approach. The models provide probabilistic predictions of beetle attack patterns that are free of assumptions required by frequentist models that are often violated in these data sets. Furthermore, we assessed seedling/sapling regeneration in response to overstory mortality and hypothesized that higher seedling/sapling establishment occurs in areas with highest overstory mortality because resources are freed from competing trees. Our results indicate that large-diameter trees were more likely to be attacked and killed by bark beetles than small-diameter trees during epidemic years for all species, but there was no shift toward preferentially attacking small-diameter trees in postepidemic years. However, probabilities of bark beetle attack and mortality increased for small diameter lodgepole pine and Engelmann spruce trees in post-epidemic years compared to epidemic years. We also show an increase in overall understory growth (graminoids, forbs, and shrubs) and seedling/sapling establishment in response to beetle-caused overstory mortality, especially in lodgepole pine dominated stands. Our observations provide evidence of the trajectories of attack and mortality as well as early forest regrowth of three common tree species during the transition from epidemic to post-epidemic stages of bark beetle populations in the field.
Keywords: bark beetles, forest mortality, Engelmann spruce -subalpine fir forests, Pinus contorta (lodgepole pine), regeneration INTRODUCTION Commencing in the latetwentieth century, the western North America bark beetle epidemic and associated forest mortality rank among the highest recorded in terms of both severity and spatial scale (Hyde et al., 2016). Comprising multiple beetle species, including mountain pine beetle (Dendroctonus ponderosae), Douglas-fir beetle (Dendroctonus pseudotsugae), western balsam bark beetle (Dryocoetes confuses), pinþon ips (Ips confusus), spruce beetles (Dendroctonus rufipennis), the outbreak has affected up to eleven million hectares across North America (Meddens et al., 2012) and accounts for ∼32% of forest mortality in the Western United States (Berner et al., 2017). The magnitude of this outbreak is attributed to forest management (e.g., fire suppression) and resulting homogenization of forest cover (Raffa et al., 2008), paired with synergistic effects of a warmer, drier climate that has the potential to benefit insects and associated fungal symbionts (Stahl et al., 2006;Bentz et al., 2010;Weed et al., 2015;Kolb et al., 2016) while exerting stress on trees (Allen et al., 2010;Lloret and Kitzberger, 2018). In addition to potential interactive effects with subsequent fire (Simard et al., 2011;Harvey et al., 2013;Hart et al., 2015;Keane et al., 2015;Sieg et al., 2017), the magnitude of this outbreak has prompted evaluation of the effects of bark beetle epidemics on ecosystem biogeochemical and water cycles and species assemblages.
Bark beetles are a native disturbance agent in North American forests, and stable populations with low beetle densities (i.e., endemic) persist in stands with suitable host trees. Due to their small number, endemic bark beetle populations rarely overcome the defenses of healthy, large-diameter trees, and attacks are generally limited to small-diameter (i.e., smaller than stand average) and weakened (e.g., stress from injury, drought, or other pathogens) host trees that are dispersed across the landscape (Waring and Pitman, 1985;Christiansen et al., 1987;Carroll et al., 2006;Safranyik and Carroll, 2006). The resulting canopy and root gaps alter resource availability to the understory (Parsons et al., 1994) and create a spatial patchwork of freshly attacked stands, un-attacked ones, and regenerating stands from previous attacks. Thus, bark beetles at endemic populations play an important role in maintaining a diversity of structure and tree age classes. Periodically, the system will shift into an alternate state (i.e., epidemic) after overcoming an outbreak threshold that is largely determined by local host tree age class (Raffa and Berryman, 1986). Population growth during epidemics is largely densitydependent (Trzcinski and Reid, 2009), although extrinsic factors such as local climate (Raffa et al., 2008;Weed et al., 2015), host tree vigor (Lewis et al., 2010), and interspecific interactions (Safranyik and Carroll, 2006) are important co-determinants of beetle population eruptions and consequent host mortality. The severity of epidemic bark beetle outbreaks is further amplified by a positive feedback loop. During epidemics, a high number of bark beetles is able to successfully overcome the defenses of large-diameter, healthy trees. These trees generally have a much greater volume of phloem, leading to increased beetle populations in the following year that in turn are able to overcome more large-diameter, healthy host trees (Boone et al., 2011). Landscapes comprised of even-aged, dense stands further facilitate swarming as the pheromones that bark beetles use to signal attacks decline exponentially with distance (Sullivan and Mori, 2009). Epidemic populations generally collapse with the advent of adverse weather conditions or with the depletion of suitable host trees (Safranyik and Carroll, 2006). Besides temporary (and often muted) alterations to water, nitrogen, and carbon cycles (Mikkelson et al., 2013;Biederman et al., 2014;Emmel et al., 2014;Frank et al., 2014;Reed et al., 2014Reed et al., , 2018Norton et al., 2015), the resulting wide-spread mortality of primarily old-growth trees has the potential to modify postdisturbance seedling/sapling establishment patterns (Pettit et al., 2019) and slow down succession (Bretfeld et al., 2015).
Empirical models based solely on indices of stand susceptibility (e.g., basal area and average tree diameter) perform reasonably well in predicting bark beetle outbreaks (Negron and Klutsch, 2017). However, the inclusion of additional ecological interactions has the potential to further strengthen these models (Nelson et al., 2008). While beetle density has been included as a risk factor in predictive models of forest damage (Lewis et al., 2010), host tree selection behavior is rarely included and in fact was suggested to "not be necessary to model explicitly" by Nelson et al. (2008). This is a reasonable assumption based on some studies showing a tight correlation between beetle behavior and population density (Wallin and Raffa, 2004;Koontz et al., 2021). However, Kausrud et al. (2011) shows that populations can shift toward epidemic behavior (i.e., targeting healthy trees) by other factors than population density, such as lowered colonization thresholds caused by drought stress or increased beetle aggressiveness, or even spontaneously. In addition, fieldcollected data on host selection behavior throughout a beetle outbreak are rare and, together with the implications of altered stand demographics on subsequent composition, remain an area of active research.
Recent research has improved our mechanistic understanding of beetle-induced tree mortality (Hubbard et al., 2013) but also highlights how the relative contribution of bark beetles to observed forest mortality varies by host tree species (Smith et al., 2015;Tai et al., 2019). In this study, we quantified signs of bark beetle attacks and resulting tree mortality in situ through annual monitoring of forest plots throughout a major bark beetle epidemic in southeastern Wyoming. To test for interspecific differences in host-insect dynamics, data were collected for three dominant forest tree species: shade-intolerant lodgepole pine (Pinus contorta) that are hosts to mountain pine beetles (D. ponderosae), shade-intolerant, co-dominant subalpine fir (Abies lasiocarpa) and Engelmann spruce (Picea Engelmannii) that are hosts to western balsam bark beetles (D. confuses) and spruce beetles (D. rufipennis), respectively. We hypothesized that during peak beetle activity (i.e., epidemic population), largediameter trees will have higher probabilities to exhibit signs of beetle attack and resulting mortality, whereas in post-epidemic years, probabilities of signs of attack and mortality will increase in smaller diameter trees as host selection shifts toward endemic behavior. We also hypothesized that opening of the canopy and root gaps from mortality of primarily large-diameter trees during the epidemic will lead to increased seedling/sapling establishment of shade-intolerant lodgepole pines whereas no such increase will be evident for shade-tolerant Engelmann spruce and subalpine fir. Our hypotheses were tested with annual vegetation data spanning 12 years coupled with Bayesian phenomenological models that produce probabilities of trees being attacked and/or successfully killed by bark beetles for each year.

Field Data Collection
All data were collected in the sub-alpine forests of the Medicine Bow-Routt National Forest in southeastern Wyoming between 2005 and 2019 ( Table 1), spanning peak epidemic bark beetle activity between 2007 and 2011 (Biederman et al., 2014;Frank et al., 2014). Accordingly, for this study we consider data collected between 2008 and 2011 as "epidemic" and between 2009 and 2019 as "post-epidemic." Primary data collection for this study occurred at the "Chimney Park" AmeriFlux site (∼2,750 m.a.s.l; US-CPk) 1 that comprises predominantly lodgepole pine (P. contorta) forest. Mature lodgepole pine forests are characterized by an even-aged canopy, sparse understory vegetation, and dry soils in late summer. Over 90% of beetleinduced mortality in these forests is caused by mountain pine beetles (D. ponderosae), with less than 10% attributed to pine engravers (Ips pini), and Pityogenes knechteli/Pityopthorus confertus (Audley et al., 2021). Each year from 2008 to 2019, at least 20 plots were sampled with the exception of 2018 when the site was partially burned in the "Badger Creek Fire" (Supplementary Table 1). Additional plots were established in nearby unburned areas to replace burned plots. In order to minimize bias due to our re-sampling representing only one site, we included data from additional sites in southeastern Wyoming (Supplementary Figure 1). These sites were part of other research projects (e.g., Tai et al., 2019) and were not sampled every year between 2008 and 2019. Plots at the "North Brush Creek, " "North Lodgepole Creek, " "Glacial Tills, " and "NoName" sites (2750-3150 m.a.s.l.) were sampled in 2011, 2013, and from 2014 to 2017, respectively. These plots comprised a mixture of lodgepole pine, Engelmann spruce (P. engelmannii), and subalpine fir (A. lasiocarpa) in the overstory with thick understory vegetation and higher soil moisture than lodgepole pine forests. Lastly, the highest elevation plots included in this 1 https://doi.org/10.17190/AMF/1246150 study are from the "GLEES" AmeriFlux site (3190 m.a.s.l.; US-GBT) 2 and comprised exclusively Engelmann spruce and subalpine fir. Plots at this site were sampled in 2005 and annually between 2009 and 2017. Spruce-fir forests are generally characterized by an uneven-age distribution with higher basal area and soil moisture compared to lodgepole pine forests. Beetleinduced mortality of Engelmann spruce is primarily caused by spruce beetles (D. rufipennis) in the Medicine Bow-Routt National Forest (Negrón and Popp, 2018), whereas mortality in subalpine fir is likely a combination of Ips spp., western balsam bark beetle (Dryocoetes confuses), drought, and fungal pathogens (Harvey et al., 2021).
Sampling protocols varied slightly between sites. In general, there were three spatial scales at which data were collected and analyzed for each site: stands, plots, and subplots. Stands comprised of clusters of up to nine proximate (tens of meters maximum distance) plots in which data were collected. Plot areas varied between different sites, ranging from 9 to 400 m 2 and were either circular or rectangular. Plot boundaries or centers were marked with semi-permanent pins and GPS referenced to facilitate annual locating of plots. All tree stems located and rooted within a plot with a diameter at breast height (DBH; measured at 1.37 m) larger than 2.5 cm were surveyed for the following characteristics: species, DBH, and condition (healthy: green, no beetle activity; beetle-green: green foliage with evidence of beetle activity; beetle-red: red foliage with evidence of beetle activity; beetle-dead: no foliage with evidence of beetle activity; and dead: no foliage with no evidence of beetle activity). In this study, we did not identify or quantify bark beetles directly but instead based evidence of bark beetle activity on the presence/absence of pitch tubes, beetle entrance holes, and boring dust (Safranyik, 1988;Safranyik and Carroll, 2006).
In all plots, additional subplots were sampled inside the plot area to assess the regenerating overstory (i.e., trees not reaching 1.37 m in height; from here on seedlings/saplings) and understory growth. Subplots to quantify seedlings/saplings were generally 1m wide transects through the plot center and spanning the largest distance across each plot. Depending on plot size, these transects ranged from 3 to 28 m 2 in area. Subplots to visually assess forest floor vegetation cover comprised of 20 × 50 cm (before 2015) and 100 × 100 cm (after 2015) Daubenmire frames (Daubenmire, 1959) at 1-5 locations inside each plot, depending on plot size.
Total understory vascular plant cover per plot was calculated as the sum of cover estimates for graminoids, forbs, and shrubs less than 1 m tall. For analysis of understory vegetation and seedling/sapling response, data were split into lodgepole pine or spruce-fir dominated stands. The cutoff for the classification of a stand was based on basal area contribution of more than 75% of either lodgepole pine or the combination of Engelmann spruce and subalpine fir trees. This separation was done to account for expected differences in seedling/sapling/sapling responses to the opening of the canopy. Lodgepole pine is a shade-intolerant species whereas both Engelmann spruce and subalpine fir are shade tolerant species that often grow in the shade of lodgepole pine (Peet, 1981). Because of different sampling protocols and resulting differences in sampling areas, all data were scaled to a common denominator (i.e., seedlings/saplings m −2 ) or converted to proportions for the analysis.

Statistical Models
To address our first hypothesis, a conditional logistic model was developed based on the phenomenological principle that the probability of trees being attacked and being killed by bark beetles each year is a logistic function of DBH. This model consisted of two parts. Equation 1 estimates the probability a tree exhibits signs of bark beetle attack, based on the tree's diameter: where φ 1 is the probability that a tree exhibits signs of beetle attack, D is the tree's diameter at breast height and B 0,1 and B 1,1 are regression coefficients. Eq. 2 estimates the probability an attacked tree is killed by the bark beetles, based on the tree's diameter: where φ 2 is the probability of beetle-induced tree mortality, D is the tree's diameter at breast height and B 0,2 and B 1,2 are regression coefficients. These equations come together to estimate the probability that any tree at a given DBH is healthy (i.e., green foliage, no signs of beetle activity in the bark; Eq. 3), beetle-green (i.e., green foliage, clear signs of beetle activity in the bark; Eq. 4), or beetle-dead (i.e., foliage red or missing, clear signs of beetle activity in the bark; Eq. 5). For lodgepole pines, the foliage of trees that have recently succumbed to beetles turn red in color (i.e., beetle-red phase) and thus are easily identified (Edburg et al., 2012). By separating beetle-dead trees between those with red foliage (recently succumbed) and those without foliage (died >1 year ago), we were able to assess year-to-year probabilities of beetle-induced mortality by DBH class for lodgepole pines. Because onset of the beetle-red phase in both Engelmann spruce and subalpine fir can be delayed (Holsten et al., 1989), we did not apply this approach to these two species.
where P(H) is the probability that any given tree is healthy, P(BG) is the probability that the tree is beetle-green, and P(BD) is the probability that the tree is beetle-dead (Figure 1). These probabilities all sum to one. All regression coefficients were fit using a Bayesian Markov Chain Monte-Carlo (MCMC) with naive priors and likelihoods; an approach that provides probabilistic predictions of beetle attack patterns that are free of typical frequentist assumptions that are commonly violated in these data sets (Olejnik and Algina, 1983;Press, 2005). Computations were performed in R Version 3.2.3 (R Core Team, 2019) using code modeled after Kruschke (2014) and the statistical package RStanARM package (Goodrich et al., 2020). The MCMC algorithm used five chains and 20,000 iterations per chain with a burn-in of 10,000 iterations. The above statistical model was fit independently for each tree species and each year across all sites, with a total of 25,412 measured trees used in the analysis (Supplementary Figure 2). Model fit (Bayes R2) was assessed following Gelman et al. (2019). In 2005, no beetle condition classes were recorded and as a result the data were not included in the statistical model. Strong selective behavior toward large trees was inferred from posterior median estimates of the slope regression coefficient when (1) the coefficient was positive and (2) the 94% credible interval did not include zero (McElreath, 2020). We then applied our model to predict the probabilities of beetle attack and post-attack mortality for small (5th percentile of sizes of sampled trees in this study) and large-diameter (95th percentile of sizes of sampled trees in this study) trees. We also predict probabilities of attack and postattack mortality of average sized, live trees that were sampled in any given sampling year to account for mortality-driven changes in average host tree sizes that potentially limit successive beetle broods to smaller sized trees.
To address our second hypothesis, a linear model within a Bayesian framework was used to assess seedling/sapling/sapling regeneration in response to overstory mortality. Similar to our model above, a Bayesian approach allows not only for the estimation of posterior means of the model parameters (i.e., slope and intercept) but also for the prediction of new observations with known uncertainty. Similar to the logistic model above, we utilized the RStanARM package (Goodrich et al., 2020) and inferred an effect of overstory mortality on seedling/sapling establishment when the 94% credible interval of the posterior slope did not include zero. Model parameter estimates and model predictions are reported as medians (lower 94% credible interval, upper 94% credible interval). Annual DBH and understory regeneration are reported as means (± standard error of the mean) to indicate uncertainty around the estimate of the population mean.

Basal Area
The inclusion of data from multiple sites that do not span the entire sampling period of this study resulted in very different year-to-year values of local basal area. Thus, we do not report absolute basal area values but instead report proportions of all Frontiers in Forests and Global Change | www.frontiersin.org FIGURE 1 | Conceptual framework of the MCMC model used to estimate probabilities of any given tree being healthy (1 − φ 1 ), beetle-green [φ 1 × (1 − φ 2 )], and beetle-dead (φ 1 × φ 2 ). Beetle-dead trees included trees in both beetle-red and beetle-gray phases. Due to the predictable beetle-red phase in lodgepole pines 1 year after succumbing to beetles, we used a separate model approach that only included beetle-red trees to assess year-to-year probabilities of mortality for that species.

Probabilities of Bark Beetle Attack and Mortality Across Tree Diameters
Posterior distributions of the conditional mean adequately represent mean proportions of observed data (Figure 3 and Supplementary Figures 3-16). Model fit ranged from less than 1% (e.g., mortality of subalpine fir 2016 and 2017) to 72% for beetle attack of lodgepole pine in 2010 (Supplementary  Tables 2, 3). Positive posterior median slope parameter estimates with 94% credible intervals not including zero indicate that increasing diameters were associated with increasing probabilities of exhibiting signs of bark beetle attack and bark beetle caused mortality in lodgepole pines throughout the entire sampling period (Figure 4 and Supplementary Table 2). Including only trees that recently succumbed to bark beetles (i.e., beetle-red phase) in the model for lodgepole pine trees lends further support that the probability of mortality increased with increasing tree diameter during peak epidemic years (2009)(2010)(2011) and also in 2012. However, in the post-epidemic years after 2012 the pattern was less clear with increasing probabilities of mortality in large-diameter trees in 2015 but decreasing probabilities in 2016, and no positive or negative trend in the other years. Probabilities of mortality stayed below 25% across all tree sizes in 2016 and 2017. This model could not be run for the years 2018 and 2019 as no trees were recorded in the beetle-red phase. The lack of beetle-red trees in those years can be interpreted as zero probability of new beetle mortality.
Similar to lodgepole pine, large-diameter subalpine fir trees were more likely to be attacked by bark beetles throughout the sampling period (Figure 5 and Supplementary Table 3). However, bark beetle caused mortality showed no relationship to tree diameters except in 2009 when larger trees were more likely to succumb to beetles (Figure 5). Probabilities of bark beetle attack were also higher in large-diameter Engelmann spruce trees in all sampling years (Figures 6A,C Supplementary Table 3). During epidemic years, up to 2011, large-diameter trees were more likely to be killed by bark

Stand Demographics and Model Predictions
Average DBH of live trees (i.e., healthy and beetle green) of all three species decreased notably during the epidemic years (Supplementary Figure 17).  (2017). This represents a decrease of 27% (lodgepole pine), 51% (subalpine fir), and 39% (Engelmann spruce) compared to the first year of sampling.
The 0.05 and 0.95 quantile of all sampled tree diameters measured in this study were 4.1 and 25.8 cm for lodgepole pine, 5.3 and 30.3 cm for subalpine fir, and 5.7 and 37.0 cm for Engelmann spruce. By parameterizing our model with the average DBH of live trees as well as hypothetical trees of 5 cm (from here on referred to as "small"), 25 cm ("medium"), and 40 cm ("large") DBH, we predict probabilities of beetle attack and post-attack mortality spanning sizes that represent over 90% of all sampled trees (Figure 7). For all three species, the predicted probabilities of being attacked by beetles are highest in the largest size class and lowest in the smallest. The predicted probability of being attacked by beetles remains under 3% for small lodgepole pine trees during the epidemic (2008-2010) but increases in the last year of the epidemic to 6.9% (94% credible interval:  clear temporal pattern is evident for Engelmann small spruce trees. An average-sized, live Engelmann spruce had a much higher predicted probability of being attacked in the epidemic years than in post-epidemic years. Similar to probabilities of attack, predicted probabilities of post-attack mortality were highest for largest and lowest for smallest lodgepole pine trees (Figure 7). Compared to epidemic years, predicted probabilities of mortality were higher for small diameter lodgepole pine trees in post-epidemic years, increasing from 1% (27.8, 0) early in 2008 to up to 60.7 (76.1, 43.1) in 2013, and ending at 30.8 (64.4, 9.3) in 2019. Predicted probabilities of mortality of medium sized trees increased during the epidemic, from 31.5 (48.6, 16.8) in 2008 to 86.5 (90.7, 81.4) in 2011, and remained high in post-epidemic years. Year-to-year mortality predictions of mortality show a similar pattern, with higher probabilities in larger trees during the epidemic, increasing probabilities in smallest trees in early post-epidemic years and a spike in 2013, when all size classes exhibit similar probabilities to succumb to beetles. Predicted probabilities and credible intervals include zero beginning in 2014, suggesting no additional trees are dying from bark beetles following that year.
Overlapping 94% credible intervals indicate no discernable size-dependence in predicted probabilities of mortality in subalpine fir. Although larger size was accompanied by greater predicted probabilities of beetle-induced mortality in epidemic years (2009)(2010)(2011), small Engelmann spruce trees were as likely as medium and large-diameter trees to succumb to bark beetles in post-epidemic years beginning in 2012 until 2017, when small trees are once again least likely and large trees most likely.
Lodgepole pine seedling/sapling density was higher in medium and high mortality than in low mortality plots in most of the sampling years ( Figure 8C). At the end of the epidemic in 2011, there were 0.4 (standard error of the mean: ± 0.2), 1.5 (±0.5), and 1.4 (±0.5) seedlings m −2 in low, medium, and high mortality plots. In 2019, seedling counts were similar to those at the end of the epidemic, with 0.3 (±0.2), 1.4 (±0.7), and 1.7 (±0.5) seedlings m −2 for low, medium, and high mortality plots. In stands dominated by Engelmann spruce and subalpine fir, seedling/sapling density also increased slightly, from 0.09 (±0.06), 0.07 (±0.07), and zero seedlings m −2 in 2009 to 0.8 (±0.3), 0.9 (±0.1), and 1 (±0.5) seedlings m −2 in 2017 for low, medium, and high mortality plots, respectively, but showed no clear response to severity of overstory mortality (Figure 8D). Seedlings/saplings in all stands comprised primarily of conspecifics of the dominant overstory species. Posterior slope parameter estimates (mean: 1.3; 94% credible interval: 0.8, 1.9; Bayes R2: 0.19) indicate a 94% probability that higher overstory mortality results in higher lodgepole pine seedling/sapling establishment in lodgepole pine dominated stands (Figure 8E and Supplementary Table 4), whereas no effect of overstory mortality on Engelmann spruce and subalpine fir seedling/sapling establishment was detected in either forest type ( Figure 8F). In spruce-fir dominated forests, subalpine fir seedlings/saplings exhibited the highest densities across the gradient of overstory mortality.

Probabilities of Bark Beetle Attack and Mortality
Our data and resulting models support the first prediction of our primary hypothesis that large trees are more likely to be attacked and killed by bark beetles during an epidemic; however, the second part of the hypothesis that states that small trees are more likely to be attacked and killed during non-epidemic years is not clearly supported. Our results suggest that largediameter trees were more likely to be attacked and killed by bark beetles than small-diameter trees during epidemic years for all species. Similar preferences for large-diameter trees have been reported for mountain pine beetles (D. ponderosae) by Negrón (2019) and spruce beetles (D. rufipennis) by Bakaj et al. (2016) in Colorado. However, there was no sustained shift toward preferentially attacking small-diameter trees in post-epidemic years (Figures 4-7). Instead, predicted probabilities of bark beetle attack (lodgepole pine) and mortality (lodgepole pine and Engelmann spruce) of small diameter trees increased and in some cases were higher or similar to those of large diameter trees only in the immediate post-epidemic years.
Host tree selection is not only a function of forest physiognomy, specifically host size and stand density, but also bark beetle behavior. For example, spruce beetles (D. rufipennis) originating from epidemic populations attacked primarily largediameter, healthy trees, and this behavior was more pronounced with a higher population density, whereas beetles from endemic populations were unlikely to attack healthy trees, even with a high beetle density (Wallin and Raffa, 2004). This shift in behavior is also implied by our data for mountain pine beetles (primary host lodgepole pine), as indicated by the more pronounced preference for large trees during the epidemic and the increasing probabilities of attack and mortality of small-diameter trees in post-epidemic years (Figure 7).
Decreases in average host tree sizes due to epidemic beetles that preferentially select large trees inevitably limit postepidemic beetle broods to smaller-size host trees. By predicting probabilities of beetle-attack and beetle-caused mortality for live trees of average DBH in a given year (Figure 7), we aim to address this potential confounding factor. Predicted probabilities of beetle impacts on average DBH trees provide further evidence of shifts in beetle behavior between epidemic and post-epidemic years for lodgepole pine (attack and mortality) and Engelmann spruce (mortality). However, it should be pointed out that we did not tag individual trees and thus newly attacked and killed individuals could not be assessed from year to year except for lodgepole pine, where our year-to-year predictions agree with results from the model approach that includes all beetle-killed trees (Figure 7). Future studies should aim toward directly measuring bark beetle densities while tagging individual trees to allow for better year-to-year tracking.
Interactions between host tree size and local environmental conditions also affect host tree selection. For example, for epidemic conditions of western pine beetles (Dendroctonus brevicomis) in the Sierra Nevada mountain range, large-diameter trees had a higher probability of mortality in hot/dry site whereas small-diameter trees had a higher probability of mortality in cool/moist sites (Koontz et al., 2021). Climate, specifically warming and drought, are critical factors affecting subalpine fir growth and mortality (Reich et al., 2016;Kelsey et al., 2018). Although western balsam bark beetles (D. confuses) accounted for 12% of subalpine fir mortality between 2008 and 2013 in the Colorado Front Range, over 72% of observed mortality was attributed to abiotic factors (Smith et al., 2015). Recent work suggests that complex, scale-dependent interactions between biotic and abiotic factors are involved in subalpine fir mortality in Colorado (Harvey et al., 2021), and that host tree selection of bark beetles was not dependent on tree diameter (Lalande et al., 2020). Of the three species sampled in our study, subalpine fir showed the weakest relationship between host tree size and probabilities of beetle attack and beetle-caused mortality (Figure 5). However, Bleiker et al. (2003) showed a clear beetle preference for large subalpine fir trees in interior British Columbia but also report other factors (resin production, recent growth) that were identified as important predictors of successful versus unsuccessful beetle attacks trees. Together with these studies, our data suggest that the relative importance of tree size on susceptibility of subalpine fir to bark beetles interacts strongly with local climate and host tree vigor. Future research needs to account for secondary factors such as pathogens and climate, and their potential interactive effects on the behavior of bark beetles and the host defenses.

Stand Demographics and Understory
Our data lend support to the hypothesis that lodgepole pine seedling/sapling establishment increases due to overstory mortality caused by bark beetles, but not subalpine fir and Engelmann spruce (Figures 8E,f). Interestingly, outbreak severity as estimated from percent of dead basal area was positively correlated with lodgepole pine seedling/sapling density, contrary to results by Kayes and Tinker (2012) for data collected in 2010 in southeastern Wyoming. The relatively broad credible interval of our predictions suggests that additional factors need to be accounted for to better predict seedling/sapling establishment. For example, soil nitrogen dynamics in experimentally logged FIGURE 8 | Time series of total understory vascular plant cover (A,B), density of total seedlings/saplings (C,D) by severity of bark beetle caused overstory mortality and density of seedlings/saplings as a function of the percentage of basal area killed by bark beetles (E,F) in lodgepole pine dominated stands (A,C,E) and in Engelmann spruce and subalpine fir dominated stands [i.e., spruce-fir stands; (B,D,F)]. Stand dominance was based on local basal area contribution larger than 75%. Data from mixed stands are shown in Supplementary Figure 18. Note that although understory vegetation cover data was collected in 2010 in lodgepole pine stands (A), no plots were classified as "high mortality" in 2010, and seedling/sapling data was not collected until 2011. Understory cover was not collected until 2011 in spruce-fir stands. Lines in panels (A-D) represent mean values, bands indicate standard errors. Lines in panels (E,F) represent the means of the posterior distributions, gray bands indicate 94% credible intervals of the predicted seedling/sapling density.
lodgepole pine forests were dependent on overstory and corresponding root gap size (Parsons et al., 1994), indicating that mortality of multiple adjacent trees rather than the overall fraction of dead trees could be a better predictor of understory regeneration. In addition, competing understory vegetation and microsite conditions can affect successful seedling/sapling establishment (Stuart et al., 1989;Kayes and Tinker, 2012;Petrie et al., 2016). Although we did not assess differences in local microclimates, our data show an increase in understory vegetation cover over time in response to medium and high overstory mortality. Increases in understory diversity and productivity were also evident in post-beetle lodgepole pine forests in western Canada (Pec et al., 2015). Ranging between 20 and 50% in 2017, our values for understory vegetation cover overlap with those reported by Collins et al. (2011) 6 years after peak tree mortality in lodgepole pine dominated forest. In their study, no new recruits were found in plots with more than 45% herbaceous plant cover, lending support to competition as an important factor to consider.
Previous studies suggest that future forest composition will be altered in response to extensive overstory mortality in lodgepole pine (Collins et al., 2010(Collins et al., , 2011Kayes and Tinker, 2012;Campbell and Antos, 2015) and other forests (Cobb et al., 2017). We did not model future stand composition but our seedling/sapling data show that subalpine fir was the second most abundant species in the understory in stands with more than 75% lodgepole pine in the overstory, and that it was by far the most common seedling/sapling in sprucefir forests.

Implications for Management
Increasing our understanding of patterns in beetle attacks and associated tree mortality enables us to improve predictions of future outbreaks and apply appropriate forest management strategies. In addition to stand density and average tree age, management operations need to take into account the age and size distributions and their interactions with beetle populations. Our field-based observations suggest that the preference to attack large trees persists independent of epidemic or endemic beetle populations. This explains why an increasing number of large trees leads to rapid population explosions whereas an increasing number of small trees does not. Paired with our observations of increased probabilities of mortality in small trees during endemic populations, our results suggest that beetles may even out host tree demographics due to increased probabilities of beetle-induced mortality of large trees during epidemic and small trees during endemic conditions. Because large trees are not likely to be killed during endemic populations while small trees exhibit higher mortality rates, long lasting endemic bark beetle populations potentially accelerate the stand aging process.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: https://doi.org/10.15786/ 14524053.v1.

AUTHOR CONTRIBUTIONS
HS, DB, and MB collected the field data with sampling design input from BE. HS performed the initial data analysis and wrote the first manuscript draft. MB performed the additional analyses with input from BE and wrote the final manuscript, with input from all authors. All authors contributed to the article and approved the submitted version.

FUNDING
Funding for this research was provided by the National Science Foundation (awards nos: 1904421 and 1208909).