Effects of Whole-Tree Harvesting on Species Composition of Tree and Understory Communities in a Northern Hardwood Forest

In the northeastern United States, whole-tree harvesting is widely used to supply fuel to biomass energy facilities, but questions remain regarding its long-term sustainability. We have previously reported findings indicating no short-term decrease in forest productivity in whole-tree harvested sites when compared with similar conventionally (stem-only) harvested sites. Here we present additional results of the same study, but focus on the effect harvest treatment has on the species composition of the regenerating forest. Within northern hardwood forests in central New Hampshire and western Maine, regeneration surveys were conducted on four (4) small clearcuts in 2010 and twenty-nine (29) small clearcuts in 2011. The species and diameter of trees > 2 m in height were recorded within 1 m or 2 m-radius plots and used to calculate the biomass fraction of each species. The 2010 study additionally measured the density of trees < 2 m in height and the diversity of understory non-tree species. Non-metric multidimensional scaling and multi-response permutation procedures were used to determine the effect of harvest treatment had on community-wide tree species composition. Potential differences were also examined on a species-byspecies basis. Both analytic methods indicated no significant differences in species composition of tree species or understory communities. Within the limits of our data, we conclude that no significant effects of residue removal on species composition are observed within our sample of northern hardwood sites at this early stage of stand development.


Introduction
The burning of woody biomass to generate electricity has the potential to supply locally produced renewable energy to the northeastern United States.With its well established forest products sector, ample forest resources and high percentage of land in private ownership to quickly meet market demands, New England represents an ideal location for a wood-based biomass industry (Benjamin et al., 2009).In forested regions such as the northeastern US, whole-tree harvesting (WTH) is commonly used to supply biomass energy plants with feedstock.During WTH, traditionally non-merchantable wood residue in the form of branches and treetops are chipped, sold, and eventually burned for energy.In conventional harvesting (CH), this residue remains on site and decomposes, returning its nutrients to the soil.WTH allows a landowner to earn additional income while supplying fuelwood to local biomass energy plants.However, the effect of residue removal during whole-tree harvesting on the future species composition and productivity of our forests is not well known.
Concerns about effects of removing additional woody residues on the future forest have long been recognized (Boyle, 1976;Kimmins, 1977;Anderson, 1985).Removal of additional nutrients contained in harvest slash may eventually leave a site unable to grow trees at its highest rate, impairing its productivity.Woody residue may provide structural benefits to wildlife, and its removal may affect animals that utilize rotting woody debris as habitat.Finally, the removal of harvest slash may change the microclimate of the forest floor, altering which seeds and seedlings are ultimately successful in regenerating and affecting the species composition of the next stand.
Branches, fine twigs and needles that make up harvesting slash contain a large proportion of a tree's nutrients and their removal could have significant effects on site fertility.The primary concern is that WTH could impair future site productivity, leaving a forest that regrows trees slowly or of poor quality.Studies examining soil nutrient content have found discrepant results; WTH harvesting was found to deplete certain soil nutrients on some sites (Olsson et al., 1996a;Vanguelova et al., 2010) but have no effect on others (Johnson et al., 1991;Olsson et al., 1996b).Several studies have measured productivity declines directly following WTH, finding trees that are smaller, have smaller diameters and lower biomass than those of CH sites (Proe et al., 1996;Egnell & Valinger, 2003;Scott & Dean, 2006;Walmsley et al., 2009).However, all of these studies were conducted in softwood plantations and additional studies show no productivity decline on similar sites (Hendrickson, 1988;Smethurst & Nambiar, 1990;Dyck et al., 1991;Sanchez et al., 2006;Tan et al., 2009;Saarsalmi et al., 2010).In a study of natural stands of northern hardwoods, we found no indication of short-term decrease in forest productivity on WTH sites when compared with similar CH sites (Roxby & Howard, 2013).The majority of previous studies do not examine ecological effects such as species composition, as most make use of artificial regeneration following harvesting.
Through removal of harvest residues, WTH likely reduces the amount of fine and coarse woody debris on a forest floor, which may have negative impacts on wildlife that require decomposing wood for habitat.Effects are likely to be both site and species specific, with one study showing decreases in richness and abundance of songbirds (Lohr et al., 2002) and another showing a decrease in small mammal populations as a result of coarse woody debris removal (Loeb, 1999).Bowman et al. (2000) found rotting logs to be important to red-backed vole populations on a landscape scale, but no evidence for a link between coarse woody debris and other small mammal populations.Coarse woody debris is also an important habitat requirement for many species of salamander, and its removal from the landscape would likely have negative effects on their populations and diversity (deMaynadier & Hunter, 1995).
We might also expect the retention or removal of woody debris to have an effect on the microclimate of the forest floor following harvest.This may regulate to some degree the tree species that are able to regenerate and recolonize the site.Proe et al. (2001) found that WTH sites experienced higher wind speeds and larger temperature fluctuations when compared with sites where harvest residues were retained.These effects quickly dissi-pated as vegetation was re-established, but would likely play a significant factor in early seedling establishment during natural regeneration.In an upland hardwoods site in Tennessee, Mann (1984) found more vigorous stump sprouting following CH, but more seedling establishment following WTH.Certain species also appeared to regenerate more strongly following a certain treatment; chestnut oak (Quercus montana), black oak (Quercus velutina) and red maple (Acer rubrum) did best following CH while black gum (Nyssa sylvatica), tulip poplar (Lireodendron tulipifera) and sourwood (Oxydendrum arboreum) responded better to WTH.Brakenhielm & Liu (1998) found that WTH significantly altered species composition, diversity and succession in a lichen-Scots pine forest in Sweden.
Many studies based in the White Mountains of New Hampshire have correlated changes in environmental factors with changes in species composition.In a study of small canopy openings, McClure and Lee (1993) showed that gap size, gap age and position within the gap were all significant predictors of tree species composition.Leak (1976) and(1978) outlined a classification scheme based on soil parent type and drainage class and correlated each category with tree species composition in late and mid-successional stands.Lee et al. (2005) also found soil substrate to be a significant predictor, and additionally showed elevation to play a major role in tree species composition.All of these studies point to the underlying effects that soil moisture and nutrient content have on regeneration.Whole-tree harvesting affects these same two factors by removing additional nutrients from the site, and likely increasing the amount of sunlight the understory receives, affecting soil moisture.This points to the potential for changes in tree species composition following WTH and natural regeneration in comparison with sites where residues are retained.We are not aware of any studies that examine the effects of residue removal or retention of the species composition of the next stand within our study region.
With whole-tree harvesting increasingly used to meet demand for woody biomass chips, research is needed to determine the ecological effects of this practice on species composition of regeneration.The vast majority of timber harvests in the northeastern United States rely on natural regeneration to repopulate stands and small changes to the microclimate of the forest floor may result in significant changes to species composition of the regeneration.Compositional changes in forest regeneration can alter the productivity and merchantability of the next stand of trees.The objective of this study is to determine the effects of residue removal from WTH on natural regeneration in northern hardwood sites in New Hampshire in comparison with conventionally harvested stands where residues were retained.Our focus is on the species composition and abundance of 1) trees > 2 m tall; 2) tree saplings under <2 m tall; and 3) non-tree understory plants.We made an additional effort to determine if any effects observed could be explained by spatial variations in radiation intensity alone.Two separate studies, different in their scope, were conducted and are summarized here.The first examines all three objectives on four sites in the Bartlett Experimental Forest in the White Mountains of New Hampshire; the second examines just the first objective on 29 sites in central New Hampshire and western Maine.

Site Description and History
Sites were measured over two summers as a part of concurrent studies examining the impacts of residue removal on site productivity.Four 12-year old small clearcuts in the Bartlett Experimental Forest (BEF) were measured during summer 2010 (Figure 1).These cuts were between 0.9 and 2.3 hectares (2.2 to 5.7 acres) in size and were harvested during November and December of 1998.Trees were whole-tree harvested in two of the cuts (residues removed) and conventionally harvested in the others (residues left on site).Sites were between 260 and 310 meters in elevation and located on glacial till parent material.Underlying bedrock was primarily igneous, with the WTH sites on Bethlehem Granodiorite and Amonoosuc Volcanic Formations.Conventionally harvested sites were located over Bethlehem Granodiorite, Granite and Littleton Formations (Lyons et al., 1997).Terrain on the sites gently sloped (6 -7 degrees) to the northeast (3 sites) or west (1 site).Although no measurements were made prior to harvest, timber sale reports reveal that the previous stands consisted primarily of northern hardwood tree species.
Twenty-nine additional clearcuts in central New Hampshire and western Maine were measured during summer 2011 (Figure 1).These cuts were not established for research purposes as were the sites at the BEF, but were selected to be as similar to the BEF sites as possible: northern hardwood stands that were commercially clearcut 10 -14 years ago, primarily for regeneration purposes, with areas that were between 0.6 and 5.7 hectares (1.6 -14.1 acres).Fourteen whole-tree harvested and 15 conventionally harvested sites that fit these criteria were located and measured.Sites were between 250 and 650 meters in elevation, located on soils with glacial till parent material, and classified as either well or moderately well drained.No data on the species composition of previous stands were available.Conventionally harvested sites were located on US Forest Service lands in the White Mountain National Forest (WMNF) while WTH sites were located on private lands.Whole-tree harvesting is not currently practiced on the WMNF due to concerns about nutrient depletion, especially on shallow or sandy sites (W.B. Leak, personal communication).

Data Collection
As our studies were conducted using different methods, they were analyzed separately; we will reference them using the year of the data collection (2010 or 2011).

2010 Study
Heights and locations of uncut trees along the edge of each clearcut were recorded during initial site visits.This information was used to accurately measure the size and shape of each site and to create a model to predict radiation intensity during the growing season.This spatial solar radiation model (described in detail in a later section) was based on simple assumptions and calculated the radiation intensity expected at any point within a site based on edge tree heights and sun position.We expected regeneration to vary based on spatial variation in sunlight levels (Finzi & Canham, 2000); areas that were shaded by nearby edge trees would be more likely to regenerate shade tolerant species than areas receiving full sunlight.
Using ArcView GIS 3.3, a 6 × 6 meter grid of potential plots was overlaid on each clearcut, creating a group of measurable plots that was dependent on cut size.For each potential plot, the radiation intensity throughout the growing season was calculated using the spatial radiation model and reported in megajoules per square meter (MJ/m 2 ).Potential plots were categorized as receiving high (>4478 MJ/m 2 ), medium (3516 -4478 MJ/m 2 ) or low (<3516 MJ/m 2 ) amounts of radiation.Actual plots to be measured were selected using stratified random sampling using radiation categories as strata.Thirty plots were allocated to each of the smaller cuts, while 50 were measured in the two larger cuts.Thus, if 60% of the potential plots were in the high radiation category and that cut was allocated 50 actual plots, 30 plots (0.6 × 50) were selected randomly from the pool of high radiation category plots to be measured.This method allowed the measured plots to be representative of the different sunlight conditions that were actually present within the cut (i.e.no radiation condition was over-or under-represented.See Table 1).
Each of the 160 plots was visited during June and July 2010.Plot centers were located using a handheld GIS unit and were generally accurate to within 5 m.Within a 2 m-radius plot, the height, diameter and species of all trees taller than 2 m were recorded.Diameter at breast height (DBH) was measured to the nearest 0.1 cm using a research grade DBH tape; small trees (<2 cm DBH) were estimated to the nearest 0.5 cm.Tree height to the topmost living leaf was recorded to the nearest 0.3 cm (0.1 inch) using a Senshin fiberglass measuring pole.Within a 1 m-radius subplot, saplings shorter than 2 m were tallied by species.We made no attempt to distinguish between saplings of seed and sprout origin.Non-tree species occurring within this subplot were also recorded.Two plant groups were difficult to separate into species and were recorded by genus (Carex.spp.and Rubus spp.); other plants we were unable to identify were excluded from our analysis.Slope and aspect of each plot were measured using a clinometer.A qualitative assessment was made as to whether or not vegetation within the plot had been browsed by moose (Alces alces) during the past year.This observation was made to categorize and isolate variation due to repeated moose browsing of regeneration.Seven plots were either located on logging roads or fell outside the actual cut area due to limitations in GPS accuracy and were excluded from our analysis.

2011 Study
Due to research resource considerations, the sampling protocol for the 2010 study was modified for the 2011 study.Plots were located according to a square grid pattern oriented north-south or east-west, with grid spacing varying by site based on the size of the cut opening.Plot spacing varied from 90 to 160 feet based on an initial goal of establishing approximately 15 plots per site.Midway through the season, data analysis revealed additional plots were needed on sites with a high coefficient of variation in tree height.We revisited these sites and measured 4 to 34 additional plots, proportional to the variation exhibited by the site.Grid spacing on sites not yet visited was decreased to accommodate about 30 plots per cut.a See text for category definitions.Number of plots measured within each site was allocated proportionally by the composition of each cut, split by radiation category.Sites were either whole-tree (WTH) or conventionally harvested (CH) and between 0.9 and 2.3 hectares in size.
With pacing difficult due to the density of vegetation, a handheld GPS unit was used to navigate to plots along each grid transect.Within each 1m-radius plot height, diameter and species were recorded for all trees > 2 m in height.Diameter at breast height was measured using a research grade DBH tape and height to the uppermost living leaf was measured using a Senshin fiberglass measuring pole.If a plot showed signs of current year moose browse damage or was composed of non-tree regeneration (i.e.field, road or small stream), this was noted.No measurements on tree saplings shorter than 2 m or non-tree understory vegetation were made due to research resource constraints.

Spatial Solar Radiation Model
Expected radiation intensity values for each plot were calculated using the GPS locations of measured plots and the positions and heights of trees along the edge of the clearcut area.The GIS programs ArcView 3.3 and Arc-GIS 10.0 were used to create this model and determine the total radiation intensity a given plot would receive during the growing season (assumed to be April 1 st through October 1 st ).First, a list of solar elevations and azimuths was generated using the Excel program SolRad v. 1.2 (Pelletier, 2011).With this list of solar positions, the program determined using trigonometry if each plot was illuminated or shaded by trees along the clearcut edge for each hour of each day during the growing season.The radiation intensity for each hour a plot received sunlight was summed and a final value, measured in megajoules per square meter (MJ/m 2 ), was generated.This summed radiation was used to categorize each plot into one of three radiation categories: low (290 -3516 MJ/m 2 ), medium (3516 -4478 MJ/m 2 ) or high (4478 -5049 MJ/m 2 ).Radiation intensity was calculated prior to data collection in 2010 and was used to proportionally randomly stratify plots within radiation category.In 2011, radiation intensity was calculated after data collection and used in analysis only.Due to GPS accuracy limitations, 16 plots in the 2011 data set appeared outside clearcut boundaries; accurate calculation of radiation intensity was therefore not possible and they were excluded from analyses requiring a radiation value.

Data Analysis and Statistics
We chose biomass fraction-the fraction of total biomass accounted for by a single species-as the metric with which to judge species abundance of stems > 2 m in height.Biomass is a straightforward measure of the amount of resources used by a given tree and a good indicator of future tree success.Since we measured stems of all size classes, using density would have biased the results toward making small, numerous species seem dominant, when in fact they were likely suppressed and may never achieve canopy dominance.Species-specific biomass equations from Jenkins et al. (2004) were used to estimate individual aboveground tree biomass from measured diameter.For some species, many potential equations exist.Whenever possible, we selected equations developed on sites near our study locations, covering the small diameter ranges of tree on our sites, and with large sample sizes.Estimates of biomass per hectare for each site were calculated by dividing the total biomass present within a plot by its area and averaging over all plots within that site.For both the 2010 and 2011 data, species dominance of stems > 2 m in height was calculated by averaging species' biomass fractions with other plots subject to the same treatment and within the same clearcut site.Since only four sites were measured during the 2010 season, the data were combined with the 29 sites from 2011 and treated as one data set for the majority of the statistical tests performed.However, since the 2010 season measured those four sites with a high sampling intensity, several statistical tests using only the 2010 data were run and are noted in our results.
For trees > 2 m in height, three treatments were tested-harvest type (WTH or CH), radiation category (low, medium or high) and edge category (within 10 m of clearcut edge or >10 m from the edge).To simplify our graphs, only the most common seven species were analyzed individually.The choice of these seven represented a natural cut-off, with additional species occurring in low enough concentrations to make drawing conclusions difficult.Lesser abundant species (nine in 2010 and eleven in 2011) were grouped and graphically displayed as "other."For trees < 2 m in height and understory plants (measured only in the 2010 season), only harvest treatment was tested.
Non-metric multidimensional scaling (a form of ordination) was used to assess the similarities in species composition between sites of trees > 2 m in height.Input variables were the biomass fraction of each of the seven most common tree species and an eighth variable containing the summed fraction of all other less common species.The values of these variables reflected the composition of an entire site and were calculated by averaging all measured plots of that site.With eight dimensions of variation-one for each measured variable-visua-lizing meaningful patterns is difficult, so ordination was used to reduce the dimensionality to make relationships easier to observe.Ordination is a common dimension reduction technique, and works to minimize the "stress" of a given configuration of replicates.The stress of a given solution is a measure of the difference between the distance of two replicates in the original dimension space and the distance between the same two in the space of reduced dimension.A low stress value shows that ordination has accurately captured the relationships between replicates present in the original higher dimensional space.Ordination works by placing points randomly in a reduced dimensional space and then iteratively moving them in a direction that most reduces the stress.A stable solution is reached when the calculated stress is low enough and stable over several iterations.This final result of a successful ordination places replicates with similar values close together and those that differ greatly far apart.Others have used non-metric multidimensional scaling with success to examine effects of environmental variables and harvest intensity on species composition (Olivero & Hix, 1998;Lee et al., 2005;McDonald et al., 2008).
Ordination was run using PC-ORD v. 4.25 (McCune & Mefford, 1999) and used the Sorensen (Bray-Curtis) distance measure to calculate distance between analysis units.To obtain the appropriate number of dimensions for ordination, preliminary ordinations were run starting at a 6-D solution, stepping down to a 1-D solution.Using a random starting configuration, 20 solutions for each dimension were generated using real data.The solutions of each dimension were examined for low stress and low instability (a measure of how much the stress varies by iteration).In each test run, the 3-dimensional solution represented a stable, low stress solution.Examinations of scree plots (which plot the final stress vs. the number of dimensions) confirmed that higher dimension solutions resulted in small decreases in final stress as compared with 3-D solutions.The 3-D solution for each test was then rerun with 500 iterations and an instability criterion of 0.0001 to be sure the lowest stress, most stable solution was found.For each ordination, the proportion of variance represented by the solution in ordination space (reduced dimension) was calculated and compared with the solution in the original space of higher dimension.
Multi-response permutation procedures (MRPP) were used to determine if tree species composition as a whole varied significantly by each of the three treatments tested-harvest treatment (WTH or CH), radiation category (high, medium or low) and edge category (edge or center).Other studies have used MRPP successfully to determine the degree that environmental variables and harvest intensity affect overall species composition (Olivero & Hix, 1998;McDonald et al., 2008).A multivariate analysis of variance (MANOVA) was used to determine whether 3-D ordination scores differed significantly by treatment.Kruskal-Wallis tests were used to determine whether individual species abundances varied significantly by treatment.For each treatment, seven Kruskal-Wallis tests were run, increasing the likelihood of falsely rejecting the null hypothesis (Type II error) that no treatment effect exists.We therefore applied the Bonferroni correction for multiple tests, which reduced our critical P value to reject the null hypothesis from 0.05 to 0.05/7 = 0.007.Species were additionally grouped by shade tolerance, and Kruskal-Wallis tests were run to determine whether each shade tolerance group varied significantly by treatment.American beech and striped maple were considered shade tolerant, red maple and yellow birch were considered mid-tolerant, and paper birch, pin cherry and bigtooth aspen were considered shade intolerant.As three Kruskal-Wallis tests were conducted for each treatment, we again applied the Bonferroni correction, reducing our critical P value for rejecting the null hypothesis from 0.05 to 0.05/3 = 0.017.Since this adjustment may be overly cautious (Gotelli & Ellison, 2004), we additionally noted non-significant Kruskal-Wallis test results with p < 0.10.Saplings 0 -1 m and 1 -2 m in height were summed separately by species by harvest treatment.This value was then divided by the total area of the plots in which they were found to obtain sapling densities.Seventeen tree species were recorded in the sapling size classes; densities of the nine most common (those with greater than 10 stems tallied total) were reported individually on our graphs.A mixed effects model was created for each species, with the number of tallied stems per plot as the dependent variable.Harvest treatment (either WTH or CH) was used as a fixed effect and a variable uniquely identifying each site was used as a random effect.
Many non-tree understory species were found in only one or two 1m-radius plots, so we focused our analyses on the most commonly occurring twelve species.In addition, any lesser abundant plants occurring in greater than three plots that showed potential differences in abundance by treatment were reported.The proportion of plots that contained each understory plant was calculated and used as a measure of that species' abundance within each treatment.Additionally, the number of species found in each plot (species richness) was calculated and used as the dependent variable in a mixed effects model.Using harvest treatment as a fixed effect and a va-riable uniquely identifying each site as a random effect, the mixed effects model determined whether observed differences in species richness could be attributed to harvest treatment.Unless otherwise noted, tests with a probability value of p < 0.05 were considered to be statistically significant.Taxonomy follows Haines (2011).

Moose Browse and Non-Forested Plots
Plots classified as containing primarily non-forest regeneration or those with current year moose browse damage (referred to hereafter as "stunted") were found to contain significantly different regeneration than those containing unbrowsed forest.F-tests showed that the sets of stunted and unbrowsed plots did not have equal variances, so Welch's t-tests (which do not assume equal variance) were used to determine whether productivity varied.Stunted plots were significantly less productive than unbrowsed forest plots (Welch's t-test; height: t = 11.35,102 d.f., p < 0.001, diameter: t = 4.99, 78 d.f.p < 0.001, biomass: t = 12.68, 129 d.f., p < 0.001).Species composition and sapling abundances were also significantly altered.Since the majority of stunted plots had reduced or no tree cover, on average they contained a greater diversity of understory plants.In an effort to reduce offtreatment variation, stunted plots (44 in the 2010 season and 277 in the 2011 season) were excluded from our subsequent analyses.While these plots are a part of actual regeneration and contribute to the composition of regenerating forests, their composition is not necessarily due to harvest treatment.We will confine our analyses to a subset of the regeneration-forest unhindered by moose browse-to reduce off-treatment variation.In the 2010 study 40 CH and 69 WTH plots meet this criterion, while 358 CH and 180 WTH plots did so in the 2011 study.

Site Variables
We measured 3922 trees in 2010 and 2951 trees in 2011 that were >2 m in height.Trees averaged 4.3 m in height, 2.4 cm in DBH and were densely spaced, with between 7400 and 30,000 trees per hectare (3000 to 12,000 trees per acre).Estimates of aboveground biomass ranged from 6 to 62 metric tons per hectare (3 to 28 tons per acre).Estimated radiation intensity did not vary significantly by harvest treatment (t-test: t-ratio = 0.98, 31 d.f., p = 0.34).Tree species composition in all height classes was dominated by a variety of hardwood species-American beech (Fagus grandifolia), yellow birch (Betula allegheniensis), paper birch (Betula papyrifera), striped maple (Acer pensylvanicum), red maple (Acer rubrum), pin cherry (Prunus pensylvanica) and bigtooth aspen (Populus grandidentata).Members of the Rubus genus (primarily blackberry and red raspberry) were abundant in the understory along with herbaceous plants such as Aralia nudicaulis, Maianthemum canadense, Dennstaedtia punctilobula, and Trillium erectum.

Species Composition Comparison by Harvest Treatment
MRPP showed that harvest type did not have a significant effect on species composition as a whole, and a MANOVA on the ordinations scores confirmed this result (Table 2).Analyses of individual species abundance (summarized in Table 3 and Figure 2) indicated that no species showed a strong affinity for a certain harvest treatment.Red maple and paper birch were slightly more abundant in CH sites, but did not meet our criteria for statistical significance.American beech was found in greater quantities in WTH sites, but this difference was again not strong enough to be deemed statistically significant.Biomass fractions of pin cherry, striped maple, yellow birch and bigtooth aspen did not vary significantly by harvest treatment.When species were grouped according to their shade tolerance, none of the three tolerance groups varied in a statistically significant manner, although all were close (Table 3).When only the intensively sampled 2010 clearcuts were analyzed, none showed patterns that was statistically significant (Table 3).
Ordination on species composition of the combined data set (both seasons) yielded a 3-D solution with a final stress value of 10.2 and a final instability value of 0.00002.The three axes combined explained 91.3% of the variation and were >99% orthogonal to each other.In Figure 3, the size of the circles in the eight graphs related to tree species are scaled by the prevalence of the species noted.Larger circles correspond to that species representing a higher proportion of the total species composition.A site's location remains constant in each graph.These eight graphs are shown to demonstrate that the ordination was successful; plots with similar species composition are, in general, grouped in ordination space.The ninth graph shows the same sites in the same locations, but without any size scaling.Instead each site is given a color based on the type of harvest conducted.
Table 2. MRPP and MANOVA tests showing no difference in overall tree species regeneration composition between whole-tree harvested and conventionally harvested clearcuts.Data from thirty-three 10 -14 year old northern hardwood sites in central New Hampshire and western Maine.

Multi-response permutation procedures (MRPP) Multivariate analysis of variance (MANOVA)
A    Red dots indicate whole-tree harvesting (WTH) while blue dots indicate conventional harvesting (CH).Little or no pattern associated with harvest treatment is evident, suggesting that harvest treatment had no observable effect on the species composition of the regenerating stand.

Species Composition Comparison by Radiation Category
The ordination for radiation intensity category yielded a 3-D solution with a final stress value of 14.2 and a final instability value of 0.006.The three axes combined accounted for explained 85.2% of the variation present originally and were >93% orthogonal.As was the case with ordination by harvest treatment category, no patterns were evident with ordination by radiation category so we have not included ordination graphs.MRPP showed that radiation intensity had no significant effect on species composition, and this was confirmed by a MANOVA on the coordinate values (Table 2).Analyses on individual species confirmed this result as well-no species showed significantly different biomass fractions between radiation categories (Table 4, Figure 4).Grouping species by shade tolerance showed similar results, with no group showing a significantly higher fraction in a radiation intensity category.When only considering data from the intensively sampled 2010 clearcuts, no species showed statistically significant patterns.Two species showed marginally significant (p < 0.10) patterns, with bigtooth aspen in highest concentrations in high radiation plots and red maple highest in plots of medium radiation intensity (Table 4).

Species Composition Comparison by Edge Category
A three dimensional ordination with a final stress value of 11.8 and a final instability value of 0.005 was generated for the 2011 dataset aggregated by edge category.The three axes explained 90.8% of the original variation and were >95% orthogonal.Similar to our earlier ordination analysis, no patterns were evident with ordination by edge category, therefore no edge ordination graphs are included.MRPP performed on the ordination scores showed that edge category had no significant effect on species composition and this was confirmed by a MANOVA on the coordinate values (Table 2).Kruskal-Wallis tests on individual species confirmed this result as well, with no species having a significantly higher biomass fraction in the center or edge regions of clearcuts (Table 5, Figure 5).Red maple, paper birch and yellow birch were nearly significant (p < 0.10), with all three having higher biomass fractions in plots in the center of clearcuts, greater than 10 meters from the edge boundary.Aggregating the species according to shade tolerance reveals that tolerant species have higher abundances in edge plots while mid-tolerant species are found more often in center plots.However, these patterns do not meet our criteria for statistical significance (Table 5).When examining just the 2010 data set, bigtooth aspen and red maple showed marginally (p < 0.10) higher abundances within "center" plots (>10 m from clearcut edge).American beech along with shade tolerant species as a whole showed marginally higher abundances within "edge" plots.

Species Composition of Saplings
Sapling regeneration (<2 m in height) in 2010 was dominated by yellow birch, red maple and American beech, which together accounted for 81% of stems 0 -1 m tall and 92% of stems 1 -2 m tall (Figure 6).Mixed effects models using harvest treatment as a fixed effect and a site identifier as a random effect showed that residue removal did not have a significant effect on sapling density of any species in either size class (Table 6).Abundances of saplings < 2 m in height were not measured in the 2011 study sites.

Understory Composition and Diversity
Twenty-five species of understory plants were identified in 40 CH plots and 30 within 69 WTH plots in the 2010 study.Eleven of the twelve most abundant understory species occurred more frequently in WTH plots than in CH plots (Figure 7).Uvularia sessilifolia was the lone exception, being found in 18% of CH plots but only 9% of WTH plots.Of lesser abundant species, Gaultheria procumbens and Monotropa uniflora were only found on CH plots (on 8% and 10% of plots respectively).Spinulum annotinum was found on 9% of WTH plots but only 3% of CH plots.A mixed effects model showed that average species richness per plot was more likely to be driven by variability between clearcut sites than harvest treatment (Table 6).

Discussion
Our data do not indicate that residue removal following whole-tree harvesting has a significant effect on the species composition or understory richness of regenerating northern hardwood forests.Ordination of both datasets using non-metric multidimensional scaling showed overall composition was not significantly different between WTH and CH treatments as demonstrated by MRPP results and confirmed by MANOVA tests on the ordination scores.We additionally could not detect any differences in overall species composition due to estimated plot radiation intensity or distance to cutting unit edge.
Comparison of how individual species varied by harvest, radiation and edge treatments revealed several patterns of note, but none met our criteria for statistical significance.Since it is possible that our use of the Bonferroni correction is overly conservative (Gotelli and Ellison, 2004) we have noted test results with an uncorrected significance value of p < 0.10 and discuss how they correlate with known ecological characteristics below.Paper birch and red maple both showed higher abundances within conventionally harvested sites, while American beech was more prevalent with WTH sites.Red maple, yellow birch and paper birch all showed slightly higher abundances in "center" plots > 10 m from the clearcut edge.No species showed even marginally significant (p < 0.10) patterns by solar radiation category.
Due to the retention of harvesting slash on CH sites, one might expect increased shading in the understory to favor trees that can regenerate in partial shade.While paper birch is thought of as shade intolerant, its seeds seem to germinate best in partial shade (Perala & Alm, 1990), with one study showing that paper birch seeds germinated best in soil exposed to 45% of full sunlight (Logan, 1965).Paper birch produce very small, winddriven seeds, and the increased shade under CH slash may provide a cool, moist environment in which to germinate.The exposed, highly variable ground temperatures in WTH sites (Proe et al., 2001) may not be as hospitable.The seeds produced by yellow birch are similar in size and germinate under similar conditions to those of paper birch (Perala & Alm, 1990).However, we found no significant patterns with respect to yellow birch abundance across harvest treatment.Paper birch and yellow birch abundances were slightly higher in the center of clearcuts of both harvest types.While this pattern was not statistically significant, it does confirm our knowledge of paper birch as a shade intolerant species that grows best in full sunlight.
Red maple showed higher abundances in conventionally harvested sites and in plots > 10 m from the clearcut edge.However, the vast majority of the red maple stems observed in both field seasons were dense clusters of stump sprouts from cut stems.Thus, observed abundances are most likely due to differences in pre-harvest red maple abundance and probably not reflective of any true harvest treatment effect on regeneration.Bigtooth aspen showed no significant patterns over any of our three treatment variables, but its abundance was likely governed by pre-harvest intensities as well.Bigtooth aspen have small, wind driven seeds, but also reproduce clonally from root suckering.Many of the aspen in our study plots were extremely tall, and likely grew from root suckers from uncut trees outside of the harvested area.In contrast to our study, Hendrickson (1988) found greater aspen regeneration following WTH as compared with CH sites in pine-aspen stands in Canada, attributing the effect to increased forest floor temperatures and sunlight due to the absence of logging slash.
American beech, a shade tolerant species that readily reproduces vegetatively from root suckering, was present in higher abundance in WTH sites.If whole-tree harvesting serves to increase the amount of sunlight present on the forest floor, we would expect it to favor shade intolerant species such as pin cherry and decrease success of shade tolerant species such as American beech.Our results show the opposite, suggesting observed patterns in American beech abundance may be due to the high degree of variability in measured regeneration.
The final two species analyzed individually, pin cherry and striped maple, showed no marginally significant patterns with respect to harvest treatment, solar radiation or proximity to clearcut edge.From knowledge of their silvicultural characteristics, we would expect pin cherry to be prominent in areas with high radiation intensity and striped maple in areas of lower radiation intensity.Our study shows no such patterns, likely an additional testament to the inherent variability present in these regenerating sites.
A mixed effects model comparing sapling abundances showed no effect of residue removal on the abundance of any species.This was likely due to low statistical power, as sapling abundances were only measured in the four clearcut sites in the 2010 field season.Our statistical power was also not strong enough to detect any meaningful differences in richness of understory plants within 1m-radius plots in our 2010 season.We may expect the removal or retention of harvest materials to affect understory non-tree regeneration in one of two ways.If the retention of harvest materials in CH increases the heterogeneity of microclimate conditions on the forest floor, we might expect a greater variety of species to colonize these niches.Alternatively, harvest material retention may shade the forest floor and decrease the amount of available growing space, lowering the abundance and richness of understory regeneration.
Our study attempted to measure the effects of residue removal, radiation intensity and distance to harvest unit edge on species composition, but unmeasured variables also likely played a role in our results.While we attempted to remove the effects of continual yearly moose browse by excluding browsed plots, other factors still remain.Species composition of the previous stand likely had significant impacts on the composition of the regeneration, especially for species with the ability to sprout from shallow roots or cut stumps.Detailed data of this kind were not available for any of the stands measured.
High site variability and low statistical power also hampered our ability to draw precise conclusions, especially when comparing sapling densities and understory richness across harvest treatment.Sapling density mixed effects models indicate that 92% -95% of the residual model variance was due to variability in sapling sizes within a given clearcut.This suggests that, to obtain a more precise estimate on sapling densities < 2 m in height, future studies would be best adopting a more intensive sampling procedure within each treatment area, with >50 plots per cut.On the other hand, 82% of the residual variance in the understory richness mixed effects model was due to variation between clearcuts.In this case, finding and measuring additional WTH and CH clearcuts would be more valuable in improving the model fit than measuring additional plots within each clearcut.

Conclusion
The data suggest residue removal has no measureable effect on the species composition of a northern hardwood regeneration 10 -14 years following clearcutting.Low statistical power and small sample sizes limited our ability to draw conclusions about sapling densities and understory richness underneath overstory trees.Our results suggest the need for additional, broader studies on the effects of residue removal on the ecology of northern hardwood forests as well as long-term monitoring of existing sites.As whole-tree harvesting becomes more common on our landscape, it becomes more vital to determine the effects of residue removal on the species composition and richness of our forests.

Figure 1 .
Figure 1.Location of the Bartlett Experimental Forest (star), site of the four patch cuts measured in 2010.Also shown are the 14 whole-tree harvested (squares) and 15 conventionally harvested (triangles) patch cuts measured in 2011 within central New Hampshire and western Maine.Location within the northeastern United States is shown on the inset map.

Figure 2 .
Figure 2. Biomass fraction in whole-tree and conventionally harvested sites for seven most common species for the 2010 (top) and 2011 (bottom) sites.Error bars show standard deviation.Data from 33 northern hardwood sites across central New Hampshire and western Maine.

Figure 3 .
Figure 3. Study sites plotted in ordination space.Each dot or circle represents a study site from the 2010 or 2011 field season.The X axis represents one of three ordination scores from the 3-dimensional ordination solution; the Y axis represents a second.The third axis is omitted due to difficulty in visualizing data in three dimensions.Data from 33 northern hardwood sites across central New Hampshire and western Maine.

Figure 4 .
Figure 4. Biomass fraction for seven most common species by plot radiation intensity category for the 2010 (top) and 2011 (bottom) sites.Error bars show standard deviations.Light intensity for each plot was estimated and categorized using the spatial light model described in the text.Data from 33 northern hardwood sites across central New Hampshire and western Maine.

Figure 5 .
Figure 5. Biomass fraction for seven most common species by plot edge category for 2010 (top) and 2011 (bottom) sites.Error bars show standard deviations.Plots within 10 meters of the clearcut boundary were considered "edge" plots; all others were considered "center".Data from 33 northern hardwood sites across central New Hampshire and western Maine.

Figure 6 .
Figure 6.Densities of saplings 0 -1 m in height (top) and 1 -2 m in height (bottom) for the 2010 sites by harvest treatment.Error bars show standard deviations.Data from four experimentally harvested northern hardwood sites in the Bartlett Experimental Forest in Bartlett, NH.

Figure 7 .
Figure 7. Proportion of plots containing the most common 12 understory plants by harvest treatment.Data from four experimentally harvested northern hardwood sites measured in 2010 in the Bartlett Experimental Forest in Bartlett, NH.

Table 1 .
Stratified random sampling of the 2010 study sample plots by radiation intensity category a .

Table 3 .
Kruskal-Wallis tests of species abundance by harvest treatment a .
a Tests were run individually on the seven most commonly occurring species as well as on groupings by shade tolerance.Results are shown for a data set containing both 2010 and 2011 sites (left) as well as on the intensively measured 2010 sites only.For entries with a p value < 0.10, the treatment with the highest biomass fraction is noted.Data collected in central New Hampshire and western Maine.

Table 4 .
Kruskal-Wallis tests of species abundance by radiation intensity category a .Tests were run individually on the seven most commonly occurring species as well as on groupings by shade tolerance.Results are shown for a data set containing both 2010 and 2011 sites (left) as well as on the intensively measured 2010 sites only.For entries with a p value < 0.10, the radiation intensity category with the highest biomass fraction is noted.Data collected in central New Hampshire and western Maine. a

Table 5 .
Kruskal-Wallis tests of species abundance by edge category.Tests were run individually on the seven most commonly occurring species as well as on groupings by shade tolerance.Results are shown for a data set containing both 2010 and 2011 sites (left) as well as on the intensively measured 2010 sites only.For entries with a p value < 0.10, the edge category with the highest biomass fraction is noted.Data from central New Hampshire and western Maine. a

Table 6 .
Mixed effects model of abundance of 0 -1 m and 1 -2 m tall saplings and understory plant richness a .
a Dashes indicate insufficient data to perform a given test for a particular species.No tests indicate a significant effect of residue removal during WTH on sapling or understory communities.Data from four twelve-year-old patch clearcuts in the Bartlett Experimental Forest in the White Mountains of New Hampshire.