Long-Term E ﬀ ects of Fuels Treatments, Overstory Structure, and Wildﬁre on Tree Regeneration in Dry Forests of Central Washington

: The long-term e ﬀ ectiveness of dry-forest fuels treatments (restoration thinning and prescribed burning) depends, in part, on the pace at which trees regenerate and recruit into the overstory. Knowledge of the factors that shape post-treatment regeneration and growth is limited by the short timeframes and simple disturbance histories of past research. Here, we present results of a 15-year fuels-reduction experiment in central Washington, including responses to planned and unplanned disturbances. We explore the changing patterns of Douglas-ﬁr regeneration in 72 permanent plots (0.1 ha) varying in overstory abundance (a function of density and basal area) and disturbance history—the latter including thinning, prescribed burning, and / or wildﬁre. Plots were measured before treatment (2000 / 2001), soon afterwards (2004 / 2005), and more than a decade later (2015). Thinning combined with burning enhanced sapling recruitment (ingrowth) into the overstory, although rates of ingrowth were consistently low and greatly exceeded by mortality. Relationships between seedling frequency (proportion of quadrats within a plot) and overstory abundance shifted from weakly negative before treatment to positive after thinning, to neutral in the longer term. However, these relationships were overshadowed by more recent, higher-severity prescribed ﬁre and wildﬁre that stimulated seedling establishment while killing advanced regeneration and overstory trees. Our results highlight the dependence of regeneration responses on the history of, and time since, fuels treatment and subsequent disturbance. Managers must be aware of this spatial and temporal complexity and plan for future disturbances that are inevitable but unpredictable in timing and severity.


Introduction
Large areas of dry coniferous forests in western North America require restoration to reduce the risk of high severity wildfire and to restore the structures and functions that characterized these forests historically [1][2][3][4][5]. Restoration treatments, including mechanical thinning and prescribed burning, have proven effective for achieving these objectives in the short term by reducing surface and ladder fuels, increasing canopy base height, and reducing canopy bulk density [6][7][8][9][10][11][12]. However, the longevity of these effects varies considerably within and among landscapes [13], and the factors influencing this variability are not well understood.
A key factor contributing to fuels-treatment longevity is the rate at which trees re-establish and grow following treatment. Understory regeneration is a primary source of ladder fuels, and as saplings Treatments were randomly assigned among units with some re-allocation [42]. 'Untreated' refers to units that were neither thinned nor prescribed burned. The stippled polygon is the area affected by the 2012 Wenatchee Complex wildfire. Each unit contains six sample plots.
The FFS experimental design included factorial combinations of thinning and prescribed-fire treatments (Figure 1). Treatments were randomly assigned to units, although the endangered northern spotted owl (Strix occidentalis caurina) and fire-control issues required some re-allocation [42]. We summarize treatments here; for details see [42]. Six units were thinned in fall 2002 and spring 2003. Units were thinned from below in a non-uniform pattern to a target basal area of 10-14 m 2 /ha such that thinning intensity varied with initial forest structure [44]. Targets reflected the historical structure, including clumped patterning of trees [42,45] and a goal of reducing fire hazard such that at least 80% of the basal area of dominant and co-dominant trees would survive if subjected to a head fire under 80 th -percentile fire-weather conditions [18]. On average, density was reduced by 60% and basal area by 50%, leaving a mean residual basal area of 17.4 m 2 /ha-somewhat larger than the target [8]. Trees were yarded by helicopter. Prescribed-fire treatments had a more complex history. Although planned in six units, only four burns were conducted in spring 2004. Fire severity was low (crown scorch on 30% of trees) and coverage was patchy (23%-51% of the forest floor) such that treatments failed to meet most surface-fuel objectives [8,46]. The two remaining burn units were Treatments were randomly assigned among units with some re-allocation [42]. 'Untreated' refers to units that were neither thinned nor prescribed burned. The stippled polygon is the area affected by the 2012 Wenatchee Complex wildfire. Each unit contains six sample plots.
Twelve experimental units were established in 2000, each~10 ha in area. A majority represented the Douglas-fir series [43]. Within each unit, six permanent plots (0.1 ha; 20 × 50 m) were subjectively placed in continuous forest vegetation to represent the dominant plant associations [42]. Plot elevations ranged from 670 to 1180 m (mean: 852 m) and slopes from 7% to 72% (mean: 41%). Prior to treatment, plots varied in relative abundance (density and basal area) of the principal tree species-from strong dominance by Douglas-fir to strong dominance by ponderosa pine [8].
The FFS experimental design included factorial combinations of thinning and prescribed-fire treatments (Figure 1). Treatments were randomly assigned to units, although the endangered northern spotted owl (Strix occidentalis caurina) and fire-control issues required some re-allocation [42]. We summarize treatments here; for details see [42]. Six units were thinned in fall 2002 and spring 2003. Units were thinned from below in a non-uniform pattern to a target basal area of 10-14 m 2 /ha such that thinning intensity varied with initial forest structure [44]. Targets reflected the historical structure, including clumped patterning of trees [42,45] and a goal of reducing fire hazard such that at least 80% of the basal area of dominant and co-dominant trees would survive if subjected to a head fire under 80th-percentile fire-weather conditions [18]. On average, density was reduced by 60% and basal area by 50%, leaving a mean residual basal area of 17.4 m 2 /ha-somewhat larger than the target [8]. Trees were yarded by helicopter. Prescribed-fire treatments had a more complex history. Although planned in six units, only four burns were conducted in spring 2004. Fire severity was low (crown scorch on 30% of trees) and coverage was patchy (23%-51% of the forest floor) such that treatments failed to meet most surface-fuel objectives [8,46]. The two remaining burn units were treated in spring 2006, when conditions were drier, severity was higher (crown scorch on 80% of trees), and coverage was more complete [8].
Following treatments, the Wenatchee Complex wildfire burned through four units in fall 2012. The units burned by the wildfire were thinned, previously burned, and untreated (Figures 1 and 2). Anecdotal observations indicated that the wildfire burned at relatively low intensity due to an inversion layer; however, it was more severe than either prescribed fire.
Forests 2020, 11, x FOR PEER REVIEW 4 of 29 treated in spring 2006, when conditions were drier, severity was higher (crown scorch on 80% of trees), and coverage was more complete [8].
Following treatments, the Wenatchee Complex wildfire burned through four units in fall 2012. The units burned by the wildfire were thinned, previously burned, and untreated (Figures 1 and 2). Anecdotal observations indicated that the wildfire burned at relatively low intensity due to an inversion layer; however, it was more severe than either prescribed fire.   Figure 2).
The trees within each plot were tagged and identified to species. At each sampling date, trees were coded as live or dead and DBH was measured on live trees. In prescribed burn or wildfire plots, presence/absence of bole char was recorded for each tree.
For each plot × sampling date, we computed tree density (number/ha), basal area (m 2 /ha), and a simple measure of overstory abundance (OA) that combined these two structural metrics ( √ (density × basal area)); Figure A1, Appendix A1), as in [44]. OA served as our primary measure of overstory structure. Because plots were thinned to a common target basal area, pre-treatment OA also serves as a proxy for thinning intensity in thinned plots. For each plot, we also calculated three metrics of overstory change during the Early-Late interval: (1) ingrowth (stems reaching 7.6 cm DBH between 2004 and 2015; number/ha); (2) mortality (trees alive in 2004 but dead in 2015; number/ha); and (3) basal area increment (BAI, difference in basal area between 2004 and 2015; m 2 /ha), capturing the net effect of ingrowth, mortality, and diameter growth of surviving trees ( Table 1). Causes of mortality were not assessed. We did not compute these metrics for the Pre-Early interval as structural changes during this interval were reported in [8].

Regeneration
Regeneration was sampled before treatment (2000 or 2001, Pre ), during the second growing season after the first prescribed burn (2005, 'Early'; third season after thinning), and in 2015 ('Late ) (Figure 2). At each sampling date, the presence or absence of conifer seedlings (<1 m tall) was recorded in each of 20 permanent quadrats (1 × 1 m) using a stratified, random design to ensure a broad distribution of sampling across each plot (see Figure A3 in [44]).
Presence/absence data were used to compute three measures of seedling response for each plot: (1) frequency at each sampling date (proportion of quadrats with seedlings), (2) colonization (number of quadrats without seedlings at the start but with seedlings at the end of the interval), and (3) persistence (number of quadrats with seedlings at both the start and end of the interval, reflecting either survival or turnover). Colonization and persistence were calculated for each interval (Pre-Early, Early-Late, Pre-Late).
In 2015, we also tallied regeneration by height class to include taller stems that were more likely to influence near-term fire behavior and forest structure. We considered three height classes: small seedlings (<0.3 m tall), large seedlings (0.3-1.37 m tall), and saplings (>1.37 m tall but <7.6 cm DBH). Small seedlings were counted in each of the 20 quadrats per plot and large seedlings and saplings in each of 40, 5 × 5 m contiguous grid cells per plot. We did not distinguish among species in our analyses because Douglas-fir accounted for most of the regeneration (~90% of stems).

Statistical Analyses
Although fuels treatments were applied to experimental units, we treated plots within units as independent samples (n = 72) for several reasons. First, a primary goal was to explore how regeneration responded to overstory structure, which varied widely at the spatial scale of plots. Plots more appropriately capture the scale at which overstory and regenerating trees interact (0.1 ha vs.~10 ha for units). Second, plots within the same unit showed wide variation in pre-treatment vegetation structure and physical environment ( Figure A2, Appendix A2). Third, not all plots within units received the intended treatments: three plots in prescribed-fire units fell outside burn perimeters and were included in analyses as unburned. We initially considered a mixed-model approach (i.e., plots nested within units, with unit designated as a random effect) but rejected it because lack of replication at the unit scale precluded testing of all disturbance histories.  We used three binary variables to represent the history of thinning and burning of each plot: 'Thin treatment (thinned or unthinned), Early fire (burned in the 2004 prescribed fires or unburned), and 'Late fire' (burned or unburned after 2004). This approach facilitated data visualization, analysis, and interpretation, despite several limitations. First, binary treatment of thinning does not capture the variability in intensity of tree removal among plots [44] ( Figure A3, Appendix A3). However, by modeling early responses as a function of pre-treatment OA, we accounted for some of this variation: thinned plots with greater pre-treatment OA experienced greater thinning intensity. Second, binary treatment of fire does not account for variation in burn severity. Nevertheless, it captures the key difference between burned (particularly late fire) and unburned plots, as illustrated by the percentage of trees charred within a plot ( Figure A4, Appendix A3). Finally, the Late fire variable does not distinguish between the 2006 prescribed fire and the 2012 wildfire, although both were more severe than the 2004 prescribed fire. As a result, Early and Late fires reflect differences in severity, as well as timing of disturbance.
We analyzed nine response variables, representing overstory dynamics, the spatial frequency of regeneration, and seedling dynamics ( Table 1). Some responses represented individual sampling dates (Pre, Early, or Late) and others represented changes over a sampling interval (Pre-Early, Early-Late, or Pre-Late). Densities of seedlings by height class were highly skewed and did not conform to model assumptions, thus we analyzed frequencies (i.e., number of quadrats or grid cells with individuals of a height class). BAI (which varied from negative to positive) was analyzed using a linear model (LM); all other responses had zero as a lower bound and were analyzed using generalized linear models (GLMs) with a negative binomial distribution.
Models represented the overstory structure (OA) of each plot, its disturbance history, and interactions between OA and disturbance history. Although pre-treatment values were highly predictive of post-treatment responses in previous analyses of understory vegetation [44], lack of pre-treatment data for several regeneration variables precluded their inclusion in these models. We added terms sequentially to each model, beginning with OA at the time of the response or at the beginning of the interval. Disturbance-related terms and associated second-order interactions were then added in chronological order. For example, responses at the Late sampling date were evaluated using the following model structure: Response~OA + Thin + Early fire + Late fire + OA × Thin + OA × Early fire + OA × Late fire + Thin × Early fire + Thin × Late fire + Early fire × Late fire (1) For responses at the Early sampling date or Pre-Early interval, we omitted terms that included Late fire , thus testing effects of disturbance until 2004. For the model of pre-treatment seedling frequency (Table 1), OA was the only model term. The significance of each term was assessed using an F-test for the LM and a chi-squared test for GLMs (i.e., Analysis of Deviance) with α = 0.05. All analyses were conducted using R (v. 4.0.0) [47] with the MASS package (v. 7.3-51.5) [48]. For each model, we present the goodness-of-fit statistic, R 2 . For the linear model, R 2 is the percentage of variance explained, and is calculated as 1 minus the ratio of the residual sum of squares to the total sum of squares. For the GLMs, R 2 is the model's discrepancy from a perfect representation of the data and is calculated as 1 minus the ratio of the deviance (fitted vs. 'perfect' model) to the null deviance ('intercept-only' vs. perfect model) [49]. We also present the percentage of each model's total explained variance or deviance attributable to each term; these percentages sum to 100% for each model. For GLMs, this percentage is the ratio of the improvement in deviance with each added term to the deviance of the full, fitted model [49]. The analysis script, data, and metadata are available online in Supplementary Materials.

Overstory Dynamics: Early-Late Interval (2004-2015)
Post-treatment ingrowth (recruitment of saplings into the overstory) was recorded in 33% of plots (24 of 72), at a low average density (10 trees/ha; range: 0-70; Figure 3). Ingrowth density was not related to OA but was related to thinning and burning history. Ingrowth density was higher in thinned than in unthinned plots. Early burning (2004) enhanced this effect (significant Thin × Early-fire interaction; Table 2), but the effect of Late fire was not significant.
analyses were conducted using R (v. 4.0.0) [47] with the MASS package (v. 7.3-51.5) [48]. For each model, we present the goodness-of-fit statistic, R 2 . For the linear model, R 2 is the percentage of variance explained, and is calculated as 1 minus the ratio of the residual sum of squares to the total sum of squares. For the GLMs, R 2 is the model's discrepancy from a perfect representation of the data and is calculated as 1 minus the ratio of the deviance (fitted vs. 'perfect' model) to the null deviance ('intercept-only' vs. ′perfect′ model) [49]. We also present the percentage of each model's total explained variance or deviance attributable to each term; these percentages sum to 100% for each model. For GLMs, this percentage is the ratio of the improvement in deviance with each added term to the deviance of the full, fitted model [49]. The analysis script, data, and metadata are available online in Supplementary Materials.  Table 1). Points represent plots with different overstory abundance (OA) values in the Early sample and different histories of thinning and fire (columns and symbols). Columns distinguish thinning and Early (2004) prescribed-fire treatments. Symbols represent the most recent fire (′p′ for prescribed fire and ′w′ for wildfire). Dashed and solid lines are modeled relationships for the Late-fire treatment based on statistically significant terms ( Table 2). For ingrowth, OA and the Late-fire treatment were not significant, thus single horizontal lines are shown.  Table 1). Points represent plots with different overstory abundance (OA) values in the Early sample and different histories of thinning and fire (columns and symbols). Columns distinguish thinning and Early (2004) prescribed-fire treatments. Symbols represent the most recent fire ( p for prescribed fire and w for wildfire). Dashed and solid lines are modeled relationships for the Late-fire treatment based on statistically significant terms ( Table 2). For ingrowth, OA and the Late-fire treatment were not significant, thus single horizontal lines are shown.
Post-treatment mortality greatly exceeded ingrowth (mean: 100 trees/ha; range: 0-510; Figure 3), especially in plots subjected to late fire. The mortality model represented the data well (R 2 = 0.64; Table 2). Mortality was unaffected by thinning but was positively related to post-treatment OA. The relationship to OA was steeper in unburned plots than in those that burned early (significant OA × Early-fire interaction). Mortality was also significantly greater in plots that experienced late fire than in those that did not.
Basal area increment (BAI), the net growth of overstory trees, was slightly positive on average (mean: 0.2 m 2 /ha) but varied widely among plots (range: −14.9-6.8; Figure 3). BAI was unaffected by thinning, positively related to post-treatment OA in plots that burned early, and negatively related to OA in plots that did not burn (significant OA × Early-fire interaction) ( Table 2). Late burning reduced BAI irrespective of previous treatment history.  Table 1 for explanations of response variables and Tables A1-A6 (Appendix A4) for full model output. For the linear model of BAI, R 2 is the proportion of the variability explained by the full model. Values for individual terms are percentages of explained variability, summing to 100%. For GLMs (all other responses), R 2 is one minus the ratio of the model's deviance from being the best (perfect) and worst (null, intercept-only) model. Values for individual terms are percentages of deviance improvement with each added term, summing to 100%. Significant terms (p < 0.05) are underlined and bold. Terms that were not tested for a particular response are denoted by a dash, -.  Prior to treatment, seedlings were present in 6% of quadrats (range: 0%-35%; Figure 4). Frequency declined with OA, although the overall model was weak (R 2 = 0.07; Table 2). especially in plots subjected to late fire. The mortality model represented the data well (R 2 = 0.64; Table 2). Mortality was unaffected by thinning but was positively related to post-treatment OA. The relationship to OA was steeper in unburned plots than in those that burned early (significant OA × Early-fire interaction). Mortality was also significantly greater in plots that experienced late fire than in those that did not.
Basal area increment (BAI), the net growth of overstory trees, was slightly positive on average (mean: 0.2 m 2 /ha) but varied widely among plots (range: −14.9-6.8; Figure 3). BAI was unaffected by thinning, positively related to post-treatment OA in plots that burned early, and negatively related to OA in plots that did not burn (significant OA × Early-fire interaction) ( Table 2). Late burning reduced BAI irrespective of previous treatment history.  Table 2 for statistical results.  Table 2 for statistical results.
In the Early sample, seedlings were present in 7% of quadrats (range: 0%-30%; Figure 4). Frequency declined with OA in unthinned plots (as it had before treatment) but increased with OA in thinned plots (significant OA × Thin interaction) ( Table 2). Early burning enhanced seedling frequency.
In the Late sample, seedlings were present in more than twice as many quadrats (16%) as the Pre or Early samples (range: 0%-55%; Figure 4). Frequency was not related to OA or thinning ( Table 2). Early burning reduced frequency (in contrast to its effect in the Early sample). Late burning increased frequency, with a greater effect in plots that were burned early (significant Early-fire × Late-fire interaction).

Regeneration Dynamics: Colonization and Persistence of Seedlings (<1 m Tall)
The vast majority of quadrats lacked conifer seedlings (<1 m tall) prior to and after treatment. Seedling colonization was often balanced by loss, and seedling persistence through or after thinning or burning was negligible ( Figure 5). During the first interval (Pre-Early), seedlings colonized 0-6 quadrats per plot ( Figure 6). As a percentage of available quadrats-those lacking seedlings at the start of the interval-colonization averaged 6% per plot (median: 5%; range: 0%-31%). Colonization declined with pre-treatment OA in unthinned plots but increased with OA in thinned plots (significant OA × Thin interaction; Table 2). Colonization was also greater in early-burned than in unburned plots.
In the Early sample, seedlings were present in 7% of quadrats (range: 0%-30%; Figure 4). Frequency declined with OA in unthinned plots (as it had before treatment) but increased with OA in thinned plots (significant OA × Thin interaction) ( Table 2). Early burning enhanced seedling frequency.
In the Late sample, seedlings were present in more than twice as many quadrats (16%) as the Pre or Early samples (range: 0%-55%; Figure 4). Frequency was not related to OA or thinning ( Table  2). Early burning reduced frequency (in contrast to its effect in the Early sample). Late burning increased frequency, with a greater effect in plots that were burned early (significant Early-fire × Latefire interaction).
(a) Pre-Early Interval (b) Early-Late and Pre-Late Intervals Figure 5. Percentage of quadrats displaying four types of seedling dynamics (absence, colonization, loss, or persistence) in the (a) Pre-Early and (b) Early-Late and Pre-Late intervals. Each quadrat was assessed for seedling presence/absence at the start and end of the interval. ′Absent′ represents quadrats that had no seedlings at either date. ′Colonization′ represents quadrats lacking seedlings at Figure 5. Percentage of quadrats displaying four types of seedling dynamics (absence, colonization, loss, or persistence) in the (a) Pre-Early and (b) Early-Late and Pre-Late intervals. Each quadrat was assessed for seedling presence/absence at the start and end of the interval. Absent represents quadrats that had no seedlings at either date. Colonization represents quadrats lacking seedlings at the start but with one or more seedlings at the end of the interval. 'Loss' represents quadrats with one or more seedlings at the start but no seedlings at the end of the interval. Persistence represents quadrats with one or more seedlings at both dates. See Section 2.3 for additional details. Note that no thinned plots experienced both Early and Late fires.
Colonization was considerably more frequent during the second interval (Early-Late; range: 0-11 quadrats per plot; Figure 6). As a percentage of available quadrats, colonization averaged 15% per plot (median: 11%; range: 0%-55%). Colonization increased with post-treatment OA but was not affected by thinning (Table 2). Colonization was reduced by early burning (in contrast to the Pre-Early period) but greatly enhanced by late burning, particularly in plots that had burned earlier (significant Early-fire × Late-fire interaction). declined with pre-treatment OA in unthinned plots but increased with OA in thinned plots (significant OA × Thin interaction; Table 2). Colonization was also greater in early-burned than in unburned plots.
Colonization was considerably more frequent during the second interval (Early-Late; range: 0-11 quadrats per plot; Figure 6). As a percentage of available quadrats, colonization averaged 15% per plot (median: 11%; range: 0%-55%). Colonization increased with post-treatment OA but was not affected by thinning (Table 2). Colonization was reduced by early burning (in contrast to the Pre-Early period) but greatly enhanced by late burning, particularly in plots that had burned earlier (significant Early-fire × Late-fire interaction).   Table 2 for statistical results.
Colonization over the combined interval (Pre-Late) was very similar to that in the second interval (range: 0-11 quadrats per plot; Figure 6). As a percentage of available quadrats, colonization averaged 15% per plot (median: 11%; range: 0%-55%). Colonization was unrelated to pre-treatment OA, but was slightly inhibited by thinning and early burning ( Table 2). In contrast, late burning enhanced colonization, promoting greater colonization of unthinned than of thinned plots, and greater colonization of plots that had burned earlier (significant Thin × Late-fire and Early-fire × Late-fire interactions; Table 2).
Seedlings persisted less often than they colonized. During the Pre-Early interval, seedlings persisted in 0-2 quadrats per plot (Figure 7). As a percentage of available quadrats-those containing seedlings at the start of the interval-persistence averaged 23% (median: 0%; range: 0%-100%; Figure 7). Persistence was similarly low for the other intervals. Models (Table 2) failed to converge due to small sample sizes (n = 32-47 plots with seedlings present at the start of an interval) and low variance in responses (0-3 quadrats with persistence). However, data distributions suggest no relationships to OA, thinning, or burning history. enhanced colonization, promoting greater colonization of unthinned than of thinned plots, and greater colonization of plots that had burned earlier (significant Thin × Late-fire and Early-fire × Latefire interactions; Table 2).
Seedlings persisted less often than they colonized. During the Pre-Early interval, seedlings persisted in 0-2 quadrats per plot (Figure 7). As a percentage of available quadrats-those containing seedlings at the start of the interval-persistence averaged 23% (median: 0%; range: 0%-100%; Figure  7). Persistence was similarly low for the other intervals. Models (Table 2) failed to converge due to small sample sizes (n = 32-47 plots with seedlings present at the start of an interval) and low variance in responses (0-3 quadrats with persistence). However, data distributions suggest no relationships to OA, thinning, or burning history.
In 2015, small seedlings were present in an average of 16% (range: 0%-80%) of the 20, 1-m 2 quadrats per plot (Figure 9). Frequency was unrelated to OA, but was slightly greater in thinned than in unthinned plots, slightly reduced by early burning, and strongly enhanced by late burning ( Table  2). Large seedlings were present in an average of 8% (range: 0%-53%) of the 40, 5 × 5 m grid cells per plot (Figure 9). Frequency was unaffected by early burning. In contrast to small seedlings, frequency was consistently lower in plots that burned late than in those that did not. Relationships with OA were complex, varying with Thin and Late-fire treatments (significant OA × Thin and OA × In 2015, small seedlings were present in an average of 16% (range: 0%-80%) of the 20, 1-m 2 quadrats per plot (Figure 9). Frequency was unrelated to OA, but was slightly greater in thinned than in unthinned plots, slightly reduced by early burning, and strongly enhanced by late burning (Table 2).
Large seedlings were present in an average of 8% (range: 0%-53%) of the 40, 5 × 5 m grid cells per plot (Figure 9). Frequency was unaffected by early burning. In contrast to small seedlings, frequency was consistently lower in plots that burned late than in those that did not. Relationships with OA were complex, varying with Thin and Late-fire treatments (significant OA × Thin and OA × Late-fire interactions; Table 2): absent late burning, frequency declined with OA in unthinned plots, but increased with OA in thinned plots. With late burning, frequency declined more steeply in unthinned than in thinned plots (Figure 9). Similar to large seedlings, saplings were present in an average of 5% (range: 0%-35%) of grid cells (Figure 9), and their frequency was consistently reduced by late burning but unaffected by early burning. Relationships to OA were strong but dependent on thinning treatment: frequency declined with OA in unthinned plots but increased with OA in thinned plots (significant OA × Thin interaction; Table 2).
Forests 2020, 11, x FOR PEER REVIEW 14 of 29 Late-fire interactions; Table 2): absent late burning, frequency declined with OA in unthinned plots, but increased with OA in thinned plots. With late burning, frequency declined more steeply in unthinned than in thinned plots ( Figure 9). Similar to large seedlings, saplings were present in an average of 5% (range: 0%-35%) of grid cells (Figure 9), and their frequency was consistently reduced by late burning but unaffected by early burning. Relationships to OA were strong but dependent on thinning treatment: frequency declined with OA in unthinned plots but increased with OA in thinned plots (significant OA × Thin interaction; Table 2).  Table  2 for statistical results. Note the differences in scale of the vertical axes.

Discussion
Few experimental studies of restoration thinning or burning in dry forests have addressed the long-term dynamics of regenerating trees [15,[50][51][52]. Long-term studies of regeneration are notably absent in forests of the Pacific Northwest (but see [53]). Our results, spanning more than decade, highlight the changing relationships of regenerating conifers to fuels treatments, residual forest structure, and wildfire. Relationships between regeneration frequency and overstory abundance (OA) shifted from weakly negative prior to treatment to strongly positive after thinning, both with and without early burning. Subsequent higher-severity burns severed these relationships, promoting  Table 2 for statistical results. Note the differences in scale of the vertical axes.

Discussion
Few experimental studies of restoration thinning or burning in dry forests have addressed the long-term dynamics of regenerating trees [15,[50][51][52]. Long-term studies of regeneration are notably absent in forests of the Pacific Northwest (but see [53]). Our results, spanning more than decade, highlight the changing relationships of regenerating conifers to fuels treatments, residual forest structure, and wildfire. Relationships between regeneration frequency and overstory abundance (OA) shifted from weakly negative prior to treatment to strongly positive after thinning, both with and without early burning. Subsequent higher-severity burns severed these relationships, promoting abundant seedling recruitment while consuming previously established seedlings and saplings. Mortality greatly outpaced ingrowth, particularly in plots with greater OA (reflecting greater initial density), resulting in little increase in overstory density or canopy infilling. However, the recent pulse of seedling establishment could signal more rapid structural change in the future.

Seedling Relationships with Overstory Structure and Disturbance
We found that thinning temporarily altered the relationship between seedling frequency and overstory structure. Plots with greater OA tended to have fewer seedlings before thinning, but were colonized more frequently afterwards. Three factors may have contributed to this changing relationship. First, greater intensity of thinning in plots with greater initial density may have created more soil disturbance, enhancing seedbed conditions for germination of Douglas-fir [33,54]. Although plausible, this explanation runs counter to observations in other systems, wherein the primary effect of thinning is to add slash or fine litter that inhibits germination and early survival [19][20][21]. Further, thinned plots were yarded by helicopter, minimizing soil disturbance. Second, production of seed, which can persist for 1 to 2 years in the seedbed [33], may have been greater in plots with greater pre-treatment OA. Finally, seedling survival may have been higher in plots with greater pre-treatment OA-a legacy of greater overstory suppression of competing understory vegetation.
Early seedling responses to fire were also positive but, unlike thinning, interactions with overstory structure did not change: these remained negative in unthinned plots and positive in thinned plots. Despite the patchiness and low severity of early fires [46], burning appeared to stimulate recruitment, likely by enhancing conditions for germination (e.g., exposure of mineral soil, mineralization of soil nitrogen, and/or consumption of competing vegetation) [21,22,25,55,56].
By the Late sample, initially positive relationships with thinning and early fire were eclipsed by effects of later, higher-severity fires, which greatly elevated seedling frequencies. Although our Late-fire term does not distinguish between the 2006 prescribed fire and 2012 wildfire, plot data indicate considerably greater colonization after wildfire. Several factors may have contributed to this result. First, prescribed fire and wildfire occurred in different seasons (spring and fall, respectively). Fall burns tend to be more severe and less patchy [57][58][59], resulting in greater exposure of mineral soil, more complete consumption of plant cover, and greater seedling establishment [21,24,25,60]. Second, Late sampling occurred sooner after the wildfire than the prescribed fire (3 vs. 9 years), reducing the temporal window for mortality. Third, annual variation in seed production, temperature, and/or precipitation during these post-fire periods may have given rise to different patterns of germination and survival [28,35]. Future measurements that extend the temporal record on wildfire sites would help elucidate the relative importance of burn severity, time since fire, and weather for seedling establishment and persistence.

Responses to Thinning and Fire History among Regeneration Height Classes
In contrast to small seedlings, taller height classes of regeneration (large seedlings and saplings, sampled in 2015) have greater potential to influence forest structure and fire behavior in the future. However, it is more difficult to interpret their development relative to the disturbance histories of plots given the absence of pre-and early post-treatment data for these height classes. Small seedlings (<0.3 m tall) likely post-date Late-fire events, but it is less clear when larger seedlings (0.3-1.37 m tall) or saplings (>1.37 m tall) established. Estimates of height growth from the literature [61][62][63][64][65] suggest that large seedlings likely established after Early fire (2004), but saplings may comprise both pre-and post-treatment recruitment.
Large-seedling and sapling responses to thinning varied with overstory abundance (Late OA) and fire. In the absence of late burning, frequency declined with OA in unthinned plots, suggesting that overstory trees continue to suppress seedling establishment and growth. Conversely, frequency increased with OA in thinned plots, suggesting that thinning disturbance benefits establishment (similar to the colonization result), but that retained trees have a facilitative effect (e.g., moderating temperature or drought, or producing seeds) [15,21,54,66,67].
Large seedlings and saplings differed distinctly from small seedlings in their responses to fire. Early fire had a negative effect on small-seedling frequency, whereas higher-severity late fire had a strongly positive effect. In contrast, early fire had little effect on larger regeneration classes, but late fire greatly reduced their frequencies. The loss of larger, established seedlings and saplings to late fire is consistent with past studies of fire effects [21,53,60,68] and is a motivation for the use of repeated burning to achieve ecological and fuels-reduction objectives during restoration [25,69,70].

Implications of Fuels Treatments and Wildfire for Structural Change and Stand Development
We found little infilling of the overstory (ingrowth) more than a decade after fuels treatments. However, our results suggest that the nature, timing, and frequency of disturbance can affect both the rate and size structure of regeneration, and therefore can influence future trajectories of structural change. Saplings recruited to the overstory at low densities and in only one-third of the plots. In fact, the maximum rate of ingrowth over the 10-year, Early-Late interval was only 70 trees/ha. Tree mortality played a distinctly greater role in structural change, with an order of magnitude more trees dying than entering the overstory. Little of this mortality was expressed as a lagged response to early fire: most death after early burning was immediate and recorded in the Early sample [8]. In contrast, late burning caused significant tree death, particularly in wildfire plots, with rates of loss proportional to pre-fire density (greatest in unthinned, previously unburned plots). Although some mortality may be attributable to factors other than, or in addition to, fire (e.g., post-fire drought or beetle attack), field observations suggest that these factors played a minor role. Basal area increment (BAI) varied widely among plots, but was unaffected by thinning, indicating that greater rates of diameter growth in lower-density thinned plots were balanced by greater densities of slower-growing trees in unthinned plots.
Overall, long-term changes in overstory structure were driven, in large part, by mortality associated with later fires. In the absence of late burning, structural changes were more subtle, reflecting the growth of existing trees rather than recruitment of new stems. Late wildfire also played a key role in the long-term dynamics of understory conifers, shifting the size structure from a relatively sparse layer of taller stems to a denser cohort of newly established seedlings. Left untreated, the near-term consequences for forest structure will hinge on the balance among seedling survival, growth, and wildfire-processes that may be difficult to predict in a rapidly changing climate [40,71,72].

Conclusions
Mechanical thinning and prescribed burning reduce fire risk in the short term, but fuels-treatment longevity is a function of subsequent forest development, including regeneration. Long-term experiments are critical to assessing these dynamics and, thus, treatment longevity. In this study, thinning and prescribed fire interacted with overstory structure and wildfire to affect the spatial and temporal dynamics of regenerating conifers. Thinning shifted the relationship with overstory structure from weakly negative (inhibitory) to positive (facilitative) to neutral, although the strength of these relationships varied with time since treatment and subsequent disturbance. Low-severity prescribed fire stimulated regeneration in the short term, presumably by enhancing seedbed conditions for germination. However, the effect was transient. More than a decade after treatment (without subsequent disturbance), seedling and ingrowth densities remained low, suggesting that in this system managers do not need to be overly concerned about rapid regeneration or canopy infilling. Instead, fuels treatments can be tailored to other ecological objectives. However, we caution that additional work is required to determine whether our results can be extrapolated to systems dominated by other regenerating tree species or experiencing more severe prescribed fires.
Wildfires are becoming more prominent as climate change exacerbates summer drought, fuel conditions, and the susceptibility of dry forests to fire. Our results highlight the importance of studying how responses to planned management activities (fuels-reduction treatments) are altered by unplanned disturbances or, conversely, how responses to wildfire are mediated by treatments. In this experiment, low-severity wildfire erased legacies of pre-treatment overstory structure and post-treatment regeneration, slowing forest growth in the short term but triggering abundant recruitment. In the absence of additional disturbance or drought-related mortality of seedlings, wildfire may set the stage for rapid structural change and elevated fire risk in the future.
Although we were able to offer plausible explanations for the relationships of regeneration to fuels treatments and wildfire, our results also raise questions that warrant further study. For example, the possibility that thinning promoted regeneration by exposing mineral soil is not consistent with past studies or with the expectation that helicopter yarding minimizes ground disturbance. Further work is needed to understand this relationship. It is also not clear to what extent responses to prescribed fires (2004 and 2006) and wildfire (2012) reflect variation in fire severity, climatic factors, or time since disturbance. Demographic studies that incorporate quantitative measures of fire severity and climate could distinguish effects of disturbance from post-disturbance processes that affect seedling establishment and survival. Similarly, it is not clear whether regeneration responses to fuels treatment are mediated indirectly through effects on competing vegetation, or through legacies of overstory influence on understory. Incorporating metrics of pre-and post-treatment understory abundance into models of regeneration response may offer insight into the roles of competing vegetation.
In sum, our results suggest that fuels-treatment longevity hinges on complex interactions in space and time among stand-scale treatments, small-scale variation in forest structure, and unplanned disturbances. Managers must plan for and respond to this complexity, including future disturbances that are inevitable but unpredictable in timing and severity. Large-scale, long-term experiments in dry, fire-dependent forests are susceptible to these confounding disturbances but, in the context of adaptive management, offer unique opportunities for learning. Acknowledgments: Richy Harrod played a key role in the design and implementation of the Mission Creek Fire and Fire Surrogate Study and in early field measurements. We appreciate the efforts of our Mission Creek field crews who assisted with data collection. We thank Rae Mecka for her contributions to field logistics, review of historical data, and database management. We appreciate the helpful comments of three anonymous reviewers.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
This appendix provides background information for the project. First, we illustrate how overstory abundance (OA), our primary measure of overstory structure, relates to the two variables that are combined in its calculation. Next, we demonstrate the variability among plots in pre-treatment vegetation structure and physical environment. Third, we demonstrate how thinning intensity and burn severity varied among plots. Finally, we provide detailed output from the linear and generalized linear models summarized in Table 2 in the main text.

Appendix A.2. Pre-Treatment Vegetation Structure and Environment
We conducted a Principal Components Analysis (PCA) to assess the variability among plots (n = 72) in pre-treatment vegetation structure and physical environment. Variables included initial tree density, basal area, OA, mean shrub cover, mean herb cover, and adjusted heat load [73]. The analysis revealed considerable spread in ordination space among plots representing the same experimental units ( Figure A2). We conducted a Principal Components Analysis (PCA) to assess the variability among plots (n = 72) in pre-treatment vegetation structure and physical environment. Variables included initial tree density, basal area, OA, mean shrub cover, mean herb cover, and adjusted heat load [73]. The analysis revealed considerable spread in ordination space among plots representing the same experimental units ( Figure A2). Figure A2. First two principal components from a PCA of an environmental variables (aHLI; an index of heat load) and five pre-treatment vegetation variables: tree density (dens.plot_pre), basal area (BA.plot_pre), OA (OA.plot_pre), mean herb cover (Herb_pre), and mean shrub cover (Shrub_pre). Plots are color-coded by experimental unit. Figure A2. First two principal components from a PCA of an environmental variables (aHLI; an index of heat load) and five pre-treatment vegetation variables: tree density (dens.plot_pre), basal area (BA.plot_pre), OA (OA.plot_pre), mean herb cover (Herb_pre), and mean shrub cover (Shrub_pre). Plots are color-coded by experimental unit.

Appendix A3. Thinning Intensity and Burn Severity
Although we modeled thinning and burning treatments as categorical variables, thinning intensity and burn severity varied widely among plots.
Thinning intensity, expressed by the percentage change in OA over the Pre-Early interval, had a bimodal distribution, with a relatively narrow range of values in unthinned plots (−7.4% to 5.2%) and a much broader range in thinned plots (−77.1% to −29.5%) ( Figure A3).   Among 21 plots experiencing early fire (i.e., 2004 prescribed fires), burn severity-expressed by the percentage of trees charred within a plot-ranged from 3% to 97% ( Figure A4a). Among 36 plots experiencing late fire (2006 prescribed fires or 2012 wildfire), burn severity ranged from 62% to 100% ( Figure A4b). Three plots in prescribed-fire units did not experience early fire and were designated as unburned. Plots are color-coded by burn status (burned/unburned). These binary designations were used for 'Early fire' and 'Late fire' explanatory variables in applicable models (see Table 2, main text). Plots are color-coded by burn status (burned/unburned). These binary designations were used for 'Early fire' and 'Late fire' explanatory variables in applicable models (see Table 2, main text).