Severity of Overstory Mortality Influences Conifer Recruitment and Growth in Mountain Pine Beetle-Affected Forests

The severity of lodgepole pine mortality from mountain pine beetle outbreaks varies with host tree diameter, density, and other structural characteristics, influencing subcanopy conditions and tree regeneration. We measured density and leader growth of shade-intolerant lodgepole pine, shade-tolerant Engelmann spruce, and very shade-tolerant subalpine fir regeneration beneath stands that experienced moderate and high overstory lodgepole pine mortality (average 40% and 85% of total basal area) a decade earlier. Lodgepole comprised >90% of the overstory basal area and mature spruce and fir were present in both mortality levels, though live basal area and disturbance history differed. Post-beetle outbreak recruitment was high in both mortality levels, but there were more lodgepole in high than moderate mortality plots (1140 stems ha−1 vs. 60 stems ha−1) and more subalpine fir in moderate than high mortality plots (4690 stems ha−1 vs. 2870 stems ha−1). Pine advance regeneration, established prior to outbreak, was more dense in high mortality than moderate mortality sites (930 stems ha−1 vs. 310 stems ha−1), but the trend was generally the opposite for the other conifers. Lodgepole recruitment increased and subalpine fir decreased with greater forest floor light availability. All species grew faster in high mortality areas than their counterparts in moderate mortality areas. However, in high mortality areas pine grew faster than the more shade tolerant species, and in moderate mortality areas spruce and fir grew faster than pine. These species-specific responses to the degree of overstory mortality will influence future stand composition and rate of forest recovery after mountain pine beetle outbreaks.


Introduction
Bark beetles killed between 10-100 percent of trees in stands affected during recent outbreaks in western North America [1] creating a range of understory conditions that will influence tree regeneration and future forest composition. Mountain pine beetles (MPB) (Dendroctonus ponderosae Hopkins) have killed lodgepole pine (Pinus contorta var. latifolia Dougl. ex. Loud.) on 10 million hectares in North America [2] and over one million hectares in Colorado and southern Wyoming since the mid-1990s [3]. Approximately 55% of beetle-affected forests in Colorado that contain lodgepole pine are mixed with shade-tolerant Engelmann spruce (Picea engelmannii Parry ex Engelm.) and subalpine fir (Abies lasiocarpa var. lasiocarpa (Hook.) Nutt.) [4]. Though it is well-established that forest canopy density regulates understory light, temperature, and other factors that determine tree recruitment and growth [5][6][7][8], the influence of levels of beetle-related forest mortality on these conditions and

Study Area
The study was conducted at the USDA Forest Service Fraser Experimental Forest (FEF) in north-central Colorado, in forests comprised of lodgepole pine (>90% of basal area before beetle-caused mortality) mixed with Engelmann spruce and subalpine fir (1%-5% of basal area each). Widespread bark beetle activity began in the area around 2002 [26] and pine mortality was generally higher in old-growth forests with larger diameter trees than in younger forests with smaller trees [27]. The study area receives 610 mm of precipitation annually (range 427-902 mm, 1981-2010), with nearly two-thirds falling as snow from October to May (USFS unpublished data 1981-2010 data). The average annual temperature is approximately 3 • C, and annual average minimum and maximum temperatures are −3.7 and 9.6 • C (PRISM Climate Group). Soils are dominated by skeletal, sandy loam Dystric and Typic Cryochrepts with 20%-30% gravel and 30%-50% cobbles [26,28].

Study Design
The study was conducted in a portion of FEF that contained high and moderate overstory mortality levels on sites with relatively consistent aspect (generally north-west), elevation (2780 m and 2930 m), slope (4%-30%), and forest composition. Old-growth lodgepole pine in the area regenerated after wildfire approximately 300 years before present [29] and second-growth stands regenerated after fire or timber harvesting that occurred between 1900 and 1940 [30,31]. Candidate stands were identified using high resolution satellite images (Google, Mountain View, CA, USA). We chose stands with the highest visible mortality from imagery as 'high mortality' and chose the lowest mortality areas nearby as 'moderate mortality' stands. Stands were then visited to confirm mortality level observed on imagery was present, pine-dominance, and spruce and fir presence. We located 10 pairs of high and moderate mortality stands. Within each pair, high and moderate mortality stands were located within 300 m of another on similar slope, aspect, and landform position to control for variation in soil and other conditions as much as possible. In each stand we established four, 40 m 2 circular plots at with a center randomly selected using ArcGIS 10.0 (ESRI, Redlands, CA, USA). Plots were clustered with one central plot and the others 30 m from the central plot at 120 degree intervals. There were seed sources of all conifer species within 15 m of each plot center.

Tree Data Collection
At each plot, we measured overstory basal area for all species using a 4.6 m 2 ha −1 basal area factor prism. Tree dbh (diameter at breast height (1.4 m)), status (live or dead), and mortality agent were recorded. Pre-beetle forest conditions were calculated by adding beetle-killed and live trees. Density of post-beetle tree recruitment and advance regeneration (trees established prior to the outbreak) was collected by species on 40 m 2 circular plots. Post-beetle recruits were identified as trees that had less than eight bud scars. Advance regeneration included all trees with eight or more bud scars that were also less than 1.4 m tall. The forest floor substrate on which each recruit had established was also recorded based on the 0.04 m 2 area at its base (or directly above the base up to the tree's height). We estimated the percent cover of the following forest floor substrates across each 40 m 2 plot: duff/litter, mineral soil, rock, woody debris, lichen/moss, forb/grass, Vaccinium spp., and other shrubs/conifer branches. Annual leader growth for 2012 was measured on one to four trees less than 1.4 m tall (recruitment and advance regeneration combined) of each species closest to plot center and up to 7.5 m from plot center. Trees with leader damage were excluded.

Canopy Density Data Collection and Processing
Hemispherical photographs were used to estimate canopy density as canopy gap fraction (CGF) using a Canon EOS 50D 15.1 Megapixel Digital SLR camera (Canon U.S.A., Inc., Melville, LA, USA) with a Sigma 4.5 mm F2.8 EX DC HSM Circular Fisheye Lens (Sigma Corporation of America, Ronkonkoma, NY, USA). Canopy gap fraction is the portion of the sky visible from a single point beneath the canopy. Photos were taken pre-dawn or post-dusk with the camera leveled at 1.4 m height at plot center. Shutter speed was 1/125th second, and exposure levels were manually set so that sky patches near the zenith were overexposed by 1-2 exposure stops [32].
Photos were processed with HemiView 2.1 SR4 Canopy Analysis Software (1998,2009 Delta-T Devices Ltd., Cambridge, UK). Because setting the threshold for canopy versus sky is subjective, photos were analyzed repeatedly until the proportion of visible sky was within 0.1 of the previous estimate [32]. We also used the following visual standard: the canopy vs. sky threshold was determined so that grouped cones on dead lodgepole pine were classified as canopy, and small diameter, solitary branches were classified as sky. To minimize topographic differences among plots we masked the hemispherical field of view at the greatest slope of the study plots (17 degrees) and determined visible sky from a 146 degree field of view [33].

Soil Temperature and Moisture
Soil temperatures at 1 cm depth below the litter/duff surface at plot center were recorded hourly from 1 July 2013 until first snowfall (21 September 2013) with Thermochron ® iButtons (Maxim Integrated™, San Jose, CA, USA). We calculated daily maximum, minimum, and mean temperatures during the sampling period.
We measured volumetric water content of the top 12 cm of soil every two weeks with a hand-held probe (Hydrosense II, Campbell Scientific, Logan, UT, USA). We sampled on rain-free days and measured soil moisture of all plots on the same day and in each high and moderate pair within 30 minutes of each other. At each plot four measurements were taken 0.3 m from plot center and composited.

Statistical Analysis
We tested the effect of mortality level (moderate vs. high pine mortality) at the stand level on preand post-outbreak basal area, recruitment, and advance regeneration density, and 2012 leader growth using PROC GLIMMIX (SAS, Version 9.3, Cary, NC, USA). We determined basal area and growth fit a normal distribution and recruitment/advance regeneration fit a Poisson distribution and used these in our models; we accounted for random effects of pair and plot within mortality level. We analyzed the effect of canopy gap fraction on species' recruitment densities and leader growth at the plot level with a linear regression model (PROC REG, SAS, Version 9.3, Cary, NC, USA).
We compared the forest floor substrates between high and moderate mortality levels using GLIMMIX, accounting for random effects of pair and plot within mortality level. For each tree species, we compared the proportion of observed trees on each substrate and the average substrate occurrence on all plots using the chi-square test of PROC FREQ (SAS, Version 9.3, Cary, NC, USA). We removed rock and down wood cover from the chi-square analysis to account for the portion of plots that were not suitable growing space. We tested effect of mortality level on temperatures and moisture with PROC GLIMMIX, accounting for random effects of pair and plot within mortality level. We analyzed the relationship between canopy gap fraction and soil surface temperatures and moistures with PROC REG.

Forest Overstory and Understory
Total pre-beetle basal area was not significantly different between high (33 m 2 ha −1 ) and moderate (32 m 2 ha −1 ) mortality areas, though there was 5 m 2 ha −1 more lodgepole pine in high mortality areas, and 2 m 2 ha −1 more subalpine fir in moderate mortality stands (Figure 1a). High mortality areas lost 85% of their basal area, while moderate areas had 40% mortality ( Figure 1b). Live basal area was 19 m 2 ha −1 in moderate mortality stands, almost four times greater than in high mortality stands. High mortality plots had 24 m 2 ha −1 dead basal area, nearly twice that of moderate mortality stands. Canopy gap fraction was much greater in high mortality stands (meaning canopy density was lower); canopy gap fraction was 0.47 in high mortality stands, 0.22 greater than in moderate mortality stands. Live basal area was a significant predictor of canopy gap fraction (R 2 = 0.38, p < 0.0001).
Forests 2018, 9, x FOR PEER REVIEW 4 of 12 used these in our models; we accounted for random effects of pair and plot within mortality level.
We analyzed the effect of canopy gap fraction on species' recruitment densities and leader growth at the plot level with a linear regression model (PROC REG, SAS, Version 9.3, Cary, NC, USA). We compared the forest floor substrates between high and moderate mortality levels using GLIMMIX, accounting for random effects of pair and plot within mortality level. For each tree species, we compared the proportion of observed trees on each substrate and the average substrate occurrence on all plots using the chi-square test of PROC FREQ (SAS, Version 9.3, Cary, NC, USA). We removed rock and down wood cover from the chi-square analysis to account for the portion of plots that were not suitable growing space. We tested effect of mortality level on temperatures and moisture with PROC GLIMMIX, accounting for random effects of pair and plot within mortality level. We analyzed the relationship between canopy gap fraction and soil surface temperatures and moistures with PROC REG.

Forest Overstory and Understory
Total pre-beetle basal area was not significantly different between high (33 m 2 ha −1 ) and moderate (32 m 2 ha −1 ) mortality areas, though there was 5 m 2 ha −1 more lodgepole pine in high mortality areas, and 2 m 2 ha −1 more subalpine fir in moderate mortality stands (Figure 1a). High mortality areas lost 85% of their basal area, while moderate areas had 40% mortality (Figure 1b). Live basal area was 19 m 2 ha −1 in moderate mortality stands, almost four times greater than in high mortality stands. High mortality plots had 24 m 2 ha −1 dead basal area, nearly twice that of moderate mortality stands. Canopy gap fraction was much greater in high mortality stands (meaning canopy density was lower); canopy gap fraction was 0.47 in high mortality stands, 0.22 greater than in moderate mortality stands. Live basal area was a significant predictor of canopy gap fraction (R 2 = 0.38, p < 0.0001).

Regeneration Density and Growth
There was more pine and spruce recruitment (trees <1.4 m tall, established after beetle), and less fir recruitment in the high mortality than moderate mortality plots (p values all <0.0001) (Figure 2). Advance regeneration density was overall lower than recruitment density. There was greater pine advance regeneration density in high mortality stands (p < 0.0001), but lower spruce regeneration density (p = 0.0342). Fir advance regeneration density did not differ between mortality levels. Regression analysis showed lodgepole recruitment increased with canopy gap fraction (R 2 = 0.22) (Figure 3a), while fir density decreased (R 2 = 0.16) (Figure 3c). Spruce showed no relationship with canopy gap fraction (Figure 3b).

Regeneration Density and Growth
There was more pine and spruce recruitment (trees <1.4 m tall, established after beetle), and less fir recruitment in the high mortality than moderate mortality plots (p values all <0.0001) (Figure 2). Advance regeneration density was overall lower than recruitment density. There was greater pine advance regeneration density in high mortality stands (p < 0.0001), but lower spruce regeneration density (p = 0.0342). Fir advance regeneration density did not differ between mortality levels. Regression analysis showed lodgepole recruitment increased with canopy gap fraction (R 2 = 0.22) (Figure 3a), while fir density decreased (R 2 = 0.16) (Figure 3c). Spruce showed no relationship with canopy gap fraction (Figure 3b).  Leader growth (of recruitment and advance regeneration) was greater in high than moderate mortality areas for all species (p values all <0.0001) (Figure 4). At high mortality levels, lodgepole pine growth was about one-third greater than spruce and fir growth. In moderate mortality, spruce and fir growth was two to three times greater than lodgepole pine growth. Regression analysis reflected the different responses of growth among species to mortality level, as measured by canopy density ( Figure 5). Growth was significantly related to canopy density for all species. This relationship was significantly stronger (slope was greater) for pine (p < 0.0001) than both spruce and fir, which did not differ. The relationship between growth and canopy gap fraction was strong for pine (R 2 = 0.48), moderate for spruce (R 2 = 0.23), and weak for fir (R 2 = 0.06). Leader growth (of recruitment and advance regeneration) was greater in high than moderate mortality areas for all species (p values all <0.0001) (Figure 4). At high mortality levels, lodgepole pine growth was about one-third greater than spruce and fir growth. In moderate mortality, spruce and fir growth was two to three times greater than lodgepole pine growth. Regression analysis reflected the different responses of growth among species to mortality level, as measured by canopy density ( Figure 5). Growth was significantly related to canopy density for all species. This relationship was significantly stronger (slope was greater) for pine (p < 0.0001) than both spruce and fir, which did not differ. The relationship between growth and canopy gap fraction was strong for pine (R 2 = 0.48), moderate for spruce (R 2 = 0.23), and weak for fir (R 2 = 0.06).

Soil Moisture and Temperature
Growing season soil temperature varied between mortality levels though moisture did not. Daily mean and maximum temperatures were greater and minimum temperatures were lower in high than low mortality plots. The overall mean temperature was 1.5 • C greater in the high (13.7 • C) than moderate (12.2 • C) mortality areas (p = 0.0004). The overall daily maximum temperature was 7.0 • C warmer in high (29.9 • C) than moderate (22.9 • C) mortality (p = 0.0002), and average minimum temperature was 0.7 • C cooler in the high (6.4 • C) versus moderate (7.1 • C) mortality plots (p = 0.0119). Canopy gap fraction explained 0.38 and 0.46 of the variation in mean (intercept: 11.0, p < 0.0001; slope: 5.4, p < 0.0001) and maximum (intercept: 16.8, p < 0.0001; slope: 25.7, p < 0.0001) soil temperatures, respectively, but a much smaller amount of the variation (R 2 = 0.16) in minimum temperatures (intercept: 7.7, p < 0.0001; slope: −2.4, p = 0.0019). Soil moisture mean was 6.7% in high mortality plots and 7.0% in moderate mortality plots. The differences in average soil moisture between pairs averaged 1.5% and all pair differences were within the ± 3% sensor accuracy margin of error. There was no relationship between canopy and soil moisture.

Forest Floor and Recruitment Distribution
Forest floor was primarily covered by duff/litter, Vaccinium spp., and wood (Table 1). Down wood cover in high mortality (16%) was nearly double that of the moderate mortality areas (9%) (p < 0.0001) and duff/needles covered 57% of plots in the moderate mortality areas, 16% more than in high mortality (p = 0.0001) ( Table 1). No other substrates differed between mortality levels. Chi-square tests showed that the distribution of recruitment across substrates was different than the observed distribution of substrates covering all plots (p = 0.0260 for pine, p = 0.0007 for spruce, and p = 0.0308 for fir). Recruits of all species were found growing on duff/needles more than expected. Seventy-one percent of lodgepole, 64% of spruce, and 69% of fir recruits were on duff/needles, though this covered only 50% of the forest floor. Eight percent of Engelmann spruce recruits were on moss/lichen, although this substrate covered only 3% of plots. Very few recruits of all species (0.7%, 0%, and 1% for pine, spruce, and fir, respectively) established on wood despite its coverage of 12% of the plots. Asterisks indicate significant differences in substrate between mortality levels (α = 0.05).

Discussion
Our results show that the composition of forests as they develop after MPB outbreaks will be influenced by the level of overstory mortality, associated changes in canopy density, and subcanopy conditions. Post-beetle-outbreak tree recruitment and advance regeneration will both contribute to development of the future forest. Lodgepole recruitment was greater in stands with high compared to moderate mortality. Pine was the fastest growing conifer in more open, high mortality stands, and spruce and fir grew faster under stands affected by moderate overstory mortality. In spite of similar pre-outbreak overstory species composition, recruitment density and leader growth of the conifer species differed between high and moderate mortality stands.
Recruitment density, composition, and growth differed with mortality level and associated changes in canopy structure and abiotic conditions at the forest floor. There were over 1100 pine recruits ha −1 in high mortality areas, where soil surface temperatures and light were high, but few pine recruits in moderate mortality areas. These patterns agree with historic reconstructions for beetle outbreaks in nearby areas [34]. Pine height growth was also about five times greater in the high mortality areas, as was expected from the high light and surface temperature adaptation of this species [35,36]. Fir recruitment was nearly twice as dense in stands with moderate mortality which is consistent with the species shade-tolerance and susceptibility to heat stress [37,38]. Somewhat surprisingly, temperature-sensitive spruce seedlings [39] were slightly denser in high mortality despite the higher soil temperatures and spruce growth increased linearly with light availability. Spruce is less shade tolerant than subalpine fir [40] and typically regenerates within partial canopy openings [41,42]. High mortality stands may have provided an ideal amount of shade for spruce recruitment, even though its density did not vary linearly with canopy density. Considering how sensitive spruce survival is to heat and drought [39,40], variation in microsite could have obscured the effect of canopy opening and further study of the interaction of these factors could improve our understanding of spruce regeneration in beetle-impacted forests.
All three conifer species established successfully on undisturbed forest floor. Litter and duff depth did not exceed 3 cm deep at these sites, similar to lodgepole pine forests elsewhere in Colorado, where pine regeneration was abundant [19,21]. These findings contrast with a regional perception that lodgepole pine requires stand-replacing disturbance and an exposed mineral soil seedbed to regenerate [18,20,23]. Herbaceous and shrub cover was low in our study area overall, unlike observations elsewhere that non-tree vegetation increases following beetle outbreaks [43]. The relative lack of understory cover and competition may have allowed high post-beetle tree regeneration densities compared to sites elsewhere with greater understory density.
In stands with high overstory mortality, pine is likely to outgrow spruce and fir and become dominant, despite the high density of advance regeneration and new recruits of the shade-tolerant conifers [36,37,40]. Conversely, in stands with moderate mortality, fir and spruce grew faster than pine and stand composition will likely shift accordingly. Our findings differ somewhat from previous studies in northern Colorado that did not distinguish between mortality levels and predict a general increase in fir abundance in beetle-killed stands in the absence of logging or further disturbance [16,21]. Our paired comparisons indicate that post-outbreak stand development and species dominance will likely differ between stands with high and moderate levels of overstory mortality. Owing to the regional abundance of forests with mortality comparable to our moderate class (average 40% of basal area mortality [2,27]), it is likely that fir will increase in a large portion of the landscape. However, our evaluation of old-growth stands with high mortality indicates that pine-dominated stands will follow bark beetle outbreaks under some conditions. The species-specific regeneration responses to different levels of overstory mortality have implications for crown fire hazard, wildlife habitat, and other ecosystem processes in the forests that regrow after bark beetle outbreaks [16,27].

Conclusions
These results show that the level of overstory mortality regulates tree regeneration and growth following bark beetle outbreaks in lodgepole pine forests in northern Colorado. Though previous research concludes that subalpine fir will become more dominant after bark beetle in the absence of additional disturbance [18][19][20][21][22][23][24], this study found that stands with high overstory mortality have dense pine recruitment and fast growth and may promote future pine dominance. Our results will help improve understanding of the canopy conditions that create distinct patterns of tree regeneration following bark beetle outbreaks and lead to differences in future forest species composition. Future work should incorporate the effect of varying mortality levels into forest growth simulations to provide more certainty about the long-term trajectory of beetle-affected forests.