Managing Wildlife Habitat: Complex Interactions With Biotic and Abiotic Disturbances

Sustainable forest management strategies include emulating historical disturbance regimes to achieve multiple objectives. Partial-harvesting strategies are used to overcome conflicts between timber production and wildlife habitat conservation; however, the potential impacts on complex disturbance interactions and ecological functions remain largely unknown. In 1984, a controlled experiment was initiated in the dry forests of central British Columbia, Canada, to test partial harvesting intended to enhance mule deer habitat while allowing timber extraction. To determine the short- and long-term impacts on complex disturbance regimes, we quantified changes in forest structure and susceptibility to western spruce budworm, Douglas-fir beetle, and wildfire. We compared structural attributes in 2014 (30 years after the first harvest) and 2015 (1 year after the second harvest) in treated forests, and contrasted them with control forests that were measured in 2015. In the short term (1 year post-harvest), partial harvesting altered forest structure by reducing total canopy cover, subcanopy tree density and basal area, and increasing the abundance of large woody surface fuels. In the long term (30 years post-harvest), the forest canopy attributes did not differ between the treatment and control areas, partly due to increased growth of subcanopy trees. Harvesting had little impact on forest susceptibility to western spruce budworm. Susceptibility to Douglas-fir beetle was lower in the short term due to fewer available mature host trees, but increased to levels similar to the control forest over the long term. Reduced canopy fuels and increased canopy base height decreased the likelihood of crown fire in favor of surface fire. In the long term, canopy fuels and likelihood of crown fire recovered, but woody fuel loads remained low after 30 years. Harvesting to enhance mule deer habitat interacts with biotic and abiotic disturbances in the short and long term. Potential cascading affects depended more on the decision to remove harvesting residuals to mitigate potential Douglas-fir beetle infestations and wildfire than on time since treatment. Provided partial harvesting occurs at intervals ≤ 30 years and residuals are immediately removed, timber extraction and mule deer habitat can be compatible with complex disturbance regimes and sustainable forest management.


INTRODUCTION
In recent decades it has been widely recognized that forest management should seek to emulate the historical range and variability of forest composition and structure to maintain biodiversity and ecological function (Keane et al., 2009;Donato et al., 2020;Čada et al., 2020). This recognition has stimulated efforts to minimize differences between managed and natural forests by modifying harvesting practices to generate spatial and temporal patterns consistent with historical disturbance regimes Harvey et al., 2002;Long, 2009;Kuuluvainen and Grenfell, 2012). However, in regions where extensive forestry is practiced, sustainable timber extraction is only one of many objectives of forest management (e.g., Sandström et al., 2011). Practices capable of addressing multiple objectives are rare (Fortin et al., 2016) and their capacity to maintain ecological function is largely unknown.
Management for multiple objectives may create conflicts such as those that arise under the simultaneous requirements to maintain wildlife habitat and sustain timber extraction (Rees, 2003;Sims et al., 2014). Harvesting to maintain or enhance wildlife habitat may impose unacceptable logistical and timber removal restrictions (e.g., Armleder et al., 1989;Sutherland et al., 2016), whereas unrestricted removal of timber may alter forest structure to the detriment of wildlife (Haggstrom and Kelleyhouse, 1996;Sullivan and Sullivan, 2001;Millington et al., 2011). Furthermore, ecological effects of forest management may not be immediately obvious and can influence long-term forest composition, structure, or function (Franklin et al., 2002;Palik et al., 2002;Hessburg et al., 2016), compromise resilience to subsequent disturbance (Radeloff et al., 2000;Hessburg et al., 2019), and inadvertently create forests prone to unprecedented disturbance severity by insects (Raffa et al., 2008) or megafires (Stephens et al., 2013(Stephens et al., , 2014. The dry interior Douglas-fir [Pseudotsuga menziesii (Mirb.) Franco] forests of western North America are prone to several biotic and abiotic disturbances including insect outbreaks (Raffa et al., 2008;Maclauchlan et al., 2018) and wildfire (Hessburg et al., 2019). These forests also provide many critical ecosystem services including the provision of timber and wildlife habitat (Curtis and Carey, 1996). In central British Columbia (BC), Canada, conservation of mule deer [Odocoileus hemionus Rafinesque] is a high priority (Lehmkuhl et al., 2001), especially in areas where populations have declined due to losses of winter habitat from clearcut harvesting (Armleder et al., 1989). To increase mule deer populations, harvesting treatments that emulate key structural aspects of critical winter habitat were devised (Dawson et al., 2007). These treatments require partial harvesting of valuable Douglas-fir and, therefore, have the potential to satisfy both objectives of wildlife habitat conservation and timber extraction. However, since partial harvesting requires removal of subcanopy trees and to a lesser degree selected canopy trees (Dawson et al., 2007;Koot et al., 2015), it has the potential to influence stand susceptibility to disturbances in the short term. Reduced tree density may lower the abundance of host trees to insects (Brookes et al., 1987;Negrón, 1998;Aukema et al., 2016) and altered vertical complexity may affect their foraging success (Wulf and Carlson, 1985;Brookes et al., 1987). Similarly, removal of ladder fuels and creation of canopy gaps may lower susceptibility to high-severity crown fires (Schoennagel et al., 2004;Agee and Skinner, 2005;Falk et al., 2007;Vaillant et al., 2015). The potential for partial harvesting to affect longer-term susceptibility to disturbances is unknown.
In this study, we quantified the impacts of partial harvesting to enhance mule deer winter habitat on forest structural attributes and susceptibility to three disturbance agents: the western spruce budworm (Choristoneura occidentalis Freeman), Douglasfir beetle (Dendroctonus pseudotsugae Hopkins), and wildfire. We evaluated the hypotheses that, in the short-term (e.g., 1 year after harvesting), the removal of subcanopy trees will reduce host availability and suitability to the western spruce budworm and Douglas-fir beetle and decrease canopy fuels and lower crown fire likelihood. In the long-term (e.g., 30 years after harvesting), recovery of understory trees within gaps will increase host availability and suitability to the western spruce budworm and increase canopy fuels and crown fire susceptibility; whereas, retention of large old trees will result in abundant hosts for the Douglas-fir beetle. Based on our results, we also present conceptual models depicting disturbance interactions in dry forests to assess trade-offs and synergies that may result from forest management to enhance mule deer winter habitat.

The Study System
Dry forests dominated by Douglas-fir in the interior of BC, Canada are the product of interactions among three primary disturbance agents: the western spruce budworm, Douglasfir beetle and fire. The western spruce budworm is a native defoliator that feeds primarily on Douglas-fir causing reduced radial growth, partial crown mortality, or tree mortality (Brookes et al., 1987;Maclauchlan et al., 2018). Episodic outbreaks occur at 30-year intervals and last from 10 to 14 years (Swetnam and Lynch, 1993), but the interval and duration vary, influencing disturbance severity (Campbell et al., 2006;Axelson et al., 2015;Maclauchlan et al., 2018). Forest susceptibility is directly related to the abundance of available hosts and is also influenced by the vertical complexity of the stand (Wulf and Carlson, 1985;Brookes et al., 1987). Dense, multilayered stands with a high proportion of host trees are very susceptible to outbreaks, especially under moisture stress (Brookes et al., 1987;Campbell et al., 2006;Flower et al., 2014a), although extreme drought can negatively impact western spruce budworm survival (Campbell, 1993).
The Douglas-fir beetle is also a native disturbance agent that feeds, reproduces, and spends most of its life cycle under the bark of Douglas-fir trees (Rudinsky, 1962;Aukema et al., 2016). Normally, low-density populations attack defensively compromised hosts, windthrown trees, large branches, and stumps. When conditions for beetle survival are good and there is an abundance of compromised hosts, as in drought, beetles can rapidly increase in number and switch to attacking large, well-defended hosts (Rudinsky, 1962;Christiansen et al., 1987). Forest susceptibility increases with the basal area and density of Douglas-fir trees (Negrón, 1998), as well as the distribution and abundance of defensively impaired hosts (Dodds et al., 2006).
In addition to biotic disturbance agents, mixed-severity fire regimes strongly influence forest dynamics in dry interior forests. Historically, frequent low-severity surface fires burned every < 10 to 40 years (Daniels, 2004;Heyerdahl et al., 2012;Harvey et al., 2017), while infrequent higher-severity fires yielded stands co-dominated by lodgepole pine (Daniels, 2004). Similar to other forests in western North America, the fire regime was disrupted in the late 19th century (Hessburg et al., 2019). Reduced fire frequency resulted in persistent dense subcanopies and accumulation of hazardous woody surface fuels within stands (Daniels, 2004), and forest homogenization at landscape scales (Perry et al., 2011;Hessburg et al., 2016Hessburg et al., , 2019. Contemporary forests are susceptible to high-severity crown fires, especially during drought (Stephens et al., 2013;Flannigan et al., 2016;Wotton et al., 2017).
Mule deer utilize multi-layered old Douglas-fir forests with moderate to high crown cover as winter habitat (Armleder and Dawson, 1992;Armleder et al., 1994). There, they feed on Douglas-fir foliage and arboreal lichen that break in the wind or due to snow loads and drop to the ground. Survival of mule deer in winter depends on this habitat because it provides snow interception, thermal cover and high-quality food (Armleder and Dawson, 1992). Selective harvesting intended to enhance mule deer habitat involves the proportional removal of Douglas-fir trees based on their diameter, with a focus on removing smalldiameter trees to create gaps, multi-layered, uneven-aged stands with a clumpy tree distribution to be sustained over the long-term (Armleder and Dawson, 1992;Dawson et al., 2007).

Study Area
Research was conducted in the Knife Creek portion of the Alex Fraser Research Forest in central BC, Canada (Figure 1). The area comprises closed-canopy forests with well-developed canopy, subcanopy, and regeneration strata dominated by interior Douglas-fir, with infrequent trembling aspen, Populus tremuloides Michx., and paper birch Betula papyrifera Marsh (Koot et al., 2015). Although present in the past, lodgepole pine is largely absent from the canopy due to a mountain pine beetle (Dendroctonus ponderosae) epidemic and salvage harvesting between 1998 and 2004 (Day, 2007). Other recent biotic disturbances include defoliation by western spruce budworm, tree mortality due to Douglas-fir beetle, and subsequent salvage harvesting of affected individuals. Modern fire records include only one 6-hectare (ha) fire at Knife Creek in September 2013 (BC Wildfire Service, unpublished data).

Mule Deer Winter Habitat Experiment
In 1983, the study area was dedicated to research on mule deer winter habitat (Armleder and Thomson, 1984). A controlled experiment was established with treatment (86 ha) and adjacent control (75 ha) areas. Similar harvesting treatments intended to improve mule deer winter habitat were applied in 1984 and 2014. In 1984, 16% of volume was removed including merchantable trees (diameter at breast height, DBH ≥ 12.5 cm) in proportion to their abundance by diameter class. Trees were harvested in clusters, leaving trees with DBH ≥ 12.5 cm, to create a multilayered stand and promote thermal and security cover (Armleder and Thomson, 1984;Armleder and Dawson, 1992). In 2014, trees from all diameter classes were harvested, with emphasis on removal of smaller diameter trees (Koot et al., 2015). Twelve to 22% of volume per treatment area was harvested to maintain a multilayered structure, create gaps to increase mule deer forage, and address declining tree vigor (Koot et al., 2015).

Research Design and Sampling
Permanent sample plots were established at 100 m intervals along 15 systematically placed transects creating a grid across the study area (Armleder and Thomson, 1984;Koot et al., 2015). Sampling in 1984 quantified stand densities using variable radius plots (basal area factor 8) of trees with DBH ≥ 12.5 cm, focussing on the treatment areas in advance of harvesting. We resampled every second plot on the grid, including 34 treatment and 33 control plots (Figure 1; Leclerc, 2017). Treatment plots were sampled in 2014 and 2015, before and after the second experimental harvest in fall 2014; control plots were sampled in 2015.
We sampled the 10 canopy and 10 sub-canopy trees closest to each plot center (Jonsson et al., 1992). In the treatment plots, different trees were measured in 2014 and 2015 due to the intervening harvesting. Tree species, state (dead or alive), DBH, height, canopy class and the distance to plot center were recorded. In every second plot, one increment core was sampled from the base of each measured tree (≤30 cm from the ground) to estimate age and assess growth history. As needed, we took multiple cores to intercept the pith or ensure the sample was within 5 years of the pith based on ring geometry.
Other biophysical attributes measured at the center of the permanent plots were location (UTM coordinates), elevation (m.a.s.l.), slope aspect and angle in degrees. To quantify canopy cover, photographs were taken c.1 m above the ground using a 35 mm digital camera fitted with a fish-eye lens with a 180 o field of view. We used the line-intercept method to quantify dead wood (hereafter "woody surface fuels") along a randomly oriented 30m transect that bisected plot center (van Wagner, 1982). For all pieces of wood (diameter ≥ 2.5 cm), we recorded location along the transect, species, and its origin either as natural or from harvesting in 2014 if it had freshly cut log ends and the bark was intact.

Forest Attributes
For each plot, the density and basal area of canopy and sub-canopy trees were calculated using a negative binomial distribution to account for the clustered spatial pattern (Lessard et al., 2002). Canopy cover was derived from the hemispherical photographs using the program Gap Analyzer Version 2 (Frazer et al., 1999). Canopy openness, the percentage of open sky after accounting for topographic shading, was subtracted from 100 to estimate canopy cover (Frazer et al., 1999). Woody fuel loads were calculated by converting log diameters to mass (kg m −2 ) (van Wagner, 1982): where, D is the wood density of 440 kg m −3 for Douglas-fir (Gonzalez, 1990), L is the length of the transect (30 m) and dia is the diameter of an individual log in cm. The mass of all logs along each transect were summed to give the plot-level large woody fuel load (kg m −2 ). Increment cores were mounted on wooden supports to dry and sanded with progressively finer-grit sand paper to reveal ring structure (Stokes and Smiley, 1968). Ring-widths were measured on digital images of each core (2400 dpi) using the program Coorecorder (Larsson, 2011a) and crossdated against an existing master chronology for Douglas-fir (Daniels, 2004) using the program CDendro (Larsson, 2011b). For cores that did not intercept the pith, the number of missed rings (7.7 ± 9.8; mean ± standard deviation) was estimated from geometric measurements of the inner-most rings (Duncan, 1989). Tree age was the difference between the calendar years of the outer-most ring and the pith, plus one, plus a correction for the number of years for trees to grow to coring height (Daniels, 2004;Daniels et al., 2017). For trees that were not cored, we developed a regression model (n = 337, R 2 = 0.308, standard error of the estimate = 1.3 years, p < 0.0001), as follows: log 10 (Age) = 1.2913 + 0.5723 * log 10 DBH + 0.5819 * crownclass − 0.3884 * log 10 DBH * crownclass where, log 10 (Age) is the logarithm of age in years; log 10 DBH is the logarithm of the DBH in cm; crownclass is a categorical variable distinguishing between canopy and sub-canopy trees; and, log 10 DBH * crownclass is the interaction between log 10 DBH and crownclass. Predicted ages were back-transformed to estimate tree age in years.

Forest Susceptibility to Disturbance
Plot-level forest susceptibility to disturbance was evaluated using empirically derived, published models that are commonly applied in western North America. For western spruce budworm, we used the index developed by Wulf and Carlson (1985). This index combines six stand structure, two climatic, and one landscape variables (Supplementary Material). The six stand structure variables for each plot were the coefficient of variation of tree DBH (CV DBH ), canopy cover (%) of all trees, host Douglas-firs only, and late-successional host trees only, as well as measures of tree vigor, and maturity. Hemispheric photographs were used to estimate total canopy cover and the cover of late-successional host trees since all plots were late-successional forest comprising 95-100% Douglas-fir trees (Leclerc, 2017). Vigor was the plot-level basal area as a percentage of the mean plot-level basal area of the upper quartile of fully-stocked control plots in our study area (50.2 m 2 ha −1 ). Maturity was the mean of the basal-area-weighted age of the Douglas-fir trees in each plot, with tree ages derived from the increment cores or predicted from the study-area specific regression model (Leclerc, 2017). All plots were classified as a warm dry Douglas-fir habitat in a dry regional climate, with the highest level of host continuity since > 76% of forests in the 400 ha surrounding the study area were dominated by Douglas-fir (Wulf and Carlson, 1985). The absolute value of each variable was converted to a weighted index value and the indices were multiplied yielding a plot-level susceptibility index from 0 to 100 (Wulf and Carlson, 1985; Supplementary Material). Plot-level index values were grouped into low (0-20), moderate (21-50), and high (51-100) susceptibility classes (Wulf and Carlson, 1985).
Susceptibility to fire was expressed as the likelihood of crown fire occurrence, calculated using the crown fire initiation and spread (CFIS) software (Cruz et al., 2004;Alexander et al., 2006). To represent conditions when wildfire risk is greatest, we simulated a wildfire burning from 13:00-15:00 h in August, using 90th percentile fire weather conditions for the Knife Creek station 225 (2008-2105, BC Wildfire Service, unpublished data). Input variables include temperature (28.4 • C), relative humidity (22%), windspeed (11.7 km hour −1 ), the initial spread (13.0) and build up (162.3) indices (van Wagner, 1987), and potential surface fuel consumption (up to 3.3 km m −2 ) (Forestry Canada Fire Danger Group, 1992). Plot-level tree densities, basal area, and mean tree height were used to derive canopy base height (m) and canopy bulk density (km m −3 ; Cruz et al., 2003), the two plot-level fuel input parameters. CFIS outputs included probability of crown fire occurrence for each plot and type of fire (e.g., surface, intermittent crown, or active crown; Alexander et al., 2006).

Changes in Forest Structure and Disturbance Susceptibility
To assess changes in forest structure due to harvesting treatments, we used a randomized complete block (e.g., control versus treatment) single-factor mixed-effects model with subsampling (e.g., plots within treatments; Leclerc, 2017). In a first set of tests, the harvesting treatment was the single factor tested and consisted of 3 levels: treatment area in 2014, treatment area in 2015, and control area in 2015. An analysis of variance (ANOVA) was conducted on each of 8 forest structure attributes, as well as western spruce budworm susceptibility, potential mortality due to Douglas-fir beetle, two canopy fuel attributes, and the probability of crown fire occurrence. In a separate test, we compared canopy and subcanopy growth after the first harvesting treatment in 1984. For each tree, the crossdated ring widths were scaled so that their cumulative values were equal to tree diameter inside the bark, which accounted for asymmetrical growth around the stem. The scaled ring widths were converted to basal area increments (BAI) using the outside-in method (Biondi and Qeadan, 2008). We compared mean BAI 30 years prior to and after the harvesting treatment in 1984. In this test, the harvesting treatment was the single factor and consisted of 4 levels: treatment and control areas for the 30 years prior to  and after  harvesting. To meet the assumptions of ANOVA, attributes were transformed using log base 10, square root and BLOM transformations (Beasley et al., 2009). When significant differences were detected (α = 0.05), pair-wise analyses of least squares means using a Bonferroni correction differentiated among treatment levels. All analyses were conducted in R 3.2.2 (R Core Team, 2015) using Type I sum of squares, and residual maximum likelihood.

Variation in Forest Structure
Five of eight stand structural attributes differed significantly before and after harvesting or between the control and treatment plots ( Table 1). The 2014 harvesting treatment significantly decreased the density and basal area of subcanopy trees and canopy cover, but the preharvest 2014 values did not differ from values in the control area. In contrast to canopy attributes, woody surface fuels increased significantly in treated plots between 2014 and 2015, but the 2015 values did not differ from values in the control area. Mean tree heights were shorter in the treatment plots in 2014 than in 2015, reflecting the harvest of subcanopy trees and subsequent taller mean tree heights. Although the density of canopy trees was greatest in the treatment plots in 2014 and decreased in 2015, pair-wise comparisons did not detect differences due to high variance among treatment levels. Similarly, canopy tree basal area and mean DBH did not differ across treatment levels.
Retrospective analysis of BAI showed subcanopy trees in the treated area grew significantly faster from 1985-2014 than the preceding 30 years but subcanopy tree growth did not change in the control area, suggesting that surviving subcanopy trees benefited from the 1984 harvesting treatment ( Table 2). The mean BAI of canopy trees did not differ between the control or treatment areas or after harvesting in 1984.

Susceptibility to Disturbance
Western spruce budworm susceptibility indices averaged 40.1 to 44.0 in the treatment and control areas, respectively, and did not differ significantly among any of the treatment levels ( Table 1).
Harvesting in 2014 decreased the average susceptibility indices by −3.9 on the relative scale from 0 to 100, but changes in individual treatment plots ranged from −23 to +20. Prior to treatment in 2014, 25 and 8 plots were ranked as moderate and high susceptibility, respectively (Figure 2). The 2014 harvesting treatment reduced susceptibility of 3 plots from high to medium but the single plot ranked low increased to medium. In the control area, 2 plots were ranked low, while 17 and 14 ranked moderate and high, respectively. Douglas-fir beetle susceptibility ranking was moderate or high in all plots since the relative basal area of Douglas-fir exceeded The BLOM transformation was applied to meet assumptions. For subcanopy trees, means (standard errors) followed by the same superscripts did not differ significantly (α = 0.05).
the 89.4% threshold (Negrón, 1998), regardless of treatment (Figure 2). Prior to harvesting in 2014, only six plots had a moderate probability of infestation across the treatment and control areas, while all other plots ranked as high. Harvesting reduced the density of Douglas-fir host trees below the 292 trees ha −1 threshold (Negrón, 1998) in 15 treatment plots, reducing their probability of infestation from high to moderate. Despite the reduction in tree density, the mean decrease in potential mortality from Douglas-fir beetle was −2.7%, although individual plots changed by −24 to +0.3%. Overall, mean potential mortality did not differ significantly between the treatment plot in 2014 or 2015 or relative to the control area ( Table 1).
For the simulated 90th percentiles of fire weather, the likelihood of crown fire occurrence was highest in the treatment plots in 2014 (61.6%), but decreased significantly following harvesting in 2014 (25.5%) to levels that were not different from the control area (23.1%, Table 1). The mean decrease in crown fire likelihood among treatment plots was −36% and ranged from 0 to −60%. These decreases in crown fire likelihood paralleled increases in canopy base height and decreases in canopy bulk density due to harvesting. Mean canopy base height was lowest in the treatment plots in 2014 (7.6 m), but increased significantly (10.5 m) and was no longer different from the mean of the control area (10.6 m). Canopy bulk density in treatment plots (0.015 kg m −3 ) decreased significantly following harvesting (0.011 kg m −3 ) to levels lower than in the control plots (0.017 kg m −3 ). Prior to treatment in 2014, 6 and 18 treatment plots had potential to support active or passive crown fire, respectively (Figure 2). Following harvesting, potential fire behavior was reduced from active to passive crown fire in 2 plots, from active or passive crown fire to surface fire in 22 plots. After harvesting in 2015, 32 of the 34 treatment plots and 32 of the 33 control areas were likely to sustain surface fire; the other three plots had the potential sustain a passive crown fire only.

DISCUSSION
Implementation of a harvesting treatment intended to enhance mule deer winter habitat, while facilitating timber extraction in a dry interior Douglas-fir forest, altered several key forest attributes over the short term, but most recovered within 30 years. The influences of the treatment on biotic and abiotic disturbances were complex, but did not necessarily preclude the potential for management to emulate the historical range and variability of the disturbance regime (e.g., Keane et al., 2009).

Mule Deer Winter Habitat Management and Key Forest Structures
During harvesting, large, old canopy Douglas-firs were retained while canopy gaps were created and subcanopy trees were targeted for removal (Armleder and Dawson, 1992;Dawson et al., 2007;Koot et al., 2015). In the short term following treatment, total canopy cover, as well as subcanopy tree density and basal area decreased, while woody surface fuels increased, a common impact of treatments that include thinning subcanopy trees (Agee and Skinner, 2005;Stephens and Moghaddas, 2005;Vaillant et al., 2015). Thirty years after the first harvesting treatment, few forest attributes differed significantly between the treatment and control areas, suggesting forest recovery. This inference is corroborated by the significant increase in basal area increments of subcanopy trees in the treatment plots, but not the control plots, after 1984. Evidently, the retention of large Douglas-fir trees, removal of < 20% of volume across all diameter classes, and harvesting at intervals of c. 30 years enhances mule deer winter habitat (Dawson et al., 2007), maintains forest structures common in the landscape (Maclauchlan and Brooks, 2009), and allows enough time for the forest to recover to sustain timber extraction and offset costs of active habitat management.

Susceptibility to Disturbance
Harvesting treatments to enhance mule deer winter habitat had mixed effects on forest susceptibility to insects and fire, both in the short (1 year) and long (30 years) terms. Contrary to our hypothesis, quantitative measures of forest susceptibility to western spruce budworm were only modestly affected in the year after treatment. The characteristics that made the Knife Creek study area suitable for mule deer winter habitat (a uniform, mature Douglas-fir dominated canopy) ensured that most stand composition and structure variables that influence susceptibility to western spruce budworm (Wulf and Carlson, 1985) were unlikely to be significantly affected in the short term by the harvesting treatment. Plots had relatively consistent cover of old Douglas-fir trees and a range of tree diameters; however, these attributes were only modestly impacted by the harvesting. Further, due to the small size of the study area, the regional climatic and landscape-scale variables known to affect susceptibility to western spruce budworm (Wulf and Carlson, 1985) were constant among plots. Over the long term, however, the harvesting treatments emulated the stand structures and frequencies expected to result from western spruce budworm outbreaks. Proportionally removing trees based on their abundances resulted in greater removal of smaller diameter trees, creating stand structures similar to defoliated stands since smalldiameter trees die at higher rates during an outbreak (Alfaro et al., 1982;Alfaro and Maclauchlan, 1992). Additionally, harvesting at 30-year intervals is within with the range of western spruce budworm outbreak frequencies in dry forests (Campbell et al., 2006;Alfaro et al., 2014;Axelson et al., 2015;Maclauchlan et al., 2018).
In support of our hypothesis, forest susceptibility to Douglasfir beetle decreased in the short term following treatment to enhance mule deer winter habitat. Where the density of Douglasfir was reduced below the critical threshold of ≤ 292 trees ha −1 (Negrón, 1998), susceptibility was reduced from high to moderate. Reducing forest density in general may be effective in preventing bark beetle infestations (Fettig et al., 2007(Fettig et al., , 2014, although the effect is enhanced when large host trees are removed (Temperli et al., 2014). During the 30 years following harvesting in 1984, fast growth of subcanopy trees suggests improved stand vigor and fewer stressed trees, which would sustain fewer beetle infestations as Douglas-fir beetle preferentially attack compromised hosts (Furniss, 1965;Negrón, 1998;Aukema et al., 2016). Further, a more vigorous stand would have a higher probability of keeping a resident beetle population in an endemic state, given that greater beetle densities are required to attack and overwhelm healthy trees (Raffa et al., 1993;Boone et al., 2011).
As hypothesized, the likelihood of crown fire occurrence was reduced in the year following treatment, reflecting the decrease in canopy bulk density and increase in canopy base height due to harvesting. Thus, harvesting to enhance mule deer winter habitat was consistent with the goal of mitigating canopy fuels to reduce fire intensity, and rate of spread, and increase forest resistance to crown fire (Stephens et al., 2009;Fulé et al., 2012) in the short term. Thinning small-diameter trees and increasing the spacing between canopy trees both reduced fuel hazards (Agee and Skinner, 2005;Stephens and Moghaddas, 2005;Stephens et al., 2009) and enhanced mule deer habitat (Armleder and Dawson, 1992;Dawson et al., 2007). Over time, the benefits of the treatments in terms of reduced canopy fuels will diminish (Stephens et al., 2009;Vaillant et al., 2015) as tree density and canopy bulk density increase and canopy base height lowers, yielding hazardous fuels. This was consistent with our observations of conditions in the long term following the initial treatment. It appears that the 30-year interval to maintain habitat while sustaining timber yield has a long-term trade-off with crown fire likelihood. Evidently, 30 years may not be the optimal treatment interval to maintain low likelihood of crown fire occurrence.
Countering the reduction in canopy fuels, surface woody fuels increased in the short term following treatment, in part because large logs were intentionally left to reduce impacts of harvesting equipment on the soil. Large woody fuels are known to contribute to spotting and smouldering combustion, which can damage tree cambium and soil (Stephens and Moghaddas, 2005). Furthermore, large surface fuel loads are common after harvesting if flammable stems, tree tops, branches, and leaf litter are not mitigated (Stephens and Moghaddas, 2005;Stephens et al., 2009). Although foliage and fine wood decays relatively quickly, residual large logs decay slowly and elevate the fuel hazard for years to decades (Campbell et al., 2016). Unfortunately, the increased woody fuel load was not reflected in the calculated likelihood of crown fire occurrence since the CFIS model (Cruz et al., 2004;Alexander et al., 2006) does not include surface fuels as an input variable, despite the welldocumented influences of surface fuels on fire spread, behavior, and effects (Hall and Burke, 2006;Collins et al., 2013). Therefore, our estimates of treatment effects on crown fire occurrence may be conservative.

Complex Disturbance Interactions
Each type of disturbance generates legacies in structure, composition, and function that interact with subsequent disturbances (Buma, 2015). Interactions among agents are complex, depend on the magnitude and time since disturbance, and vary across spatial and temporal scales (Jenkins et al., 2008;Hicke et al., 2012). Figure 3 depicts our conceptual model of baseline historical conditions and disturbance interactions in the Douglas-fir-dominated forests at Knife Creek. Individual disturbance agents create self-reinforcing structures that perpetuate their regimes (Brookes et al., 1987;Hessburg et al., 2007), except the Douglas-fir beetle that kills susceptible host trees (Aukema et al., 2016). Frequent surface fires limit woody surface fuel loads, forest density, canopy fuel loads, and crown fire (Schoennagel et al., 2004;Hessburg et al., 2007), but perpetuate multilayered forest canopies optimal for western spruce budworm feeding and inter-tree dispersal (Brookes et al., 1987;Hood and Bentz, 2007). Large, old trees commonly survive surface fire, providing high-quality hosts for Douglas-fir beetles (Negrón, 1998;Aukema et al., 2016). In contrast, periodic crown fires consume surface and canopy fuels, limiting subsequent fire (Schoennagel et al., 2004). Crown fire also kills host trees for insects (McCullough et al., 1998), although scorched trees are preferentially colonized by Douglas-fir beetle immediately following fire (Furniss, 1965;Hood and Bentz, 2007). Western spruce budworm defoliation initially reduces canopy and ladder fuels (Alfaro et al., 1982;Lynch and Moorcroft, 2008;Sturtevant et al., 2012), limiting fire propagation through the canopy (Cohn et al., 2014). Tree mortality caused by western spruce budworm or Douglas-fir beetle generates fine surface fuels, canopy gaps, and multilayered canopies prone to surface fires in the short term (Lynch and Moorcroft, 2008;Hicke et al., 2012;Donato et al., 2013), although impacts on fire occurrence are ambiguous (Flower et al., 2014b). In the mid-term, snags fall so large woody surface fuels accumulate, increasing the intensity of surface fire and potential for crown fire (Hummel and Agee, 2003;Jenkins et al., 2008;Vaillant et al., 2015). Ultimately, a discontinuous canopy cannot sustain crown fire (Jenkins et al., 2008;Donato et al., 2013). Between biotic agents, moderate defoliation by the western spruce budworm stresses and predisposes trees to Douglas-fir beetle (Hadley and Veblen, 1993), while Douglas-fir FIGURE 3 | Conceptual model of historical disturbance interactions in dry Douglas-fir-dominated forests. Solid color arrows and "+" symbols indicate net positive interactions between disturbances. Dashed color arrows and "−" symbols indicate net negative effects. Arrow colors correspond to the causal disturbance agent and point toward the affected disturbance; labels identify key mechanisms of interaction. beetle removes host trees for the western spruce budworm (McCullough et al., 1998).

Harvesting Impacts on Disturbance Interactions
The implementation of harvesting treatments to enhance mule deer winter habitat affected disturbance interactions in a multitude of ways over time (Figure 4) and could potentially unbalance the complex interactions shown in Figure 3. In the short term (e.g., 1 year after harvesting), tree harvesting decreased the likelihood of crown fire and susceptibility to Douglas-fir beetle by temporarily disrupting fuel dynamics and the generation of susceptible host trees (Figures 4A,C). However, a well-documented but unsustainable interaction is tree mortality due to infestation by Douglas-fir beetles following harvesting (Aukema et al., 2016). Management options to counter this effect include preventing infestation (e.g., removing logs and other harvesting residuals; Figures 4A,B; Armleder and Thomson, 1984;Dawson et al., 2007) or managing infestation (e.g., destroying infested trees; Figures 4C,D; Aukema et al., 2016); each option impacts subsequent disturbances. Management by removing harvesting residuals also mitigates the surface fuel hazard, enhancing the short-term shift from crown to surface fire ( Figure 4A). In the long term (e.g., 30 years after harvesting), canopy fuels generally recover, but large surface woody fuels may not. Where harvesting residuals were mitigated, lower intensity surface fires perpetuate complex vertical and horizontal canopy FIGURE 4 | Conceptual models of disturbance interactions following harvesting to enhance mule deer winter habitat in Douglas-fir-dominated forests. Four scenarios (A-D) contrast time since harvesting (1 year or 30 years, columns) and options to mitigate Douglas-fir beetle infestations after harvesting (preventing or managing infestatios, rows). Solid rectangles indicate self-reinforcing disturbances; open rectangles indicate disrupted disturbances. Black arrows indicate management impacts on disturbance: "+" symbols indicate net positive effects, "−" symbols indicate net negative effects, and "o" indicates no net effect or a neutral interaction. Solid color arrows and "+" symbols indicate net positive interactions between disturbances. Dashed color arrows and "−" symbols indicate net negative effects. Arrow colors correspond to the causal disturbance agent and point toward the affected disturbance; labels identify key mechanisms of interaction.
structures that are susceptible to the western spruce budworm and scorch trees making them suitable hosts for the Douglasfir beetle (Figure 4B). The combined effects of sustained low woody fuel loads, surface fires, defoliation by western spruce budworm, and infestations by Douglas-fir beetle perpetuate stands with a more open canopy, inhibiting potential crown fire spread (Jenkins et al., 2008;Donato et al., 2013;Cohn et al., 2014). In contrast, managing infestations but leaving harvesting residuals in situ increases the potential intensity of surface fire and propensity for crown fire, despite the reductions in canopy fuels following harvesting ( Figure 4C). Crown fires exacerbate infestations by Douglas-fir beetle in the short term. In the long term, this system is more likely to be dominated by crown fires (Figure 4D). With more crown fires, susceptibility to western spruce budworm and Douglas-fir beetle lower as available hosts are killed and stand structure simplifies (McCullough et al., 1998).
Harvesting to enhance mule deer winter habitat evidently disrupts historical disturbance interactions and shifts susceptibility to various disturbance agents depending on time since treatment and whether residual logs are removed (Figure 4). Fire is omnipresent. Whether surface or crown fires burn depends on mitigation of harvesting residuals and the subsequent cascade of interactions. Surface fires generate complex forest structures and interactions with insects, while crown fires simplify forests and exclude insects. Further, the likely occurrence of surface versus crown fire depends on time since treatment and whether harvesting residuals are removed. Of the four scenarios based on current mule deer winter habitat management (Figure 4), the combination of harvesting with mitigation of harvesting residuals was closest to our conceptual model of historical disturbance regime interactions (Figure 3). Given the potential for prolonged fire seasons creating conditions more conducive to lightning ignitions and intense fires (Flannigan et al., 2016;Wotton et al., 2017), we recommend managing mule deer habitat to avoid future crown fires by mitigating surface woody fuels.

Management Recommendations
Failure to consider the complexity of ecosystems in their management is a frequent cause of failures in sustainable resource management efforts (Ascher, 2001;Allen and Gunderson, 2011). Single-objective management programs reduce complexity by placing greater emphasis on a particular value while failing to consider potential consequences on other values. Traditionally, extensive management for timber production was prioritized in the dry forests of British Columbia by using clearcut harvesting, regenerating dense forests, and suppressing wildfires. Cumulatively, these strategies degraded mule deer habitat (Armleder et al., 1989), requiring regulation. Designed to counter these negative impacts, the harvesting treatment to enhance mule deer habitat emphasizes partial harvesting to remove subcanopy trees and create canopy gaps, while sustaining timber extraction during periodic maintenance treatments (Dawson et al., 2007). Based on our assessments 1 and 30 years after harvesting, we conclude that timely removal of logging residuals immediately following partial harvesting is necessary to simultaneously address short-term susceptibility to Douglas-fir beetle and mitigate against intense surface fire. Although forest recovery after 30 years generated structures and disturbance interactions that were consistent with our best understanding of historical regimes, forest susceptibility to insects and crown fire were relatively high. Thus, in dry Douglas-fir forests, we recommend maintenance treatments at intervals ≤ 30 years, accompanied by immediate removal of harvesting residuals, to maintain both timber production and critical mule deer winter habitat and achieve sustainable forest management.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
M-AL organized the data collection and analyzed the data with guidance from LD and AC. M-AL led the writing of the manuscript. All authors contributed to discussing and revising the manuscript.