Epidemic Spruce Beetle Outbreak Changes Drivers of Engelmann Spruce Regeneration

Climate-mediated disturbances outside the range of historical variability can have severe consequences on vital, post-disturbance regeneration processes. High-elevation forests of the Rocky Mountains that are dominated by Engelmann spruce (Picea engelmannii) and subalpine fir (Abies lasiocarpa) are expected to be sensitive to climate change. Additionally, these forests have experienced recent epidemic spruce beetle (Dendroctonus rufipennis) outbreaks that have often resulted in >95% mortality of overstory Engelmann spruce. Therefore, the future distribution of Engelmann spruce forests depends largely on natural regeneration processes. We examined Engelmann spruce seedlings across gradients in soil moisture and stand structural conditions 20 yr post-disturbance on the Markagunt Plateau in southern Utah. All Engelmann spruce seedlings were mapped, measured, and aged, and aspects of stand structure and the microclimate were measured. The goal of our research was to infer processes affecting Engelmann spruce establishment by determining if patterns of advance regeneration that established before the outbreak (~60% of individuals) differed from seedlings that established during and immediately following the outbreak (combined into one group, ~40% of individuals). A generalized linear multi-model approach identified that the density of advance regeneration (seedlings/saplings) was negatively influenced by historical competition with overstory trees. In contrast, post-outbreak regeneration was related to microclimate conditions, including positive relationships with climatic moisture deficit and July soil water content. All seedlings were not significantly clustered around Engelmann spruce snags; however, there was evidence of facilitation of post-outbreak seedlings by pre-outbreak seedlings at higher elevation sites with lower moisture deficit. Together, these findings suggest post-outbreak seedlings were not moisture-limited at lower elevations but instead encouraged by higher evapotranspiration. Moreover, facilitation at higher elevations likely resulted from how pre-outbreak seedlings modify snowpack and associated seedbed environments. Our study provides insight for managing Engelmann spruce after a beetle outbreak. In these forests, preand post-outbreak regeneration can increase resilience to climate–disturbance interactions, but are patchy and structured at different scales. Therefore, the presence of advance regeneration and the likelihood of post-outbreak seedlings depend on local environment (soil moisture and stand structure) and could be taken into account to most effectively plan post-disturbance planting activities.


INTRODUCTION
Climate change is expected to directly alter the geographic distributions of tree species by affecting physiological processes that govern establishment, growth, and mortality, and indirectly by increasing the frequency and severity of disturbances (Dale et al. 2001, Rehfeldt et al. 2006, Lenoir et al. 2008, Turner 2010, Walter et al. 2017. Climate-mediated disturbances that are outside the range of historical variability may have profound impacts on species distributions through their effects on regeneration processes (Johnstone et al. 2016, Hansen et al. 2018. Regenerating tree seedlings and saplings may be particularly sensitive to changes in climatic conditions (Davis et al. 2019). Thus, disturbancerelated increases in radiation and temperature at the ground level may compound effects of changes in the macroclimate, triggering shifts in species distributions (De Frenne et al. 2013). Understanding the effects of climate on postdisturbance tree regeneration is important given the implications for future species distributions, biodiversity, and ecosystem structure and function.
Forest regeneration relies on complex interactions among biotic and abiotic processes that vary across different scales of space and time (Dodson et al. 2014, Kern et al. 2017. Following high-severity, large-scale disturbances, there is a marked decline in host species seed production yet an increase in available resources for competitors. Most tree species do not maintain a long-lived seed bank, and abundant seed production does not begin until individuals are large Shepperd 1984, Burns andHonkala 1990). Moreover, competition with associated non-host trees can limit post-outbreak regeneration (Veblen et al. 1991, DeRose and Long 2007, 2010, Dodson et al. 2014, Temperli et al. 2015. In addition to these biotic effects, the removal of closed, continuous canopies also connects the macroclimate to the microclimate affecting understory environmental conditions (Carlson and Groot 1997, Gray et al. 2002, Lazarus et al. 2018. For example, openings in the canopy increase solar radiation in the understory which in turn increases daytime soil temperatures during the growing season (Carlson and Groot 1997, Raymond et al. 2006, Schatz et al. 2012). However, reduction in belowground biomass (e.g., root system) due to canopy tree mortality and associated competition for soil water in gaps may result in higher levels of soil moisture (Gray et al. 2002, Ritter et al. 2005, Burton et al. 2014a. Moreover, increases in snowpack in canopy gaps can further provide a soil moisture subsidy for regeneration (Hardy et al. 1997, Mahat and Tarboton 2014, but see Biederman et al. 2012, Pugh and Gordon 2013. Therefore, while increases in resources following disturbance may boost growth and survival of established seedlings, low seed production, competition, and the relatively severe microclimate may discourage the establishment of new seedlings and increase mortality rates of existing seedlings (De Frenne et al. 2013, Frey et al. 2016. Following large-scale or severe disturbances, altered microclimates can increase the relative importance of facilitation compared to competition (Bertness and Callaway 1994, Callaway 1995, Callaway and Walker 1997, Holmgren et al. 1997. Stand structural characteristics (live canopy, snags, and/or downed wood) can modify local environmental conditions-thereby influencing seedling establishment patterns through buffering of soil and air temperature (Germino et al. 2002) and/or microsite heterogeneity (Wild et al. 2014). Therefore, postdisturbance regeneration could be mediated by interactions with snags, seedlings, and/or downed wood (Marzano et al. 2013). These interactions may be particularly strong at sites with more extreme climatic conditions (Callaway et al. 2002), which makes it important to understand these processes across bark beetle outbreaks and other widespread disturbances that may span the entire range of species like Engelmann spruce (Picea engelmannii).
Warmer temperatures can hasten the life cycle of some bark beetles (Dendroctonus spp.), thereby increasing the frequency and severity of large-scale bark beetle mortality events in response to global warming (Paine et al. 1997, Hansen et al. 2001, Hebertson and Jenkins 2008, Raffa et al. 2008, Bentz et al. 2010. Bark beetle epidemics have affected over 47 million ha in North America between 1998 and 2008 (Raffa et al. 2008), with that area continuing to ❖ www.esajournals.org 2 November 2019 ❖ Volume 10(11) ❖ Article e02912 increase. Recent spruce beetle (Dendroctonus rufipennis) epidemics have devastated over 2.5 million ha of spruce forests in the western United States and Alaska from 1997 to 2013 (USDA 2009(USDA , 2010(USDA , 2012(USDA , 2015. An epidemic spruce beetle outbreak on the Markagunt Plateau in southern Utah that occurred from~1990 to 2000 provides an excellent opportunity to study post-disturbance regeneration processes and interactions with climate. The Markagunt Plateau is located near the southern range limit for Engelmann spruce. High-elevation forests of spruce and subalpine fir (Abies lasiocarpa) occupy sites spanning a broad gradient in elevation and climate conditions (Appendix S1: The main goal of this study was to determine how large-scale, climate-mediated disturbances alter regeneration patterns across gradients in moisture availability and vegetation conditions as indicated by changes in (1) drivers of seedling/sapling density and (2) interactions between seedlings and snags, conspecific seedlings, and (3) downed wood. To address these questions, we compare the density and spatial patterns of Engelmann spruce regeneration that established prior to the outbreak (pre-outbreak, i.e., advance regeneration) to seedlings that established during and after the outbreak (combined into one group, hereby referred to as post-outbreak seedlings). We expected the distribution of post-outbreak seedlings to be more sensitive to recent climatic factors, such as soil moisture and temperature, than advance regeneration resulting in relatively low post-outbreak seedling densities, particularly on more stressful sites (Obj. 1). However, at these more stressful sites, we expected to find an effect of facilitation on post-outbreak seedling establishment. Specifically, we expected post-outbreak seedlings to be clustered around Engelmann spruce snags and seedlings (Obj. 2), and more prevalent on the north side of downed wood (Obj. 3), an area that could buffer climate extremes after losses in canopy cover.

Study area
The Markagunt Plateau is located in the Dixie National Forest in southwestern Utah (Fig. 1). Part of the Colorado Plateau, this high-elevation area experienced an epidemic spruce beetle outbreak that started in~1987 and over the next decade led to over 95% mortality of mature Engelmann spruce across the landscape (DeRose and Long 2007, Pettit 2018. Prior to the outbreak, forests ranging from about 2700 to 3400 m asl (meters above sea level) were dominated by Engelmann spruce with common associates that included subalpine fir and aspen (Populus tremuloides), and less common associates Douglas fir (Pseudotsuga menziesii) and limber pine (Pinus flexilis).
The Markagunt Plateau is near the southwestern edge of the geographical range of Engelmann spruce: an area where future projections indicate the loss of suitable habitat for the species based on bioclimatic envelopes (Rehfeldt et al. 2006). Climate on the Markagunt Plateau is characterized as humid continental (Gillies and Ramsey 2009) with mean summer (JJA) temperatures ranging from 8.5°to 14.8°C and winters (DJF) from À3.0°t o À8.0°C from 1986 to 2016 across all study sites (Wang et al. 2016). Average annual precipitation on the plateau across all study sites from 1986 to 2016 was 760 mm with an average of 54% falling as snow (Wang et al. 2016). Precipitation on the plateau was bimodal with monsoonal rain from the southwest in the summer and snowfall from the Pacific in the winter (Mock 1996). Climatic moisture deficit (CMD), defined as potential evapotranspiration minus precipitation, ranged from 179 mm (low deficit) to 416 mm (high deficit) across our sites (Wang et al. 2016).

Data collection
Sixteen plots were sampled across the Markagunt Plateau in the summer of 2017 (Fig. 1). Nine plots were located at sites used during previous studies on the plateau by DeRose and Long (2007, 2010) that were still accessible and had not been disturbed by human activity (e.g., two of their original 11 plots were excluded due to logging and recent fire activity). Seven additional plots were added to expand the range of climate variability captured (i.e., more sites at the extreme high and low elevations). To maintain continuity in sampling design, we followed similar plot selection criteria as described in DeRose and Long (2012), by randomly locating new stands in areas that were dominant in Engelmann spruce composition prior to the outbreak, within a systematic grid that covered the entire plateau. We collected vegetation data at all 16 plots following the design of the Forest Inventory and Analysis (FIA; USDA Forest Service 2016). Each plot consisted of four circular subplots (168 m², 7.31 m radius) where all trees (≥12.7 cm diameter at breast height, dbh) were measured for species, status, dbh, and height, and mapped to the subplot center. Across each subplot, we also counted the number of saplings (>2.54 cm and <12.7 cm dbh) per species. Additionally, at each subplot the slope, aspect, and elevation were recorded. Six supplementary plots were added at the plots with the most extreme environments (three highest and three lowest CMD) at least 50 m from the original plot footprint to increase the sample depth for the spatial analysis (see Seedling interactions with snags and conspecific seedlings; Fig. 1). On these supplementary plots, all vegetation measurements were taken but no microclimate data were collected. The study plots spanned an elevation gradient of 2708-3307 m, which corresponded to a CMD gradient of 179-416 mm (Appendix S1: Table S1).
In addition to the FIA protocol, we measured and mapped Engelmann spruce seedlings (<2.54 cm dbh) on the subplot. We recorded height and basal diameter of each seedling and estimated age to the decade by counting bud scars (Niklasson 2002, Zielonka 2006, Ba ce et al. 2011. From the decadal age estimates, we retrospectively divided seedlings into two categories based on timing of establishment. Pre-outbreak seedlings were defined as those that established prior to the start of the epidemic outbreak which included seedlings in decade classes above thirty (i.e., established prior to 1987; for the density models, ❖ www.esajournals.org 4 November 2019 ❖ Volume 10(11) ❖ Article e02912 saplings were included with pre-outbreak seedlings and thus referred to as advance regeneration, see Density models of seedlings/saplings). Postoutbreak seedlings were established during and/or after the start of the outbreak and thus were in the decade classes up to thirty (i.e., established 1987present). Bud scar counts are most accurate for younger-aged seedlings, including the time period that encompassed where we distinguished between pre-and post-outbreak seedlings. Putatively less reliable age estimates for seedlings of advanced ages (greater than three decades) in the understory are of less concern because these seedlings were established long before the outbreak and therefore were correctly placed in the pre-outbreak cohort. We also recorded the substrate type where the seedling was rooted based on four major groups (duff, log, mineral soil, and rock). An additional microplot (13.5 m², 2.07 m radius) was established around each seedling, where we mapped and measured all saplings, downed wood, and any trees that fell outside of the subplot, as well as tallied all seedlings by species. Decay classes were recorded for each log following FIA protocol (USDA [U.S. Department of Agriculture] Forest Service 2016), and percent cover by vegetation, early-stage decay wood, and late-stage decay wood was estimated (percent categories: 0, 1, 5, 15, 25, 50, 75, 85, 95, 99, and 100).
To quantify microclimate variability among subplots, we measured temperature and volumetric soil water content (SWC) using Thermochron iButtons (Maxim Integrated, California, USA) and a FieldScout portable time-domain reflectometer (Spectrum Technologies, Illinois, USA), respectively. iButtons were deployed in the center of each subplot 30 cm above the ground on a PVC pipe using reflective tape and insulation, and were left on site for three months (July-September). SWC was measured in nine places in each subplot (in the subplot center and at 4 m in each cardinal and intercardinal directions). SWC measurements were repeated every month for three months (July-September). Measurements were taken during a one-week window across all plots in each month of sampling with the notable exception for the month of July at plot BHF which was inaccessible due to the Brian Head Fire until late in July (round 1, 6-24 July; round 2, 31 July-3 August; round 3, 20-22 September).

Soil moisture and microclimate
To explore fine-scale climatic drivers of patterns in seedling density, we calculated subplotspecific values of soil moisture and microclimate CMD to be used in our seedling density models (see Density models of seedlings/saplings). In many forests of the western United States, SWC typically declines over the course of the growing season as water accumulated during winter is depleted by transpiring vegetation without being replenished by precipitation and/or snowmelt (e.g., Gray and Spies 1996). However, in the southwestern United States, summer precipitation in the form of monsoon rains makes this pattern less predictable. Additionally, SWC can vary spatially with gap size, within-gap position, and substrate type (Gray et al. 2002, Raymond et al. 2006, Burton et al. 2014a). To account for spatiotemporal trends in SWC, we created a linear model of SWC with the unique site/subplot combination, Julian day, time since last rain, and time of day as predictors (Appendix S2: Fig. S1). We then characterized subplot-to-subplot variation in SWC with model estimates at each subplot at the beginning and end of the measurement period when differences across sites were the largest (day 187 and day 265; Appendix S2: Fig. S1) while holding time since last rain and time-ofday at their mean (24 h and 12:00, respectively; for complete model selection process and model results, see Appendix S2).
From the iButton temperature data, we calculated microclimate CMD, by subtracting precipitation from potential evapotranspiration, for July, August, and September at each subplot (Appendix S1: Fig. S1). Potential evapotranspiration was estimated following the Hargreaves method (Yates and Strzepek 1994) using downscaled monthly radiation values in addition to the iButton temperature data measured on each subplot. Precipitation values were obtained from ClimateWNA for each plot center (Wang et al. 2016).

Density models of seedlings/saplings
To assess how the beetle outbreak has altered the drivers of seedling/sapling density (Obj 1), we developed models of Engelmann spruce understory density (seedlings and saplings) measured across sixteen stands in five steps that accounted for a hierarchy of controls (Table 1; ❖ www.esajournals.org 5 November 2019 ❖ Volume 10(11) ❖ Article e02912 Burton et al. 2014b). For the density analysis, we included saplings in addition to seedlings in order to better characterize pre-outbreak densities, that is, advance regeneration, but due to sampling constraints, all other analyses included only seedlings. In the first step, we assessed alternative models of fine-scale abiotic variables including the microclimate conditions of each subplot, followed by the subplot average of biotic measurements recorded in the 2.07 m radius microplots around each seedling (Table 1). In step two, we expanded the scope to include understory and overstory subplot-level metrics characterizing forest structure (Table 1). Finally, we included landscape-level variables accounting for the impact of the macroclimate (subplotlevel slope, aspect, and elevation, and plot-level CMD from ClimateWNA; Wang et al. 2016; Table 1). In each step, we considered alternative models with all plausible variables and biological interactions within that level. We selected the best model from each step based on Akaike's information criterion (AIC) before proceeding to the next step where alternative models with additional variables and interactions were added to the model selected in the previous step (Burnham and Anderson 2002, Burton et al. 2014b).
Interactions across levels were only tested if the main effect from the previous step was already present in the best-selected model. Models were built from the micro-to macroscale because our system comprises a single vegetation type, and as such, we expect microscale variables to generally be more important than macro. Landscapescale climate variables may only become important once controlling for these microscale variables. Potential fixed effects between levels were poorly correlated (Appendix S1: Fig. S2), and thus, the mechanistic interpretation of fine-scale fixed effects was likely not due to landscape-level processes. We developed models for (1) the density of Engelmann spruce seedlings and saplings that established prior to the beetle outbreak, hereby referred to as advance regeneration (saplings were assumed to have originated prior to the outbreak) and (2) the density of Engelmann spruce seedlings that established during or after the beetle outbreak (post-outbreak seedlings; for full model selection steps and AIC values, see Pettit 2018). Although saplings did not have their own microplot for the fine-scale biotic measurements, the subplot average of pre-outbreak seedling microplots was used to characterize the entire subplot. For both models, we used a negative binomial generalized linear model due to detected overdispersion when fitting a Poisson distribution model. Since subplots were separated by a distance greater than the on-site tree height potential (35 m), we used the subplot as the sampling unit (n = 64). When we tested the random effect of plot, it did not improve the model performance (AIC increased by ≤2) and model results did not change (see Appendix S1: Table S2 for model results with random effect of plot); thus, final models do not include random effects (Zuur et al. 2009). After final models were selected, we plotted the marginal effect of each predictor while holding all other terms fixed at their mean. Negative binomial models were

Seedling interactions with snags and conspecific seedlings
To assess how large-scale, climate-mediated disturbances alter regeneration patterns as indicated by changes in interactions (Obj. 2), we analyzed the spatial relationship of seedlings with snags and conspecific seedlings. Two main mechanisms responsible for driving spatial patterns in community structure are environmental heterogeneity (e.g., soil characteristics) and species interactions (e.g., facilitation/competition). The dominant process can be difficult to distinguish when analyzing the spatial patterns of seedlings to other conspecific seedlings; however, the spatial patterns of seedlings to other structural components (snags and downed wood) are suggestive of facilitation through amelioration of environmental conditions (Wild et al. 2014). In harsh environments, like those present postdisturbance, spatial clumping of seedlings to other conspecific seedlings has been interpreted as facilitation (Casalini and Bisigato 2018). Thus, the importance of facilitation can be inferred from the spatial clumping of seedlings and their spatial associations with biological legacies such as snags. Snags, if killed during the outbreak, could have functioned as competitors for preoutbreak seedlings, whereas for post-outbreak seedlings, they may have moderated the harsh microclimate following the death of the canopy trees. To examine where seedlings were spatially random, clustered, or uniform, we analyzed univariate and bivariate versions of the pair correlation function, g(r), on the point patterns of seedlings and Engelmann spruce snags across twenty-two plots (16 main plots and 6 supplementary plots; Fig. 1) together, and also in three CMD categories based on natural breaks in the distribution of CMD values across plots: high (304-418 mm); moderate (243-275 mm); and low (179-205 mm). The pair correlation function computes the average number of points lying within a certain distance (r) of each data point compared to a null model of complete spatial randomness, where g(r) = 1 (Baddeley et al. 2016). Values larger than one indicate more points at that distance than random (i.e., clumped), and values less than one indicate a more uniform distribution. The pair correlation function is similar to Ripley's K but distinguishes between specific interpoint distances instead of computing a cumulative average, making it a more robust analysis especially at close distances (Wiegand and Moloney 2004, Perry et al. 2006, Law et al. 2009, Wild et al. 2014. We first analyzed the univariate version of the pair correlation function to determine the relationships of seedlings. To compare the spatial relationship of seedlings to Engelmann spruce snags, we used a bivariate version of the pair correlation function. Both univariate and bivariate point pattern analyses were analyzed using the spatstat package (Baddeley and Turner 2005) in R (R Core Team 2016) with confidence envelopes based on 99 simulations of complete spatial randomness. Because of the small area encompassed by each FIA subplot and the sensitivity of spatial analyses to plot size, all four subplots at one site were combined into one larger continuous spatial plot by truncating the edges of each circular plot and combining the four smaller square plots into one larger plot 20.85 9 20.85 m. This approach has been found to maintain the forest spatial structure in stands with random and uniform tree distributions (Woodall and Graham 2004). Seedlings established pre-and postoutbreak were modeled together and also independently to test for differences associated with the outbreak. Bivariate and univariate pair correlation analyses were pooled for all sites and for sites by CMD categories. Pooling allows for comparison of replicated spatial point patterns created as a result of the same spatial process. It takes the weighted average of each individual analysis and accounts for within-group variation by calculating confidence intervals (Baddeley et al. 2016).

Seedling interactions with downed wood
To assess how the beetle outbreak has changed interactions between seedlings and downed wood (Obj. 3), we tested for differences in the proportion of (1) pre-and (2) post-outbreak seedlings per subplot in categorical directions (centering on N, S, E, and W) from downed wood in two separate nonparametric tests. Using a Kruskal-Wallis test followed by a post hoc ❖ www.esajournals.org 7 November 2019 ❖ Volume 10(11) ❖ Article e02912 Dunn's test, we determined which parameter pairs were significant for pre-and post-outbreak seedlings separately. For this analysis, three criteria were used: First, only seedlings that had a piece of downed wood in their 2.07 m radius microplot were included. Second, that downed wood must have been down before the seedling established, as indicated by log decay class and age of individual seedlings. Third, if multiple logs met the second condition, only the direction from the closest log was included when calculating proportions. To visualize the spatial relationships, we created circular histograms (rose diagrams) matching the cardinal directions for pre-outbreak and post-outbreak seedlings.

RESULTS
A total of 505 Engelmann spruce seedlings/ saplings were documented across all 16 plots; 302 from before the outbreak and 203 during/ after the outbreak (60% and 40%, respectively; Appendix S1: Table S3). The greatest number of spruce seedlings/saplings was found on the SNO site (144; 2143 stems/ha) and the least at SPC (3; 45 stems/ha; Appendix S1: Table S3). Most seedlings were found growing on mineral soil; however, when subdivided by timing of establishment, a larger proportion of post-outbreak seedlings were found on downed wood compared to pre-outbreak (Appendix S1: Table S3).

Density models of seedlings/saplings
The final model of advance regeneration (D-squared = 52%, Table 2) indicated a negative relationship of seedling/sapling density to basal area of snags (Fig. 2a) and percent cover by early-stage decay downed wood (Fig. 2b). The advance regeneration model also included a positive relationship of seedling/sapling density to percent understory cover (Fig. 2c), and SWC on Julian day 265 (September 22nd; Fig. 2d) which marks the transition from rain to snow on the Markagunt Plateau (Appendix S1: Fig. S3). The final model of post-outbreak seedlings (D-squared = 27%, Table 2) indicated a positive relationship of seedling density to microclimate CMD (Fig. 3a), percent cover by late-stage decay downed wood (Fig. 3b), and SWC on Julian day 187 (July 6th; Fig. 3c).

Seedling interactions with snags and conspecific seedlings
The univariate g(r) for seedlings indicated clumping regardless of timing of establishment (Appendix S1 : Fig. S4). The bivariate g(r) for all seedlings in relation to snags indicated uniformity at close distances (1-2 m), mainly driven by pre-outbreak seedlings (Fig. 4a, b). The pre-outbreak seedlings had a uniform arrangement from Engelmann spruce snags at close distances (1-2 m) with slight clumping at larger distances (3-4 m; Fig. 4b). Post-outbreak seedlings were randomly spaced from Engelmann spruce snags at all distances (Fig. 4c). When the bivariate relationship of post-outbreak seedlings to Engelmann spruce snags was grouped by CMD, there was a uniform arrangement at high and mid CMD sites, but random at low CMD sites (Fig. 5a). The bivariate relationship between postoutbreak seedlings and pre-outbreak seedlings indicated clumping at low CMD sites (Fig. 5b).

Seedling interactions with downed wood
The proportion of seedlings in different directions from downed wood only varied for pre-outbreak seedlings (P < 0.1) not for Notes: Soil water content in September and July (SWC.265 and SWC.187), percent understory cover as measured on the seedling microplot, percent cover by early-stage decay downed wood (% Early-stage DW), percent cover by latestage decay downed wood (% Late-stage DW), snag basal area (Dead BA), and microclimate climatic moisture deficit (CMD). For full model selection process, see Pettit (2018).
❖ www.esajournals.org 8 November 2019 ❖ Volume 10(11) ❖ Article e02912 post-outbreak seedlings (P = 0.71). Following a post hoc Dunn's test, the largest difference (P = 0.0071) was the increased prevalence of preoutbreak seedlings on the east compared to west side of downed wood (Appendix S1: Fig. S5). However, the rose diagram for combined preand post-outbreak seedlings in relation to downed wood indicated pre-outbreak seedlings distributed on all sides of downed wood, whereas post-outbreak seedlings were found predominately on the west side of logs (Fig. 6).

DISCUSSION
Our results highlight how climate-mediated disturbances can alter abiotic and biotic conditions in the understory, which in turn affect forest regeneration. Post-outbreak forest structure consisted of a combination of seedlings that established prior to and following the disturbance. Before the disturbance, live canopy trees exerted a negative control on seedling abundance. In contrast, seedlings that established after the disturbance were not associated with snags and instead were more abundant at warmer sites with greater climatic moisture deficit, given adequate soil moisture during the growing season. Post-outbreak seedlings were positively influenced by downed wood and seedlings of older cohorts indicating that conditions following canopy removal were perhaps detrimental to seedlings. The high-severity beetle outbreak across our sites altered regeneration patterns for Engelmann spruce compared to pre-disturbance ❖ www.esajournals.org 9 November 2019 ❖ Volume 10(11) ❖ Article e02912 conditions, demonstrating how disturbance interactions with climate and stand structure may determine the future tree species composition of these subalpine forests.

Patterns of advance regeneration are primarily related to historical competition
Undisturbed high-elevation spruce-fir forests are characterized by attenuated light reaching the forest floor and also presumably by increased competition for soil moisture and nitrogen that would inhibit regenerating seedlings (Ritter et al. 2005, Hill et al. 2018. The negative relationships of advance regeneration to snag basal area and early-stage decay downed wood indicated regeneration was largely limited by overstory competition prior to the outbreak (Knapp and Smith 1982, Gray and Spies 1996, Feller 1998. For example, we encountered fewer seedlings/saplings at sites that prior to the outbreak had more standing live trees, particularly Engelmann spruce (Fig. 2). Although tree cover can facilitate Engelmann spruce seedlings at the alpine treeline ecotone (Maher et al. 2005), closed-canopy stands below treeline are more likely to result in competition for resources rather than continued facilitation (Hill et al. 2018). The forests we studied were below treeline, and the spatial distribution of pre-outbreak seedlings we found was more consistent with the patterns expected in undisturbed stands with closed canopies-through competition with extant Engelmann spruce (Svoboda et al. 2010 ; Fig. 4b).
The positive relationship of the density of advance regeneration to late September soil moisture and the spatial distribution of seedlings around downed wood is consistent with recent literature that suggests spruce regeneration is moisture-limited in closed-canopy stands (Andrus et al. 2018, Hill et al. 2018. Seedlings are often rooted in shallow soil layers that can dry out quickly; thus, elevated soil moisture levels late in September, after the summer monsoon, can extend the growing season and enable continued growth and survival of tree seedlings (Padilla and Pugnaire 2007). Similarly, the east side of downed wood may provide a more favorable microsite, not only by shielding seedlings from direct radiation during the hottest part of the day (Kluber et al. 2009), but also by slowing snowmelt and thus elevating soil moisture in the spring before the monsoon when temperatures are generally high and soil moisture levels are low (Cunningham et al. 2006, Hu et al. 2010. While soil moisture levels and other environmental conditions may have differed when seedlings established, relationships observed here may reflect the importance of the balance of moisture and microclimate for post-outbreak seedling survival. Advance regeneration density was also positively related with understory cover, suggesting understory vegetation may function to facilitate spruce regeneration. However, understory vegetation likely responded to disturbance more quickly than seedlings, and cover values may not directly indicate what was present prior to the outbreak. Instead, the relationship may reflect greater local-scale site productivity and/or availability of nutrients, which can provide a net benefit to both seedlings and understory vegetation. Rather than co-occurrence, competition between the two may be manifested as reduced seedling growth (Loranger et al. 2017).

Post-outbreak regeneration is less predictable
Following the gradual removal of the canopy during disturbance-related spruce mortality, light availability generally increases as do extreme climatic conditions (i.e., variation in snowpack, temperature, and CMD; Carlson and Groot 1997, Hardy et al. 1997, Gray et al. 2002, Raymond et al. 2006. Increases in light availability can benefit post-outbreak seedlings given sufficient growing season soil moisture. The density of post-outbreak seedlings increased with soil moisture provided from monsoonal rainfall consistent with the known sensitivity of Engelmann spruce seedlings to drought (Knapp andSmith 1982, Lazarus et al. 2018). However, contrary to moisture limitation, post-outbreak seedling density was positively related to CMD (Fig. 3b). It is important to note this relationship accounts for the positive relationship with soil moisture. Together, this suggests that these high CMD sites contained enough water for seedlings to establish (our model controlled for variation in %SWC), and the positive relationship with high CMD Fig. 4. Bivariate point pattern analysis for all (a), pre-(b), and post-outbreak (c) to Engelmann spruce snags. Black solid line is the bivariate (seedling, tree) pooled isotropic-corrected estimate of g(r). Red dashed line represents the bivariate (seedling, tree) theoretical Poisson g(r) representing complete spatial randomness. Gray bands indicate the bivariate (seedling, tree) upper and lower limits of two-sigma CI based on isotropic-corrected estimate of g(r). Values of g(r) > 1 indicate clustered patterns, whereas values < 1 indicate spatial uniformity.
❖ www.esajournals.org 11 November 2019 ❖ Volume 10(11) ❖ Article e02912 reflects high potential evapotranspiration, or associated temperatures, rather than a moisture limitation. It is also conceivable that increases in temperature and CMD were associated with increased solar radiation in openings, which could also be an important factor affecting spruce seedling establishment and survival (Schatz et al. 2012). With sufficient soil moisture to prevent seedling desiccation, higher light levels coupled with warmer temperature and/or a longer growing season at lower elevations can benefit seedling establishment and survival, consistent with recent findings (Hill et al. 2018).
Whether the CMD relationships that we observed reflect hysteresis in response to strengthening droughts, or simply characterizes the longer-term trends of establishment response, is impossible to evaluate without longer-term climate or establishment data. Increases in climate extremes or resource heterogeneity following the outbreak appear to have caused intraspecific interactions to have been increasingly important for post-outbreak seedlings. Although post-outbreak seedlings were found in clusters (Appendix S1: Fig. S4), these clusters were not found around Engelmann spruce snags (Figs. 5c and 6a) as has been found for P. abies in central Europe (Pomeroy et al. 2009, Ba ce et al. 2011, Wild et al. 2014). Compared to monospecific stands of P. abies, after the Fig. 5. Bivariate point pattern analysis for post-outbreak seedlings to Engelmann spruce snags (a-c) and postoutbreak to pre-outbreak seedlings (d-f). Spatial patterns are divided by categorical macroclimate CMD values: low (179-205 mm; a, d); moderate (243-275 mm; b, e); and high (304-418 mm; c, f). Black solid line is the bivariate (seedling and tree and pre-and post-outbreak seedlings) pooled isotropic-corrected estimate of g(r). Red dashed line represents the bivariate (seedling and tree and pre-and post-outbreak seedlings) theoretical Poisson g(r) representing complete spatial randomness. Gray bands indicate the bivariate (seedling and tree and pre-and post-outbreak seedlings) upper and lower limits of two-sigma CI based on isotropic-corrected estimate of g(r). Values of g(r) > 1 indicate clustered patterns, whereas values < 1 indicate spatial uniformity. Note the difference in y-axis scale from top row to bottom. ❖ www.esajournals.org 12 November 2019 ❖ Volume 10(11) ❖ Article e02912 Markagunt Plateau beetle outbreak subalpine fir was by far the most abundant tree species across our study sites. We also regularly observed regeneration of subalpine fir in tight clusters around Engelmann spruce snags, indicating that seeds of this species are overwhelmingly abundant (Appendix S1: Fig. S6; DeRose and Long 2010). In contrast, post-outbreak Engelmann spruce seedlings relied on interactions with cohorts of older Engelmann spruce seedlings and downed wood. Engelmann spruce seedlings established preferentially on downed wood (Appendix S1: Table S3). Its use as a substrate could offset any soil moisture limitations given the increased water-holding capacity of decayed logs (Knapp and Smith 1982, Zielonka 2006, Macek et al. 2017. While this is a well-documented strategy, it is important to note that survival on bare mineral soil may be higher over the long term (Noble andAlexander 1977, Windmuller-Campione et al. 2017) due to a lack of successful rooting to acquire soil moisture on downed wood (Grossnickle 2005). The positive relationship between post-outbreak seedling density and microclimate CMD suggests that post-disturbance spruce regeneration is temperature, or energy, limited. This may be due not only to temperature, but also associated effects on snow depth and retention. In contrast to the closed-canopy forest understory, where increased snowpack provides a moisture subsidy for spruce seedlings (Andrus et al. 2018, Hill et al. 2018, in post-disturbance canopy openings snow may inhibit seedlings from becoming established and effectively shorten the growing season (Hardy et al. 1997). The increased abundance of post-outbreak seedlings at sites with high CMD (warmer temps), on the W/SW side of downed wood, and clustering around pre-outbreak seedlings at high-elevation sites where snowpack would be deepest and persist longest are all indications of the importance Fig. 6. Spatial histogram of pre-outbreak compared to post-outbreak seedlings from closest piece of downed wood that would have been down at the time of seedling establishment. Histogram bin width is 5 degrees, and y-axis scale represents 2, 4, 6, 8, and 10 seedlings per bin.
❖ www.esajournals.org 13 November 2019 ❖ Volume 10(11) ❖ Article e02912 of extended growing season due to lower albedo and greater local heating. Similar to Wild et al. (2014), this explanation also invokes the same mechanistic process to seed capture via early snowmelt and subsequent establishment on moist mineral soil free of competition, but implicates pre-outbreak seedlings and/or downed wood instead of snags. Although the spatial patterns of seedlings can be influenced by multiple processes (seed dispersal limitations, soil properties, and facilitation/competition), we argue that the combined spatial pattern of post-outbreak seedlings to themselves, to older cohorts, and downed wood all point to the microclimate amelioration provided by these various structures in a more severe post-disturbance environment. In order to better understand the mechanistic drivers of seedling regeneration patterns, a more direct measurement of spatial patterns of snow depth/melt throughout the winter/spring and its effects on seedling photosynthesis, and/or more high-resolution information about edaphic factors are needed.

Differences between pre-and post-outbreak seedlings likely do not reflect life-history stages
Our study analyzed the density of extant Engelmann spruce seedlings and saplings, and, therefore, has life-history stage limitations. In the absence of a spruce beetle outbreak, the density of advance regeneration would be indicative of future forest composition. This would be accompanied by predictable changes in conditions that occur through the establishment and growth of seedlings (Schupp 1995, Raymond et al. 2006. However, given the disturbance, the difference in factors driving density patterns that we observed between the advance regeneration and postoutbreak seedlings suggests that establishment patterns dictating future spruce forests were strongly modified by the epidemic beetle outbreak. As a result, we do not expect advance regeneration to be the sole indicator of future spruce forests; rather, post-disturbance changes in microsite that influence seedling establishment will need to be considered. In combination, the differences in seedling establishment that were created as a result of the epidemic beetle outbreak, and the continued impact of climate change on disturbance regimes, are likely to alter future Engelmann spruce regeneration ecology.

Management implications
Previous research has predicted that Engelmann spruce in ecosystems similar to the Markagunt Plateau may be at risk of extirpation based on climate predictions (Rehfeldt et al. 2006). However, our results indicated some level of spruce establishment across a range of environmental conditions. While it appears that this system was not resistant to an epidemic outbreak (sensu DeRose and Long 2014), a combination of advance regeneration and somewhat limited post-outbreak spruce establishment imparted some level of resilience~20 yr after the outbreak. If planning goals include maintaining Engelmann spruce in areas that have not yet been affected by an outbreak, management practices could be focused on creating resilience in terms of promoting advance regeneration-easily achieved through silvicultural treatments (Windmuller-Campione and Long 2015, Windmuller-Campione et al. 2017). However, if Engelmann spruce is desired in the wake of large outbreaks, where altered regeneration processes have the potential to influence future species composition, post-outbreak microclimate conditions should be considered when planning salvage harvests or planting operations. In other climatically limited systems, projected increases in temperature are expected to reduce low-elevation habitat (Koo et al. 2015). However, given adequate soil moisture in monsoon-influenced regions, even lower elevations with high moisture deficit could be suitable habitat for Engelmann spruce seedlings whether regenerated naturally or planted. Regardless of the approach, careful consideration of the microsite during management should allow for persistence of spruce at lower elevations despite warming from climate change (Kueppers et al. 2017).