Pattern and Drivers of White-Tailed Deer (Odocoileus virginianus) Herbivory on Tree Saplings across a Plateau Landscape

White-tailed deer (Odocoileus virginianus) populations are impacting long-term regeneration across eastern United States forests. Deer distribution and resulting herbivory patterns are variable across a landscape due to habitat patchiness and topography. It is poorly understood how features associated with topography control deer herbivory. We examined the heterogeneity of deer herbivory as it affects sapling densities across a single forest-type landscape on the Cumberland Plateau. The 1242 hectare site represented a peninsula of tableland that transitioned from developed land to forest and was surrounded on three sides by a bluff, irregularly punctuated by drainages. We examined the spatial variability of deer impacts on sapling density and modeled the relative importance of plateau accessibility features related to topography, proximity to edge, and deer culling as predictors of sapling variation. We used a stratified random design to sample sapling density across the landscape in 2012 and 2015. The intensity of deer herbivory on saplings varied, with the fewest saplings in forests surrounded by residential development. Our model predicted that plateau accessibility measures best determined sapling densities, followed by distance from edge and deer culling measures. Our results suggest that herbivory impacts may not be homogeneous in a contiguous uniform landscape if there are topographic barriers.


Introduction
Across the eastern United States, increasing white-tailed deer (Odocoileus virginianus Zimmermann) populations represent a challenge to natural area managers due to their effects on long-term forest regeneration and consequences for forest structure and composition [1][2][3][4][5][6]. Forest overstory composition is largely determined by the successful recruitment of younger sapling stages, and the differential growth and survival of saplings are good predictors of canopy layer composition [7][8][9]. White-tailed deer have been shown to reduce the growth and survival of seedling and small sapling size classes for a variety of tree species in the eastern deciduous forest [2,[10][11][12][13]. At a landscape scale, when deer population densities are high, deer herbivory can cause saplings to be generally underrepresented in plant communities, allowing sapling density to be used as a measure of deer impacts on forest regeneration within natural areas [2,14,15].
Deer herbivory patterns are variable across a landscape due to nonrandom deer distribution and movement patterns in response to habitat patchiness and topographic features, coupled with selective herbivory preferences [16][17][18][19][20][21][22]. In a landscape representing a mosaic of different habitat patches (such as early/late succession, different plant communities), deer utilize patches with varying intensity as mediated by patch connectivity [21,22]. It has also been well established that deer herbivory is increased around forest edge [23]. Forest edge represents a sudden shift from forest to open habitat and is associated with human activity and land use including roads, residential clearings, farm fields, utility corridors and fencerows [10,[24][25][26][27]. Deer herbivory may also be heterogeneous at the landscape level as a result of landscape features that restrict movement and foraging behavior [18]. Deer typically follow less steep pathways, such as streams [28][29][30], and steep bluff topography may discourage deer movement and herbivory [31]. The relative importance of landscape features associated with topography in controlling the heterogeneity of deer browse is poorly understood [32].
In order to assess the role of topography in driving patterns of deer herbivory across the landscape, it is necessary to control for habitat variation, so that deer impacts can be considered directly reflective of deer density patterns and not a function of browse preference. In this study, we examine the heterogeneity of deer herbivory as it affects sapling densities across a plateau landscape characterized by a single forest type. The objective of this study is to contrast the relative importance of topographic features that may restrict deer movement to that of habitat edge as landscape-level drivers determining the heterogeneity of deer herbivory on tree saplings.

Study Area
Our study took place on a 1242 hectare tract of upland oak-hickory forest located on the surface of the Cumberland Plateau on the campus of The University of the South (UoS) in Sewanee, TN, USA. The upland topography gently varies from convex to concave surfaces with only subtle changes in vegetation that map these gradients across the landscape. The study site is bounded on three sides by a steep bluff transition to a cove community perforated in places by stream drainages. On the fourth side, the forest rapidly transitions to agricultural and residential land use. Elevations of the study area vary from 580-625 m above mean sea level [33].
The forest canopy is generally dominated by oaks (Quercus alba L., Q. coccinea Muenchh., Q. montana Willd., and Q. velutina Lam.) and hickories (Carya glabra Miller, C. pallida (Ashe) Engl. & Graebn., and C. tomentosa (Lam.) Nutt.) along with Acer rubrum L., Nyssa sylvatica Marshall, Oxydendrum arboreum (L.) DC., and Prunus serotina Ehrh. with a distinct understory layer of shrubs and small trees including: Sassafras albidum L., Vaccinium corymbosum L., Vaccinium stamineum L., and Viburnum acerifolium L. White-tailed deer (Odocoileus virginianus) were largely eliminated from the Sewanee area due to severe hunting pressure during the 1930s, but were reintroduced in the 1950s, along with a predator trapping program which targeted bobcats and feral dogs [34]. By the late 1990s, the decline in sapling densities within forested areas of the UoS was attributed to deer herbivory [35]. Deer culls have been performed in selected areas of various sizes within the study area since 2010 for the purpose of management of deer populations in the vicinity of agricultural and residential land use.

Field Methods
To examine sapling density across the landscape, we established 46 stratified random 20ˆ50 m box transects between June and August 2012. Within each transect, there were 10 randomly distributed circular plots (5 m diameter). Within each circular plot, saplings (0.25-3 m height) were recorded by species and counted. The number of saplings per circular plot was added together, so a single value that represented the total number of saplings was generated for each transect. Plot locations were revisited in July 2015, and tree saplings were resampled within a box transect with the same area as the initial survey. The entire 20ˆ50 m box transect was surveyed for sapling species and number. While sampling methodology differed between the two years, the area sampled and plot location remained the same.
Four fenced exclosures within the study area, established by the UoS prior to the study, were used to assess the plant community condition in the absence of deer. Two 16 m 2 exclosures were established in 1996, and two 24 m 2 exclosures were established in 1998 ( Figure 1). Exclosures were fully sampled in 2012 and 2015 after excluding a 0.5 m buffer zone near the fence.

Sapling Density Patterns
The range of variation in sapling density across all plots and exclosures was characterized for 2012 and 2015 with box plots (Figure 2) [36]. The spatial variation of sapling density across the study area was interpolated from sampling points using ordinary kriging in ArcGIS 10.1 [37], with a spherical semivariogram model and a variable search radius that used the nearest 12 points for the interpolation. If kriging interpolated to values outside of the possible range (e.g., negative values), the sapling density is assumed to be zero.

A Model of Sapling Density Based on Landscape Predictors
We evaluated several predictors of sapling density on the landscape that were grouped into three classes: edge predictors (7), deer cull predictors (3), and plateau accessibility predictors (4) ( Table 1).  Edge predictors included distance (m) of a sampling point from the nearest edge in each forest edge category. We defined forest edges as belonging to one of the following categories: residential development, road, cleared road, recent (present in 2012; there was no change in the landscape between 2012 and 2015) clearcut, livestock pasture/agriculture, and lake. The nearest edge predictor denotes the distance (m) to the nearest type of edge in these categories.
Deer cull predictors were generated using deer cull area locations and the number of deer culled (per km 2 ) within each area. In 2014, 14 plots fell within deer cull areas, and 31 were not within deer cull areas. For all plots outside of the deer cull areas, we assumed that the deer cull was equivalent to zero deer. Deer cull area predictors include the number of deer culled per km 2 in 2014, 5-year deer cull, and 3-year deer cull.
Plateau accessibility predictors related to topography include latitude and derived measures which describe the steepness and inaccessibility of movement from the adjacent coves to the top of the Cumberland Plateau. Using the digital elevation model (DEM), we derived potential plateau access zones by determining the minimum elevation around each individual pixel (10 m cell size) along the bluff edge within a 50 m radius and then subtracting this minimum from the actual elevation of the pixel. This result gave a measure of potential plateau access locations along the bluff edge which had low differences between the neighborhood minimum and the local elevation. If the range of difference between the bluff highs and local lows in elevation were 4 m or greater, then the location was defined as a plateau access zone, and its width was evaluated. To characterize the number of nearby access zones in relation to sampling points, we identified the closest bluff edge to each sampling point and evaluated the number and average width of potential access zones 500 m of either side of the closest bluff edge. We evaluated the following predictors of plateau accessibility: distance to plateau access zone, number of nearby (within 500 m of bluff edge) plateau access zones, and average nearby (within 500 m of bluff edge) plateau access zone width.
We used an information theoretic approach [38] to evaluate the association of each of these parameters to our observed landscape patterns of sapling density. To avoid multicollinearity, we evaluated correlation using Pearson's Correlation Coefficient and p-values to determine the likelihood of correlation. If the p-value was 0.05 or less, we did not include predictors within the same candidate model. Using a global model with uncorrelated predictors representing each major category, we evaluated model fit assuming an underlying negative binomial distribution using a chi-square test. We chose the following predictors for the global model that were highly correlated with the other predictors in their group: latitude (correlated with number of nearby plateau access zones (p-value < 0.0001), average nearby plateau access zone width (p-value = 0.0007)), distance to nearest edge (correlated with cleared road (p-value < 0.0001) and clearcut (p-value = 0.027)), 5-year deer cull (correlated with 3-year deer cull (p-value < 0.0001) and average deer culled per km 2 in 2014 (p-value < 0.0001)), and sampling year (Table A1). We also evaluated if variation in sampling protocols in 2012 and 2015 altered model likelihood and found that the global model including time did not perform as well as the global model excluding time (DAIC = 1.6; RMSE = 1.057). To reflect different sampling methodologies, we included time as a random variable in all subsequent models. To separately examine the effect of treatment (exclosure or sampling plot) and year, we used the Kruskal-Wallis test [39].
Models were fit in R using the MASS package [40]. To determine the relative plausibility of each model given our data, we used Akaike's Information Criteria (AIC) [41] adjusted for small sample sizes using the small-sample bias adjustment (AICc) [42]. Model weights were calculated for each candidate model and were ranked from most likely to least likely, based on their weight [38]. Models were only included in the confidence set if they had Akaike weights greater than 10% of the best fitting model's Akaike weight [43]. Model fit was also evaluated using the root mean square error (RMSE). We performed a step-wise protocol by assessing the best-fitting predictor within each category (edge predictors, deer cull predictors, plateau accessibility predictors) before assessing the fit of the best predictors relative to one another. To evaluate the importance and influence of each predictive parameter, we calculated importance weights and scaled estimates to assess the biological relevance of each parameter.

Results
At a landscape level, the sapling density in forested areas varies widely, with the fewest saplings in forest surrounded by residential development and associated edge habitat (Figure 1). The sapling densities in the exclosures were within the range or higher than the highest sapling densities found in the surveyed plots (Figure 1 In comparing the sapling densities within plots for the two years, we found the range in variation increased between 2012 and 2015 (Bartlett's test K 2 = 19.14, p < 0.001) and the median number of saplings per plot decreased (Figure 2). Using the Kruskal-Wallis test, we found that exclosure had a significant effect on sapling density (X 2 = 16.64, df = 1; p-value < 0.001). Year had no effect on sapling density (X 2 = 0.39, df = 1; p-value = 0.53).
The global model for sapling density adequately described the data (X 2 = 39352, df = 89; p-value = 1). The best fitting predictor from each category included latitude (plateau accessibility predictors), distance from residential edge (edge predictors), and 5-year deer cull (deer cull predictors). Akaike model weights suggested that the best-fitting model of sapling density incorporated the variable latitude. Many of these predictors were correlated, so they could not be included in the same model (Table A1). Because of this, we compare the individual effects of each predictor in our models. Including the predictor latitude resulted in a model likelihood that was 52.3 times greater than the next best fitting model which included distance from residential edge, and 65.8 times more likely than the model which included the 5-year deer cull average ( Table 2). Sapling density was positively related to increasing latitude, marginally positively related to distance to residential, and negatively related to 5-year deer cull average ( Table 2). An increase of 0.01 degrees latitude corresponded with an increase of 1.42 saplings. With a 1000 m increase in distance from residential areas, sapling density increased by 1.22. An increase of 0.1 deer culled per km 2 resulted in a decrease of 0.64 saplings (Table 3). Table 2. Akaike importance weights corrected for small sample sizes (AICc) for parameters in confidence sets of linear regression models for prediction of sapling density. The difference in AICc (DAICc) is pictured along with mean square error (MSE) and root mean square error (RMSE).

Model
Parameter

Discussion
Our study showed that within a continuous forest community on the Cumberland Plateau, the effect of deer herbivory on forest regeneration was heterogeneous at the landscape level. Our model results indicated that there are a number of variables that may be driving this pattern, with the strongest predictor being plateau accessibility; latitude was the plateau accessibility measure that produced the best-fitting model of sapling density. There was a clear north-south gradient of sapling density across the landscape, with the lowest density in the southern areas and the highest density in the north. This pattern is not reflected in the exclosure data, which suggesting that this trend is a function of browse and not some other underlying environmental factor. There are a number of factors that could generate this latitudinal difference: (1) hunting pressure and food sources off the plateau; or (2) topographic variation restricting plateau accessibility to deer that is best captured by differences between northern and southern bluffs in the study area. Latitude was highly correlated with most of our more direct measures of plateau accessibility (2/3 of our other plateau accessibility predictors). It is apparent that the north side of the plateau differs from the south side of the plateau in terms of plateau access zone concentration and spacing, which may affect deer movement patterns (Figure 1).
Studies have shown that edge and patch quality generate heterogeneously browsed landscapes [17,32,44]. Our study also finds this to be case, with proximity to residential areas being the second best predictor of sapling density. We found that fixed water availability in the form of five impoundments in the study area did not have an effect on sapling density. Despite repeated deer culls in the areas of lowest sapling density, we found that sapling density in these areas remained low. There are two possibilities here: either the hunting pressure is not sufficient to decrease the number of deer necessary to increase sapling densities, or not enough time has transpired to allow for change in sapling recovery. Studies have shown that it can take decades for tree populations to recover from intense deer browse [4,45,46]. Hardwood species may require more than 10 years of deer reduction/exclusion in order for seedling and sapling densities to increase. This may be in part due to the reduction in seed sources and the past elimination of persistent root-sprouting individuals [47]. The reductions in saplings will eventually lead to reductions in understory and canopy trees [46]. Being on top of the plateau, reductions in densities of canopy trees that are not being replaced make them more prone to windthrow, further affecting seed sources and regeneration potential within tree populations [46].
Our results suggest that when managing habitat that is contiguous and uniform, one cannot assume that herbivory impacts will be homogeneous, particularly if the habitat is constrained by topographic barriers. The southern Cumberland Plateau has been designated as one of the more resilient sites in the southeast in the face of climate change [48], however, managing forest regeneration in the face of climate change may be complicated by deer browse effects on sapling growth and survival. The resiliency inherent in landscapes [48] is dependent on the ability of plant communities and populations to change in regard to microclimatic variation. Restrictions on sapling distribution caused by deer herbivory may preclude this ability for tree species to track microclimatic change over time [9].

Conclusions
At a landscape level, the intensity of deer herbivory varies widely, with the fewest saplings in forested areas surrounded by suburban development and associated edge habitat. On the surface of the Cumberland Plateau, the degree of accessibility between plateau and adjacent cove habitat appears to play a strong role in determining browse heterogeneity across the landscape. This finding suggests that access to the plateau is a primary driver of spatial patterns of deer impacts. Topographic features associated with the bluff line that may alter the foraging behavior of deer are the primary predictors of sapling density. Studies have shown that edge and patch quality clearly generate heterogeneously browsed landscapes, but on the Cumberland Plateau topography plays a critical role in restricting foraging behavior of deer. Our results suggest that it cannot be assumed that herbivory impacts will be homogeneous if the habitat is constrained by topographic barriers such as found on the edge of the plateau.
Accurately quantifying the underlying spatial heterogeneity of deer herbivory and identifying landscape features driving movement patterns will provide critical information for developing effective management strategies to protect biodiversity and increase landscape level resiliency. As deer move according to season, reproductive status, and availability of food and water, their movement patterns are influenced by landscape features, which ultimately dictate the distribution of their herbivory impacts [32,49]. Spatial mapping of deer herbivory has historically been used to quantify population distributions [50,51], but can also represent a management tool to assess the distribution and severity of deer impacts and predict where the impacts will be the most detrimental on the landscape. Understanding where deer impacts are most concentrated can help managers reduce human-wildlife conflicts, such as vehicle collisions, increasing Lyme disease prevalence, and economic losses, while improving forest regeneration and diversity.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

DAICc
Difference in Akaike's information criteria with small sample size correction AIC Akaike