The Impact of Variable Retention Harvesting on Growth and Carbon Sequestration of a Red Pine (Pinus resinosa Ait.) Plantation Forest in Southern Ontario, Canada

As atmospheric carbon dioxide concentrations continue to rise and global temperatures increase, there is growing concern about the sustainability, health, and carbon sequestration potential of forest ecosystems. Variable retention harvesting (VRH) has been suggested to be a potential method to increase forest biodiversity, growth, and carbon (C) sequestration. A field trial was established in an 88-year-old red pine (Pinus resinosa Ait.) plantation in southern Ontario, Canada, using a completely randomized design to examine the response of tree productivity and other forest values to five harvesting treatments: 33% aggregate retention (33A), 55% aggregate retention (55A), 33% dispersed retention (33D), and 55% dispersed retention (55D) in comparison to an unharvested control (CN). In this study, we explored the impacts of VRH on aboveground stem radial growth and annual C increment. Standard dendrochronological methods and allometric equations were used to quantify tree- and stand-level treatment effects during a five-year pre-harvest (2009–2013) and post-harvest (2014–2018) period. Tree-level growth and C increment were increased by the dispersed retention pattern regardless of retention level. At the stand level, the total C increment was highest at greater retention levels and did not vary with retention pattern. These results suggest that the choice of retention level and pattern can have a large influence on management objectives as they relate to timber production, climate change adaptation, and/or climate change mitigation.


INTRODUCTION
Atmospheric carbon dioxide (CO 2 ) concentrations have increased from approximately 280 ppm in the mid-1700s to >411 ppm in 2020 due to fossil fuel combustion and other anthropogenic activities (Dlugokencky and Tans, 2018;Lindsey, 2018). Depending on global efforts to reduce emissions in the coming decades, atmospheric CO 2 concentrations are predicted to further increase to as high as 1000 ppm by 2100 (IPCC, 2014;Lindsey, 2018). Forests play an important role in the global carbon (C) cycle and can be managed to offset a proportion of these emissions by sequestering atmospheric C via photosynthesis into biomolecular C sinks (Unwin and Kriedemann, 2000;Pan et al., 2011). However, climate change and the associated projected increase in frequency, extent, and severity of natural disturbances represent a serious threat to the ecological functioning of forest ecosystems and their capacity to act as significant C sinks (Gauthier et al., 2014;Peterson et al., 2014).
A variety of silvicultural approaches can be used to manage forests to minimize the adverse effects of climate change stresses on productivity (Puettmann, 2011;Peterson et al., 2014;Williamson et al., 2019) and maintain or enhance the rate of C sequestration (Colombo et al., 2005;Peterson et al., 2014). Reduction of stand density by thinning and use of partial harvesting systems can reduce resource competition and increase the growth and C sequestration rates of residual trees and stands (Davis et al., 2009;Dwyer et al., 2010;Nunery and Keeton, 2010), as well as improve soil moisture availability, tree water status, and resilience to periodic drought (D'Amato et al., 2013;Sohn et al., 2016;Bradford and Bell, 2017). Variable retention harvesting (VRH) is a relatively new approach to partial harvesting used to balance economic with ecological objectives by creating a post-harvest forest structure that emulates that resulting from the natural disturbance regime (Franklin et al., 2007;Palik and D'Amato, 2019). Specifically, this approach retains a proportion of canopy trees and other biological structures at various densities and spatial patterns to enhance the structural complexity, biodiversity, and resilience of the resulting stand, while maintaining acceptable levels of productivity (Franklin et al., 2007;Palik and D'Amato, 2019). The primary considerations in application of VRH systems are the trees and structures to be retained, and the level, spatial pattern, and the period of retention best suited to the forest type being managed (Palik and D'Amato, 2019). Generally, VRH is used to create stands containing areas of both uniformly dispersed trees of different densities and aggregated patches of residual trees differing in size and shape (Palik and D'Amato, 2019;Franklin and Donato, 2020).
The effect of VRH on productivity and other forest values has been the focus of several, operational scale experiments established in North America over the past two decades (Palik et al., 2002(Palik et al., , 2014Aubry et al., 2009;Puettmann et al., 2016;Xing et al., 2018;Beese et al., 2019). The impact of VRH on productivity varies with forest type, tree species, the level and spatial pattern of tree retention, time after treatment, and other factors (Roberts and Harrington, 2008;Xing et al., 2018). In general, VRH has been shown to increase stem growth rates of individual trees as compared with unharvested areas (Roberts and Harrington, 2008;Powers et al., 2010;Xing et al., 2018). Tree growth rates and mortality are negatively related to retention level, but this relationship varies with species (Maguire et al., 2006;Xing et al., 2018). As well, trees located near the perimeter of retained patches have shown higher (Powers et al., 2009) or no difference in growth rate (Roberts and Harrington, 2008) compared to those in the interior. At the stand level, unharvested control plots have greater total biomass or volume growth compared to VRH treatments regardless of the level or pattern of retention (Maguire et al., 2006;Palik et al., 2014). Given the initial emphasis on VRH effects on tree regeneration and biodiversity, it remains unclear which level and pattern of tree retention best balances timber production with climate change adaptation and mitigation objectives for different forest types and site conditions in natural stands and plantations (Xing et al., 2018;Franklin and Donato, 2020).
Red pine (Pinus resinosa Ait.) plantations were established in many areas in the Great Lakes region of North America beginning in the early 1900s to control soil erosion on abandoned agricultural land and to rehabilitate large cutover and burned forests (Cayford and Bickerstaff, 1968;Johnson, 1995;Buckman et al., 2006). Red pine was favored for these efforts due to its commercial value, genetic uniformity, suitability to evenaged management, and ability to grow well in drought-prone, nutrient-poor soils (Stiell, 1959;Buckman et al., 2006). Many of the red pine plantations established on degraded agricultural land are located south of its natural range in highly fragmented, urbanized landscapes where the area of forested land has been dramatically reduced (Borczon, 1982;Elliott, 1998). These plantations are typically being managed to restore them to native forest types through traditional row or selection thinning methods coupled with regeneration of desired species by planting or release of advance reproduction (Parker et al., 2008;Abella, 2010).
In this study, we examined the impact of VRH on residual tree-and stand-level growth and C increment in an 88-yearold red plantation established beyond its natural range limit in southern Ontario, Canada. We used dendrochronological methods to examine VRH effects because tree rings provide fine spatial and temporal resolution in radial stem growth (Davis et al., 2009;Dye et al., 2016;Pompa-García et al., 2018) and are reliable proxy measures of tree-and stand-level C sequestration rates (Babst et al., 2014). The specific objectives of this study were to determine the effect of: (1) the level and spatial pattern of retention on post-harvest tree-and stand-level stem growth and C increment during the 5-year post-harvest period; and (2) tree location within aggregate patches of residual trees on tree stem growth and C increment.

Site Description
The study site is located within the Turkey Point tract of the St. Williams Conservation Reserve (SWCR) located in Norfolk County, in southern Ontario, Canada (42 • 42 16 N, 80 • 21 29 W) ( Figure 1A). The SWCR consists of approximately 1034 ha of Crown Land managed in collaboration with the Ontario Ministry of Northern Development, Mines, Natural Resources and Forestry (NDMNRF). The primary management objective of SWCR is to increase biodiversity, protect speciesat-risk, and restore these lands to their native ecological communities (SWCRMP (St. Williams Conservation Reserve Management Plan), 2007). This study was conducted in a 20.2 ha FIGURE 1 | (A) Location of the field site in relation to the Great Lakes in southern Ontario, Canada. The blue dashed line represents the approximate southern limit of red pine trees in Ontario as defined by Rudolf (1990). (B) Aerial view and location of the 20 treatment plots with five treatments: unharvested control (CN), 33% dispersed retention (33D), 55% dispersed retention (55D), 33% aggregate retention (33A), and 55% aggregate retention (55A). red pine plantation established in 1931. The area was most likely occupied by mixed oak and pine forests prior to settlement of the region. This plantation is part of the Turkey Point Flux Station and is associated with the Global Fluxnet and Global Water Futures research network (Peichl et al., 2010). The plantation was thinned around 1959-1960 when every fourth row of trees was removed. The plantation was thinned a second time in February 2014 to demonstrate and test the VRH approach with the longterm objective of restoring the area to a native forest type while enhancing biodiversity through the creation of a structurally complex residual tree canopy. A typical operational approach to the management of this plantation would be to reduce stand density according to guidelines for uniform shelterwood with the objective of initiating regeneration (Ontario Ministry of Natural Resources and Forestry (OMNRF), 2015).
A pre-harvest survey of the area conducted in 2011 using a series (n = 60) of 12.62 m radius plots indicated the canopy was almost entirely composed of dominant and codominant red pine, with a few scattered white pine (Pinus strobus L.) and black oak (Quercus velutina Lam.). The plantation was overstocked with an average total basal area (BA) of 38.4 ± 5.5 m 2 /ha and a density of 641.3 ± 103.2 trees/ha. Red pine averaged 28.6 ± 4.3 cm in diameter at 1.3 m (DBH) and ranged from 17.2 to 47.2 cm. Height of red pine averaged 23.8 ± 2.8 m. The subcanopy consisted of naturally established Prunus, Acer, and Quercus spp. of intermediate crown class and averaged about 70 trees/ha and 13.2 cm in DBH across the study area. The VRH study uses a completely randomized experimental design, with four replicates of five treatments ( Figure 1B). Treatments were applied to 20 plots, each averaging 0.89 ha in area and separated by fire access roads ( Figure 1B). The treatments include: 33% dispersed retention (33D), 55% dispersed retention (55D), 33% aggregate retention (33A), 55% aggregate retention (55A), and an unharvested control (CN). The dispersed treatments were created according to provincial guidelines for using the shelterwood silvicultural system for regeneration of red and white pine forests (Ontario Ministry of Natural Resources and Forestry (OMNRF), 2015). Tree marking of the dispersed treatments was aimed at uniform spacing of residual trees with overstory crown closure of 33 or 55%. The aggregate treatments were harvested to retain four (two large and two small) uncut, roughly circular patches over 33 or 55% of the plot area. Over the entire study area, 4,737 red pine with DBH of 16-48 cm were harvested, with ca. 74% trees removed having DBH of 24 to 32 cm. Pre-and post-harvest structural features and photographs of the VRH treatments are presented in Table 1 and Figures 2A-F, respectively.

Dendrochronological Methods
Cores were collected from at least 30 red pine trees within each of the 20 treatment plots at ∼1.3 m above ground (DBH) during  -harvest (2009-2013) and post-harvest (2014-2018) mean and standard deviation (SD) of stand density, basal area (BA), and diameter at 1.3 m (DBH) of living overstory trees ≥10 cm DBH for each treatment.

Treatment
Pre-harvest Post-harvest  the spring and summer of 2019 using a 5 mm increment borer (Haglöf, Sweden). For each tree, two cores were collected at approximately 90 • from each other to account for variability within individual trees. In the 33A and 55A treatments, cores were collected from an equal number of trees from the interior and exterior of the patch to account for potential spatial variation in radial growth with location. As well, exterior trees sampled were evenly distributed around the patch perimeter to account for possible effects of aspect on radial growth (Roberts and Harrington, 2008;Powers et al., 2009). The DBH was measured for each tree sampled for dendrochronological assessment of annual C increment. Tree cores were prepared for ring width measurement according to standardized dendrochronological methods (Stokes and Smiley, 1968;Speer, 2010). Cores were sanded with increasingly finer grades of sandpaper until the tree rings were clearly visible (Stokes and Smiley, 1968) and visually crossdated using the list method (Yamaguchi, 1991). Ring widths were measured to the nearest 0.001 mm using CooRecorder and CDendro (Larsson, 2018). Visual crossdating was verified using the COFECHA software program (Holmes, 1983(Holmes, , 1994 before assigning calendar dates to each tree ring (Speer, 2010).

Climate
The climate in southern Ontario is classified as Warm-Summer Humid Continental Climate (Kottek et al., 2006) and is characterized by cool winters and hot and humid summers (Crins et al., 2009). Climate data for the study area for the period of the tree-ring chronologies were obtained from Environment Canada climate stations located in Delhi, ON and Tillsonburg, ON located within 34 km from TP31 (Supplementary Table 1 and Supplementary Figure 1) 1 . From 1937 to 2018, the mean annual temperature was 8.0 • C, the mean annual maximum temperature was 13.2 • C, and the mean annual minimum temperature was 2.9 • C. Precipitation was evenly distributed throughout the year and the total annual precipitation averaged 962 mm.

Percent Growth Change
Percent growth change (%GC) over two discrete 5-year periods for each sampled tree was calculated using ring widths and the radial growth average technique (Nowacki and Abrams, 1997) as: where M 1 is the mean pre-harvest radial growth (2009)(2010)(2011)(2012)(2013) and M 2 is the mean post-harvest radial growth (2014-2018).

Annual Stem Biomass and C Increment
Tree ring derived DBH values and allometric equations were used to estimate annual stemwood biomass and C increment for each sample tree. It is important to note that estimating aboveground biomass from annual growth rings is a source of uncertainty (Martin-Benito et al., 2021). The allometric equations developed by Lambert et al. (2005) were selected for use because they were constructed from a large number of red pine trees (n = 371), primarily from Ontario, with an average DBH and height representative of the study plantation. Allometric equations using tree height and DBH measurements provide more accurate estimates of aboveground biomass, but annual height measurements of our sample trees were not available (Lambert et al., 2005;Cienciala et al., 2008). The allometric equation used to estimate the biomass of the individual tree stemwood is represented by the general form: where y wood is the dry biomass of the stemwood (kg), β wood1 and β wood2 are model parameters with coefficient estimates, D is DBH (cm), and e wood is the error term (Lambert et al., 2005). For a mean prediction of biomass, e wood was assumed to be 0, as it is assumed e wood follows a normal distribution with a mean of 0. Based on the β wood1 and β wood2 estimates for red pine, the allometric equation used for this study is as follows (Lambert et al., 2005):

Calculating Tree-Level Stemwood C Increment
The DBH of each tree in each study year (2009-2018) was calculated as: where DBH k is the value for any k year (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017), j is the final year of our assessment period (2018) and the upper limit of the summation, and RWi is the sum of the ring widths from year i to year j (Dai et al., 2013). The estimated DBH for each year and the allometric equation were used to estimate the total stemwood biomass increment for each year as: Annual biomass year k = total biomass year k − total biomass year k − 1 (5) where k represents the year of interest. Although C concentration may vary among species and plant parts (Lamlon and Savidge, 2005;King et al., 2007), C increment was estimated by multiplying, annual stemwood biomass by 0.5, assuming that approximately 50% of red pine stem dry mass is C (King et al., 2007;Dai et al., 2013). Using the generic value of 50% may not accurately estimate the absolute amount of C sequestered, but it is suitable for comparison of VRH treatment effects on C increment.

Stand-Level C Estimates
Stem density (trees/ha) of living red pine canopy trees for all treatment plots during the pre-harvest period was estimated using 2011 pre-harvest inventory data. The total number of red pine trees removed from each treatment plot was subtracted from its pre-harvest stand density to obtain initial (2014) postharvest density. Annual changes in living red pine stand density during the post-harvest period were estimated using a subsample (n = 10-38) of canopy trees in each plot that were monitored for mortality (Supplementary Table 2). Because of the relatively small number of trees monitored, stand density for each plot was estimated using the mean percent mortality for a given treatment (Supplementary Table 2). Mean annual tree-level C increment (kg C/tree) for each treatment plot was multiplied by stem density to calculate pre-and post-harvest stand-level estimates of mean annual C increment (kg C/ha).

STATISTICAL ANALYSIS
One-way analysis of variance was used to quantify VRH treatment effects on %GC and tree-and stand-level annual C increment. Where significant treatment effects were exhibited, a post hoc Tukey test was performed to quantify differences among treatment means (computed with https://astatsa.com/OneWay_ Anova_with_TukeyHSD/). Differences among treatments in annual C increment between pre-and post-harvest periods, and between interior and exterior trees within and among the aggregated treatments, were quantified using t-tests.

Treatment Effects on Tree-Ring Growth Responses
The tree-ring chronologies for all 20 treatment plots spanned the period from 1936 to 2018 and exhibited comparable descriptive statistics ( Table 2). The mean series intercorrelation was 0.646 and ranged from 0.551 to 0.736, indicating that ring growth of all sampled red pine trees was influenced by a common growth factor. The mean annual ring width growth increment was 1.64 mm and ranged from 1.47 to 1.89 mm. Autocorrelation was high (<0.871) for all treatment plots, with a mean value of 0.903 ( Table 2), indicating that tree growth in a given year is influenced by growth in the previous year. The DBH of the 600  trees sampled averaged 31.1 ± 4.4 cm and ranged from 21.1 to 45.8 cm ( Table 2). Prior to harvest, ring-width index (RWI) was similar among treatments (Figure 3). Following harvest, RWI of all VRH treatments exhibited similar year-to-year high-frequency variability as the CN treatment. Beginning in 2016, the third growing season after harvest, the 33D and 55D treatments exhibited a trend toward higher RWI as compared with the 33A, 55A, and CN treatments (Figure 3). Treatment type had a significant (p < 0.05) effect on mean annual %GC over the postharvest period, with values being significantly higher in the 33D than the 55A treatment ( Figure 4A). In the aggregate treatments, mean annual %GC of the post-harvest period of exterior trees was significantly (p < 0.05) higher than interior trees in the 33A treatment but did not differ (p = 0.139) with tree position in the 55A treatment ( Figure 4B). For both aggregate treatments combined, mean annual %GC was significantly (p ≤ 0.001) higher in the exterior trees.

Treatment Effects on C Increment
Treatment effects on pre-and post-harvest rates of mean annual C increment depended on the spatial scale at which this variable was expressed (Figure 5). At the individual tree level, mean annual C increment in the post-harvest period was significantly lower in the CN (p ≤ 0.001), 33A (p = 0.026), and 55A (p ≤ 0.001) treatments compared to the pre-harvest period ( Figure 5A). Preand post-harvest mean annual C increment per tree did not differ in the 33D (p = 0.724) and 55D (p = 0.528) treatments. Mean annual C increment in the exterior trees in the postharvest period was significantly higher than the interior trees in the 33A [t(38) = 4.623, p < 0.001], 55A [t(38) = 4.793, p ≤ 0.001] and both treatments combined [t(38) = 6.677, p < 0.001]. At the stand scale, mean annual C increment was significantly (p ≤ 0.001) lower in the post-harvest period for all five treatments (Figure 5B).
Treatments also differed in tree-and stand-level mean annual C increment within the post-harvest period (Figure 6). At the individual tree scale, the two dispersed retention treatments had significantly (p < 0.05) higher mean annual C increment than the CN and the two aggregate treatments (Figure 6A). For the aggregate treatments, tree location had a significant effect on growth, with exterior trees having significantly (p < 0.001) higher mean annual C increment than the interior trees ( Figure 6B). By comparison, at the stand scale, mean annual C increment was significantly (p < 0.01) higher in the CN treatment than all VRH treatments except the 55D treatment ( Figure 6C). Mean annual C increment did not differ between dispersed and aggregate treatments of the same retention level, but mean annual C increment was significantly (p < 0.05) higher in the 55D than the 33D and 33A treatments (p < 0.01).
Patterns of inter-annual variability in mean annual C increment rates at the individual tree-and stand-level during the 10-year pre-and post-harvest period suggest that the two dispersed retention treatments began to exhibit greater C uptake compared to the three other treatments beginning in 2016, the third growing season after harvest (Figure 7).

Impact of Retention Level and Pattern on Stem Growth
Stand density reduction through thinning or partial harvest systems is employed in part to release residual canopy trees from neighborhood competition and increase tree growth rates (Gilmore and Palik, 2006). The pattern and intensity of tree removal determine the degree of release from competition. Annual variation in RWI from residual chronologies indicates FIGURE 5 | Box and whisker plots comparing pre-harvest (2009-2013) and post-harvest (2014-2018) treatment effects on (A) mean annual tree C uptake (kg C/tree) and (B) mean annual stand C uptake (kg C/ha). The * indicates a significance level of p < 0.05 and the * * indicates a significance level of p < 0.001. The X in the middle of each box represents the treatment mean.
FIGURE 6 | Box and whisker plots comparing treatments during the post-harvest (2014-2018) for (A) mean annual tree C increment (kg C/tree), (B) mean annual tree C increment (kg C/tree) between interior and exterior trees in combined aggregate plots, and (C) mean annual stand C increment (kg C/ha). Means sharing the same letters are not significantly different. For (B), the * * indicates a significance level of p < 0.001. The X in the middle of each box represents the treatment mean.
that red pine radial growth exhibited a strong response to both interannual climatic variability and stand density changes over our 10-year study period. The shared annual variability in RWI is associated with climate-growth relationships established for red pine. Cooler, wetter growing season conditions favor higher growth in red pine growing in natural stands (Graumlich, 1993;Girardin et al., 2006) and plantations (Magruder et al., 2013;Ashiq and Anand, 2016). The gradual decline in RWI exhibited prior to harvest for all treatments is likely due to high inter-tree competition and growth suppression when stand density and BA were comparatively high. Post-harvest tree RWI varied with VRH treatments and indicated that stem radial growth rate of residual canopy trees depended primarily on the pattern of retention as it influenced the degree of reduction in neighborhood competition. The trend toward higher RWI in the dispersed treatments is most likely due to near-complete release from competition by adjacent trees in the uniformly spaced trees. In contrast, only the perimeter trees in the aggregate treatments experience a partial release from competition and more limited increase in resources (Chen et al., 1992;Sherich et al., 2007;Powers et al., 2009). The positive mean %GC values for the dispersed treatments and negative mean %GC for the aggregate treatments are consistent with trends in RWI. Further, although exterior trees had higher mean %GC than interior trees, mean values of %GC for exterior trees were negative and largely similar to that of the CN treatment. Interestingly, %GC of exterior trees was higher than interior trees in the 33A treatment and both aggregate treatments combined, but not the 55A treatment, suggesting that exterior trees in smaller patches with higher perimeter to area ratios may experience greater resource availability than in larger patches.

Impact of Retention Level and Pattern on Tree-and Stand-Level C Increment
Not surprisingly, the pattern of annual changes in C increment over our 10-year study period was consistent with that for RWI, with only trees in dispersed treatments exhibiting a gradual increase in C increment after harvest. The significantly lower post-harvest C increment for the control and aggregate treatments indicate that trees in residual patches respond similarly to the unharvested treatment plots, despite the comparatively higher radial growth observed for exterior trees that experienced a partial reduction in competition. These results are consistent with the findings for natural red pine stands (Powers et al., 2009(Powers et al., , 2010 and other tree species (Deal and Tappeiner, 2002;McDonald and Urban, 2004;Sherich et al., 2007;Roberts and Harrington, 2008). In contrast to our findings, reduced growth of Douglas-fir in aggregate patches was attributed to poor growth response of exterior trees (Maguire et al., 2006).
The similarity in post-harvest mean C increment for trees in the dispersed treatments is likely a reflection of the delayed, gradual increase in radial growth following harvest that was not expressed until the third growing season after treatment. This response is probably due to so-called "thinning shock" that results from the abrupt exposure of trees to higher radiation, air temperature, wind speed, and a more evaporative atmospheric environment in thinned stands (Youngblood, 1991;Latham and Tappeiner, 2002;Bebber et al., 2004). This gradual, growth response is thought to be related to a brief period of physiological and morphological acclimation of trees to the new post-harvest microclimate when carbohydrate allocation to root growth may be favored over aboveground growth (Bladon et al., 2006;Powers et al., 2010). The comparatively greater variation observed in mean C increment in the dispersed treatments is probably associated with both heterogeneity in growth resources and differences in acclimation capacity among residual trees of different size and crown class (Roberts and Harrington, 2008). Trees of dominant and co-dominant crown classes with higher live crown ratio and leaf area show more immediate, positive growth response to thinning than subordinate, smaller diameter, less vigorous canopy trees of intermediate crown class (Horton and Bedell, 1960;Gilmore and Palik, 2006).
The two dispersed treatments exhibited significantly higher mean C increment in the post-harvest period than the control and aggregate treatments, while the C increment of the aggregate treatments did not differ from the unharvested control treatment, consistent with the findings of Roberts and Harrington (2008). This indicates that in our study, mean C increment is dependent on the pattern of retention, and not the level of retention. Maguire et al. (2006) also found higher stem growth in the dispersed compared to aggregate treatments due to increased access to resources and minimal crown overlap.
Expressed at the stand level, mean C increment was higher during the pre-harvest than the post-harvest period for all five treatments. Decreasing mean annual C increment in the control treatment reflects increasing stand basal area and tree growth suppression in unharvested plots, which is the rationale behind the use of frequent, periodic thinning to maximize growth and reduce mortality in managed red pine stands (Gilmore and Palik, 2006). Despite increased growth rates due to complete or partial release of residual canopy trees from competition, the greatly reduced stem density due to the thinning process reduced the total stand-level C increment in most VRH treatments (Maguire et al., 2006). This is not surprising given that thinning increases individual tree growth rates but rarely increases stand-level rates (Gilmore and Palik, 2006;Davis et al., 2009), since it takes time for these faster individual rates to compensate for lower densities.
Mean stand-level C increment in the post-harvest period was significantly higher in the unharvested control treatment than all other treatments, regardless of level or pattern of retention. Again, this is an expression of the much higher stand density in the control treatment despite lower tree C increment. In contrast to findings with tree-level C increment, stand-level C increment is driven by retention level. Standlevel C increment was significantly lower in the 33% compared with the 55% retention level regardless of pattern. In natural red pine stands, stand-level growth was higher in unharvested versus VRH treatments but was only marginally higher in dispersed treatments compared to aggregate treatments of the same retention level (Palik et al., 2014).

CONCLUSION
In this study of VRH effects on tree-and stand-level growth and C increment of red pine, we used tree rings and allometric equations to estimate stem radial growth and C increment, respectively. While using dendrochronological methods to estimate annual productivity has potential limitations (Martin-Benito et al., 2021), our comparatively large sample of randomly selected residual trees and our focus on a short, 10-year growth period, coupled with pre-and post-harvest tree and stand measurements, eliminated or reduced these sources of error (Babst et al., 2014;Dye et al., 2016). The results of this study indicate that VRH has a significant effect on tree-and stand-level stem growth and C increment of an 88-year-old red pine plantation that was significantly related to the level and spatial pattern of retention. Dispersed VRH treatments at both levels of retention were more effective in promoting post-harvest tree-level growth and C increment than aggregate treatments. In contrast, standlevel growth and C increment were highest for greater levels of retention, regardless of pattern. We recognize that stand C sequestration depends on above and below ground growth response of all forms of vegetation as well as treatment effects on C contained in dead biomass and mineral soil. Our focus on living red pine canopy trees clearly does not quantify the potential effects of VRH on total C accumulation.
The results of this study are applicable to management of red pine plantations in southern Ontario and other parts of the Great Lakes region where balancing timber production with climate change adaptation and/or mitigation is an objective (D'Amato et al., 2011). Where improved tree growth rates are desired, dispersed retention was the superior treatment. Where mitigation through enhanced stand-level C sequestration is a primary objective, a higher level of retention is desirable. In stands where adaptation and mitigation are of high priority, a mixture of dispersed and aggregate treatments of different retention levels may best meet both objectives.

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

AUTHOR CONTRIBUTIONS
JZ, MP, and MA devised the research objectives and formulated the research questions. JZ conducted the field data collection with assistance from SM and performed all laboratory analysis. MP contributed to the statistical analysis. WP and KE designed and established the variable retention harvest experiment. JZ wrote the manuscript with assistance, critical feedback, and revisions from MP, SM, WP, KE, and MA. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We also acknowledge the continued in-kind support of the Ontario Ministry of Natural Resources and Forestry (OMNRF), Environment and Climate Change Canada, and the St. Williams Conservation Reserve Community Council (SWCRCC). We thank Audrey Heagy of the SWCRCC for her continued support and collaboration. Finally, we gratefully acknowledge the contributions of Steve Williams (retired) of OMNRF. Without his collaborative spirit, professional curiosity, and desire to "try something new, " this experiment would not have been possible.