Interactive effects of chronic deer browsing and canopy gap disturbance on forest successional dynamics

The interaction of browsing by white-tailed deer (Odocoileus virginianus) and canopy gap disturbances may affect long-term tree composition and lead to significant changes in forest structure. We used an individual-based forest gap model (ZELIG) to better understand the aggregate and interactive impacts of these processes on the long-term (200 years) successional dynamics of a mesic deciduous forest. We parameterized ZELIG to: (1) simulate successional dynamics within a temperate deciduous secondary forest typical of eastern North America; (2) simulate browsing impacts by white-tailed deer; and (3) simulate gap-scale disturbance of variable size and frequency. Our estimates of browsing impacts by species were derived from a 20-year, four-hectare deer exclusion study. Model calibration matched observed tree species composition, density by size class, and total basal area (39.92 m ha 1 vs. 37.13 m ha ). Simulated deer browsing had little impact on total basal area over two centuries. However, deer browsing had substantial impacts on community composition, creating a less diverse understory, lower species richness, and decreased abundance of Quercus species, while retaining the dominance of Liriodendron tulipifera. Simulated gap disturbances exacerbated the impacts of chronic deer browsing and these impacts became stronger over time. Our analyses suggest that recent increases in white-tailed deer density within many forests of eastern North America will result in altered community dynamics that persist beyond the sapling level, and that any increases in overstory disturbance frequency will exacerbate


INTRODUCTION
In the secondary forests that cover much of eastern North America, overstory gap-scale disturbance is a dominant mechanism driving succession.Over time, canopy gap formation and replacement creates forests with complex age and size structures, and patchy species compo-sition in the canopy (Runkle 1982, Whitmore 1989).The role of biological drivers in successional processes is not as well understood, but is also potentially important.Herbivory by whitetailed deer (Odocoileus virginianus), for example, can reduce the survival and growth of several woody species and change the dominance rank of tree species at the sapling stage (Rooney and Dress 1997, Rooney 2001, Russell et al. 2001, Côte et al. 2004, McShea 2012).Exactly how deer browsing, in and of itself, may affect long-term successional processes is not well understood.However, it has been suggested that through preferential browsing deer alter sapling establishment rates, and may shift forests onto alternative succession pathways (Augustine et al. 1998).The potential for deer herbivory to affect succession has increased dramatically in recent decades as white-tailed deer populations have increased-in some cases to densities .50deer km À2 -throughout many forests of eastern North America (McShea et al. 1997, McShea 2012).
Forest succession is an inherently slow process that is difficult to observe directly, and even more difficult to manipulate experimentally.Consequently, our understanding of long-term successional dynamics and the consequences of alternative disturbances on succession has benefited from the use of mechanistic, individuallybased forest models.These ''gap models'' have been used worldwide to reliably simulate changes in tree composition and growth (Botkin et al. 1972, Shugart 1984, Shugart 1998, Bugmann 2001, Bugmann et al. 2001, Holm et al. 2012).The primary emphasis of these models is on the response of individual trees to their prevailing abiotic conditions.In particular, these models focus on an individual tree's response to light availability at height intervals from the canopy to the regeneration layer.The secondary emphasis is on how trees modify abiotic conditions for other trees, with an emphasis on estimating the extinction of light through the canopy as a function of leaf area.Other abiotic factors, such as soil attributes, temperature, and disturbances, are variably incorporated depending on the model (Shugart 2002).Typically, regeneration is treated as a stochastic process, with the probability of establishment constrained by the same abiotic factors that constrain growth.Gap models have been used extensively to forecast forest change due to varying types and levels of disturbances, such as windstorms (Mailly et al. 2000), fungal outbreaks (Shugart and West 1977), and deer browsing (Seagle andLiang 2001, Kramer et al. 2003).
Both overstory disturbance and deer browsing have been implicated in shifting successional pathways for forests in eastern North America, but there has been little consideration of how these factors might interact.Any impacts of browsing on succession are likely to be larger in the presence of overstory gap disturbance, as these events tend to accelerate the process of species turnover (Connell and Slatyer 1977).To examine issues related to deer browsing and gap disturbance in a Mid-Atlantic forest, we to used ZELIG, a popular simulator developed in the original gap-model paradigm.ZELIG has been used to simulate forest succession dynamics throughout North America (Cumming and Burton 1993, Lauenroth et al. 1993, Urban et al. 1993, Busing and Solomon 2004), and locally in eastern U.S. and Mid-Atlantic forests (Coffin andUrban 1993, Seagle andLiang 2001).
We incorporated tree regeneration data from a 20-year deer exclusion study into the parameterization of ZELIG to explore the potential successional consequences of chronic over-browsing with and without gap disturbances.This approach allowed us to examine interactions between press and pulse disturbances (sensu Bender et al. 1984) in a way that would be next to impossible in the field.Our deer exclusion study is among the longest running in the region; nonetheless, the observed impacts of deer browse so far exist only within the seedling and sapling age classes.In particular, sapling density is 4.1 times greater where deer have been excluded compared to where deer browse is unrestricted (McGarvey et al. 2013).Given that few individual sapling stems ever reach a forest canopy, the long-term consequences of selective browsing by deer remain uncertain.However, because it is advanced regeneration, and not new establishment, that tends to claim any newly available growing space (McClure et al. 2000, Barker Plotkin et al. 2012), the differences we see in the understory tree community attributable to preferential deer browse have the potential to eventually manifest in the overstory.
In the only similar study of which we are aware, Seagle and Liang (2001) adapted an earlier version of ZELIG to evaluate the longterm impacts of deer browsing (at multiple intensity levels) in an eastern North American riparian forest.This study did not include gap disturbances.Also, because they did not have empirical data on the impacts of browse on tree regeneration, they simulated browsing impacts using a generalized deer preference ranking for each tree species.Their results suggested that deer browsing would reduce total forest basal area, mainly by reducing regeneration and growth of early successional Liriodendron tulipifera.Additionally, greater browsing intensities altered species composition (toward a late successional Fagus grandifolia-dominated forest), thus accelerating the rate of forest succession.However, long term trajectory of succession was not affected (Seagle and Liang 2001).
We hypothesize that browsing pressure will interact with gap disturbances to create alternative successional pathways.Following the theory of ''press-pulse'' events (Bender et al. 1984, Arens andWest 2008), the concurrence of the press disturbance of chronic deer browsing and the pulse disturbance of sudden forest gap creation will exacerbate the effects resulting from browsing alone.By integrating the field data into the gap model's regeneration sub-routine, we addressed two questions: (1) What are the longterm impacts of chronic deer browse on a deciduous forest and its successional dynamics?
(2) How are these impacts affected by the creation of overstory gaps of various sizes and frequencies?

Study site and field sampling procedure
Field data for this study came from the 25.6 ha Smithsonian Institute Global Earth Observatory (SIGEO; www.si.sigeo.edu)forest dynamics plot at the Smithsonian Conservation Biology Institute (SCBI), a 1,295-ha research facility located 3 km southeast of Front Royal, VA, USA (38854 0 N, 78809 0 W).The SCBI-SIGEO plot is located in a mesic, mature secondary mixed deciduous forest, with overstory tree ages ranging from 84 to 124 years (Bourg et al., in press;McGarvey et al. 2013).The forest canopy is predominantly closed, but small gaps frequently open portions of the canopy.The overstory is strongly dominated by Liriodendron tulipifera (tulip poplar), followed by Quercus alba L., Q. rubra L., Q. prinus L. and Q. velutina Lam.(white, northern red, chestnut, and black oak, respectively), Carya glabra Mill.and C. tomentosa L. (pignut and mockernut hickory), Nyssa sylvatica Marsh.(black gum), and Fraxinus americana (white ash).Prominent understory tree components include Asimina triloba L. (pawpaw), Carpinus caroliniana Walter (American hornbeam/ ironwood), Cercis canadensis L. (eastern redbud), and Cornus florida L. (flowering dogwood).
Within the SCBI-SIGEO plot, two 4-ha subplots were used to calibrate the ZELIG forest simulator.The first subplot has been enclosed by a 2.4m wire fence since 1990 in order to exclude deer (i.e., deer exclusion plot).The second subplot is a control plot that was delineated based on similar overstory species composition to the deer exclusion plot (see McGarvey et al. 2013 for details on selection of control plot).In 2008 all woody stems !1 cm DBH in the 25.6-ha SCBI-SIGEO plot were identified to species, measured for DBH, tagged and mapped using the methodology of Condit (1998).For the purposes of this study, we classified stems between 1 and 5 cm DBH as saplings, which matches the sapling size range in ZELIG.In order to answer our larger questions (what are the long-term impacts of chronic deer browse on a deciduous forest and its successional dynamics), the field data from the two 4-ha subplots was crucial to calibrating the forest demographic model ZELIG for browsing and nobrowsing simulations.Results from the 20-year deer exclusion plot and the control plot explain the effect of chronic deer browsing on the understory.In order evaluate if these effects progress into the overstory, ZELIG was then used to simulate long-term projections of species composition, basal area, and successional trajectories for a deciduous forest under deer browsing.

ZELIG simulation model
ZELIG is an individual-based forest gap model that simulates the growth and development of each individual tree (Urban 1990, 2000, Urban et al. 1991, 1993, Holm et al. 2012).ZELIG is based on the original principles of the JABOWA (Botkin et al. 1972) and FORET gap models (Shugart and West 1977).As in many gap dynamic and individual-based models, the main routines include growth, mortality, regeneration, and tracking of environmental conditions.The model simulates forest stands by tracking all trees as they grow, die, and regenerate across many plots (here: 400-m 2 plots, replicated 50 times).The ZELIG framework assumes the maximum po-tential behavior for all forest processes and then constrains these behaviors depending on the resources available.Potential tree regeneration, growth, and survival are decreased depending on the following environmental constraints: light conditions, soil moisture, level of soil fertility resources, and temperature.Specific details on methodical approaches used in the model can be found in Urban (1990Urban ( , 2000)), Urban et al. (1991Urban et al. ( , 1993)), and Cumming and Burton (1993).ZELIG was chosen for this study because of its use and validation in similar forest types (Urban et al. 1991, Coffin and Urban 1993, Larocque et al. 2006, 2011, Elton 2011), and the availability of the versatile disturbance routine added into ZELIG-TROP (Holm et al. 2012).

Model parameterization
We parameterized ZELIG for the SCBI-SIGEO forest plot by updating environmental and climatic parameters specific for the region and species parameters for 26 tree species (Tables 1  and 2).Of the 26 species, six of the common understory species were grouped and labeled as ''understory'' for analysis purposes: Amelanchier arborea (common serviceberry), Asimina triloba (pawpaw), Carpinus caroliniana (American hornbeam/ironwood), Cercis canadensis (eastern redbud), Cornus florida (flowering dogwood), and Hamamelis virginiana (witch hazel).The majority of the species parameters came from the validated FORCLIM gap model of the Eastern North American forests (Bugmann 1994, Bugmann andSolomon 1995).The following seven species were not present in the FORCLIM model: Acer negundo, Ailanthus altissima, Amelanchier arborea, Asimina triloba, Cercis canadensis, Celtis occidentalis, and Hamamelis virginiana.Parameter values for these seven species came from a combination of the following sources: FORET (Shugart and West 1977), TREGROW combined with ZELIG (Elton 2011), USDA/USFS species fact sheets Note: Agemax, maximum age for the species; DBHmax, maximum diameter at breast height (cm); HTmax, maximum height (cm); shape adjustment factor for height; G, growth rate scaling coefficient; DegDMin and DegDMax, minimum and maximum growing degree days; (full parameter explanation found in original ZELIG paper: Urban 1990).
ZELIG required additional species parameters that were not present in the FORCLIM model, and therefore were determined from a variety of sources.The crown shape parameter, used to help determine leaf distribution based on the shape of the crown, was estimated from Gilman and Watson (1993) from USDA/USFS species fact sheets.All species were given the stress mortality probability factor of P S ¼ 0.369, consistent with stress mortality estimates found by Shugart (1984) and Van Daalen and Shugart (1989).Age-related or natural death of tree species was estimated using the natural mortality equation found by (Botkin et al. 1972, Shugart 1984), but modified for lifespan characteristics typically observed in the SCBI-SIGEO forest plot, which was estimated as 10% of individuals reach their species-specific maximum age (AGEMAX), i.e., m i ¼ 2.303/AGEMAX.
The two species-specific parameters relating to sapling regeneration (relative sapling establishment rate, RSER, and Stocking) were estimated using site-specific measurements from the two 4 ha subplots located in the SCBI-SIGEO forest dynamics plot.Stocking was the proportion of area occupied by each species, or species abundance in a given area.Relative sapling establishment rate (RSER) represented the germination potential, in terms of ''potential'' number of saplings that were likely to germinate and reach breast height under optimal conditions.RSER was defined as the potential number of saplings that can be established per square meter, and was used to rank regeneration potential for each species.We calculated two datasets of RSER and stocking values for all 26 species in order to capture the difference in sapling regeneration v www.esajournals.organd abundance with respect to 20þ years of deer exclusion.The first dataset used sapling (1-5 cm DBH) census data from inside the 4-ha deer exclusion subplot and the second from sapling census data inside the 4-ha control subplot.Multi-stemmed understory tree saplings exhibiting a shrubby growth form were excluded from the parameterization (e.g., C. caroliniana).However, these species were included when they grew as single stems.We further parameterized ZELIG for the Mid-Atlantic deciduous forest region by modifying the environmental parameters (Table 3).We calculated mean monthly temperature (0.29-23.828C) and precipitation (3.49-27.62cm) from an on-site weather station, which were used to simulate weather conditions in ZELIG.Additional environmental parameters modified for the study area included latitude, longitude, mean elevation, solar radiation, soil water carrying capacity, and soil wilting point.

Model calibration
ZELIG was calibrated for the Mid-Atlantic region (specifically the SCBI-SIGEO plot) by comparing the observed SCBI-SIGEO data to model outputs, which used plot data for relative sapling establishment rates and stocking percentages.Specifically, calibration was confirmed by comparing ZELIG to census plot data including total basal area (m 2 ha À1 ), species composition and individual species basal area (m 2 ha À1 ), total stem density (N ha À1 ), size class distribution (DBH, cm size class), and relative sapling abundance (%).We calibrated ZELIG for sitespecific growth by initiating the model with a species composition and tree size (DBH, cm), which was an average of the combined control and deer exclusion 4-ha subplots.Census data from the combined 8-ha area were reduced to the ZELIG plot scale (400 m 2 ) such that DBH for each species was taken from regular intervals (deciles) in order to portray a clear representation of the species DBH size range.All calibration simulations were run for 320 years and replicated for 50 independent 400-m 2 plots.The time at which ZELIG began to reach a stable state, and coincided with the observed mature mixed deciduous forest, was seen at year 120 (total basal area and species composition used to determine stable state).All calibration results from ZELIG (i.e., basal area, stem density) were averaged over 100 years (stand age of 120 to 220 years old), and from an average of 50 plots.

Browsing and disturbance treatments
To better understand the long-term effect of deer herbivory on forest succession and growth, we crafted two suites of simulations.The first suite incorporated the effects of deer browsing on tree establishment and the second suite had establishment rates based on twenty years of deer exclusion.ZELIG, like most gap models, ''regenerates'' stems at 2.5-cm DBH.Accordingly, we used the sapling class stage from the browsed and un-browsed census data to simulate the impacts of deer herbivory and deer absence, respectively.Both simulations started with a 120year spin-up period using the deer free establishment rate to replicate the first 120 years of the site, which had low levels of deer herbivory until the mid 1970s.The ''with deer'' simulations were followed by current levels of deer browsing for the last 200 years of simulation (year 120-320).We introduced effects of deer herbivory at simulation year 120 by setting the sapling establishment rate and sapling stocking density to values recorded from the unfenced subplot.The ''without deer'' simulations did not include any deer effects throughout the 320-year simulation period.For these simulations, we set the sapling establishment rate and stocking density to values recorded the deer exclusion subplot.
To better understand the long-term interactive effects of deer browsing and overstory gap disturbances, we crafted ten gap disturbance treatments and applied them to the ''with deer'' and ''without deer'' scenarios described above.The first treatment had zero gap disturbances throughout the 320-year simulation.This treatment simulated a continuous, intact forest canopy, with only natural senescence of large trees creating canopy gaps.Each of the remaining nine treatments began after the initial spin-up period of 120 years.The gap treatments consisted of a disturbance intensity of removing one, two, or four dominant trees at the plot scale (400 m 2 ), at a disturbance frequency of every 5 years (high disturbance), 30 years (moderate disturbance), or 50 years (low disturbance).The trees were removed by randomly selecting a tree with a DBH in the top 10% of standing trees, regardless of species.The range of variability in gap size and gap frequency captured by our simulations was designed to emulate the wide diversity of chronic (e.g., insect or pathogen) and stochastic (e.g., wind burst) disturbances that result in tree death in these forests (Runkle 1985, Yamamoto 2000).

Analyses
We calculated the mean basal area, and standard deviation, by species across the fifty replicates and plotted changes throughout the 320-year simulation.We then compared changes in the relative basal area for each species.We also tested the difference in total basal area with and without deer browsing pressure, calculated as a percent difference between browse (B) and no browse (NB) for both the mean and the standard deviation between replicates.To examine the relationship between levels of gap disturbance and deer browsing, all forest species data (average basal area for all species) were ordinated using nonmetric multidimensional scaling (NMS) in the PCOrd 6.0 software (McCune and Mefford 2011).Therefore in order to examine the changes in community composition from the last 100 years of simulation, we constructed an ordination using the following parameters and settings: relative Sorensen distance measure, varimax rotation, instability criterion of 0.00001, 250 iterations, 50 random starts, and a starting maximum of six possible axes which stepped down dimensionally to an axes solution.The appropriate number of dimensions (axes) was determined by plotting final stress against the number of dimensions using a scree plot.We used a nonparametric multivariate analysis of variance (PerMANOVA) to assess the relative influence of deer browsing and gap disturbance

ZELIG calibration for SCBI-SIGEO
The ZELIG calibration produced quantitatively comparable forest metrics (within 7% difference) to the observed field data, in terms of diameter distribution, stem densities, and basal area (Figs.2-4).We compared total stem density (N ha À1 ) between the ZELIG simulation and the observed SCBI-SIGEO forest averaged from both the control and deer exclosure subplots (Fig. 2).The observed stem density was 1,298 stems ha À1 (dashed black line in Fig. 2).The average ZELIG stem density, after the spin-up period, (last 200 years and over 50 replicates) was 1,391 stems ha À1 , with a standard deviation (grey area) of approximately 676 stems.Towards the end of the simulation, the average observed stem density fell within one standard deviation of the simulated stem density without exception.The model produced stem frequencies that followed a similar pattern in size class distribution to the observed temperate forest (Fig. 3).Once the simulation reached a mature population (around year 120) the size class distribution was a classic reverse-J shape.The size class with the most frequent number of stems was 1-5 cm DBH and decreased for each subsequently increasing size class.ZELIG slightly underestimated the number of larger DBH stems within the 20-30 cm, 30-40 cm, and 40-50 cm size classes, and overestimated stems in the lower size classes (5-10 cm and 10-15 cm).
Model calibration successfully simulated observed total basal area (Fig. 2) and basal area contribution by individual species (Fig. 4).The observed average total basal area (m 2 ha À1 ) between both the control and fence subplot was 37.13 m 2 ha À1 .The average ZELIG basal area (averaged for the last 200 years after the spin-up Fig. 2. ZELIG simulated total stem density (N ha À1 ) 6 SD and simulated total basal area (m 2 ha À1 ) 6 SD, compared to the observed stem density and total basal area from the SCBI-SIGEO forest dynamics plot inventoried in 2008 (dashed lines), averaged from both the control and deer exclosure plot.Reported after forest reached stable state at year 120 and continued until year 320.
v www.esajournals.orgperiod) was 39.92 m 2 ha À1 (SD: 1.75).ZELIG was successful at predicting the top three common species by basal area in the SCBI-SIGEO temperate forest plot: L. tulipifera, Q. velutina, and Q. rubra.In addition to ZELIG successfully predicting the common and rare species, the model was also in close agreement with simulating the observed basal area value for each of the 26 species (Fig. 4).
Simulated forest stands were subjected to browsing based on total sapling density recorded in the SCBI-SIGEO forest control plot.The calibration testing produced sapling abundances that were 47.6% lower than those on the SCBI-SIGEO plot.The observed sapling density from the control plot was 168 saplings ha À1 and mean Fig. 3. ZELIG simulated stem density distribution (N ha À1 ) for 11 DBH (cm) size classes 6 SD, compared to the observed size class distribution (N ha À1 ) of the SCBI-SIGEO forest dynamics plot averaged from both the control and deer exclosure plot.The average simulated size class distribution was reported after reaching the stable state, at 120 years.Fig. 4. Relative basal area (%) 6 SD for 26 tree species from the SCBI-SIGEO forest dynamics plot, averaged from both the control and deer exclosure plot, as well as from the ZELIG simulation at year 120.v www.esajournals.orgtotal sapling density predicted by ZELIG was 88 saplings ha À1 (SD: 15).However, the relative sapling density was more accurate when evaluating each of the 26 species individually, with no significant difference between the observed and ZELIG values per species (two-sample t test, t (49,2.01)¼ 1.27, p ¼ 0.211); we considered this test a more important metric to evaluate community composition.

Browsing treatments
Our simulation results suggested that deer browsing does not notably impact overall basal area.Average simulated basal area over the last 200 years with no browsing was 41.59 m 2 ha À1 (SD: 13.87) versus 39.92 m 2 ha À1 (SD: 13.74) with browsing.However, deer herbivory did play a strong role in determining forest composition, as seen by altered species composition compared to conditions prior to implementing browsing.The PerMANOVA of the three treatment effects  ) for 26 species reported in terms of percent contribution of total basal area for a simulation with no deer browsing occurring during the full simulation (using the deer exclusion input data), and no gap disturbances.Species abbreviations can be found in Table 1.The group ''understory'' species consist of the following species: Amelanchier arborea (common serviceberry), Asimina triloba (pawpaw), Carpinus caroliniana (American hornbeam/ironwood), Cercis canadensis (eastern redbud), Cornus florida (flowering dogwood), and Hamamelis virginiana (witch hazel).v www.esajournals.orgFig. 6. (A) Percent contribution of total basal area for 26 species reported by ZELIG, for a simulation with deer browsing occurring in the last 200 years, after an initial spin-up period void of browsing (up to year 120, dashed white line), and gap disturbances of removing 1 tree at the plot scale (400 m 2 ), occurring at a high frequency of every 5 years, (B) occurring at a moderate frequency of every 30 years, and (C) occurring at a low frequency of every 50 years.(D) Percent contribution of total basal area for 26 species reported by ZELIG for a simulation with no deer browsing occurring during the full simulation (using the deer exclusion input data), and gap disturbances of removing 1 tree at the plot scale (400 m 2 ), occurring at a high frequency of every 5 years, (E) occurring at a moderate frequency of every 30 years, and (F) occurring at a low frequency of every 50 years.Species abbreviations can be found in Table 1.v www.esajournals.orgFig. 7. (A) Percent contribution of total basal area for 26 species reported by ZELIG, for a simulation with deer browsing occurring in the last 200 years, after an initial spin-up period void of browsing (up to year 120, dashed white line), and gap disturbances of removing 2 trees at the plot scale (400 m 2 ), occurring at a high frequency of every 5 years, (B) occurring at a moderate frequency of every 30 years, and (C) occurring at a low frequency of every 50 years.(D) Percent contribution of total basal area for 26 species reported by ZELIG for a simulation with no deer browsing occurring during the full simulation (using the deer exclusion input data), and gap disturbances of removing 2 trees at the plot scale (400 m 2 ), occurring at a high frequency of every 5 years, (E) occurring at a moderate frequency of every 30 years, and (F) occurring at a low frequency of every 50 years.Species abbreviations can be found in Table 1.
indicated that deer browsing explained the most variation in tree community composition, followed by gap frequency, and finally by gap intensity (Table 4).This finding became more pronounced by the last 100 years of simulation.Deer browsing maintained the dominant species, L. tulipifera, increased the relative abundance of all three hickory species, decreased three out of the four Quercus species, and substantially inhibited the six common understory species (Fig. 5).Over the 200-year period, Q. rubra was the only Quercus species that increased in relative abundance with deer browsing (from 2.9 to 3.3 m 2 ha À1 ).Each of these results was measured by change in relative basal area (m 2 ha À1 ).

Overstory species response to disturbance and browsing treatments
Implementing gap disturbances coupled with deer herbivory resulted in a significant impact to overstory species composition and individual species response (Figs.6-8).However, similar to simulations without gap disturbance, total basal area was not notably affected by deer browsing (over a period of 200 years) for any of the gap treatments (Tables 5 and 6).We therefore examine relative basal area throughout the rest of the paper.
Disturbance frequency (every 5, 30, or 50 years) exacerbated the effects of deer browsing more than changes to disturbance intensity (removing one, two, or four trees per plot).Overall, browsing impacts on community composition were greater in a high frequency, high intensity disturbance regime (i.e., every five years, four tree gap) (Fig. 8A, D), than in a low frequency, low intensity disturbance regime (i.e., every fifty years, one tree gap) (Fig. 6C, F).For example, browsing treatments decreased the relative abundance of all Quercus species, with a larger decrease seen during frequent disturbances.During low frequency disturbance scenarios, the effects of deer herbivory on changes to species relative basal area were less severe in three ways: (1) Quercus decreased, but at a slower rate; (2) the decrease in understory species was smaller compared to high disturbance frequencies; and (3) hickories increased in a pattern similar to those seen in high and medium disturbance frequencies, with no major difference connected to browsing (Fig. 6A, D vs. C, F).
Increasing the disturbance intensity above one tree amplified the browsing effects, but did not lead the forest to diverge into new successional pathways .For example, when two trees were removed from the canopy, the decrease in understory species due to browsing was offset due to larger canopy gaps and more light and resources reaching the forest floor.Therefore, understory species experienced a greater increase (Fig. 7) than in the low intensity disturbance scenarios.After removing four trees in each disturbance event, understory basal area increased more significantly (Fig. 8).The presence of deer assisted in retaining the dominance of L. tulipifera with an increase in basal area ranging from 3.0% to 52.2% in the first hundred years of simulation.This included a 7.3-26.8%difference in standard deviation (SD) between the replicates (Tables 5 and 6).The increase in L. tulipifera was more pronounced in the second 100 years of simulation and with increased frequency of disturbance.
Unlike L. tulipifera, deer browsing and disturbance negatively affected the four Quercus species (Tables 5 and 6).Quercus basal area increased from 11.5 to 14.2 m 2 ha À1 in the absence of deer browsing.The addition of deer browsing slightly reduced Quercus basal area from 11.5 to 10.9 m 2 ha À1 over the 200-year simulation.Under these conditions, Q. velutina experienced the largest reductions, followed by Q. prinus, then Q. alba, then Q. rubra.The introduction of gaps lowered Quercus basal area irrespective of browse impacts-but this was because the gaps were often created by removing overstory Quercus trees.Nonetheless, in all disturbance scenarios, the introduction of deer browsing lowered Quercus basal area.The relative difference in basal area ranged from 9% in the small infrequent gap scenario to 277% in the large moderately frequent gap scenario.
When grouped together, the six understory species showed the largest response to deer browsing, with a significant decrease in understory basal area (Tables 5 and 6).Without browsing, the understory basal area increased from 1.7 to 2.2 m 2 ha À1 .The addition of deer browsing reduced understory basal area from 1.7 to 1.4 m 2 ha À1 .Once both browsing and gap effects were included, the relative basal area was reduced by 20% to 130%.Furthermore, the basal ), occurring at a high frequency of every 5 years, (B) occurring at a moderate frequency of every 30 years, and (C) occurring at a low frequency of every 50 years.(D) Percent contribution of total basal area for 26 species reported by ZELIG for a simulation with no deer browsing occurring during the full simulation (using the deer exclusion input data), and gap disturbances of removing 4 trees at the plot scale (400 m 2 ), occurring at a high frequency of every 5 years, (E) occurring at a moderate frequency of every 30 years, and (F) occurring at a low frequency of every 50 years.Species abbreviations can be found in Table 1.area change in understory species was significantly exacerbated with high disturbance intensity (four tree gaps).During these higher intensity and more frequent disturbances, two counteracting processes occurred; increased light to the forest floor allowed the understory species to increase, but deer presence then significantly decreased these species.For the understory species examined, deer browsing most dramatically reduced C. caroliniana (Tables 5 and 6), which would have made up the largest portion of the understory, ranging from 22% to 73%, if deer browsing was prevented.Additionally, at the highest disturbance levels and in the absence of deer browsing pressure, C. caroliniana became the dominant species in terms of basal area, even over L. tulipifera (Table 6).
Variation in overstory disturbance was largely captured by Axis 2. The distance in ordination space between browsed and un-browsed plots increased with the presence of more frequent and more intense gap disturbance.

Species response to gap disturbance and browsing treatments
Within our simulations, gap disturbances exacerbated the impacts of long-term deer browsing and affected overstory species composition.Therefore, predicted increases in disturbances (e.g., storms, invasive pests, fire) due to climate change (Dale et al. 2001, IPCC 2007, Fig. 9. Non-metric multidimensional scaling results for 18 gap disturbance and deer browsing treatments (labeled as disturbance intensity_disturbance frequency_browsing type).Main separation between the nonbrowsing (NB) and browsing (B) treatments, with the second group distinction between gap disturbance frequencies of no disturbance present (ND, black), 5 years (green), 30 years (blue), and 50 year (orange) frequency intervals.
v www.esajournals.orgAdam et al. 2009) will exacerbate deer browsing impacts.The ''press-pulse'' theory of ecological perturbations (Bender et al. 1984) was demonstrated in this study, in that the concurrence of press disturbances (stress from browsing) and drastic pulse events (gap disturbances) significantly altered forest community composition.Moreover, disturbance frequency played a stronger role than disturbance intensity in these alterations.The highest disturbance frequency (i.e., every 5 years) produced the greatest species response, and was a stronger driver of species change.
Even when gap disturbances were added to natural senescence, we found that total basal area was not affected in our forest simulations.A constant basal area was maintained over the 200year simulations with the help of advanced regeneration, increased light levels reaching the forest floor (which increased the survival of seedlings), and species being released from suppression by canopy tree mortality.The forest canopy is predominantly closed for Mid-Atlantic mesic deciduous forests, but small gaps frequently open portions of the canopy (Abrams and Copenheaver 1999).Our model results support advanced regeneration as a main factor allowing the forest to resist dramatic shifts in ecosystem processes and invasion of new species-very similar to the results reported by Barker Plotkin et al. (2012).One of the exotic species present in our study site and included in our simulations is Ailanthus altissima.A. altissima did not represent a large portion of the relative basal area in the observed data or in the modeled data before implementing the browsing treatments and gap disturbances (Fig. 4, AIAL).Upon simulation of gap disturbances and browsing, A. altissima did marginally increase, but not enough to become a dominant species.For example, the largest increase in relative basal area of A. altissima went from 0.06% of the total basal area during the nodisturbance treatment to 0.16% during a frequent disturbance treatment (removing one tree, every five years), suggesting that the forest can resist invasion by exotic tree species.
As previously discussed, the effects from gap disturbances coupled with deer herbivory did produce changes to community composition.For the purposes of this study we only simulated single tree gap disturbance, as these are the most common canopy perturbations within this region (Runkle 1982).Effects of larger disturbances such as clear cuts, or full removal of all species due to disease or pest outbreaks were not examined.If either of these scenarios were coupled with deer herbivory, different successional pathways may have resulted.For example, if deer browsing were prevented in conjunction with gap disturbances, C. caroliniana (ironwood) would become a major species in the understory, both in numbers and basal area.Along with habitat type, density of C. caroliniana was a main driver determining gap regeneration abundance and composition (Millington et al. 2013).With the likelihood of canopy disturbances increasing due to climate change (Dale et al. 2001), C. caroliniana could become an important species with a strong role in regeneration dynamics within the SCBI-SIGEO plot and similar forests.
The increases in hickories (Carya spp.) produced in all of our simulation scenarios have implications for both forest carbon storage and wildlife food resources.Greater Carya basal area may result in higher levels of carbon sequestration due to the greater wood densities of hickories.Iverson and Prasad (2001) showed that oak-hickory forests would expand on average by 34% in area in the eastern U.S. under five modeled climate change scenarios, and Scheller and Mladenoff (2005) also found increases in this forest type with climate change in northern Wisconsin, albeit delayed due to seed dispersal limitations.Given these findings, the larger growth of hickories at the expense of Quercus displayed in our modeling should be considered when assessing forest carbon storage potential in the future.
The shift away from Quercus species in the forest canopy could impact the populations of hard mast-consuming wildlife species.Longterm data on Quercus acorn and Carya nut production within the SCBI-SIGEO plot exclosure and two nearby unenclosed replicates indicate that current mean annual acorn production has been 2.6 times the level of Carya nut production over a 25-year period, measured in kg/ha (Bourg et al. 2013;W. J. McShea, unpublished data).Carya seeds have a more limited suite of dependent species than do Quercus seeds (McShea and Healy 2002).Any reduction is mast production by Quercus would have serious repercussions to wildlife, especially under conditions of continued high deer densities (McShea and Schwede 1993, McShea and Healy 2002, McShea et al. 2007).It is unclear whether and how a drastic reduction in acorn mast here, or in other similar eastern forests, could be compensated for by greater hickory mast as a food resource, and is an issue which needs further study.

Management purposes
This study adds knowledge on the effects of current deer densities at the forest stand level, a topic that has management as well as ecological implications (Russell at al. 2001).Our modeling effort points to two factors shaping succession pathways within eastern deciduous forest; light availability and browsing preferences from deer.Therefore, management of forests and deer are tightly coupled, and any forest management that alters sub-canopy light levels will alter the impact of deer on forest succession patterns.Opening forest canopy through activities such as selective harvest or exurban development will increase the forcing factor of deer populations on the succession process.Therefore, all forest management activities should consider the deer density in the forest planning stage and maintain low densities throughout the early stages of succession.Deer density must be considered for multiple years, as it is the transition from seedling to sapling, not simply seedling establishment, where deer have their strongest selection impact (McGarvey et al. 2013).
The impact of deer will be most obvious in forests with the goal of maintaining or restoring Quercus as a dominant canopy component.Selective harvest within Quercus stands, without implementing deer management plans, will decrease the proportion of Quercus in the subsequent canopy.The maintenance of Quercus forest communities is a major concern due to the decrease in Quercus regeneration from fire suppression, intensive deer browsing, and increases in exotic species, pests and pathogens (Dey 2002, Abrams 2005, McShea et al. 2007, McEwan et al. 2011).Our modeling results show Quercus regeneration failures will be exacerbated by increases in disturbance frequency.Central and southern Appalachian public forests are increasingly managed under a multi-use frame-work that includes timber management and recreation, as well as shifting to private ownership (Alvarez 2007, Drummond and Loveland 2010, McShea 2012).Disturbance regimes differ across each of these land uses but parcelization and ownership transfers of private lands generally lead to increased disturbance rates (Zheng et al. 2010) and thereby increased impacts by deer.Parcelization itself has been correlated with increased deer density, possibly due to limiting the movement of hunters (Lovely et al. 2013).In addition to disturbances associated with land use, natural disturbances are likely to increase with climate change (Dale et al. 2001, IPCC 2007, Anderegg et al. 2013), warranting the need to include deer management plans in all plans to buffer forest systems from climate change impacts (Thompson et al. 2009).
Within mature closed canopy forests, such as those on many National Park Service units or public lands that emphasize forest preservation, the impact of deer on forest succession is not as severe a forcing agent on forest succession as light limitations due to the low disturbance levels.However, in our models, deer impacted understory forest components and understorydependent species will be impacted by deer through their reductions of the volume and diversity of this forest component.Multiple studies have demonstrated a link between forest understory characteristics and understory-dependent birds (McShea andRappole 2000, Martin et al. 2011) or invertebrates (Martin et al. 2010).Land managers concerned with these species groups should maintain low deer densities regardless of disturbance regimes.
Although not directly covered in our analysis, the linkage of disturbance and deer browsing, dictates that land managers consider issues of scale.With the exception of broad insect outbreaks or hurricane damage, the scale of disturbance in eastern forests is far smaller than the range and distribution of deer populations (Runkle 1982, Frelich 2002).States manage and measure deer density at the coarse scale of a county or a region, because the mobility of deer allows them to exploit shifting resources, especially in a landscape that involves agriculture or timber extraction.Management of deer populations should continue to be at the landscape level, to compensate for the mobility of deer and the shifting arrangement of canopy gap disturbances.Most individual landowners will not be able to remove sufficient deer at the stand level to reduce deer impacts on forest succession.

CONCLUSION
While our deer exclusion plot is among the longest running in the region, the observed impacts of deer browse are confined to the seedling and sapling age classes (McGarvey et al. 2013).Through the use of an individual-based demographic forest model, our results showed that the differences we see in the understory tree community do manifest into the canopy.Our study yielded three main outcomes: (1) total basal area was not affected by deer herbivory or gap disturbance; (2) forest tree species composition was altered, and the forest became less diverse, in the presence of deer and common rates of gap disturbance; and (3) gap disturbances exacerbated deer browsing impacts and affected overstory species composition.
Compared to previous deer exclusion studies, we were able to include data from a 20þ year deer exclusion plot into a forest simulation model that operates at the individual tree level (ZELIG).This allowed us to successfully compare the longterm effects of deer browsing to its absence.Additionally, the use of a gap model allowed us to manipulate the gap disturbance regime in relation to deer herbivory on a mesic deciduous forest.Based on our results from a 200-year simulation period, we conclude that in the absence of disturbance, and under current deer population densities, total basal area and ecosystem structure will be stable, new successional trajectories will not occur, and a shift in ecosystem state will not be triggered.When gap disturbances common to eastern deciduous forests are added to natural senescence rates and current deer browsing levels, changes in individual species abundance and community composition will occur.Under these simulated conditions, understory species will be suppressed by deer browsing, but increase due to canopy disturbances (specifically C. caroliniana), L. tulipifera will become more dominant, and Quercus species will undergo further decreases in basal area.Additionally, the strength of these results will increase over time, as exemplified by the larger percent difference in basal area between browsing and no browsing treatments in the final 100 years of the 200-year simulation.We have shown that the effects of deer herbivory are not limited to the understory layer, but also affect overstory community composition over long time periods.
Fig. 1 reports all ZELIG RSER parameter values for saplings in the deer exclusion (no browsing) and control subplot (browsing), and are consistent with values found in McGarvey et al. (2013).

Fig. 1 .
Fig. 1.Relative sapling establishment rate (RSER) for saplings inside the 4-ha fence exclosure (No Browsing) and for saplings inside the 4-ha control plot (Browsing), of the SCBI-SIGEO forest dynamics plot used as a parameterization input variable in ZELIG.
treatments on community composition (measured by basal area).To perform the PerMANO-VA we used the ''Adonis'' function and a Bray Curtis distance matrix within the Vegan Community Ecology Package(Oksanen et al. 2009) in the R statistical language (R Development Core Team 2006).We performed two PerMANOVA tests, one on the average basal area by species for 50 years near the beginning of treatment (year 170-220), and for the last 50 years of treatment (year 270-320).

Fig. 5 .
Fig. 5. (A) ZELIG simulated basal area for 26 species reported in terms of percent contribution of total basal area, for a simulation with deer browsing occuring in the last 200 years, after an initial spin-up period void of browsing (up to year 120), and no gap disturbances.(B) ZELIG simulated basal area (m 2 ha À1) for 26 species reported in terms of percent contribution of total basal area for a simulation with no deer browsing occurring during the full simulation (using the deer exclusion input data), and no gap disturbances.Species abbreviations can be found in Table1.The group ''understory'' species consist of the following species: Amelanchier arborea (common serviceberry), Asimina triloba (pawpaw), Carpinus caroliniana (American hornbeam/ironwood), Cercis canadensis (eastern redbud), Cornus florida (flowering dogwood), and Hamamelis virginiana (witch hazel).

Fig. 8 .
Fig. 8. (A) Percent contribution of total basal area for 26 species reported by ZELIG, for a simulation with deer browsing occurring in the last 200 years, after an initial spin-up period void of browsing (up to year 120, dashed white line), and gap disturbances of removing 4 trees at the plot scale (400 m 2), occurring at a high frequency of every 5 years, (B) occurring at a moderate frequency of every 30 years, and (C) occurring at a low frequency of every 50 years.(D) Percent contribution of total basal area for 26 species reported by ZELIG for a simulation with no deer browsing occurring during the full simulation (using the deer exclusion input data), and gap disturbances of removing 4 trees at the plot scale (400 m 2 ), occurring at a high frequency of every 5 years, (E) occurring at a moderate frequency of every 30 years, and (F) occurring at a low frequency of every 50 years.Species abbreviations can be found in Table1. NB)

Table 1 .
Species-specific allometric and ecological parameters for the 26 tree species used in ZELIG, and found in the SCBI-SIGEO forest dynamics plot.All species were assigned a probability factor of stress mortality of 0.369, probability factor of natural mortality of 2.303, and zone of seed influence of 200.Parameters found from the literature and previous gap models are in normal text, while parameters found from the mid-Atlantic SIGEO plot are in boldface (including the natural mortality probability factor).

Table 2 .
Continued species-specific allometric and ecological parameters for the 26 tree species used in ZELIG, and found in the SCBI-SIGEO forest dynamics plot.Parameters found from the literature and previous gap models are in normal text, while parameters found from the mid-Atlantic SIGEO plot are in bold (including the natural mortality probability factor).The ''control'' represents the 4 ha control subplot with deer browsing present, and the ''fence'' represents the 4 ha deer exclusion subplot.

Table 3 .
Environmental or site parameters used in the ZELIG model for the mid-Atlantic SIGEO forest dynamic plot.

Table 4 .
Summary of results from a permutation-based multivariate analysis of variance (PerMANOVA) which examined the effects of deer browsing, gap frequency, and disturbance intensity by tree removal on tree species community composition (as measured by change in basal area).NS ¼ not significant at a ¼ 0.05.