Forest Ecology and Management 389 (2017) 149–157 Contents lists available at ScienceDirect Forest Ecology and Management journal homepage: www.elsevier.com/locate/foreco Reserve tree mortality in two expanding-gap silvicultural systems 20 years after establishment in the Acadian forest of Maine, USA David R. Carter ⇑, Robert S. Seymour, Shawn Fraver, Aaron Weiskittel School of Forest Resources, University of Maine, 5755 Nutting Hall, Orono, ME 04469, USA a r t i c l e i n f o Article history: Received 23 August 2016 Received in revised form 21 December 2016 Accepted 28 December 2016 Keywords: Retention forestry Expanding-gap Retention tree Femelschlag Windthrow Late-successional a b s t r a c t Land managers are increasingly called upon to maintain ecosystem function and restore late-successional forest structures by retaining trees (‘‘reserves”) in harvest prescriptions. Such retention practices often result in mortality of reserve trees owing to increased exposure to wind and ‘gap shock’, thereby compromising management objectives. This study investigated the mortality of reserve trees (n = 787) retained over 20 years in the Acadian Forest Ecosystem Research Project (AFERP), a long-term ecological forestry experiment in central Maine, USA. Cumulative mortality across 18 species was very low (8.4% ± 3.3%, mean ± SD) relative to similar studies throughout the world. Annualized mortality of reserve trees of the two silvicultural systems was 1.1% for the large-gap treatment and 0.4% for the small-gap treatment. Cumulative mortality of Tsuga canadensis was lowest (2%) among species with >50 individuals, while the cumulative mortality of Thuja occidentalis was highest at 19%. Over 59% of cumulative mortality was wind-related. Annualized wind-related mortality of the reserve trees was 0.6% in the large-gap treatment and 0.2% in the small-gap treatment. More vigorous trees showed lower mortality rates. Cumulative mortality and wind-related deaths were somewhat higher in large gaps, and were influenced by spatial position within gaps. The non-traditional, expanding-gap, multi-aged silvicultural systems in this study may explain the lower incidence of windthrow and gap shock experienced by reserve trees relative to the annualized mortality rates from studies of similar, though predominantly single-cohort, silvicultural systems. Ó 2016 Elsevier B.V. All rights reserved. 1. Introduction Centuries of intensive, commodity-focused forest management have simplified forest structure worldwide, likely causing reductions in biodiversity and ecosystem services (FAO, 2010; Gustafsson et al., 2012; Lindenmayer et al. ,2012a). Recent decades have witnessed an increased interest in improving these management practices, intending to halt or reverse this simplification (Franklin, 1989; Seymour and Hunter, 1992; Franklin et al., 2007; Bauhus et al., 2009). Today, silvicultural practices in many regions are guided by our knowledge of natural disturbances, assuming that the complexity of post-disturbance forest structures will benefit both biodiversity and ecosystem services. For example, retaining particular forest structures at harvest reduces the often striking contrast between intensively managed forests and those allowed to develop under natural conditions (Hunter, 1990; Lindenmayer et al., 2012b). ⇑ Corresponding author at: Department of Forest Resources, University of Minnesota, Green Hall, 2005 Upper Buford Circle, St. Paul, MN 55108, USA. E-mail address: david.r.j.carter@gmail.com (D.R. Carter). http://dx.doi.org/10.1016/j.foreco.2016.12.031 0378-1127/Ó 2016 Elsevier B.V. All rights reserved. The structures targeted for retention -- be it in a dispersed or aggregated configuration (Franklin et al., 1997) -- are often large or old trees, or trees with particular attributes such as cavities, that may promote biodiversity (Hunter, 1990; Gustafsson et al., 2010). In silviculture, trees retained for purposes other than seed propagation and light regime mitigation are termed ‘‘reserves” (or reserve trees), defined as any tree, pole-sized or larger, that is left unharvested after the conclusion of the regeneration period, typically to enhance biodiversity, improve aesthetics, or grow to biological maturity (Helms, 1998, p. 153). Though reserve trees produce seed and effect the light regime in the systems in which they are retained, this is not the justification for their retention in these systems. Instead, reserve trees are retained to provide continuity from the pre- to post-harvest stand and thus serve as refugia or ‘lifeboats’ for associated organisms as the stand recovers and develops following harvest (Gustafsson et al., 2010). Maintaining these features on managed landscapes serves as a ‘‘coarse-filter” approach (Hunter, 1990; Manning et al., 2006) for sustaining critical habitat for many species, supplying ecosystem services, and maintaining structural diversity (Lorimer and White, 2003; Lindenmayer et al., 2012b; Fedrowitz et al., 2014; Hämäläinen 150 D.R. Carter et al. / Forest Ecology and Management 389 (2017) 149–157 et al., 2016). These features also create the potential for restoration of environmental conditions associated with structurally complex forests (Franklin et al., 1997), and they contribute significantly to carbon storage and sequestration (Stephenson et al., 2014). Realizing many of these ecological benefits that retained trees provide necessitates their survival while the system recovers post-harvest. Despite the growing interest and widespread adoption of retaining biologically important trees at the time of harvest (often called retention forestry, sensu Lindenmayer, 2012b), few empirical studies have evaluated the long-term fate of the reserve trees themselves (but see Buermeyer and Harrington, 2002; Busby et al., 2006; Jönsson et al., 2007; Walter and Maguire, 2007; Urgenson et al., 2013). Of particular concern for forest managers is (1) the risk of windthrow for recently isolated retained trees that have not yet developed wind-firmness (Issac, 1940; Lohmander and Helles, 1987) and (2) the risk of crown dieback and poor vigor sometimes experienced by recently-isolated trees (Kozlowski et al., 1991; Jönsson et al., 2007), that is, the ‘thinning shock’ or ‘gap shock’ well known to field foresters. Studies that address these risks over extended periods are uncommon simply because retention forestry has not been in practice long enough to provide such data. Previous research on windthrow following partial harvest shows that retained trees are initially at increased risk of wind damage as a result of exposure for a short period of 2 to 5 years, after which vulnerability decreases (Scott and Mitchell, 2005; Busby et al., 2006; Gibbons et al., 2008; Rosenvald et al., 2008). Elevated rates of mortality can persist beyond a decade, however (Ruel et al., 2000; Caspersen, 2006). Informed selection of retained trees could aid in minimizing losses of retained trees to windthrow and gap-shock. Some species (e.g. Abies balsamea (Ruel, 2000) and Picea rubens (Canham et al., 2001)) and trees with particular attributes have shown to be more vulnerable to wind- and exposure-related deaths than others. Some trees may become more wind firm over time in harvest openings by increasing diameter growth (Holgen et al., 2003), strengthening root systems (Peterson, 2004), or modifying canopy shape (Canham et al., 2001). Even within a given region, species’ relative resistance to uprooting can vary with soil type (Élie and Ruel, 2005). Patterns in the direction of windthrow can also provide management insight and inform actions that may reduce the risk of wind damage by mitigating the exposure of individual trees (Rosenvald et al., 2008). In this study, we capitalize on a long-term, well-replicated experiment, namely the Acadian Forest Ecosystem Research Program (AFERP), designed to test ecological forestry concepts. The experiment – established in 1994 – includes two expanding-gap regeneration systems with reserve trees (Seymour et al., 2006), and it conveniently includes a range of tree species, tree sizes, and gap sizes. The experiment is located on the Penobscot Experimental Forest (PEF) in east-central Maine, USA. Here we explore both the spatial and temporal aspects of dispersed retention tree mortality over a 20-year period. Our specific objectives were to (i) assess treatment effects on reserve tree mortality; (ii) determine meaningful tree-level predictors of mortality; and (iii) investigate the influence of exposure on reserve tree mortality. 2. Methods 2.1. Study site AFERP is located in the Penobscot Experimental Forest (PEF), of central Maine, USA (44_510 N, 68_370 N). The PEF is 1618 ha and part of the Acadian Forest Region. The Acadian Forest Region has a cool, humid climate; it is characterized as a transitional ecotone between the northern boreal, spruce-fir forest type, and southern broadleaf forests. The mean annual temperature is 7.06 C°, with approximately half of the 106 cm of precipitation falling in 156 days from May through October (Arsenault et al., 2011). The soils are derived from glacial till and range from well-drained loams, stony loams and sandy loam ridges, to poorly drained loams and silt loam flat areas (Saunders and Wagner, 2008; Arsenault et al., 2011). Softwoods in the region are red, white, and black spruce (Picea rubens, P. glauca, and P. mariana, respectively), balsam fir (Abies balsamea), eastern white pine (Pinus strobus), eastern hemlock (Tsuga canadensis), and northern white cedar (Thuja occidentalis). Common hardwoods found in this region include red maple (Acer rubrum), quaking and big-tooth aspen (Populus tremuloides and P. grandidentata), and paper and yellow birch (Betula papyrifera and B. alleghanienses). 2.2. Project background The silvicultural prescription in AFERP is loosely based on the German ‘‘Femelschlag” with two expanding-gap silvicultural regeneration methods designed as multi-aged, naturaldisturbance-based silvicultural systems. The area harvested per cutting cycle, by design, emulates the 1% annual disturbance frequency common to the region (Seymour et al., 2002, 2006; Arsenault et al., 2011). There are nine, approximately 10 ha research blocks (i.e. research areas; hereafter ‘‘RAs”) to which treatments were randomly assigned; six are treated while three serve as unmanaged controls. Three RAs are treated with an irregular group shelterwood with reserves in which 20% of the total area in the RA is harvested every 10 years (hereafter ‘‘large-gap treatment”) for 5 cutting cycles. The remaining three blocks were treated with a group selection treatment in which 10% of the total area -- half that of the large-gap treatment -- in the RA is harvested every 10 years (hereafter ‘‘small-gap treatment”) for 10 cutting cycles. The large-gap treatment was designed to encourage natural regeneration of tree species of intermediate shade tolerance and to maintain stands of mid-successional status. The small-gap treatment was designed to encourage shade-tolerant species and accelerate development of late-successional stands. When possible, initial gap harvests were centered on areas of sufficient understory stocking of desired species. Occasionally these areas were the product of small, pre-existing natural gaps formed by single-tree mortalities. Both the large- and small-gap treatments are applied using a 10-year cutting cycle. After their respective 10- and 20-year regeneration period (i.e. elapsed time before harvest re-entry), the gaps are asymmetrically expanded. At the time of this study, both treatments had undergone two entries. Gaps in the large-gap treatment had received a single expansion and the small-gap treatment had two entries of gap-creation without expansion. Two RAs are harvested in a year, one large- and one small-gap treatment, requiring three successive years to complete a project-wide cutting cycle. AFERP was established in 1994 and the initial treatments began in 1995 and ended in 1997, while the second entry spanned from 2005 to 2007. These harvests were conducted in the winter when conditions were optimal to reduce logging damage. Trees were felled manually, delimbed at the stump, then skidded tree-length to the roadside with small cable skidders. Prior to treatment, the stands that comprise the research areas were fairly homogeneous. There were no significant differences in density (2404 ± 138 trees ha 1, mean ± SD), basal area (BA) (37.6 ± 1.1 m2 ha 1), or volume (283.6 ± 10.4 m3 ha 1) among RAs prior to treatment and, compositionally, balsam fir, red maple, and hemlock were present at the greatest densities (Arsenault et al., 2011). D.R. Carter et al. / Forest Ecology and Management 389 (2017) 149–157 In these large- and small-gap treatments, retention exists in three conditions: unharvested (i.e. yet-to-be harvested) matrix, permanent reserve trees in gaps, and temporary ‘‘overwood” trees (defined later) in gaps that can be harvested in future entries. The latter two conditions are retained in a dispersed configuration throughout the harvested gaps. The permanent reserve trees represent the only retention of the three aforementioned conditions that are to remain after the entirety of the given RA has been harvested. These trees are therefore the focus of this study. Within the treated RAs, the permanent, reserve tree retention ranged from 11 to 18% of the initial stand BA. The original study plan (ca. 1994) specified permanent retention to be 10% of the stand BA prior to treatment, but actual retention has been slightly higher. These trees were intended to be retained permanently as dispersed reserve trees, distributed throughout the gaps as they are created. Reserve trees were selected to have the following characteristics: (1) possessing cavities and wildlife value, (2) locally uncommon tree species, (3) relatively large size, and (4) the potential for high timber value under a long rotation. A total of 18 species was retained as permanent reserves. Softwoods retained were balsam fir, red spruce, red pine (Pinus resinosa), white pine, northern white-cedar, and hemlock. Hardwoods were red maple, sugar maple (A. saccharum), eastern shadbush (Amelanchier canadensis), yellow birch, paper birch, American beech (Fagus grandifolia), white ash (Fraxinus americana), black ash (F. nigra), ironwood (Ostrya virginiana), big-tooth aspen, trembling aspen, and red oak (Quercus rubra). All of the reserve trees were permanently tagged and georeferenced using a Trimle GeoXH with sub-meter precision. Within the treated RAs, an additional 1–8% of the initial stand BA was retained as temporary ‘‘overwood”. In contrast to reserve trees, overwood trees are those that are retained either to grow into merchantable size-classes (e.g. immature growing stock) or as additional overstory cover in areas of low seedling density (e.g. low advance-regeneration stocking). The presence of overwood in areas of low advance-regeneration stocking can promote the establishment of desired species by mitigating the understory light regime to deter the growth of unwanted, shade-intolerant species. planar break; fungal-related). Uprooted trees were assumed to have died from wind-related causes. The height, DBH, and position within the gap were collected for each down tree and the azimuth of the direction of fall (degrees) was determined for all windrelated deaths. Some data were derived from previous inventories. When height or LLC could not be collected for down trees because decay or damage prevented accurate reconstruction of the formerly living stem or branches, these data were obtained from the most recent inventory of the tree while it was still living. Previous DBH measurements were used to calculate the annual diameter increment (ADI; growth rate). All overwood trees over 25 cm in DBH were inventoried and used in part to calculate retention density within gaps, a predictor variable. Although a smaller DBH threshold could have influenced some of the trends observed in the analysis, we believe the 25-cm DBH threshold accurately defines trees that were likely to be in the main canopy prior to initial gap creation. The initial and expanded gap areas were calculated, as were metrics related to reserve tree exposure. These included (1) nearest distance to original gap edge (i.e. reserve trees in unexpanded gaps (small-gap system only) and gaps prior to expansion (large-gap treatment only)), (2) nearest distance to expanded gap edge (i.e. reserve trees in the expanded gaps (large-gap treatment only)), and (3) nearest distance to original gap edge from within the gap expansion (large-gap treatment only) (Fig. 1). This final metric measures the degree of partial exposure experienced by reserve trees prior to gap expansions and any potential influence partial exposure may have had on decreasing mortality (i.e. ‘‘preadaptation” – a partial exposure that may gradually prepare a tree for full exposure; Rosenvald et al. (2008)). These metrics were computed in ArcGIS, ver. 10.2 for each reserve tree based on its history of exposure through the study period. If a tree died prior to the expansion of the gap it was in, based on a previous inventory, this was accounted for when describing its exposure history. Mortality rates can be calculated in many ways (Sheil et al., 1995). We used the common instantaneous rate measure, k, calculated using the following formula: k¼ 2.3. Data collection Since 1994, the RAs in AFERP have been inventoried on a 5-year interval that corresponds to each RA’s respective harvesting schedule. Using 0.05 ha fixed-area circular sampling plots, overstory trees and saplings are monitored and stem-mapped (i.e. distance and azimuth from plot center). These plot data were used to locate ‘‘paired analogues” (defined later) of the reserve trees in the unharvested control. The complete set of reserve trees is inventoried on the same 5-year interval. Overstory trees and saplings in the fixedarea plots and the reserve trees each have the following attributes recorded and monitored: species, diameter at breast height (DBH; cm), total height (m), lowest-live-crown (height to lowest live whorl; LLC; m), and condition (live, declining, or dead). In the summer of 2014, all reserve trees (n = 787) within the six treated RAs were re-inventoried. For each stem, we recorded species, DBH, total height, height to lower live whorl (LLC), condition (live or dead), and position within the gap (x-y coordinates used to relocate stems and compute distance to gap edge). If a tree was found dead, decay-class, using a five-class system, (USDA, 2012) was noted. We recorded the condition of each tree with structural failure as one of three classes: trunk broken, brash, or uprooted (based on the protocol of the International Tree Failure Database (Smiley et al., 2006)). If the trunk was broken, the stem was visually inspected for evidence of wind-related damage (vertical shear of wood fibers; wind-related) or brash failure (horizontal, 151 lnðn0 Þ lnðnt Þ t Fig. 1. An aerial perspective of an expanded gap. Green represents the entirety of the gap in 2014 with purple representing uncut, forest matrix. The line within the gap is the original gap boundary prior to expansion (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.). 152 D.R. Carter et al. / Forest Ecology and Management 389 (2017) 149–157 where n0 is the number of stems in the initial population of reserve trees and nt is the number of stems surviving to time, t. Time, t, is the number of years between two censuses. The annualized mortality was calculated for each treatment by averaging across the rates for each entry year. The same approach was used to calculate the wind-related mortality, separately. When data were available, annualized wind-related mortality rates were also calculated from previous studies and included for comparison in the Supplementary Table S.4. 2.4. Statistical analysis 2.4.1. Treatment effects The experiment is a replicated design with research areas (RAs) treated as random, dictating a multi-level mixed-effects analysis. The replicates of the treatments occur at the RA-level. To assess treatment effects on reserve tree mortality (objective (i)), windrelated and overall mortality of the small- and large-gap treatments were analyzed using a mixed-effects ANOVA in R 3.2.1 (R Project, 2015). Treatment was considered a fixed effect and RA, a random effect. Overall mortality rates in the control were compared to those in the treated RAs for conspecific trees of comparable size. This comparison was done by stratifying the DBH data from the treated RAs and the control into 5 cm size-classes to match the species composition and size-class distribution of each other (hereafter referred to as ‘‘paired analogues”). Diameter distributions did not differ between treatments (KS test; p > 0.05). Direction of fall data from reserve trees within the harvested gaps were used in a circular statistical software package, Oriana 2.0 (Kovach Computing), to determine directional uniformity within treatments. For all tests, we used an alpha of 0.05 as the basis for statistical significance. 2.4.2. Tree-census and exposure variables To address objectives (ii) and (iii), we developed two types of models, one based on tree-census variables, the second based on tree-exposure variables. This division reflects the two distinct stages in implementing this silvicultural prescription: (1) reserve tree selection based on tree attributes and (2) spatial aspects of the harvest-gap layout. Within these separate analyses, generalized linear mixed- models (RA as the random variable) with a binomial error distribution and logit link were constructed to predict response variables (1) wind-related mortality and (2) overall mortality. In analyzing the reserve tree (i.e. control data excluded) treecensus and exposure variables, seven reserve trees of three species were omitted due to their low abundances: (Abies balsamea (n = 4), Amelanchier canadensis (n = 1), Ostrya virginiana (n = 2)). Similarly, due to small sample sizes, Fraxinus spp., and Populus spp. were collapsed into respective genera. Tree-census variables consisted of species, crown ratio, ADI (annual diameter increment), slenderness, DBH, height, and cavity presence. Crown ratio was calculated by subtracting LLC (lowest live crown) from height and dividing the total by height. ADI was calculated by using previous diameter measurements and averaging the growth rates between inventories for each reserve tree after the initial harvest. Slenderness was calculated by dividing height by DBH. Tree-census variable means can be found in Supplementary Table S.1. Exposure variables consisted of treatment, retention BA, retention stem density, gap area (ha), duration of exposure (years), distance to the nearest original gap edge (m), distance to the nearest expanded gap edge (m), and distance to the nearest original gap edge from within the gap expansion (m) (Fig. 1). Exposure variable means can be found in Supplementary Table S.2. Model construction was conducted in three steps: (1) predictors were assessed for significance in ANCOVAs per response variable because of the large number of variables in this study. (2) Predictors with low p-values were then incorporated into generalized linear mixed models. Two-way interactions were considered. ‘‘Species” and ‘‘Treatment” interactions with tree-census and exposure variables, respectively, were tested. Higher-level interactions were not considered due to the observational approach used here. In such an approach, higher-level interactions are often uninterpretable. Post-hoc pairwise comparisons of overall mortality between the treated RAs and the control among species were conducted using the Tukey’s Honest Significant Difference (Tukey’s test) function in the multcomp package (Hothorn et al., 2008). Due to the relatively small sample sizes for some species in this study, models predicting the same response, i.e. wind-related or overall mortality, were compared using the Akaike information criterion with a small-sample bias adjustment (AICc) using the AICc function in the ‘AICcmodavg’ library in R (Burnham and Anderson, 2002; Mazerolle, 2015). The AICc was used to determine the model of ‘best fit’ and parsimony for a given response variable within a set of models. Models with the lowest AICc were considered the highest supported model. (3) Model performances were finally assessed using the ‘area under the curve’ (AUC) output from the pROC function in the ‘pROC’ library in R (Robin et al., 2011). The AUC value reflects the accuracy of a model. It is the percentage of correctly classified random pairs in a binary model and is an indicator of the ability of the model, based on the predictive strength of its independent variables, to correctly classify binary outcomes. 3. Results 3.1. Mortality of reserve trees Overall, the cumulative mortality rate of reserve trees was 8.4% ± 3.3%. The annualized mortality rate for the reserve trees in the large-gap was 1.1% and 0.4% in the small-gap. In total, 66 out of 787 reserve trees died since the onset of the experiment in 1994. Wind-related mortality accounted for 59.1% of the total mortality (n = 39); 33.4% by trunk breakage and 25.7% by uprooting (Table 1). The annualized wind-related mortality rate was 0.2% in the first entry and 0.6% in the second entry. Other conditions contributing to the overall total included 18.2% standing dead and 22.7% fallen as the result of brash failure (Table 1). For species with 20 or more individuals, Thuja occidentalis (19.0%), B. papyrifera (16.7%), and Picea rubens (11.1%) had the highest cumulative mortality rates (Table 1). 3.2. Treatment effects Wind-related and overall mortality, tested separately, were significantly greater (p = 0.004 and 0.004, respectively) in the largegap treatment (cumulative: 6.3% and 10.4% (Fig. 2); annualized: 0.6% and 1.1%) than in the small-gap treatment (cumulative: 2.3% and 4.2% (Fig. 2); annualized: 0.2% and 0.4%). There was no significant pattern in the direction of fall for either the large- or small-gap treatment (Rayliegh test p = 0.161 and p = 0.33, respectively; Fig. 3). The large-gap treatment appeared to maintain a diametrically bimodal direction of fall pattern, but after using the standard angle-doubling procedure (Mahan, 1991), this relationship proved to not be significant (p > 0.05). The greater rate of mortality among Abies balsamea (p = 0.06) and Picea rubens (p = 0.0002) paired analogues in the control relative to conspecific reserve trees in the treated RAs resulted in a significantly lower overall mortality in the treated RAs (p < 0.001; 8.1% and 17.8%, respectively; Table 2; Supplementary Table S.3.). 153 D.R. Carter et al. / Forest Ecology and Management 389 (2017) 149–157 Table 1 Number of reserve trees in harvested RAs and mortality data for different species covering the duration of the study. Tree species Softwoods Abies balsamea Picea rubens Pinus resinosa Pinus strobus Thuja occidentalis Tsuga canadensis Hardwoods Acer rubrum Acer saccharum Amelanchier canadensis Betula alleghaniensis Betula papyrifera Fagus grandifolia Fraxinus americana Fraxinus nigra Ostrya virginiana Populus grandidentata Populus tremuloides Quercus rubra Total Trees that died during the study No. Mortality (%) Brash Died standing Trunk broken Uprooted Wind-related deaths (%) 4 153 12 112 84 146 0 11.1 8.3 8.9 19 2.1 0 4 0 3 1 0 0 1 0 2 2 2 0 6 0 4 6 0 0 6 1 1 7 1 0 70.6 100 50 81.3 33.3 119 33 1 12 24 17 22 1 2 9 18 18 3.4 3 100 8.3 16.7 11.8 0 0 0 0 33.3 0 2 0 0 0 3 1 0 0 0 0 1 0 0 0 1 1 1 0 0 0 0 0 2 0 2 0 0 0 0 1 0 0 0 0 3 0 0 1 0 0 0 0 0 0 0 0 0 0 50 100 0 0 0 50 0 0 0 0 50 0 787 8.4 15 12 22 17 59.1 Fig. 2. The overall and wind-related (species pooled) mortality in the large-gap and small-gap treatments. These figures are intended to convey differences in mortality related to exposure resulting from treatment. No other within-species comparisons between the control and treated RAs were significant. 3.3. Tree-census variables An ANCOVA of the tree-census variables predicting windrelated mortality showed that species (p = 0.04), crown ratio (p = 0.04), ADI (p = 0.005), and slenderness (p = 0.02) were the only significant variables (Table 3). An ANCOVA of the tree-census variables predicting overall mortality showed that species (p < 0.001) and ADI (p < 0.001) were significant and crown ratio (p = 0.06) was marginally significant (Table 3), each addressed below. The highest supported generalized linear model of tree-census variables predicting wind-related mortality included ADI (coefficient: 10.53) and slenderness (coefficient: 4.15; Table 4; AUC = 0.91). The highest supported generalized linear model of tree-census variables predicting overall mortality included ADI (coefficient: 7.25), only (Table 4; AUC = 0.85). 3.4. Exposure variables An ANCOVA of the exposure variables predicting wind-related mortality showed that treatment (p = 0.02), duration of exposure (p < 0.001) and distance to expanded gap edge (p = 0.04) were significant predictors (Table 5). Similarly, treatment (p = 0.004) and duration of exposure (p < 0.001) were the only significant variables and distance to original gap edge (p = 0.06) was marginally significant in predicting overall mortality (Table 5). The highest supported generalized linear model of exposure variables predicting wind-related mortality included treatment, duration of exposure (coefficient: 0.20), and distance to expanded gap edge (coefficient: 0.04; AUC: 0.71; Table 6). The large-gap treatment had significantly more cases of wind-related mortality than the small-gap treatment. Wind-related mortality showed a significant negative relationship with duration of exposure and a marginally significant positive relationship with distance to expanded gap edge (Fig. 4). In these data is the presence 154 D.R. Carter et al. / Forest Ecology and Management 389 (2017) 149–157 Fig. 3. Circular histograms of direction of fall in degrees for wind-related mortality of trees by treatment. Table 2 Significance of treatment (harvested vs. control), species, and their interaction from the ANOVA model of overall mortality. Variable F-ratio p-value Treatment (Levels: Harvested, Control) Species (Levels: 13) Treatment * Species 48.45 35.75 5.03 <0.001 <0.001 <0.001 Table 3 Result of ANCOVA models testing tree-census predictors for wind-related mortality and overall mortality. Tree-census variables Species Crown ratio ADI (cm yr 1) Slenderness (HT/DBH) DBH (cm) HT (m) Cavity (binary) Mean Wind-related mortality Overall mortality p-value p-value Levels = 13 0.53 0.46 0.60 0.04 0.04 0.005 0.02 <0.001 0.06 <0.001 0.30 35.5 20.0 n = 17 0.29 0.78 0.37 0.64 0.99 0.28 of a data-point with high leverage (distance to expanded gap edge = 47.5 m). With this data-point removed, the marginal significance of the predictor is lost but the direction of relationship remains the same. This point was not an outlier so it remained in the analysis but its influence should be noted. Table 5 Summary table of tested exposure predictors from ANCOVA in explaining the variance in wind-related mortality and overall mortality. Exposure variables Mean Wind-related mortality (n = 39) p-value Overall mortality (n = 66) p-value Treatment (gap size; m2) Large-gap: 6325 Smallgap: 1157 4628.8 9.6 0.02 0.004 0.13 0.51 0.51 0.65 98.0 0.97 0.87 11.3 <0.001 <0.001 6.6 0.60 0.06 12.71 0.04 0.30 10.78 0.65 0.60 Gap area (m2) Retention density BA (m2 ha 1) Retention density (stems ha 1) Duration of exposure (years) Distance to original gap edge (m) Distance to expanded gap edge (m) Distance to original gap edge from expansion (m) The highest supported model of exposure variables in predicting overall mortality included treatment, duration of exposure (coefficient: 0.28), and distance to original gap edge (coefficient: 0.11; AUC: 0.72; Table 6; Fig. 5). The large-gap treatment had significantly more incidences of mortality than the small-gap treatment. Table 4 Summary of the most highly supported generalized linear model assessment of tree-census variables in predicting wind-related mortality and overall mortality (AUC = Area Under Curve; AICc = Akaike’s corrected Information Criteria). Wind-related mortality model Wind-related mortality AUC Wind-related mortality AICc ADI + Slenderness ADI + Species * Slenderness Species * ADI + Slenderness Species + ADI + Slenderness Overall mortality model ADI Species * crown ratio Crown ratio Species 0.91 0.96 0.97 0.97 Overall Mortality AUC 0.85 0.82 0.66 0.73 96.2 99.4 118.9 119.8 Overall Mortality AICc 162.3 330.4 333.2 431.9 Wind-related mortality DAICc 3.2 22.7 23.6 Overall Mortality DAICc 168.1 170.9 269.6 D.R. Carter et al. / Forest Ecology and Management 389 (2017) 149–157 155 Table 6 Summary of the most highly supported generalized linear model assessment of exposure variables in predicting wind-related mortality and overall mortality (AUC = Area Under Curve; AICc = Akaike’s corrected Information Criteria). Windrelated mortality DAICc Wind-related mortality model Windrelated mortality AUC Windrelated mortality AICc Treatment + duration of exposure + distance to expanded gap edge Treatment + duration of exposure Duration of exposure + distance to expanded gap edge Duration of exposure Treatment Treatment + distance to expanded gap edge Distance to expanded gap edge Overall mortality model 0.71 289.1 0.70 290.3 1.2 0.59 293.2 4.1 0.66 0.58 0.68 297.4 307.3 308.8 8.3 18.2 19.7 0.62 Overall Mortality AUC 0.72 311.4 Overall Mortality AICc 404.2 22.3 Overall Mortality DAICc 0.73 405.4 1.2 0.72 407.3 3.1 0.69 415.1 10.8 0.69 419.1 14.9 Treatment + duration of exposure + distance to original entry edge Treatment + duration of exposure + distance to original entry edge + distance to expanded gap edge Treatment + duration of exposure Duration of exposure + distance to expanded gap edge Duration of Exposure Fig. 5. Overall mortality risk as distance from first entry edge increases. ence of dispersed, emergent Pinus strobus canopies and the relatively small gap sizes in AFERP may have created an aerodynamically unique system that is less vulnerable to catastrophic winds. Alternatively, Ruel (1995) argued that the most influential factor regarding mature tree mortality in partial cuttings and variable retention harvests is the condition of adjacent stands. The adjacent stands in AFERP are mostly fully stocked forest matrix or relatively light partial harvests. The greater mortality of conspecific trees in the uncut control blocks is evidence that the beneficial effects of releasing these species from competition outweighed risks from exposure as reserve trees. The lack of any directional pattern of windthrow could be explained by storm-generated winds. Taylor and MacLean (2007) found that as wind-stress caused by the prevailing wind direction stimulates growth that resists windthrow in that direction, the tree becomes more asymmetrical and vulnerable to windthrow by nonprevailing winds generated by powerful storms that gust in directions in which the tree is not yet adapted. 4.1. Tree-Census variables Fig. 4. Wind-related mortality risk of reserve trees as distance increases from expanded gap edge. 4. Discussion Our documented overall mortality rates (8.4% ± 3.3% cumulative; 1.1% and 0.4% annualized overall mortality for the largeand small-gap, respectively) are much lower than most other studies, with surprisingly low wind-related mortality rate over the entire 20-year period (4.8% cumulative; 0.6% and 0.2% annualized for the large- and small-gap, respectively; Table S.4.). Most other studies were conducted in single-cohort silvicultural systems in which reserve trees are isolated in large, heavily cut blocks, whereas in this study, reserves trees were left in much smaller gaps (most well under one hectare) surrounded by unharvested matrix conditions. O’Hara (2014) postulates that complex, multiaged structures like those created in this study may have greater resistance to windthrow because individual trees are preconditioned to wind-stress prior to emerging into the canopy. The pres- Metrics associated with vigor – crown ratio and high growth rate – were negatively associated with mortality, which is in agreement with the results of Gibbons et al. (2008). Our findings concur with several mechanistic theories in the literature that assert that vigorous trees are less vulnerable to wind-related mortality. Other studies have demonstrated that, following a harvest, certain trees may maintain greater adaptive capacities that decrease mortality (e.g. Canham et al., 2001; Holgen et al., 2003; Peterson, 2004). From a land manager perspective, these results and others similar to it would support retaining trees with large crowns when growth rate data are unavailable as large crowns are associated with high growth rate and vigor. This contradicts findings from other studies in which large crowns were often associated with increased windthrow as a result of a presumed ‘‘sail-effect” creating increased wind-stress (Zeng et al., 2004). Crown ratio was negatively associated with wind-related deaths in this study, perhaps due to its associations with tree vigor. In our data, crown ratio was significant when predicting ADI (coefficient: 0.41, p = 0.0002). Cavity presence, surprisingly, was not a significant predictor of windrelated death or overall mortality. Our finding that increased slenderness (high ratio of height to diameter) is associated with decreased wind-related deaths was likely an artifact of the relatively high occurrence of wind-related mortality in Thuja occidentalis -- a species with highly tapered boles and the lowest average slenderness (0.53) value across 156 D.R. Carter et al. / Forest Ecology and Management 389 (2017) 149–157 treatments. The slenderness-species interaction was significant when tested, though not a component of the highest supported model predicting wind-related mortality (Table 4). When the data were subset by the five predominant species and tested, slenderness was only significant in predicting wind-related mortality for Thuja occidentalis (coefficient: 3.2; p = 0.04). Species showed marked variation in mortality rates. Mortality of Thuja occidentalis was greatest of any species with a large population size (n > 50). This species is commonly found on wet sites with shallow rooting zones, and thus naturally more prone to windthrow (Élie and Ruel, 2005). Relative to our silvicultural objectives, however, mortality of only 19% is clearly a positive outcome. Tsuga canadensis had the lowest mortality of all species with a large population size. This species is commonly retained as a reserve tree owing to its relatively low economic but high ecological value as a long-lived, wildlife species. Picea rubens, a shallow-rooted, highvalue species often considered prone to windthrow (Cary, 1896; Canham et al., 2001; Fraver and White, 2005), maintained a low mortality rate in this study (11.1%; Table 1) and had a significantly lower mortality rate than its paired analogues in the control (Table S.3). This finding suggests that concerns commonly held by land managers over the poor windfirmness of Picea rubens in this region may be exaggerated or perhaps influenced by edaphic or topographical factors not common in this study area. Conversely, this species could simply be responding favorably to these treatments. toring of these systems as gaps are increasingly expanded will be critical in affirming this hypothesis. 4.3. Management implications Based on these findings, to realize the benefits of live-tree retention, silvicultural prescriptions should gradually increase exposure of retained trees and target the retention of trees with high growth rates and high crown ratios in order to decrease mortality of reserve trees. Retaining Tsuga canadensis and A. rubrum over B. papyrifera, Populus spp. and Thuja occidentalis would also decrease overall mortality rates. Adherence to these guidelines, however, would compromise the objective of retaining a diversity of tree species of varying age and size. Selecting the individual of each species with high vigor may thus be the optimum strategy. It is important to note that modest levels of reserve-tree mortality, as observed here, are arguably desirable, as they serve to enrich the pool of large downed woody material that is often deficient on the forest floor of managed forests. Future studies on reserve trees should focus on metrics defining species-specific vigor and exposure risk. Ideal metrics are those that can be quickly assessed in the field. Studies may also consider exploring the effects of tree retention on maintaining predisturbance environmental conditions post-harvest. Overall, the results highlight the viability of non-traditional, expanding-gap, multi-aged silvicultural systems in promoting reserve tree survival in mixed softwood and hardwood forests such as the one examined in this analysis. 4.2. Exposure variables Overall mortality was negatively associated with distance to original gap edge, while distance to expanded gap edge was positively associated with wind-related mortality. These phenomena could be explained by (1) damage from logging equipment as heavily used trails tended to be on the edge of gaps (i.e. negative relationship between mortality and distance to original gap edge/closer to where equipment trails are located) (Thorpe et al., 2008), and (2) wind-stress as the gaps expanded (i.e. positive relationship between windthrow and distance to expanded gap edge/ greater exposure). Contrary to the findings in the literature, we found retention density, although not varied experimentally, was not a significant predictor of wind-related deaths or overall mortality. Similarly, incidence of wind-related death and overall mortality showed negative relationships with duration of exposure -- meaning most trees that died, did so shortly after treatment. This suggests that the dispersed retention configuration in this study provides an adequate spacing (3.75 m2 BA/ha 1 to 11.25 m2 BA/ha 1 dispersed) to minimize competition for the resources necessary to respond positively to treatment. Likewise, these multi-aged systems appear to provide adequate collective resistance to wind-stress via retention densities and gap sizes. Another factor leading to the low rate of mortality found in this study may be the prescription itself and how it developed over time. Many gaps were originally created in pre-existing natural openings with adequate advance regeneration stocking. These gaps were then subsequently expanded incrementally over successive entries. Furthermore, in areas of low advance regeneration stocking, additional overwood was retained to promote regeneration of shade-tolerant species, thereby increasing the retention density in those gaps. These harvesting practices may ‘pre-adapt’ (Rosenvald et al., 2008) or gradually increase exposure of the mature trees retained allowing for the trees to develop resistance to windthrow. Given the comparably low rate of wind-related death, this project may have refined a prescription that allows for uniquely low mortality rates of reserve trees. Continued moni- Acknowledgements We would like to thank our collaborators at the PEF, as well as, field and lab assistants Paul Szwedo, Danae Shurn, Emily Anderson, and consultants Alissa Smith, Louis Morin, and Alan White for their hard work and intellectual investment in this study. We would also like to thank McIntire-Stennis for providing the funding for this research. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.foreco.2016.12. 031. References Arseneault, J.E., Saunders, M.R., Seymour, R.S., Wagner, R.G., 2011. First decadal response to treatment in a disturbance-based silviculture experiment in Maine. For. Ecol. Manage. 262, 404–412. Bauhus, J., Puettmann, K., Messier, C., 2009. Silviculture for old-growth attributes. For. Ecol. Manage. 258, 525–537. Beurmeyer, K.R., Harrington, C.A., 2002. Fate of overstory retention trees and patterns of regeneration 12 years after clearcutting with reserve trees in southwest Washington. Western J. Appl. For. 17, 78–85. Burnham, K.P., Anderson, D.R., 2002. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. Springer-Verlag, New York, NY, USA. Busby, P.E., Adler, P., Warren, T.L., Swanson, F.J., 2006. Fates of live trees retained in forest cutting units, western Cascade Range, Oregon. Can. J. For. Res. 36, 2550– 2560. Canham, C.D., Papaik, M.J., Latty, E.F., 2001. Interspecific variation in susceptibility to wind-throw as a function of tree size and storm severity for northern temperate tree species. Can. J. For. Res. 31, 1–10. Cary, A., 1896. 3rd Annual Report of the Maine Forest Commissioner. Department of Conservation, Augusta, ME. Caspersen, J., 2006. Elevated mortality of residual trees following single-tree felling in the northern hardwood forests. Can. J. For. Res. 36, 1255–1265. Élie, J.G., Ruel, J.C., 2005. Windthrow hazard modelling in boreal forests of black spruce and jack pine. Can. J. For. Res. 35, 2655–2663. Food and Agricultural Organization of the United Nations (FAO), 2010. Global Forest Resources Assessment 2010: Main Report. FAO Forestry Paper no. 163. Franklin, J.F., 1989. Toward a new forestry. Am. For. (Nov–Dec issue), 37–44. D.R. Carter et al. / Forest Ecology and Management 389 (2017) 149–157 Franklin, J.F., Berg, D.R., Thornburgh, D.A., Tappeiner, J.C., 1997. Alternative silvicultural approaches to timber harvesting: variable retention harvest systems. In: Kohm, K.A., Franklin, J.F. (Eds.), Creating a Forestry for the 21st Century. Island Press, Washington D.C, pp. 111–139. Franklin, J.F., Mitchell, R.J., Palik, B.J., 2007. Natural Disturbance and Stand Development Principles for Ecological Forestry. Gen. Tech. Rep. NRS-19. U.S. Department of Agriculture, Forest Service, Northern Research Station, Newtown Square, PA. Fraver, S., White, A.S., 2005. Identifying growth releases in dendrochronological studies of forest disturbance. Can. J. For. Res. 35, 1648–1656. Fedrowitz, K., Koricheva, J., Baker, S.C., Lindenmayer, D.B., Palik, B., Rosenvald, R., Beese, W., Franklin, J.F., Kouki, J., Macdonald, E., Messier, C., Sverdrup-Thygeson, A., Gustafsson, L., 2014. Can retention forestry help conserve biodiversity? A meta- analysis. J. Appl. Ecol. 51, 1669–1679. Gibbons, P., Cunningham, R.B., Lindenmayer, D.B., 2008. What factors influence the collapse of trees retained on logged sites? A case-control study. For. Ecol. Manage. 255, 62–67. Gustafsson, L., Baker, S.C., Bauhus, J., Beese, W.J., Brodie, A., Kouki, J., Lindenmayer, D.B., Lõhmus, A., Pastur, G.M., Messier, C., Neyland, M., Palik, B., SverdrupThygeson, A., Volney, W.J.A., Wayne, A., Franklin, J.F., 2012. Retention forestry to maintain multifunctional forests: a world perspective. Bioscience 62, 633–645. Gustafsson, L., Kouki, J., Sverdrop-Thygeson, A., 2010. Tree retention as a conservation measure in clear-cut forests of northern Europe: a review of ecological consequences. Scandanavian J. For. Res. 25, 295–308. Hämäläinen, A., Hujo, M., Heikkala, Osmo., Junninen, K., Kouki, J., 2016. Retention tree characteristics have major influence on the post-harvest tree mortality and availability of coarse woody debris in clear-cut areas. For. Ecol. Manage. 369, 66–73. Helms, J.A., 1998. The Dictionary of Forestry. Society of American Foresters, Bethesda, MD, p. 210. Holgen, P., Søderberg, U., Hånell, B., 2003. Diameter increment in Picea abies shelterwood stands in northern Sweden. Scand. J. For. Res. 18, 163–167. Hothorn, T., Bretz, F., Westfall, P., 2008. Simultaneous inference in general parametric models. Biomet. J. 50, 356–363. Hunter Jr., M.L., 1990. Wildlife, Forests, and Forestry: Principles for Managing Forests for Biodiversity. Prentice Hall, Englewood Cliffs. Isaac, L.A., 1940. Mortality of Douglas-Fir Seed Trees on Cut-Over Lands. USDA Forest Service Research Note PNW-31, pp. 11–12. Jönsson, M.T., Fraver, S., Jonsson, B.G., Dynesius, M., Rydgård, Esseen, P-A., 2007. Eighteen years of tree mortality and structural change in an experimentally fragmented Norway spruce forest. For. Ecol. Manage. 242, 306–313. Kozlowski, T.T., Kramer, P.J., Pallardy, S.G., 1991. The Physiological Ecology of Woody Plants. Academic Press, San Diego. Lindenmayer, D.B., Franklin, J.F., Lõhmus, A., Baker, S.C., Bauhus, J., Beese, W., Brodie, A., Kiehl, B., Kouki, J., Pastur, G.M., Messier, C., Neyland, M., Palik, B., Sverdrup Thygeson, A., Volney, J., Wayne, A., Gustafsson, L., 2012b. A major shift to the retention approach for forestry can help resolve some global forest sustainability issues: retention forestry for sustainable forests. Conserv. Lett. 5, 421–431. Lindenmayer, D.B., Laurence, W.F., Franklin, J.F., 2012a. Global decline in large old trees. Science 338, 1305–1306. Lohmander, P., Helles, F., 1987. Windthrow probability as a function of stand characteristics and shelter. Scand. J. For. Res. 2, 227–238. Lorimer, C.G., White, A.S., 2003. Scale and frequency of natural disturbances in the northeastern US: implications for early successional forest habitats and regional age distributions. For. Ecol. Manage. 185, 41–64. Mahan, R.P., 1991. Circular Statistical Methods: Applications in Spatial and Temporal Performance Analysis. United States Army Research Institute for the Behavioral and Social Sciences Special Report 16. Manning, A.D., Fischer, J., Lindenmayer, D.B., 2006. Scattered trees are keystone structures– implications for conservation. Biol. Conserv. 85, 9–19. 157 Mazerolle, M.J., 2015. AICcmodavg: Model selection and multimodel inference based on (Q)AIC(c). R package version 2.0-3 http://CRAN.R-project. org/package=AICcmodavg. O’Hara, K.L., 2014. Multiaged Silviculture: Managing for Complex Structures. Oxford University Press, United Kingdom. Peterson, C.J., 2004. Within-stand variation in windthrow in southern boreal forests of Minnesota: is it predictable? Can. J. For. Res. 34, 364–375. R Core Team, 2015. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, ISBN 3-900051-07-0. URL http://www.R-project.org/. Robin, X., Natacha, T., Hainard, A., Tiberti, N., Lisacek, F., Sanchez, J.-C., Miller, M., 2011. PROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinform. 12, 77. Rosenvald, R., Lõhmus, A., Kiviste, A., 2008. Preadaptation and spatial effects on retention-tree survival in cut areas in Estonia. Can. J. For. Res. 38, 2616–2625. Ruel, J.-C., 1995. Understanding windthrow: silvicultural implications. For. Chronicle 71, 434–445. Ruel, J.-C., 2000. Factors influencing windthrow in balsam fir forests: from landscape studies to individual tree studies. For. Ecol. Manage. 135, 169–178. Saunders, M.R., Wagner, R.G., 2008. Long-term spatial and structural dynamics in Acadian mixed-woods stands managed under various silvicultural systems. Can. J. For. Res. 38, 498–517. Scott, R.E., Mitchell, S.J., 2005. Empirical modelling of windthrow risk in partially harvested stands using tree, neighborhood, and stand attributes. For. Ecol. Manage. 218, 193–209. Seymour, R.S., Hunter, M.L., 1992. New Forestry in Eastern Spruce-Fir Forests: Principles and Applications to Maine. Maine Agricultural and Forestry Experiment Station Miscellaneous Report 382. Seymour, R.S., Guldin, J.M., Marshall, D., Palik, B., 2006. Large-scale, long-term silvicultural experiments in the United States: historical overview and contemporary examples. Allgemeine Forst und Jagdzeitung 177, 104–112. Seymour, R.S., White, A.S., deMaynadier, P.G., 2002. Natural disturbance regimes in northeastern North America—evaluating silvicultural systems using natural scales and frequencies. For. Ecol. Manage. 155, 357–367. Sheil, D., Burslem, D.F.R.P., Alder, D., 1995. The interpretation and misinterpretation of mortality rate measures. J. Ecol. 83, 331–333. Smiley, E.T., Matheny, N., Clark, J., 2006. International Tree Failure Database User Manual. USDA Forest Service Version 1.1. Stephenson, N.L., Das, A.J., Condit, R., Russo, S.E., Baker, P.J., Beckman, N.G., Coomes, D.A., Lines, E.R., Morris, W.K., Rüger, N., Álvarez, E., Blundo, C., Bunyavejchewin, S., Chuyong, G., Davies, S.J., Duque, Á., Ewango, C.N., Flores, O., Franklin, J.F., Grau, H.R., Hao, Z., Harmon, M.E., Hubbell, S.P., Kenfack, D., Lin, Y., Makana, J.-R., Malizia, A., Malizia, L.R., Pabst, R.J., Pongpattananurak, N., Su, S.-H., Sun, I.-F., Tan, S., Thomas, D., van Mantgem, P.J., Wang, X., Wiser, S.K., Zavala, M.A., 2014. Rate of tree carbon accumulation increases continuously with tree size. Nature 507, 90–93. Taylor, S.L., MacLean, D.A., 2007. Spatiotemporal patterns of mortality in declining balsam fir and spruce stands. For. Ecol. Manage. 253, 188–201. Thorpe, H.C., Thomas, S.C., Caspersen, J.P., 2008. Tree mortality following partial harvests is determined by skidding proximity. Ecol. Appl. 18, 1652–1683. Urgenson, L.S., Halpern, C.B., Anderson, P.D., 2013. Level and pattern of overstory retention influence rate and forms of tree mortality in mature, coniferous forests of the Pacific Northwest, USA. For. Ecol. Manage. 308, 116–127. USDA, 2012. National core field guide. Version 6.0. Walter, S.T., Maguire, C.C., 2007. Conifer response to three silvicultural treatments in the Oregon Coast Range foothills. Can. J. For. Res. 34, 1967–1978. Zeng, H., Peltola, H., Talkkari, A., Venäläinen, A., Strandman, H., Kellomaki, S., Wang, K., 2004. Influence of clear-cutting on the risk of wind damage at forest edges. For. Ecol. Manage. 203, 77–88.