Complex interactions of deer herbivory, soil chemistry, and competing vegetation explain oak-hickory forest tree regeneration in central Pennsylvania, USA

The root causes of forest tree regeneration failure are difficult to resolve, although numerous studies show ungulate herbivory, soil conditions, and competition from undesirable vegetation as likely contributors. To better understand the relative importance of each issue, we conducted a 7-year manipulative experiment to assess the interactive effects of white-tailed deer (Odocoileus virginianus) herbivory, soil acidity, and competing vegetation on tree regeneration in oak-hickory forests of central Pennsylvania, USA. Outcomes depended on initial tree seedling abundance, and all three factors had significant interactions. At low initial seedling abundance, fencing resulted in the greatest increase, but all treatments had a positive effect on seedling growth and abundance. At higher initial seedling abundance, abundance failed to recover 7 years after herbicide treatment and soil pH was an important predictor. When soil pH was > 4.6 from lime application, seedling growth and abundance in unfenced controls with high initial abundance was comparable to the fenced-only treatment. Competing vegetation, assumed to be a symptom of excessive, long-term deer herbivory, does not seem to be the primary factor limiting tree regeneration in our study area. Ameliorating acid deposition warrants greater consideration as a management action because it could provide long-lasting benefits compared to short-term fence installations.

Methods
2.1 Study area
The study area was located on portions of the Rothrock and Bald Eagle state forests in central Pennsylvania in the Ridge and Valley Physiographic Province. The study area was primarily a contiguous forest dominated by hardwood stands of red oak (Quercus rubra), chestnut oak (Q. montana), black oak (Q. velutina), and scarlet oak (Q. coccinea) at higher elevations and red oak and yellow-poplar (Liriodendron tulipifera) at lower elevations. Red maple (Acer rubrum) was common in the understory and black birch (Betula nigra) was oftentimes an abundant species in early successional stands. The mean tree basal area across the study area from 2014-2021 was 22.40 m2/ha (SD = 11.97, n = 261). The most common plant community was the oak-heath association (Zimmerman et al., 2012) with an oak overstory, and an understory dominated by mountain laurel, huckleberry (Gaylussacia spp.), and blueberry (Vaccinium spp.), especially at higher elevations on dry soils. At the time of full plot establishment in 2014, at least one of these three taxa was present at 85% of plots (n = 261); mountain laurel was present at 73% of plots, blueberry species at 71% of plots, and huckleberry species at 58% of plots.  Furthermore, mountain laurel volume was highly correlated with the summed volume of all three taxa combined (R2 = 0.72). In 2018–2019, deer density was greater in Rothrock State Forest (7.9–9.2 deer/km2) than in Bald Eagle State Forest (2.8–3.3 deer/km2). Sampling sites were 400–700 m above sea level and all plots within a site were no more than 36.6m from the next nearest plot, laid out in a network bounded by a 110m by 52m perimeter (Fig. 2). Study area soils are derived from a number of sandstone, siltstone, and shale and include Typic Hapludults, Typic Fragiudults, Aquic Fragiudults, and Typic Dystrudepts. This study did not require ethical approval and permits to conduct fieldwork were obtained from the Pennsylvania Department of Conservation and Natural Resources (DCNR).

2.2 Experimental design
From 100 potential sampling locations established across the two state forests, we randomly selected 24 vegetation sampling sites. One of these sites occurred within a deer exclosure previously erected by managers, but evidence of deer presence (pellets) and herbivory was noted during our site visit in 2015. At each of the 24 sampling sites, we established 11 circular 14.62-m diameter plots (plots 1–5 in 2013 and plots 6-11 in 2014). At the same time, we nested a single circular 3.30-m diameter microplot 3.66 m due east of each plot center. Microplots were numbered the same as the plot in which they were nested.

We assigned one of seven treatment combinations to microplots 5–11 at each location. The treatments were fence (always at microplot 5), and then herbicide, lime, combinations of fence + herbicide, fence + lime, lime + herbicide, and a three-way combination of fence + lime + herbicide were randomly assigned to microplots 6-11. Plots 1–4 served as controls throughout the study and their microplots did not have any treatment applied. We installed fencing and applied liming treatments in the summer of 2014 (27 May–14 Aug) immediately following a full vegetation inventory (described below). Fences were constructed of a 1.8-m tall polyvinyl deer fence attached to posts and overstory trees encompassing the entire microplot plus a minimum 0.75 m buffer. We maintained fences yearly at the time of data collection, except in 2020. We applied to each microplot pelletized dolomitic limestone (CaMg(CO2)3) (consistent with other studies) at a rate of 5 Mg ha-1 (or 6.7 kg per microplot). We conducted herbicide treatments following the 2015 summer vegetation monitoring season over an 11-day period from 17–28 Aug 2015. We applied 1.9 L of a non-targeted herbicide, broadcast oil-water emulsification tank mix of 1.5% v/v ester triclopyr and 3.0% v/v amine glyphosate combined with a 1.5% v/v emulsifier and 0.05% v/v non-ionic surfactant. All vegetation was treated within the microplot due to the non-selective nature of the herbicide, promoting the competitive release of newly emerging vegetation following application.

2.3 Vegetation and Soil Monitoring
We sampled vegetation at each microplot across all 24 sites before treatment application in 2014 and 2015 and after treatment application in 2016, 2018, and 2021. We considered any live, arborescent, woody species that was ≤2.54 cm diameter at breast height (DBH) a tree seedling and counted individuals of each species by height class (class categories are: 1 = < 0.15 m, 2 = 0.15 – 0.3 m, 3 = 0.3 – 0.9 m, 4 = 0.9 – 1.5 m, 5 = > 1.5 m). DBH was measured 1.37m from the ground. For each taxon not considered a tree species, we recorded a cover in 10% increments, ranging from 0-10% (class 0) to 90-100% (class 9) by the same height class categories used for seedlings. We converted each cover class estimate to its midpoint. Because taxa can overlap at different heights, the amount of cover in each microplot could exceed 100%.

We estimated the volume (m3) of mountain laurel by multiplying the midpoint of the total estimated horizontal cover in each microplot by the upper bound of the height class category. For analysis, we used the mean volume of mountain laurel across survey years after herbicide treatment (2016-2021). We calculated basal area as the total cross-sectional area (m2) of measured tree boles >12.7 cm DBH within the plot area (168 m2) and converted values to m2 per ha. Mean BA (m2 ha-1) across all survey years (2014-2021) was used as a covariate for analysis. For canopy openness, we used the logit-transformed proportion of open sky at 1 m above ground level assessed for each microplot with fisheye lens photos and Gap Light Analyzer.

We collected soil samples by hand from a 0.6 m diameter soil pit to a depth of 40 cm at microplots 1, 5, 6, 7, 8, 9, 10, and 11 at each sampling site in both 2014 and 2016 and at microplots 2–4 when time allowed. We pooled the Oa and A horizons together for analysis because the A horizon was thin and indistinguishable from the Oa horizon. We dried each sample on kraft paper at ~23 °C, ground it using mortar and pestle, and sieved it to 2 mm to remove coarse debris. We measured pH using 0.01 M CaCl2 with a 1:5 soil-to-solution ratio (Ag Analytical Services Laboratory in University Park, PA).

2.4 Statistical analysis
Our experiment with treatments of herbicide, lime, and fencing was identified and executed a priori to any statistical analysis. The statistical analysis was not exploratory but tested the three treatments and interactions. We conducted our analyses using two different response variables because the seedling response to treatments could occur via mortality, recruitment, and growth. Counts of seedlings (Count) only assess abundance via mortality and recruitment, while aggregate height (AggHt) is an index that incorporates both seedling abundance and growth. Aggregate height was calculated as the number of seedlings in each height class multiplied by the midpoint of the height class in meters; for the tallest height class (> 1.5m), we set 1.8 m as the midpoint. The response variables, based on Count and AggHt, were calculated as the change between 2014 and 2021, where a positive change indicated an increase over time.

We used a generalized linear mixed model (GLMM) where the treatments of fencing (FENCE), herbicide (HERB), and lime (LIME) were binary fixed effects. We included Count or AggHt in 2014 (COUNT14 or AggHt14) as a covariate in all models because declines in the response variable could not be greater than what existed on the microplot in 2014. Site was included as a random intercept effect to account for potential similarities among microplots within a sampling site and each observation at the microplot was treated as a random effect to account for overdispersion.

We fit models that included all main effects and interactions, all main effects and 2- and 3-way interactions, all main effects and 2-way interactions, and only main effects. We used Akaike’s information criterion (AIC) to identify the best model, which among the suite of models developed identifies the model that minimizes the number of parameters and maximizes the amount of variation explained by the model, rather than a hypothesis testing approach based on P-values (or confidence intervals) and likelihood-ratio tests. We ranked the top models by AIC to show how well models explained patterns in our data, and in the results, we explained the patterns that existed in the best model.

To further explore the seedling growth responses to liming we used soil pH in 2016 (PH16; after lime and herbicide were applied) as a continuous covariate, because soil pH varied among microplots both before and after treatment. We used AggHt as the response variable in a GLMM with FENCE, HERB, PH16, and AggHt14 as main effects, and included all interaction terms. The site was used as a random intercept effect and each observation at the microplot was used as a random effect to account for overdispersion. We created additional models that included overstory basal area (BA), volume of mountain laurel (KALMIA; after herbicide treatment), or proportion of canopy openness (OPEN) to incorporate the effects of light availability and competition for additional resources on seedling growth. We conducted all analyses in R using the package glmmTMB.