Effect magnitudes of operational-scale partial harvesting on residual tree growth and mortality of ten major tree species in Maine USA Forest and Management

Silvicultural systems based on partial harvesting (PH) have become increasingly common in areas historically dominated by clearcut-based systems in response to ecological and social concerns. Current knowledge about the effects of PH is based primarily on stand responses from controlled experiments from limited number of sites. A broader scope of inference is needed to fully understand the effects of PH when applied operationally at a landscape scale. We used 835 permanent monitoring plots from the US Forest Service Forest Inventory and Analysis network to quantify the effect magnitude (EM) of low- (5 – 40% basal area removal) and high-intensity (41 – 80% basal area removal) PH on the periodic diameter growth and mortality of residual trees of ten major tree species in Maine, USA. Our results showed that the EM and timing of statistically significant effect (p < 0.05) of the two intensities of PH varies across species. Tree diameter growth responses to high-intensity PH was rapid and prolonged for sugar maple ( Acer saccharum ) and American beech ( Fagus grandifolia ), while other species did not show significant responses for at least five years. High-intensity PH increased mortality of American beech, balsam fir ( Abies balsamea ), red spruce ( Picea rubens ), northern white cedar ( Thuja occidentalis ), and red maple ( Acer rubrum ), but not for other tree species. No relationship was found between species shade tolerance and species-specific response magnitude to PH in terms of diameter growth and tree- or stand-level (absolute or relative) mortality. This result suggests that species shade tolerance does not always determine the type or magnitude of response that a given species has to increased resource availability following PH and that other functional traits may be more important. Overall, the findings provide strong evidence that subsequent tree responses to PH are not consistent among contrasting species for at least 15 years following harvest and longer-term regional assessments are needed.


Introduction
Forest management approaches using clear cutting can often simplify vegetation composition and produce structural homogeneity, especially when followed by artificial regeneration practices (Paquette and Messier, 2010;Puettmann et al., 2015). Over the past several decades, ecosystem-based approaches have emerged as a model for forest management in parts of North America (Gauthier et al., 2009;Keeton, 2006) and elsewhere in the world (Gustafsson et al., 2012;Kuuluvainen, 2009;Lindenmayer et al., 2012). These approaches generally aim to ensure forest resilience and productivity by maintaining natural ecosystem processes and biodiversity (D'Amato and Palik, in press;Grumbine, 1994;Palik et al., 2020). In the forests of the north-eastern US and southeastern Canada, where small-scale natural disturbances generally modulate stand structural and compositional development (Angers et al., 2005;Fraver et al., 2009;Frelich, 2002), partial harvesting (PH) of various intensities and spatial configurations have been suggested to emulate natural dynamics, thus achieving ecological objectives with acceptable economic outcomes (Raymond et al., 2009;Seymour et al., 2002).
Our current understanding about the effects of PH has generally been based on limited number of controlled experiments where trees are carefully selected to create theoretical post-harvest stand patterns (Bataineh et al., 2013;Leak, 2005;Prévost and Gauthier, 2013;Raymond and Bédard, 2017). This may create uncertainties for wider application of our current knowledge, especially when PH is applied operationally at a landscape scale. In such harvest operations, the target volume is generally extracted from a larger area compared to controlled experiments under more varied operational constraints and conditions. Although some studies have quantified the effects of operational-scale PH on natural regeneration (Bose et al., 2020;Bose et al., 2016;Guay-Picard et al., 2015;Montoro Girona et al., 2018), there are few studies on growth and mortality responses of residual trees and stands following PH. The knowledge gap is important considering the large-scale application of PH in this region and the lack of experiments that encompass the degree of complexities associated with operational-scale PH. Further, the existing knowledge on growth and mortality responses to PH is limited to just a few tree species even at the experimental level (e. g., Baral et al., 2016;Carter et al., 2017a;Forget et al., 2007;Jones and Thomas, 2004).
Forested landscapes across north-eastern USA are dominated by mixed-species ecosystems (Seymour, 1995), where species-specific differences in survival and growth play a key role in community composition (Gravel et al., 2010;Kobe et al., 1995) and the trade-off between high-light growth and low-light survival has often been used to explain species-specific responses in these ecosystems (Canham, 1989;Pacala et al., 1996). The resource conditions and exposure generated by PH are expected to differentially affect tree species based on light requirements and growth forms. Shade-intolerant species are generally expected to show greater levels of increased growth compared to shade-tolerant species in response to increased light levels (Valladares and Niinemets, 2008). Growth and mortality responses may also vary across tree morphology in which shade-tolerant conifers may have higher crown depth compared to shade-intolerant deciduous species (Cole and Lorimer, 1994;Weiskittel et al., 2009). In addition, species associated with shallow root systems are generally more vulnerable to windthrow induced tree mortality compared to species with deeper root systems (Canham et al., 2001) and an increased occurrence of windthrow is often reported after PH (Bose et al., 2014c). This suggests a range of dynamic yet interactive factors may control species response to PH.
In addition to species-specific responses, growth and mortality responses of residual trees after PH are often associated with neighbourhood conditions created by the harvest operations (Cortini and Comeau, 2020), pre-harvest growing conditions (Bose et al., 2014a;Montoro Girona et al., 2017;Thorpe et al., 2007), tree and stand characteristics prior to harvest (Anning and McCarthy, 2013;Baral et al., 2016), and post-harvest competitive status (Bose et al., 2014a;Jones et al., 2009). The actual effect magnitude (EM) of PH can be difficult to quantify due to lags between the occurrence of PH and subsequent tree responses (Jones et al., 2009;Thorpe et al., 2007). This response is exacerbated by differences in physiological acclimation potential among tree species to cope with altered growing conditions such as changes in growth dominance (Lemire et al., 2020), inter-tree competition (Anning and McCarthy, 2013), wind penetration in the residual stands (Ruel, 2000), and damage from harvesting operations (Bose et al., 2014c). Variation in harvest intensity, method, and timing make the situation even more complex, making it difficult to compare effects across PH treatments in time and space. When addressing and trying to understand speciesspecific responses to PH, it is crucial to quantify the actual EM of PH and how it varies with harvest intensity in space and time, including any potential lag effects.
In this study, we used 835 permanent sample plots from the Forest Inventory and Analysis (FIA) network of the US Forest Service across the State of Maine, USA to: (1) quantify the EM of low-intensity (5-40% of overstory basal area removal) and high-intensity (41-80% of basal area removal) PH on the subsequent diameter growth and mortality responses of ten major tree species, and (2) assess whether EM of preharvest tree size and stand basal area (an indicator of suppression history) were significant across species and across unharvested and two types of partially harvested stands. We hypothesized that the EM of two types of PH on tree growth and mortality will decrease with species shade tolerance and will follow a gradient suggested by Niinemets and Valladares (2006): balsam fir > eastern hemlock > sugar maple > American beech > red spruce > northern white cedar > red maple > white pine > yellow birch > paper birch. In addition, we expected a negative effect of pre-harvest tree size and of pre-harvest overstory basal area on subsequent tree growth, but a positive effect on tree mortality. Finally, it was expected that the EM will vary among tree species based on timing and intensity of the PH, as well as the shade tolerance ranking outlined above.

Study sites
Our study sites are in the Maine, USA covering a large and diverse geographic region (43 • 08 ′ N-47 • 43 ′ N, 66 • 99 ′ W-71 • 07 ′ W; Fig. 1). Maine has a cool humid continental climate with a mean annual precipitation of 901-1501 mm and mean annual temperature of 1.2-8.4 • C. The region has an approximately 6-month growing season extending from April to September, but the length of the growing season varies from north to south. Soils are generally classified as podzols, and soil fertility varies between forests with softwoods or hardwoods, where the latter has higher productivity than the former (Bose et al., 2016;Seymour, 1995).
Tree species composition in these forests is dominated by shadetolerant species because of long fire-return intervals and absence of clearcutting (Seymour, 1995). American beech (Fagus grandifolia Ehrh.), Fig. 1. The sample plot locations across Maine, USA. In the two partially harvested plots, harvesting was conducted during 1999-2003 and not reharvested during the examined period. The background map was downloaded using the function map_data from ggmap package in R (Kahle and Wickham, 2016), while the ggplot function of the package ggplot2 in R was used to create the map (Wickham, 2016). sugar maple (Acer saccharum Marsh.), and red maple (Acer rubrum L.) are the dominant shade-tolerant deciduous species, while balsam fir (Abies balsamea (L.) Mill.), red spruce (Picea rubens Sarg.), eastern hemlock (Tsuga canadensis (L.) Carriere), and northern white cedar (Thuja occidentalis L) are common shade-tolerant conifers. Among shade-intolerant or mid-tolerant tree species eastern white pine (Pinus strobus L.) is the dominant conifer, while yellow birch (Betula alleghaniensis Britton) and paper birch (Betula papyrifera Marsh.) are the most frequently found deciduous species (Bigelow and Canham, 2002;Olson and Wagner, 2011).

Data
FIA is a nationwide forest inventory that uses a consistent sampling protocol across the country (http://www.fia.fs.fed.us). FIA sampling protocol consists of four points arranged in a cluster, where point 1 is the center, and points 2-4 are located 36.58 m from point 1 at azimuths of 0, 120, and 240 • , respectively. Each point serves as the center of a 7.3 m fixed-radius subplot where all trees with a diameter at breast height (dbh; 1.37 m) ≥ 12.7 cm are measured for dbh by species and characterized as "live" of "dead". (https://www.fia.fs.fed.us/library/databas e-documentation/).
In the context of this study, we characterized PH as any harvest that removed 5-80% of overstory (trees ≥ 12.7 cm at DBH) basal area in a forest stand (Bose et al., 2020). The harvest intensity was quantified from the ratio of preharvest to the postharvest overstory basal area of a stand. We then selected all FIA plots that received PH since 1999 and had a preharvest measurement, at least three postharvest measurements, and were not re-harvested during 1999-2016. A total of 424 plots met these criteria. We then randomly selected an equal number plots (n = 424) that were not harvested from 1999 to 2016. We excluded 13 of the 424 unharvested plots due to very low overstory basal area (<1 m 2 ha − 1 ). As we were assessing the post-harvest diameter growth and mortality responses to PH treatment, it was important to consider only those unharvested plots whose stand structure and composition was not significantly different compared to the pre-harvest stand structure and composition of partially harvested stands (Bose et al., 2020). We characterized 5-40% and 41-80% basal area removal as low-intensity and high-intensity PH, respectively, based on the median of harvesting intensity distribution received by our sample plots (see details in Bose et al., 2020). From the 424 original partially harvested plots, we had 171 and 253 low-intensity and high-intensity PH plots available for study, respectively.

Data analysis
We quantified EM of two types of PH on tree-level diameter (DBH) growth and tree-and stand-level mortality of ten commonly occurring tree species in Maine, USA. These include yellow birch, paper birch, eastern white pine, balsam fir, red spruce, eastern hemlock, northern white cedar, American beech, sugar maple, and red maple. For diameter growth, we quantified the absolute diameter growth (ADG) and relative diameter growth (RDG). The ADG and RDG was defined by the absolute and relative diameter (i.e., diameter at breast height) growth performance between two successive measurement periods, and was quantified using the following equations: where i = measurement period 1 ……… n, and j = measurement period i + 1.…… n + 1. Preliminary analysis showed similar trends with ADG and RDG values. We therefore, retained RDG in the analysis to better assess the relative differences among various tree species and PH treatments being examined, as well as to capture the dynamics following PH (Bose et al., 2018a).
For mortality of each species, we quantified (i) probability of tree mortality (i.e., live or dead), (ii) absolute stand-level tree mortality: total number of trees that were living in the first measurement but dead in the following measurement, and (iii) relative stand-level tree mortality: percentage of dead trees relative to the total number of trees (i.e., live and dead) across three measurement periods. Tree-level and stand-level (absolute and relative) mortality were quantified separately for each species. We considered three post-harvest measurement periods, 1-5 years, 6-10 years, and 11-15 years since the onset of the treatment.
For examining whether the EM of pre-harvest tree-and stand-level factors on diameter growth and mortality were statistically significant across unharvested and partially harvested stands, we considered preharvest tree dbh and total stand basal area of live trees. For standlevel mortality (absolute and relative) analysis, we considered average tree dbh as a proxy of tree size. In addition, we quantified the ratio between conspecific basal area and total stand basal area (RCT), which measures the relative importance of intra-specific against inter-specific competition for explaining tree growth and mortality of each tree species. However, preliminary analysis showed insignificant effects (p > 0.05) of RCT on tree growth and mortality of all species and the variable was excluded from the analysis.

Statistical analysis
We first tested the differences in pre-treatment stand-and tree-level characteristics between unharvested and PH plots. For this analysis, we performed linear mixed-effect modelling using the function lme in the nlme package in R (version 3.2.5) (Pinheiro et al., 2014;R Development Core Team, 2018), where harvesting treatments (i.e., unharvested and partially harvested) were treated as fixed effects variable. For standlevel variable, we considered subplots nested within plots, plots nested within counties, and counties nested within management units as random effects, while trees nested within subplots, subplots nested within plots, plots nested within counties, and counties nested within management units as random effects for tree-level variable.
For the two research questions, we measured the EM by the slope of the mixed-effect model. This is a common statistical approach for measuring EM of a treatment (e.g., Forrester, 2019;Komatsu et al., 2019;Vitasse et al., 2019). For the first research question, each response variable (i.e., relative diameter growth, probability of tree mortality, absolute stand-level mortality, and relative stand-level mortality) was modelled as a function of treatments (two levels: unharvested vs. lowintensity PH or unharvested vs. high-intensity PH). In addition to fixed effects (i.e., effects from PH), the associated random effects were incorporated into the model (see above). The modelling was performed separately for each period and for each species. The linear mixed-effect modelling (Zuur et al., 2009) was performed for continuous data (i.e., diameter growth) using the lme function of the R package nlme (Pinheiro et al., 2014;R Development Core Team, 2018). The generalized linearmixed effect modelling was targeted at count (i.e., absolute mortality), percentage (i.e., relative mortality), and presence-absence (i.e., ADG = Tree diameter at measurement period j − Tree diameter at measurement period i (1) RDG = Tree diameter at measurement period j − Tree diameter at measurement period i Tree diameter at measurement period i probability of tree mortality) data, which were performed using the glmer function of the R package lme4 (Bates et al., 2017). In glmer modelling, the binomial family was used for percentage and presenceabsence data, while Poisson family was considered for the count data. We visually verified the assumptions of normality and variance homogeneity of the residuals. We used square root and/or log transformation when needed for continuous data. The relative mortality analysis was not performed for sugar maple and eastern hemlock to avoid heteroscedasticity.
For the second research question, the response variables (i.e., diameter growth, probability of tree mortality, and relative stand-level mortality) were modelled as a function of pre-harvest tree size and total stand basal area. The modelling was performed separately for each treatment (i.e., unharvested, low-intensity PH, and high-intensity PH), species, and period. For stand-level mortality, we only considered the relative mortality as the results from first research question showed that the absolute mortality did not capture the actual differences between unharvested and partially harvested stands. Similar to the analysis for the first research question, we used linear mixed-effect modelling for continuous data (i.e., relative diameter growth) and generalized linearmixed effect modelling for percentage (i.e., relative stand-level mortality) and presence-absence (i.e., probability of tree mortality) data.

Pre-treatment characteristics of studied trees and stands
The mean dbh of live trees was similar in unharvested and in the PH stands before treatment (p > 0.590) except for red spruce and eastern hemlock ( Table 1). Trees of those two species were slightly larger in high-intensity PH stands (mean ± SD, red spruce: 22.2 ± 6.9 cm and eastern hemlock: 25.9 ± 9.5 cm) compared to unharvested stands (mean ± SD, red spruce: 21.0 ± 7.0 cm and eastern hemlock: 24.2 ± 9.7 cm) ( Table 1). Conspecific overstory basal area (CBA) was similar in unharvested and PH stands for balsam fir and yellow birch, but not for other species where CBA was slightly greater in high-intensity PH stands for American beech and red maple yet lower for sugar maple compared to unharvested stands (Table 1). In low-intensity PH stands, the CBA was slightly higher for American beech, red maple, paper birch, red spruce, and white pine, but lower for sugar maple, eastern hemlock, and northern white cedar compared to unharvested stands (Table 1).

Tree relative diameter growth (RDG)
For the entire 15-year post-harvest period, our descriptive analyses showed that high-intensity PH resulted in the highest RDG in sugar maple and lowest in northern white cedar when compared to unharvested stands, while low-intensity PH resulted in the highest RDG in sugar maple and lowest in red spruce (Table 2).
Our mixed-effect model analyses showed that all species displayed positive RDG responses to PH, however, the response magnitude and timing of response varied among species (Fig. 2). RDG responses to highintensity PH was rapid from sugar maple and American beech, while other species took more than five years to detect a statistically significant response. The RDG response of balsam fir, paper birch, yellow birch, eastern hemlock, red maple, white pine, and red spruce to highintensity PH was significant from 6-10 to 11-15 years since the onset of the treatment, however, statistically significant growth response was only visible at 11-15 years since the onset of the treatment for northern white cedar (Fig. 2). Low-intensity PH did not produce any change in RDG of American beech, northern white cedar, paper birch, yellow birch and red spruce, but observed in balsam fir, eastern hemlock, red maple, white pine, and sugar maple. Sugar maple is the only species that displayed a rapid and persistent RDG response to the low-intensity PH treatment (Fig. 2).
The effect of time was not statistically tested. Our analysis showed whether the EM was significant or not across measurement periods and the size of the EM (i.e., slope) and associated standard errors indicated how the EM evolved over time since harvesting (Fig. 2). Based on these results, the response magnitude of RDG to low-intensity PH was smaller compared to the response magnitude for high-intensity PH irrespective of species and time since treatment (Fig. 2). The response magnitude to high-intensity PH increased over time since treatment for white pine, red spruce, and northern white cedar, while it decreased for sugar maple. The patterns of response magnitude with time since treatment was curvilinear for American beech, balsam fir, eastern hemlock, paper birch, yellow birch, and red maple (Fig. 2).

Table 1
Total number of sample plots and trees and pre-harvest stand and tree characteristics of ten species across three partial harvesting treatments. The underlined values indicate a significant difference (p < 0.05) between unharvested and that partial harvesting treatment. Note. Low: 5-40% basal area removal and High: 41-80% basal area removal.

Tree-level mortality
High-intensity PH increased the probability of individual-level tree mortality of balsam fir, paper birch, and red spruce, but not of other species. These significantly higher mortality probabilities in highintensity PH stands occurred during 1-5 as well as 6-10 years since harvesting for balsam fir and 6-10 years since harvesting for paper birch and red spruce. Low-intensity PH increased the probability of tree mortality for paper birch and yellow birch, but not for other species (Fig. 3).

Stand-level mortality
For the entire 15-year post-harvest period, our descriptive analyses showed that when compared to unharvested stands high-intensity as well as low-intensity PH resulted the highest absolute mortality of red spruce and northern white cedar yet was lowest for American beech and sugar maple ( Table 2). The highest relative mortality was observed in northern white cedar and red spruce in both high-intensity and lowintensity PH stands, while the lowest relative mortality was observed in red maple species (Table 2).
Our mixed-effect modelling analysis indicated that PH did not make a statistical difference compared to unharvested stands in terms of absolute tree mortality of balsam fir, northern white cedar, eastern hemlock, paper birch, yellow birch, sugar maple, and white pine irrespective of harvesting intensities and time since treatment (Fig. 4). Total number of dead trees of American beech decreased in high-intensity and lowintensity treatments at 6-10 as well as at 11-15 years since the onset of the treatment. In addition, low-intensity PH decreased the number of dead trees of red maple but increased the number of dead trees of red spruce at 11-15 years since the onset of the treatment (Fig. 4).
Although high-intensity PH decreased the absolute mortality (i.e., total number of dead trees) of American beech, it increased the relative mortality (i.e., percentage of dead beech trees to total number of beech trees) during the initial 1-5 years after harvest (Fig. 5). High-intensity PH also increased the relative mortality of balsam fir, red spruce, northern white cedar, and red maple, but not the other species. In contrast to high-intensity PH, low-intensity PH did not make any significant difference in terms of relative mortality across all species and measurement periods (Fig. 5).

Effects of pre-harvest tree size and stand basal area on growth and mortality
Pre-harvest tree size (dbh) and overstory basal area were negatively correlated with post-harvest diameter growth across species. However, based on the size of the slope and associated standard errors, our results showed that the magnitude of the negative effect was higher in PH stands compared to unharvested stands ( Fig. 6 and Table SM1). Unlike unharvested stands, the magnitude of the negative effect of tree dbh increased over time in high-intensity PH stands for balsam fir, red spruce, eastern hemlock, and white pine (not statistically tested). Overall, tree dbh had a stronger negative effect on post-harvest relative diameter growth compared to overstory basal area irrespective of PH treatments and species (Fig. 6 and Table SM1).
Pre-harvest tree size and overstory basal area were mostly found insignificant in explaining the tree-level mortality (Table SM2). Among the few significant effects, pre-harvest overstory basal area was positively related to the tree-level mortality of balsam fir in unharvested and PH stands, where the EM of pre-harvest overstory basal area was higher in high-intensity PH when compared to unharvested stands (Table SM2).
Pre-harvest tree size and overstory basal area had a minor effect (often not statistically significant; p > 0.05) on relative mortality across PH stands and species (Table SM3). We detected no effect of either preharvest factors on relative mortality of northern white cedar or white pine. However, both factors had a significant positive effect on the relative mortality of American beech in high-intensity PH stands at 6-10 years since the treatment (Table SM3). In contrast to American beech, overstory basal area (but not tree size) had a significant positive effect on relative mortality of balsam fir, paper birch, red maple, and red spruce in unharvested stands but not in high-intensity PH stands (Table SM3). Our descriptive statistics showed a greater difference across species in terms of their mortality responses to overstory basal area when compared to tree size irrespective of PH treatments ( Fig  SM3).

Discussion
Our analysis revealed interspecific variation in overstory tree diameter growth and mortality among the major tree species in Maine USA with and without partial harvesting (PH). We showed direction Table 2 Magnitude of change (%) in high-intensity and low-intensity partial harvested (PH) stands compared to unharvested stands in terms of relative diameter growth and stand-level absolute and relative tree mortality of ten tree species of Maine, USA. Absolute mortality: total number of dead trees per hectare and relative mortality: total number of dead trees relative to total number of dead and live trees per hectare. See equation (2) in Method section for relative diameter growth. Growth and mortality data represent the average of three measurements which were conducted at five-year interval over a 15-year post-harvest period. (positive or negative), magnitude, and temporal trajectories of responses in individual tree diameter growth as well as tree-and stand-level mortality to both low-(5-40% basal area removal) and high-intensity (41-80% basal area removal) PH. We had expected that by increasing overall light availability PH would increase diameter growth across all tree species, but the magnitude of the growth enhancement would be negatively correlated with the shade tolerance for the tree species. However, we detected no relationship between the EM of PH intensity and shade tolerance (Fig. 7). We also expected that PH will increase the mortality of residual trees because of potential post-harvest shock and increased wind penetration into residual stands. Our results showed that only high-intensity PH caused increased mortality of residual trees Effect magnitudes (i.e., slope of the linear mixed-effect models) of two partial harvesting treatments on tree-level relative diameter growth (RDG) over a 15year post-treatment period. The effect sizes represent the slope of the mixed-effect models where RDG is modelled as a function of unharvested vs low intensity and unharvested vs high intensity partial harvestings. The error bars represent the mean ± standard error (sample size for each species is presented in Table 1) and the * sign indicates the statistically significant (p < 0.05) difference between unharvested and that partial harvesting treatment.
compared to unharvested stands, but that greater mortality in highintensity PH stands did not occur for all species and was mostly limited to initial 5-10 years since harvesting operations.

Diameter growth
We found that the diameter growth of shade-tolerant sugar maple and American beech increased immediately following the high-intensity PH, while mid-tolerant white pine and yellow birch and intolerant paper birch needed more than five years to display statistically significant diameter growth increases following high-intensity PH (Fig. 2). In addition, the EM of high-intensity as well as low-intensity PH was larger for shade-tolerant sugar maple compared to shade-intolerant paper birch or mid-tolerant yellow birch and white pine. These findings may Fig. 3. Effect magnitudes (i.e., slope of the generalized linear mixed-effect models) of two partial harvesting treatments on tree-level mortality over a 15-year posttreatment period. The error bars represent the mean ± standard error (sample size for each species is presented in Table 1) and the * sign indicates the statistically significant (p < 0.05) difference between unharvested and that partial harvesting treatment.
have a relationship with trees position within the canopy prior to harvesting (Wright et al., 1998). Shade-intolerant species may have a greater stress from long-term suppression compared to shade-tolerant species (Kobe et al., 1995), which may negatively affect their response to increased light availability (Baral et al., 2016). Another possible explanation is that early successional shade-intolerant species may already be in the canopy with high crown exposure (Smith et al., 2016), therefore, their growth might already be at maximum potential and thus, further harvesting did not increase their diameter growth. In addition, PH may not sufficiently increase the light needed for shadeintolerant or mid-tolerant species to produce a substantial diameter growth response. For example, our results showed significant negative effects of tree dbh on diameter growth in two types of PH stands for most of the shade-tolerant species (e.g., beech, sugar maple, balsam fir) but not for the shade-intolerant species (e.g., white pine or paper birch) (Table SM1). This indicates that the diameter growth of shade-tolerant species was more related to their size compared to shade-intolerant species consistent with findings of Kuehne et al. (2020). This may also Fig. 4. Effect magnitudes (i.e., slope of the generalized linear mixed-effect models) of two partial harvesting treatments on stand-level absolute tree mortality over a 15-year post-treatment period. The error bars represent the mean ± standard error (sample size for each species is presented in Table 1) and the * sign indicates the statistically significant (p < 0.05) difference between unharvested and that partial harvesting treatment. The missing data (e.g., eastern hemlock responses in 11-15 years since harvesting) indicates that the data was not suitable (e.g., too many zeros) for a statistical analysis.
indicate a greater responsiveness of shade-tolerant species to PH treatments compared to shade-intolerant species.
Our results showed that the EM of high-intensity PH on diameter growth was larger than low-intensity PH irrespective of species. This indicates that higher light availability and lower post-harvest inter-tree competition resulted from high-intensity PH was crucial for improved diameter growth of all tree species. The positive effect of PH intensity on tree diameter growth has also been reported by other studies conducted in this region (Carter et al., 2017a;Prevost et al., 2010) as well as in other neighbouring regions (Bose et al., 2014a;Thorpe et al., 2007;Xing et al., 2018).
Existing literature suggests a strong role of pre-harvest tree size and overstory basal area on the post-harvest tree growth responses where both factors are often found to be negatively correlated with tree growth (Baral et al., 2016;Carter et al., 2017a;Montoro Girona et al., 2016;Moussaoui et al., 2020). In this study, we anticipated negative effects from those two factors on tree growth and compared the EM of the negative effects across ten species as well as across unharvested and two Effect magnitudes (i.e., slope of the generalized linear mixed-effect models) of two partial harvesting treatments on stand-level relative tree mortality over a 15-year post-treatment period. The error bars represent the mean ± standard error (sample size for each species is presented in Table 1) and the * sign indicates the statistically significant (p < 0.05) difference between unharvested and that partial harvesting treatment. The missing data (e.g., red spruce responses in 11-15 years since harvesting) indicates that the data was not suitable (e.g., too many zeros) for a statistical analysis.
types of PH stands. Although both of those factors had a significant negative effect on tree growth of a number of species in unharvested stands, the differences among species as well as over the 15-year monitoring period was lower in unharvested stands compared to the two types of PH stands (Fig. 6). We also observed stronger negative effects on tree growth of both of those factors in PH stands compared to unharvested stands (Table SM1), which indicates a greater role of preharvest tree and stand characteristics for determining post-harvest tree growth in PH stands compared to unharvested stands.
Our descriptive statistics showed that the EM of pre-harvest tree size varies largely among species compared to the EM of pre-harvest overstory basal area (Fig. 6). We did not observe statistically significant (p < 0.05) effect of pre-harvest overstory basal area (an indicator of overall levels of resource competition) on diameter growth of most of our tree species in high-intensity PH stands with the exception of northern white cedar and paper birch (Table SM1). However, pre-harvest tree size was strongly related to diameter growth of almost all species in highintensity PH stands, where the strongest negative effect was observed in American beech and balsam fir yet weakest in red spruce (Table SM1 and Fig. 6).

Mortality
Post-harvest tree mortality is a major concern associated with PH (e. g., Bose et al., 2014c;Montoro Girona et al., 2019;Thorpe et al., 2008). Several studies reported that stand density reduction by PH can increase wind penetration and generate harvesting-shock in residual stands and can physically damage residual trees through uprooting or stem breakage (e.g., Carter et al., 2017b;Coates, 1997;Ruel and Gardiner, 2019;Ruel et al., 2001). We found high-intensity PH (41-80% basal area removal) increased probability of individual-level mortality of balsam fir, red spruce, and paper birch but not of the other species examined Fig. 6. Effect magnitudes (i.e., slope of the linear mixed-effect models) of pre-harvest overstory basal area and tree dbh on diameter growth of ten tree species in unharvested, low-intensity partially harvested, and high-intensity partially harvested stands over a 15-year post treatment period. The error bars represent the mean ± standard error (sample size for each species is presented in Table 1). See appendix-1 for statistical significance. (Fig. 3). However, the statistically significant increase of tree-level mortality in high-intensity PH stands compared to unharvested stands of those three species did not persist during the later post-harvest periods (i.e., 6-10 and 11-15 years since harvesting) (Fig. 3) indicating residual trees acclimated to the new conditions relatively rapidly.
Our stand-level analyses (absolute and relative stand-level mortality) provided further insights on post-harvest tree mortality. We showed that although total number of dead trees per hectare (i.e., absolute mortality) did not increase (rather decreased) after high or low-intensity PH, the relative mortality (i.e., percentage of dead trees to the total number of trees) increased after high-intensity PH of several species including American beech, balsam fir, northern white cedar, red maple, and red spruce (Fig. 5). Compared to unharvested stands, we observed higher stand-level mortality in high-intensity PH stands than in low-intensity PH stands (Table 2). Several other studies that separated post-harvest tree mortality by high-intensity and low-intensity PH also observed greater mortality in high-intensity PH stands (see synthesis in Bose et al., 2014c). These studies characterized that high-intensity PH generally Fig. 7. The relationship between species shade tolerance and effect magnitude of two types of partial harvesting treatments (low-intensity and high-intensity) in terms of growth, tree-level mortality, absolute stand-level mortality, and relative stand-level mortality for ten tree species. The effect magnitudes represent the average of three values (i.e., three post-harvest measurement periods). produced greater harvesting shock and logging damage to residual trees and created greater opening for wind penetration into residual stands, therefore, higher windthrow induced tree mortality (e.g., Harvey and Brais, 2007;Lavoie et al., 2012;Montoro Girona et al., 2019;Ruel and Gardiner, 2019;Ruel et al., 2003;Solarik et al., 2012). Another explanation of greater mortality in our high-intensity PH stands could be related to removal of higher percentage of merchantable or healthy trees (Harvey and Brais, 2007), which resulted in the dominance of smallsized or weaken trees in post-harvest stands that are more vulnerable to post-harvest wind or insect damages (Bose et al., 2014b;Solarik et al., 2012).
We observed strong differences among tree species in terms of their post-harvest tree-and stand-level mortality (Table 2 and Figs. 3-5). Overall, conifer species had higher mortality in high-intensity PH stands compared to deciduous species (Table 2). However, among conifer species, eastern hemlock and white pine were not significantly affected by either of the two types of PH treatments (Figs. 3-5). We observed significantly higher mortality of balsam fir, red spruce, and northern white cedar in high-intensity PH stands, where most of the dead trees of northern white cedar and balsam fir were smaller than 20 cm at dbh (Fig. SM2). These smaller trees of extremely shade tolerant species such as balsam fir and northern white cedar may have suffered from thinning shock, which could have been triggered by high-intensity PH. The higher mortality of red spruce in high-intensity PH stands is probably due to windthrow damages because of their shallower root systems, which was also observed by other studies conducted in this region (Canham et al., 2001). Moreover, balsam fir and red spruce can also be affected by spruce budworm (Choristoneura fumiferana Clem.) outbreaks which is often found to be more detrimental in post-windthrow stands (Girard et al., 2014). Among all species, sugar maple had the lowest mortality in the two PH stands (Table 2 and  . Sugar maple displayed the highest diameter growth responses among the ten species we examined (Table 2 and Fig. 2) indicating that the species did not suffer at all in the two types of PH stands examined in this analysis. It is important to mention that sugar maple is primarily found on welldrained, nutrient rich sites, which are less prone to windthrow and contain high levels of resources to support growth.
Based on the findings of earlier studies, we anticipated that the preharvest overstory basal area would have a positive relationship with tree mortality because of suppression effect (Bose et al., 2018b;Das et al., 2016), while pre-harvest tree size would have a negative relationship with mortality because of general increase in self-thinning after PH (Bose et al., 2014b;Montoro Girona et al., 2019). We therefore wanted to compare the EM of pre-harvest overstory basal area and tree size on mortality of all species. Our tree-and stand-level analysis detected almost similar mortality responses to pre-harvest tree size and overstory basal area across species (Table SM2 and SM3). Our results showed comparatively minor influences of pre-harvest tree size and overstory basal area on mortality than the influences they had on diameter growth (Table SM1, SM2, and SM3). The direction of the mortality responses (i. e., positive or negative) are not consistent to our expectation (i.e., positive effect of overstory basal area and negative effect of tree size). This indicates a greater complexity in tree mortality when compared to diameter growth responses to PH treatments. Similar findings have also been reported by other studies (e.g., Bose et al., 2018b;Carter et al., 2017b;Montoro Girona et al., 2019;Ruel and Gardiner, 2019;Thorpe et al., 2008).
Generally, post-harvest tree mortality is difficult to predict because the probability of mortality can be related to a large number of factors in addition to tree size and overstory basal area. This may include the trees' location within a stand and its proximity to logging trails (Thorpe et al., 2008), while stand neighbourhood conditions such as a clearcut areas can intensify wind penetration into residual stands (Ruel and Gardiner, 2019), and harvesting methods such as group retention can result in lower mortality compared to dispersed retention of residual trees (Solarik et al., 2012). Despite the robust and systematic design of FIA experimental plots, we cannot exclude the potential effects of logging trails and plot neighbourhood conditions (Thorpe et al., 2008), which will be a subject of subsequent analysis. Our study also lacks direct measurements of climate (precipitation, temperature, and solar radiation) and geographic factors (soil, exposure to wind, and depth of water table), which generally modulate wind penetration into residual stands and determine levels of windthrow-induced tree mortality (Ruel and Gardiner, 2019). In addition, the FIA plots were measured periodically every five years, a finer resolution of the data such as those generated by dendrochronological approaches could provide more insightful understanding on growth and mortality response to PH. Our study provides quantitative knowledge related to species specific responses to low-and high-intensity PH. Using this framework, future studies could examine if PH shift the growth dominance from one species to another or from one size class to another (Forrester, 2019;Lemire et al., 2020). In addition, PH could change the competition dynamics in residual stands. Therefore, incorporating neighbourhood competition (Canham et al., 2004) in the analysis could improve the strength of the growth and mortality analysis.

Management implications
The current understanding on the effect of PH on tree growth and mortality is based primarily on experimental studies with relatively consistent site conditions and where trees are carefully selected to create homogeneous harvesting patterns (e.g., Bose et al., 2014b;Carter et al., 2017a;Leak, 2005;Ruel and Gardiner, 2019;Thorpe et al., 2008). This limits the wider application of the knowledge from those experimental studies, especially when PH is applied operationally at the landscape scale. Using 835 permanent monitoring plots of USDA Forest Service FIA, our study provides valuable quantitative knowledge on postharvest residual tree growth and mortality of ten major yet ecologically contrasting tree species of north-eastern USA and south-eastern Canada. Our analyses showed that the EM of the low-intensity PH (5-40% of basal area removal) and high-intensity PH (41-80% of basal area removal) on tree growth and mortality varied strongly across species. Among the ten species examined in this study, sugar maple showed the highest diameter growth and lowest mortality responses to lowintensity as well as high-intensity PH. This result may provide valuable feedback to forest managers, researchers, and ecologists who have been working on PH based silvicultural activities over the past several decades for sugar maple forests.
We observed lower growth and higher mortality responses in red spruce and northern white cedar, which may indicate that these two species are not well-suited for PH treatments irrespective of PH intensity. We found no relationship between species shade tolerance and species-specific response magnitude to PH in terms of diameter growth or tree-and stand-level (absolute or relative) mortality. Our results might suggest that species shade tolerance should not always be the primary attribute for determining the response magnitude of a species to any PH treatment as the other factors related to species, tree, and site conditions dynamically modulate the response patterns. Our results clearly indicate that the residual tree mortality in two types of PH stands was fairly limited across species (Table 2) and therefore, foresters should not avoid using PH for the potential that it would increase tree mortality. Overall, our analysis highlights the complexity of tree-and stand-level response to PH that suggests the need for continued long-term evaluation of these treatments, particularly at the regional scale given the broad application of these methods currently.

CRediT authorship contribution statement
A.K.B. Conceptualized the ideas, designed Methodology, analysed the Data and led the Writing of the manuscript; R.G.W, A.R.W., and A.W. D contributed in Data analysis and Writing of the manuscript.