DOI: 10.5061/dryad.z8w9ghxnk
title: Complex interactions of deer herbivory, soil chemistry, and competing vegetation explain oak-hickory forest tree regeneration in central Pennsylvania, USA

description: 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 &gt; 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.
description: 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 (<em>Quercus rubra</em>), chestnut oak (<em>Q.
 montana</em>), black oak (<em>Q. velutina</em>), and
 scarlet oak (<em>Q. coccinea</em>) at higher elevations and
 red oak and yellow-poplar (<em>Liriodendron tulipifera</em>)
 at lower elevations. Red maple (<em>Acer rubrum</em>) was
 common in the understory and black birch (<em>Betula
 nigra</em>) was oftentimes an abundant species in early successional
 stands. The mean tree basal area across the study area from 2014-2021 was
 22.40 m<sup>2</sup>/ha (SD = 11.97, <em>n</em> =
 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
 (<em>Gaylussacia</em> spp.), and blueberry
 (<em>Vaccinium</em> 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 (<em>n</em> =
 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 (R<sup>2</sup> = 0.72). In 2018–2019, deer
 density was greater in Rothrock State Forest (7.9–9.2
 deer/km<sup>2</sup>) than in Bald Eagle State Forest (2.8–3.3
 deer/km<sup>2</sup>). 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). <em>2.2 Experimental
 design</em> 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(CO<sub>2</sub>)<sub>3</sub>)
 (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.
 <em>2.3 Vegetation and Soil
 Monitoring</em> 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 = &lt; 0.15 m, 2 =
 0.15 – 0.3 m, 3 = 0.3 – 0.9 m, 4 = 0.9 – 1.5 m, 5 = &gt; 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
 (m<sup>3</sup>) 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
 (m<sup>2</sup>) of measured tree boles &gt;12.7 cm DBH
 within the plot area (168 m<sup>2</sup>) and converted values
 to m<sup>2</sup> per ha. Mean BA (m<sup>2
 </sup>ha<sup>-1</sup>) 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 CaCl<sub>2</sub> with a 1:5
 soil-to-solution ratio (Ag Analytical Services Laboratory in University
 Park, PA). <em>2.4 Statistical
 analysis</em> 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 (&gt; 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.
description: # Complex interactions of deer herbivory, soil chemistry, and competing
 vegetation explain oak-hickory forest tree regeneration in central
 Pennsylvania, USA
 [https://doi.org/10.5061/dryad.z8w9ghxnk](https://doi.org/10.5061/dryad.z8w9ghxnk) Summarized seedling, competing vegetation, and soil data by vegetation treatment type (lime, herbicide, fence, combination, or control). Data collected for the Deer-Forest Study in Rothrock and Bald Eagle State Forests in central, PA, USA.  ## Description of the data and file structure This is a single Excel file with 21 columns and 262 rows (first-row headers).  COLUMN VARIABLE Followed by description. * STUDY_AREA: State Forest Name where plots are located (Rothrock or Bald Eagle). * UID: Unique ID is a combination of state forest, plot, and subplot ID that denotes where data are collected. The first digit is state forest (1=Rothrock, 2=Bald Eagle), the second two digits are plot ID (01-50), and the last two digits are subplot ID (01-11). * PLOT: Plot number where data are collected (1-50). * SUBPLOT: Subplot (or microplot for cover and seedling data) where data were collected * TRT: Type of treatment applied to the microplot (C = control (no treatment), F = fence, H = herbicide, L = lime). Combination treatments are denoted by multiple letters. * FENCE: Yes or No whether the microplot received a fencing treatment * LIME: Yes or No whether the microplot received lime treatment * HERB: Yes or No whether the microplot received and herbicide treatment * delta.COUNT_6: Change in the aggregate height of seedlings across all height classes from 2014 to 2021 (2021 aggregate height under 6" - 2014 aggregate height under 6"). * AGGHT14: 2014 aggregate height of seedlings summed across all height classes greater than 6" (6-12", 1-3', 3-5', and &gt;5'). Calculated as the midpoint of each height class multiplied by the total count within a microplot, then summed and converted to m. Midpoint for &gt;5 set at 71 inches. * AGGHT21: 2021 aggregate height of seedlings summed across all height classes greater than 6" (6-12", 1-3', 3-5', and &gt;5'). Calculated as the midpoint of each height class multiplied by the total count within a microplot, then summed and converted to m. Midpoint for &gt;5 set at 71 inches. * delta.AGGHT: Change in the aggregate height of seedlings across all height classes greater than 6" from 2014 to 2021 (2021 aggregate height over 6" - 2014 aggregate height over 6"). * COUNT14: 2014 count of seedlings summed across all height classes greater than 6" (6-12", 1-3', 3-5', and &gt;5'). * COUNT21: 2021 count of seedlings summed across all height classes greater than 6" (6-12", 1-3', 3-5', and &gt;5'). * delta.COUNT: Change in the summed count of seedlings across all height classes greater than 6" from 2014 to 2021 (2021 count over 6"- 2014 count over 6"). * PH14: 2014 pH values at each microplot extracted via CaCl. * MN: 2014 KCl Extractable Mn (cmol(+)/kg) from soil analysis at each microplot. * CANOPY: Logit transformed canopy openness from fisheye lens photography in 2016 at each microplot. * BA: Average basal area in m^2/ha at each subplot based on tree DBH data collected in 2014, 2015, 2016, 2018, and 2021. Calculated as the mean of the sums of individual tree basal area in m2 ((3.14159*(DBH_cm/2)^2))/(100^2)) by year multiplied by 59.3 to determine BA per ha (mean(2014_sum, 2015_sum ...)*59.3). * PRE_TRT_KALA_VOL: Total volume of mountain laurel (*Kalmia latifolia*) in 2014 prior to application of any microplot treatments. Calculated as the upper bound of the height class in inches (6", 12", 36", 48", 60", 71") converted to m and multiplied by the midpoint of the percent cover class (5, 15, 25, 35, 45, 55, 65, 75, 85, 95). Values reported as m^3. * PRE_TRT_ALL_COMP_VOL: Total volume of all competing vegetation (mountain laurel (*Kalmia latifolia*), blueberry (*Vaccinium* spp.), and huckleberry (*Gaylussacia* spp.) in 2014 prior to application of any microplot treatments. Calculated as the upper bound of the height class in inches (6", 12", 36", 48", 60", 71") converted to m and multiplied by the midpoint of the percent cover class (5, 15, 25, 35, 45, 55, 65, 75, 85, 95). Values reported as m^3. * POST_TRT_KALA_VOL: Average total volume of mountain laurel (*Kalmia latifolia*) from 2015, 2016, 2018, and 2021 after the application of microplot treatments. Values are the mean values from all years, calculated from the upper bound of the height class in inches (6", 12", 36", 48", 60", 71") converted to m and multiplied by the midpoint of the percent cover class (5, 15, 25, 35, 45, 55, 65, 75, 85, 95). Values reported as m^3. NAs indicate data not collected. ## Sharing/Access information NA ## Code/Software R Core Team. 2022. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. glmmTMB package Brooks M. E., K. Kristensen, K.J. van Benthem., A. Magnusson, C.W. Berg, A. Nielsen, H.J. Skaug, M. Maechler, and B.M. Bolker. 2017. glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. The R Journal, 9(2), 378–400. (https://journal.r-project.org/archive/2017/RJ-2017-066/index.html).