Journal of Ecology
BRITISH
ECOLOGICAL
SOCIETY
Journal of Ecology 2017, 105, 1386-1398 doi: 10.1111/1365-2745.12748
Deer-mediated changes in environment compound the direct impacts of herbivory on understorey plant communities
Autumn E. Sabo*'1, Katie L. Frerker2, Donald M. Waller3 and Eric L. Kruger1
1 Department of Forest & Wildlife Ecology, UW-Madison, 1630 Linden Dr., Madison, Wl 53706, USA; 2US Forest Service, Superior National Forest, 8901 Grand Ave. Place, Duluth, MN 55808, USA; and 3Botany Department, UW-Madison, 430 Lincoln Dr., Madison, Wl 53706, USA
Summary
1. In forests of eastern North America, white-tailed deer (Odocoileus virginianus) can directly affect, via herbivory, the presence, abundance and reproductive success of many plant species. In addition, deer indirectly influence understorey communities by altering environmental conditions.
2. To examine how deer indirectly influence understorey plants via environmental modification, we sampled vegetation and environmental variables in- and outside deer exclosures (10-20 years old) located in temperate forests in northern Wisconsin and the Upper Peninsula of Michigan, USA. We assessed how excluding deer affected understorey community composition and structure, the soil and light environment, and relationships between direct and indirect effects, using non-metric multidimensional scaling (NMDS), mixed linear models and nonparametric multiplicative regression (NPMR).
3. Excluding deer altered sapling communities and several aspects of the understorey environment. Excluding deer from plots with lower overstory basal area increased sapling abundance, decreasing the amount of light available to groundlayer plants. Exclusion also reduced soil compaction and the thickness of the soil E horizon.
4. The composition of understorey communities covaried in apparent response to the environmental factors affected by exclusion. In several common species and groups, E horizon thickness, compaction, openness, and/or total (sapling and overstory) basal area were significant predictors of plant frequency.
5. Complementary analyses revealed that deer exclusion also altered the frequency distributions of several species and groups across environmental space.
6. Synthesis. Deer alter many facets of the understorey environment, such as light availability, soil compaction and thickness of the soil E horizon, which, in turn, appear to mediate variation in plant communities. Those environmental modifications likely compound direct impacts of herbivory as drivers of understorey community change. Our results provide evidence that deer effects on the environment have important implications for forest composition. Thus, we suggest a re-examination of the common assumption that understorey community shifts stem primarily from tissue removal.
Key-words: community composition, non-metric multidimensional scaling, nonparametric multiplicative regression, openness, soil compaction, soil E horizon, soil fertility, tree regeneration, understorey vegetation, white-tailed deer (Odocoileus virgnianus)
Introduction
Ungulates are widely regarded as keystone herbivores and ecosystem engineers in wetland, grassland and forest systems around the world (e.g. reviewed by Hobbs 1996;
""Correspondence author. E-mail: aesabo@wisc.edu
Waller & Alverson 1997; Knapp et al. 1999; Cote et al. 2004; Barrios-Garcia & Ballari 2012). Native and introduced ungulates shape plant community composition and structure directly via frugivory, grazing of grasses, and browsing of woody plants and forbs (Bodmer 1990). Selective browsing has profound influences on vegetation communities (reviewed by Hanley 1997), with deer (members
© 2017 The Authors. Journal of Ecology © 2017 British Ecological Society
Deer effects on forest understoreys 1387
of the Cervidae family) being common browsers in temperate forests (Putman 1988).
In temperate deciduous forests of eastern North America, white-tailed deer (Odocoileus virginianus) overabundance can markedly alter forest composition, structure and function (Rooney 2001; Russell, Zippin & Fowler 2001; Rooney & Waller 2003; Cote et al. 2004). In particular, large increases in deer abundance since the mid-1900s have led to tree regeneration failures (e.g. Rooney & Waller 2003; Cote et al. 2004; White 2012), pronounced declines in the relative abundance of palatable vs. less palatable plants (e.g. Frelich & Lorimer 1985; Tilghman 1989; Miller, Bratton & Hadidian 1992; Anderson 1994; Balgooyen & Waller 1995; Fletcher et al. 2001; Horsley, Stout & DeCalesta 2003; Frerker, Sabo & Waller 2014; Nuttle, Ristau & Royo 2014; Bradshaw & Waller 2016), losses of plant species richness and community diversity (e.g. Rooney & Dress 1997; Horsley, Stout & DeCalesta 2003; Bressette, Beck & Beauchamp 2012; Shelton et al. 2014; Habeck & Schultz 2015), changes in community developmental trajectories (Frelich & Lorimer 1985; Augustine, Frelich & Jordan 1998; Nuttle, Ristau & Royo 2014) and enhanced ecosystem vulnerability to biotic invasions (e.g. Baiser et al. 2008; Knight et al. 2009; Kalisz, Spigler & Horvitz 2014; Davalos, Nuzzo & Blossey 2015; Dobson & Blossey 2015).
A number of these changes result directly from deer her-bivory, which leads to declines in plant vigour, competitive ability and fecundity (e.g. Balgooyen & Waller 1995; Shelton & Inouye 1995; Augustine & Frelich 1998; Fletcher et al. 2001; Webster, Jenkins & Rock 2005), increased susceptibility to other biotic or abiotic stressors (e.g. Long, Pendergast & Carson 2007; Davalos, Nuzzo & Blossey 2014), and in many instances, decreased survival (e.g. Tripler et al. 2005; Krueger et al. 2009). Moreover, the direct effects of her-bivory, along with the impacts of associated deer activities (e.g. defecation, bedding and trampling), may perturb the forest understorey environment, e.g. by altering microclimate (Yamada & Takatsuki 2015), soil biotic and physical properties (e.g. Heckel et al. 2010; Bressette, Beck & Beauchamp 2012; Chips et al. 2014; Shelton et al. 2014), and the overall availability as well as spatial patterning of understorey light (Murray, Webster & Bump 2013) and soil resources (e.g. Seagle 2003; Jensen et al. 2011; Murray, Webster & Bump 2013; Tahtinen et al. 2014).
Results of several studies have indicated that deer-mediated environmental modifications may drive pervasive changes in the demography of woody and herbaceous plant species (Dufresne et al. 2009; Knight et al. 2009; Heckel et al. 2010; Holmes & Webster 2011; Nuttle, Ristau & Royo 2014; Tahtinen et al. 2014; Dobson & Blossey 2015). Here, we adopt a comprehensive approach, asking how such environmental perturbations interact with direct consequences of deer herbivory to reshape woody and herbaceous communities. We take advantage of the comparisons possible between conditions in-and outside deer exclosures located in hardwood-dominated forests in northern Wisconsin and the Upper Peninsula of Michigan, USA (Fig. 1). Previous work with these exclosures
< Ottawa
* Kemp
0 1	100 mi I
1 0	I 100 km
Fig. 1. Map showing exclosure locations in Wisconsin and the Upper Peninsula of Michigan, USA. Exclosures and paired deer-access plots were located in the Kemp Natural Resources Station ('Kemp'), three State Parks in Door County, Plum Creek industrial timberland ('Plum Creek') and Ottawa National Forest ('Ottawa').
(Frerker, Sabo & Waller 2014) revealed that deer exclusion affects populations of common understorey species, increasing tree regeneration and shrub and forb cover, and decreasing the abundance of graminoids. To extend this study, we have coupled vegetation and environmental data from the same exclosures to test these hypotheses:
Hi. Deer exclusion alters the structure and composition of understorey plant communities;
H2: Changes resulting from deer exclusion alter key facets of the understorey environment, such as light availability and soil properties;
H3a: Changes in understorey communities following deer exclusion stem, in part, from changes in the understorey environment identified in H2;
H3b: Alternatively, or additionally, deer exclusion modifies the distribution of understorey plant species across environmental space by allowing palatable taxa to occupy a larger portion of their fundamental niche.
Materials and methods
STUDY LOCATIONS
In Summer 2011, we visited 17 exclosures in upland forests of northern Wisconsin and Michigan's Upper Peninsula (UP) in stands dominated by Acer saccharum with variable,  but often substantial,
© 2017 The Authors. Journal of Ecology © 2017 British Ecological Society, Journal of Ecology, 105, 1386-1398
1388  A. E. Sabo et al.
amounts of Tsuga canadensis in the overstory (Fig. 1, Table 1, Table SI, Supporting Information). In each stand, we sampled vegetation and measured key environmental variables inside the fences and in paired 'deer-access' plots outside the fences. The latter were in close proximity to the exclosures (always within 100 m), and were placed in locations with similar management histories and overstory structure and composition. Additional details about the study locations are provided in Frerker, Sabo & Waller (2014).
VEGETATION SAMPLING
We employed three sampling designs to accommodate the 100 m2, 182 m2 and 2-8 ha exclosures (Fig. 2). Depending on the size of each exclosure, we sampled herbaceous and woody understorey vegetation in 16 to 84 1-m2 quadrats arrayed along transects that varied in arrangement, length and number (Table 1, Fig. 2). We measured two subplots (each 2500 m2) within the 2 ha and larger exclosures and paired deer-access plots. Within every quadrat, we recorded all vascular plant species present. Across larger areas within each plot, we identified the species and measured the diameter at breast height (DBH, 1-4 m) of saplings and trees. In 2 ha and larger plots (Ottawa and Plum Creek locations), we tallied small saplings (0-6 cm < DBH <3-2 cm) within 1-5 m of the entire length of the transects. We also measured DBH of all large saplings (3-2 cm < DBH <10 cm) and trees (DBH >10 cm) within 3 m of the transects (Fig. 2). In the smaller plots (100-182 m2, in Door Co. and Kemp), we tallied all small saplings and recorded DBHs of all larger saplings within each plot. Additionally, we measured the diameter of trees >10 cm DBH that were tallied in variable-radius plots sampled with a prism (basal area factor = 2-3 m2 ha-1) held over the centre of each smaller plot. For plant species authorities, see USDA, NRCS (2016).
Table 1. Site locations (latitude and longitude in decimal degrees), along with soil taxonomie class, number of exclosures (#) per location, installation date, size of exclosures and population density estimated by state natural resources agencies based on harvest data collected in deer management units during 2010-2011 (Michigan Department of Natural Resources and Environment 2010; Wisconsin Department of Natural Resources 2010). Also included are the overstory condition or manipulation just prior to exclosure installation and the total number of quadrats (Quad) sampled in each exclosure or paired deer-access plot. The three state parks (Peninsula SP, Whitefish Dunes SP, and Potawatomi SP) are located in Door County
							Deer		
							density	Overstory	
Location	Lat, Lon	g	Soil Taxonomic Class	#	Date	Size	knT2	condition	Quad
Kemp NRS	45-84, -	89-67	Sandy, mixed, frigid Entic Haplorthods	10	2001	100 m2	2-11	Undisturbed, blowdown	16
Peninsula SP	45-15, -	87-22	Loamy, mixed, active, frigid Lithic Hapludolls; Coarse-loamy, mixed, superactive, frigid Typic Haplorthods; Rock outcrop	2	1991, 1992	182 m2	4-8	Undisturbed	21
Whitefish	44-93, -	8718	Loamy, mixed, active, frigid Lithic	1	1992	182 m2	4-8	Undisturbed	21
Dunes SP			Hapludolls						
Potawatomi SP	44-87, -	87-43	Coarse-loamy, mixed, superactive, frigid Typic Haplorthods	1	1992	182 m2	4-8	Undisturbed	21
Ottawa NF	46-40, -	88-90	Sandy, isotic, frigid Typic Haplorthods; Sandy, mixed, frigid Entic Haplorthods -Sandy, mixed, frigid Alfic Haplorthods -Coarse-loamy, mixed, superactive, frigid Alfic Haplorthods	2	1998, 2002	2 ha	4-8	Gaps created	84
Plum Creek	45-83, -	87-51	Fine-loamy, mixed, active, frigid Inceptic	1	1996	8 ha	12	Thinned	84
timberland			Hapludalfs - Coarse-loamy, mixed, superactive, frigid Aquic Argiudolls						
© 2017 The Authors. Journal of Ecology © 2017 British Ecological Society, Journal of Ecology, 105, 1386-1398
100 m2
2-5 m
2+ ha
5 m
■ ■ I
4 m
25 m
■ Groundlayer quadrat (3 additional consecutive quadrats, not pictured because randomly located, were added to each transect in 2+ ha plots at 2 5 m spacing)
Fig. 2. Vegetation sampling design for the 100 m2 (Kemp), 182 m2 (Door Co.) and one of two subplots in the 2+ ha (Ottawa and Plum Creek) exclosures and paired deer-access plots.
CHARACTERIZING THE UNDERSTOREY ENVIRONMENT
In all exclosure and paired deer-access plots, we measured biotic and abiotic variables to characterize the understorey environment. We estimated openness (a proxy for light availability) at the groundlayer, using fisheye lens photographs obtained at 0-5 m height (Promis et al. 2011), once in the centre of the 100-m2 exclosure/deer-access plots, once in the centre of the first, third and fifth transects of the 182-m2 paired plots, and every 15 m along each of the three transects within the large exclosure/deer-access subplots (Fig. 2). We recorded the diameter of coarse woody debris (CWD) that intersected each transect, and estimated the cross-sectional area per km of transect
Deer effects on forest understoreys 1389
(assuming a circular cross-section). Additionally, we used a rapid assessment protocol to rank exotic earthworm invasion (Loss et at 2013) in one 1-m2 quadrat near the centre of each 100-m2 plot, once at a random location close to each transect for 182-m2 plots, and twice at random locations between transects in each of the two subplots of 2 ha and larger plots. At each earthworm assessment site, we measured the thickness of leaf litter and surface soil horizons in a 20-cm core. We measured soil compaction once per quadrat with a soil compaction meter (Lang penetrometer, Inc., Gulf Shores, AL, USA). We also collected three soil cores (10-cm depth) per plot in areas with representative vegetation, and submitted mixed composite samples (one per plot) to the Soil Testing Laboratory at the University of Wisconsin-Madison. Samples were analysed for total N using the Kjeldahl method, P and K using the Bray and Kurtz P-l procedure, Ca and Mg via ammonium acetate extraction, organic matter (OM) by loss-on-ignition and pH in a 1:1 soil:water slurry. Finally, for every plot, we extracted coarse-scale soil type and texture data from soil survey maps (USDA, NRCS 2013).
STATISTICAL ANALYSES
For every plant species encountered in one or more quadrats across our 17 exclosure/deer-access pairs, we estimated plot-level frequencies as the proportion of plot quadrats containing that species. We also summed the proportional frequencies of species to form groups by plant family, growth form (e.g. vascular cryptogams, graminoids, forbs, etc.), ecological status (e.g. exotic species and late-successional shrubs - Dirca palustris, Taxus canadensis, Acer spicatum, Cornus rugosa, Diervilla lonicera, Symphoricarpos albus, Viburnum aceri-folium and Mitchella repens), and for species found in at least 20% of the plots, by dispersal mode. We determined the dispersal mode of our 26 most common quadrat species by referring to a plant traits database compiled from literature and field specimens (Amatangelo et al. 2014). The sapling (0-6 cm < DBH <10 cm) and tree tallies were converted to plot-level averages for hardwood and conifer sapling and overstory basal area (BA, m2 ha-1). We also generated plot averages for all the environmental data sampled at the quadrat, transect or subplot levels.
To test Hi, that deer exclusion alters community composition, we first employed non-metric multidimensional scaling (NMDS) to characterize variation in sapling, shrub plus tree seedling and herbaceous communities. Prior to ordinating the vegetation communities, we deleted species that occurred in fewer than 5% of our plots, leaving 131 herb, shrub and tree species for the frequency data and 20 species for the sapling BA dataset. Sapling abundance data were square-root transformed and Wisconsin double standardized (McCune, Grace & Urban 2002). We used Bray-Curtis distances to calculate all the species distance matrices. Using the same dataset, we then assessed the significance of differences in treatment centroids, after accounting for variation among location and plot pairs, for sapling, woody under-storey and herbaceous communities using permutational multivariate analysis of variance (PERMANOVA; Anderson 2001). In addition, we tested for treatment differences in variance using permutational analysis of multivariate dispersions (PERMDISP; Anderson 2004). In an extension of the herb and shrub frequency analyses conducted in Frerker, Sabo & Waller (2014), we used mixed models, with plot pair and location as random effects, to determine the significance of deer exclusion on the sapling BA of tree species or genera that occurred in at least 20% of the plots. We include location in these models primarily to account for the potentially confounding influence of variation among locations in the number of plot pairs.
To test H2, that deer exclusion altered key biotic and/or abiotic aspects of the understorey environment, we again used NMDS. This time, rather than ordinating species data, we characterized variation in the collective suite of environmental variables across all exclosure and deer-access plots. Separate NMDS ordinations were generated for the environments of the groundlayer communities (shrubs plus tree seedlings, and herbs) and sapling layer. All 20 environmental variables (Table SI) were included in the groundlayer environment NMDS, while 17 of the 20 variables were included in the sapling environment NMDS (excluding openness, the BA of conifer saplings, and the BA of hardwood saplings because they were confounded with the corresponding vegetation data). Prior to NMDS, we normalized environmental variables by the maximum so that each ranged between 0 and 1, and then calculated Euclidean dissimilarities. We assessed the influence of individual environmental variables on the NMDS axes using Kendall's Tau. Finally, we tested the significance of treatment differences in (i) environmental NMDS scores, after accounting for variation among location and plot pairs, using PERMANOVA and PERMDISP; and (ii) individual environmental variables, using mixed models with plot pair and location as random effects.
We adopted a multifaceted approach to test our two competing hypotheses (H3a and H3b) that contrasts in species composition between exclosures and deer-access plots correspond with differences in the environment (H3a), or that relationships between species composition and environmental gradients differ between exclosure and deer-access plots (H3b). First, we fit environmental vectors to the NMDS ordinations of sapling, shrub plus tree seedling and herbaceous communities. Next, we employed several tactics to examine relationships between the environment and species abundances. We used Kendall's Tau to determine which environmental NMDS axes were correlated with plot-level vegetation NMDS axis scores. We employed mixed models, with location as a random effect, to assess the influences of deer exclusion or environmental variables (as fixed effects) on abundances (frequencies or BA) of individual species that occurred in at least 20% of the plots and species groups. We also calculated divergence in groundlayer species abundance, and divergence in environment, between paired plots by subtracting exclosure values from corresponding deer-access values. Finally, we again employed mixed models (with location as the random effect) to examine relationships between the species and environmental divergences. As detailed in tables in the results, we usually transformed abundances (e.g. with square-root) to reduce heteroscedasticity and/ or improve model fit. Sapling abundances were not tested against deer-mediated environmental variables because they were frequently confounded.
We used nonparametric multiplicative regression (NPMR) to detect whether deer exclusion altered the manner in which individual species or species groups were distributed in environmental space. Unlike linear mixed models, NPMR does not require assumptions about model form and it simultaneously tests a set of predictors and their interactions (McCune 2006). We conducted NPMR using the stepwise free search function in Hyperniche 2.0 (McCune & Mef-ford 2011), with the model type set to local mean quantitative, Gaussian weighting and medium overfitting. The dependent variables again were vegetation NMDS axis scores and plot-level abundances (frequencies or BA) of individual species or groups, while the predictor variables were axis 1 and 2 scores from the corresponding environmental NMDS and a categorical variable indicating deer treatment. Among the resulting models, we deleted all but those with the highest cross-calibration R2 (xR2), which is calculated based on   agreement   between   observed   values   and corresponding
© 2017 The Authors. Journal of Ecology © 2017 British Ecological Society, Journal of Ecology, 105, 1386-1398
1390  A. E. Sabo et al.
predictions using a leave-one-out procedure during each of n model runs, where n is sample size. For species with best-fit models that included the deer treatment term and had a positive xR2 value, we compared the best-fit model against one that included only our environmental NMDS axes 1 and/or 2 as predictors. When the increase in xR2 resulting from addition of the deer treatment to the model was greater than the 3% convention used in Arkle, Pilliod & Welty (2012), we report the difference, recognizing that a small improvement in the cross-validated R2 is more meaningful than comparable changes in traditional R2 (Bruce McCune, personal communication). Tolerance values, model smoothing parameters that increase as species' distributions broaden across the range in predictor variables, and sensitivities, which increase with predictor importance (e.g. sensitivity equals 10 when a 10% change in the predictor produces a 10% change in the response and a sensitivity of zero means no change in the response), were provided for NMDS axis 1 and/or 2 when either or both appeared in the best-fit model. Tolerance and sensitivity were not calculated for treatment, as it was a categorical variable.
Parametric statistics and Kendall's Tau correlations were conducted in JMP Pro 11, while NMDS ordination, PERMANOVA and PERM-DISP were performed using the metaMDS, adonis and betadisper functions, respectively, in R v. 3.2.2 and Vegan package 2.0-10 (Oksanen et al. 2013). The sample sizes for analyses were 34 when examining data across all plots, and 17 for tests of divergence between deer-access and exclosure treatments in plot pairs. Regarding results of mixed models, the significance of treatment effects, relationships between community composition and environment, and their interactions, was evaluated using the sequential Bonferroni procedure (Rice 1989) to account for an increased risk of type I error that accompanies repeated tests of the same hypothesis on many species.
Results
EXCLOSURE EFFECTS ON THE UNDERSTOREY PLANT COMMUNITY
Composition varied considerably (Sabo et al. 2017), both within and across study locations, in communities of tree saplings (Fig. 3a), shrubs plus tree seedlings (Fig. 3b), and herbs (Fig. 3c). Scores of common species along vegetation NMDS axes are presented in Table S2. In support of Hl5 results of PERMANOVA indicated that exclosures differed from deer-access plots in the composition of saplings (P = 0-003) but not groundlayer woody plants (shrubs plus tree seedlings, P = 011) or herbs (P = 0-34). According to results of PERMDISP, dispersion of ordination data about their respective centroids did not differ between treatments for any under-storey community (P > 0-31, data not shown).
For saplings, the treatment difference in the average score along NMDS axis 2 coincided with increases in the BA of Betula papyrifera, Populus spp. and Tsuga canadensis inside exclosures (P < 0-05 after square-root transformation of BA, although none of these increases were significant according to the sequential Bonferroni test). Correspondingly, under less dense overstories, sapling BA was higher inside than outside of exclosures (Fig. 4a). This divergence was driven largely by treatment differences in the abundance of hardwood rather than conifer saplings (Table SI).
0-5 -
-0-5 -
(a)
+
A *
AA AA
A       A á
CM in X 03
CO Q
10 -
00 -
-10 -
(b)
A
A
A A4
+
10 -
00 -
-10 -
(C) A *A
Location
• Door Co. A Kemp
■ Ottawa +  Plum Creek
Deer
• Access
• Exclosure
A ++• \    A +
—r
-1
1
NMDS axis 1
Fig. 3. Plot positions in ordinations of sapling basal area (a, stress = 0-24), tree seedling plus shrub frequency (b, stress = 0-13) and herb frequency (c, stress = 019). Based on results of PERMANOVA, sapling non-metric multidimensional scaling (NMDS) scores differed between exclosure and deer-access plots (P < 001). [Colour figure can be viewed at wileyonlinelibrary.com]
EXCLOSURE EFFECTS ON THE UNDERSTOREY ENVIRONMENT
NMDS ordinations of the sapling and groundlayer environments (Fig. 5) revealed substantial overlap between exclosure and deer-access plots in environmental space. In a test of H2, results of PERMANOVA revealed that, after accounting for variation among location and plot pairs, deer exclusion modified the groundlayer environment (P = 0-02) but not the sapling environment (P = 0-11). Relationships between environmental NMDS axes and individual environmental variables are provided in Table S3.
Openness was higher outside than inside exclosures, on average (P < 0-05, 12-4% v. 7-5%, Table SI), and the contrast became larger as overstory BA (S0) decreased (Fig. 4b). This trend was closely coupled with the aforementioned effect of deer exclusion on sapling BA (Ss), which was amplified under sparse overstories (Fig. 4a). Consequently, while openness declined with increasing B0 (P = 0-009, n = 34), it was more strongly and negatively related to variation in total BA (ST), calculated as the sum of BQ and Bs (P < 0-0001,
© 2017 The Authors. Journal of Ecology © 2017 British Ecological Society, Journal of Ecology, 105, 1386-1398
Deer effects on forest understoreys 1391
15
10 -
(a)
Deer
• Access Exclosure
30
2,1
A		
A		
A	A	(b)
A		Location
		•  Door Co.
Q. O
10
A Kemp NRS ■ Ottawa NF +  Plum Creek
0 15 30 45
Overstory basal area (m2 ha-1)
Fig. 4. Plots of (a) sapling basal area (Bs, m2 ha-1) and (b) openness (O) against overstory BA (B0). Bs was negatively related with Ba in the exclosures (solid line), but not in the deer-access plots (dashed line). Conversely, O was negatively related with Ba in the deer-access plots (solid line), but not in the exclosures (dashed line). Mixed models, with location as a random effect, and Ba, exclosure treatment (D), and their interaction as fixed effects, were based on data from all 34 plots, and response variables were square-root transformed: (a) Bs05 = 2-87 - 0-048ßo - 1-59D + 0-039BcA (b) O05 = 2-61 - 0003ßo + 1-88D - 0-054BoD. All fixed effects and their interactions were significant (P < 005) except for B0 (P = 0-81) in the openness model. [Colour figure can be viewed at wileyonlinelibrary.-com]
n = 34, data not shown). Correspondingly, after accounting for variation due to SQ, Sx increased significantly in exclosures. Specifically, in a mixed model with location as a random effect and B0 and exclosure treatment (D) as fixed effects: ST05 = 2-61 + 0-57Bo°5 - 1-91D + 0-33Bo°5A where D = 0 and 1 for exclosures and deer-access areas, respectively (P < 005 for both fixed effects and their interaction, n = 34, data not shown).
Additionally, based on mixed models with location and plot pair as random effects and treatment as the fixed effect, soil compaction was 9% lower (P = 0-003, Table SI) and soil E horizon was 27% thinner (P = 0-03, Table SI) in exclosures compared with deer-access plots. Soil organic matter and nutrient concentrations, pH, A and O horizon thickness, volume of CWD and earthworm invasion scores did not differ between treatments (all P > 010 in mixed models, Table SI). We did not test for treatment differences in soil texture (i.e. % sand, silt or clay) because, for these data, we relied on coarse-level mapping rather than direct field sampling in our plots.
RELATIONSHIPS BETWEEN UNDERSTOREY COMMUNITIES AND ENVIRONMENT
Patterns of community composition correlated with several aspects of the understorey environment (Table 2). Environmental variables were more closely related to woody and herbaceous groundlayer communities than to the sapling community. Of the environmental variables that differed by treatment, thickness of the soil E horizon was correlated with NMDS axis scores for all three communities and openness was correlated with scores for the two groundlayer communities (Table 2). Additionally, axis 1 of the environmental NMDS was positively correlated with axis 2 of the sapling vegetation NMDS (P = 0-03) as well as axis 1 of both the woody and herbaceous groundlayer NMDS (P < 0-0001),
0-5
CM
co CO Q
-0-5
(a)
A*. +
(b)
Deer
• Access Exclosure
Location
• Door Co. •     4 Kemp
■ Ottawa + Plum Creek
-1
1 -1 0
NMDS axis 1
Fig. 5. Ordinations of environmental conditions across all deer-access and exclosure plots. To depict the understorey environment for tree saplings (a), we conducted non-metric multidimensional scaling (NMDS) of 17 plot-level abiotic (soil surface horizons, chemistry, texture and compaction) and biotic (earthworm invasion rank, CWD and basal area of trees) environmental variables. Stress = 014. To depict the groundlayer environment for herbs, shrubs and tree seedlings, we generated an NMDS of 20 plot-level environmental variables (b), which included all variables from (a) with the addition of openness, the basal area of hardwood saplings and the basal area of conifer saplings. Stress = 012. Based on results of PERMANOVA, NMDS scores of the groundlayer environment differed between exclosure and deer-access plots (P < 005). [Colour figure can be viewed at wileyonlinelibrary.com]
© 2017 The Authors. Journal of Ecology © 2017 British Ecological Society, Journal of Ecology, 105, 1386-1398
1392  A. E. Sabo et al.
Table 2. Environmental loadings that were correlated to non-metric multidimensional scaling (NMDS) ordinations of sapling, seedling and shrub and herb communities (Fig. 3) by significance category
Significance of environmental correlates
Community	P < 0001	P < 001	P < 005	P < 010
Saplings	Hardwood overstory basal area		N, % organic matter, Ca, Mg, % sand, % silt, thickness of E and A horizons	O horizon thickness
Seedlings and	N, pH, Ca, Mg, K, % organic	Openness, % silt		Conifer sapling basal
shrubs	matter, % sand, thickness of E and A horizons, hardwood overstory basal area			area, O horizon thickness
Herbs	N, pH, Mg, % organic matter, % sand, thickness of E horizon	Ca, openness, hardwood overstory basal area, thickness of A horizon	% silt	K
while axis 2 of the environmental NMDS was positively correlated with axis 2 of both the woody and herbaceous ground-layer NMDS (P < 0-04) (data not shown).
In support of H3a, results of mixed models reveal that the abundances of numerous species and species groups (46 of 55 tested) covaried with one or more environmental variables that were modified by the exclosure treatment (Tables S4-S7). For example, graminoid frequency increased in more compacted soils (P < 00001, Fig. 6a), while Liliaceae frequency declined as the soil E horizon thickened (P = 0-003, Fig. 6b). Moreover, the frequencies of 15 species and groups (n = 34, Table 3) were explained by mixed models that included additive combinations of compaction, E horizon, openness and/or total BA. When included as a fixed effect in these additive models, deer exclusion was never significant (P > 0-16). Divergence between the paired plots (deer-access minus exclosure) in the frequency of several common species and groups correlated with the corresponding paired plots' divergence in environmental conditions (Fig. 7, Table S8). In several cases (e.g. Carex arctata and common herbs with unassisted dispersal), these correlations matched the significant overall relationships we observed between species/ group frequency and environmental conditions (Tables 3 and S4-S7).
Despite indications that shifts in composition are linked to changes in key environmental conditions, we also found evidence, albeit limited, that excluding deer altered relationships between understorey composition and environment, supporting our final hypothesis (H3b). Across environmental space, the abundance of four species or groups differed across the fence, as was evidenced in the NPMR results for Betula papyrifera saplings, several shrubs and one native forb (Trientalis borealis) (Table S9). The genus Lonicera, with both native and exotic species, was generally more common in deer-access plots (Fig. 8a). Higher abundances of B. papyrifera, late-successional shrubs and T. borealis, were observed across a broader environmental space in exclosures compared to deer-access plots (Fig. 8b-d, Table S9). For most species, however, distributions across environmental space were best explained with models containing only the two environmental NMDS axes (Table S10).
Discussion
Although our study focused on temperate forests dominated mostly by northern hardwoods, there was considerable variation, both within and across exclosure sites, in key ecosystem properties, including edaphic conditions, forest structure and
Deer
• Access Exclosure
A. 'A*
(b)
A
Location
•  Door Co.
A Kemp ■ Ottawa + Plum Creek
r - _,
0-15 0-2 0-25 0-3
Soil compaction (MPa)
3 6 9
E horizon thickness (cm)
12
Fig. 6. Plots of (a) graminoid frequency against soil compaction, and (b) Liliaceae frequency against thickness of the soil E horizon, across all 17 exclosures (orange symbols) and deer-access plots (black symbols). Mixed models, in which frequency was square-root transformed, yielded significant trends (dashed line) for graminoids (P < 0001) and lilies (P = 0003). In neither case was the effect of deer treatment, or its interaction with the continuous independent variable, significant (P > 01). Thus, a single trend was fit to all 34 data points. [Colour figure can be viewed at wileyonlinelibrary.com]
© 2017 The Authors. Journal of Ecology © 2017 British Ecological Society, Journal of Ecology, 105, 1386-1398
Deer effects on forest understoreys 1393
Table 3. In additive models, environmental factors affected by deer exclusion, including thickness of the soil E horizon (E, cm), soil compaction (C, MPa), openness (O, %) and total BA (BT = overstory BA + sapling BA, m2 ha-1) each explained a significant amount of variation, across exclosure and deer-access plots, in the frequencies of several groundlayer species and species groups. Mixed models, which included location as a random effect, and soil compaction and one of the other three environmental variables as fixed effects, were based on data from all 34 plots. Collinearity precluded models with more than two environmental variables and certain variable combinations, E + BT, E + O, or BT + O. Frequencies were square-root transformed to decrease heteroscedasticity. Below, we provide intercepts and fixed effect coefficients for models with significant effects and, for species with more than one significant model, corrected Akaike's information criterion (AICc)
Species/group
Model parameter coefficients
AICc
Tree seedlings
Abies balsamea
Acer rubrum
Pinus strobus Tilia americana Herbs
Carex arctata
Dryopteris carthusiana Maianthemum canadense
Oryzopsis asperifolia
Taraxacum officinale Trientalis borealis
Veronica officinalis Abiotic dispersal Ingestion dispersal Exotics Graminoids
-0-21 -0-59** -0-57** -0-04 -0-22 0-22
-0-59***
-0-54***
-0-74*
019
005 -0-98** -0-38 -0-41** -0-27 -0-36 -0-84** -0-44** -0-36 -0-39 -0-27 -0-85*
-0018**
-0013* -0013**
-0011* 0010*
-0-032**
-0-02**
-0015**
-0-030***
-0011*
-0-035** -0-023*
1-53* 2.03**
1- 45*
-0-98**
3.32**** 3-63**** 3-05**** -1-90*** 3-65***
2- 38* 2-81*
2-84*
2-73**
2- 99**
3- 39**
2- 59***
3- 67** 2-97**
4- 36** 7.47**** 8.44****
-0 009***
-0-017**
-0016**
0006** 0-008****
0003**
0006***
0-017****
0-01***
0009**
0010* 0-017****
0010**
0-7 4-0 -4-1 3-4
-12-6 -120 -9-2
200 30-4 34-1 15-3 17-7
21-3
24- 0
25- 2
36- 9
37- 3
Asterisks denote significance of individual model coefficients: *P < 010; **P < 005; ***P < 001; ****p < 0001.
disturbance histories. Thus, detection of a consistent exclusion-mediated shift towards a more abundant and diverse sapling community, supporting IJ4, is notable. Regarding the absence of a discernible shift in the groundlayer, studies demonstrating an unambiguous effect of deer on the overall composition of herbaceous communities are particularly uncommon (Mudrak, Johnson & Waller 2009; Jenkins et al. 2014; Pendergast et al. 2016). Signatures are often confined to variation in individual or population-level attributes (e.g. Anderson 1994; Kirschbaum & Anacker 2005; McGraw & Furedi 2005; Frerker, Sabo & Waller 2014; Dobson & Blos-sey 2015). Habeck & Schultz (2015) identified several potential causes for the paucity of reported impacts on community-level metrics of herbaceous species composition, including legacy effects of chronic deer overabundance leading to slow groundlayer responses. A growing body of evidence provides compelling support for this tenet (Royo et al. 2010; Tanent-zap, Kirby & Goldberg 2012; Nuttle, Ristau & Royo 2014; Pendergast et al. 2016). We propose an additional possibility:
deer effects on community composition are often partially obscured by complex species responses to the array of influences deer exert on the understorey environment.
Thus, we encourage further exploration of deer impacts on individual environmental variables, and their subsequent indirect effects on community composition. In our study, exclo-sures clearly affected both biotic and abiotic environmental factors (H2), including sapling abundance and light availability. Because saplings compete for soil resources as well as cast shade, increases in their abundance within exclosures (e.g. Kuijper et al. 2010; Kain et al. 2011; Tanentzap et al. 2011; Bressette, Beck & Beauchamp 2012) could have strong cascading effects on groundlayer vegetation. Notably, while it was not a focus of our hypotheses, we also observed, as have others (e.g. Kuijper et al. 2010), that variation in local environment can exert an important influence on understorey responses to deer. In particular, deer exclusion impacts on the sapling community varied substantially across the locations we studied, and were especially pronounced in stands with
© 2017 The Authors. Journal of Ecology © 2017 British Ecological Society, Journal of Ecology, 105, 1386-1398
1394  A. E. Sabo et al.
0-7
-0-7
(a)	(b)	
		
	• V. •	•
		• •
•		
-5 0 5
A E horizon thickness (cm)
0 20     -20 0 20
A Openness (%) A Total BA (m2 ha-1)
Fig. 7. Treatment divergences (A) in species or group frequency plotted against corresponding divergences in environment: (a) A Viola spp. frequency against A E horizon thickness, (b) A Carex arctata frequency against A openness, and (c) A ingestion-dispersed herbs against A BT (total basal area). To calculate divergence (A), exclosure values were subtracted from deer-access values for each of the 17 plot pairs. Trends were significant in all cases (P < 0 04 in mixed models with location as a random effect, Table S8). [Colour figure can be viewed at wileyonlinelibrary.-com]
10-C
Fig. 8. Variation in species abundances (basal area [m ha~ ] or frequencies) across our ordination-based characterization of environmental space (NMDS axes 1 and 2 from Fig. 5). Data are values from each of the 34 exclosure and deer-access plots. Circle size denotes basal area or frequency and '+' marks the location of abundance-weighted centroids. The genus Lonicera (native and invasive shrubs) was more frequent in deer-access plots (a), while Betula papyrifera saplings (b), late-successional shrubs (c) and Trientalis borealis (d) were more abundant with broader distributions across environmental space in exclosures compared to deer-access areas. [Colour figure can be viewed at wileyonlinelibrary.com]
lower canopy coverage (less overstory BA), as was also reported by Wright et al. (2012) in New Zealand forests with introduced ungulates. These results highlight the context dependency of deer impacts on forest systems (Horsley, Stout & DeCalesta 2003; Tripler et al. 2005; Tremblay, Huot & Potvin 2007; Dufresne et al. 2009; Collard et al. 2010).
Several other studies have also assessed the effects of deer exclusion on understorey light environments. Similar to our study, light availability was greater in deer-access plots than 20-year-old exclosures corresponding with exclusion-mediated increases in understorey tree density (Tanentzap et al. 2011). In another case (Nuttle et al. 2011), leaf area index was still lower in controls due to species composition differences, despite tree basal areas converging after fences were removed. In contrast, two studies concluded that light availability was not significantly different between exclosures and deer-access areas (Holmes & Webster 2011; Murray, Webster & Bump
2013). Such results could reflect the fact that these studies differed from ours in measurement protocols, exclosure ages, overstory composition and structure, and stand disturbance histories, any of which could affect light levels or their relation to deer exclusion.
Our study corroborates findings by Heckel et al. (2010) and Shelton et al. (2014) that soil compaction is lower in exclosures than in adjacent deer-access areas. This may reflect the direct effect of deer trampling, as we found no differences in leaf litter thickness or fine root abundance (data not shown). Heckel et al. (2010) surmised that recovery from trampling could account for the similar effect they saw in 15-year-old exclosures but Shelton et al. (2014) questioned the importance of trampling in driving differences observed just 2 years after exclusion, suggesting instead that deer exert indirect effects like reducing mycorrhizal activity (or fine root abundance).
© 2017 The Authors. Journal of Ecology © 2017 British Ecological Society, Journal of Ecology, 105, 1386-1398
Deer effects on forest understoreys 1395
We also observed that excluding deer reduced the thickness of the E horizon in forest soils. This appears to be the first report of this soil effect. Characteristics of the E horizon and soil podsolization are shaped by a complex suite of climato-logical, geological and biological factors (Sanborn, Lamon-tagne & Hendershot 2011). Thus, exclusion impacts could result from the effects deer have on vegetation composition and nutrient inputs. Thinner E horizons may be associated with increases in hardwood sapling BA, as podsolization tends to be negatively related to the abundance of hardwood species relative to conifers (N0rnberg, Sloth & Nielsen 1993; Sanborn, Lamontagne & Hendershot 2011). Without deer depositing urea and faeces, soils within exclosures may also experience less leaching of soil dissolved organic carbon (Chantigny 2003). However, total N concentrations were similar in- vs. outside exclosures, as has been noted by others (e.g. Hobbs 1996; Pastor et al. 1998; Persson, Danell & Bergstrom 2000; Seagle 2003).
Several lines of evidence support our basic premise that responses of understorey species to deer exclusion often reflect the effects deer have on environmental factors such as light availability, soil compaction and E horizon thickness (H3a). The frequencies of several species covaried with these key environmental variables. This covariation is also consistent with the ecology of these species, such as the shade tolerance of forest herbs, Maianthemum canadense and Trientalis borealis, the intolerance of Rubus idaeus, and known demographic responses to these environmental factors. Godefroid & Koedam (2004) also reported that the abundance of several congeners of the graminoid taxa present in our study increased with more soil compaction at their sites. Furthermore, increased soil compaction is implicated as an important non-consumptive consequence of deer contributing to declines in growth and fecundity of several forest groundlayer species (Heckel et al. 2010).
While few studies have explored how the E horizon affects understorey communities, Wilde & Leaf (1955) found strong relationships between the abundances of many species (including several studied here - Table S4) and the extent of podsolization in temperate hardwood-conifer forests. The thickness of the E horizon may reflect the extent to which organic matter, base cations and other nutrients in upper soil horizons are depleted via leaching of organic acids (Lund-strom, van Breemen & Bain 2000). Indeed, across our sites, E horizon thickness was inversely correlated with nearly every measure associated with soil fertility (e.g. OM, N, P, Ca, Mg, pH, R2 for linear relationships ranging from 0-13 to 0-28, results not shown). Therefore, E horizon thickness may provide an integrative index of fertility in the rhizosphere that, in turn, appears to drive understorey community composition (Fig. 3; Hutchinson et al. 1999; Burton et al. 2011; McEwan & Muller 2011).
Excluding deer had fairly modest effects on the environmental variables, we measured (9-18% of the observed range for each variable, Table SI), which could reflect the relatively young age of these exclosures (i.e. that deer effects on soils persist long after exclusion). Nevertheless,  there was a
discernible treatment difference in the environmental ordination space for the groundlayer community (Fig. 5) and excluding deer had important effects on community and environmental properties known to affect the distribution and abundance of species (Figs 6-8, Tables 3 and S4-S8). In combination over several years, subtle environmental changes appear to have substantial cumulative impacts on understorey communities.
Lastly, we found some support for the hypothesis (H3b) that deer can modify differences between the realized and fundamental niches of plant species (Hutchinson 1959). This perturbation, manifested by shifts in species' relative abundances along environmental gradients, might have resulted from the removal of key stressors (e.g. herbivory) and/or the amplification of others (e.g. competition for light and other resources) within exclosures. In one of the few studies addressing niche shifts due to ungulates, two native herb species expanded their spatial distributions to exploit a wider range of environmental conditions when ungulate pressure was removed (Gomez 2005). The opposite effect has also been found in exotic herbs, where abundances can drop rapidly following deer exclusion (Kalisz, Spigler & Horvitz 2014; Davalos, Nuzzo & Blossey 2015). However, the best-fit models for few species/groups included both deer and environmental predictors. Thus, responses of understorey species to deer exclusion may actually reflect the effects that deer have on environmental factors to a greater degree than the direct impacts of deer on understorey species via herbivory.
Variation in ecosystem properties among our study locations precluded certain potentially informative data analyses. For example, three of the four oldest exclosures lie on fertile soils derived from dolomitic parent material in Door County (Table SI) and have also had the highest historic deer population densities (15 deer km~2, or potentially much higher, reflecting historical hunting restrictions in state parks - Wisconsin Department of Natural Resources 2010). Because location-to-location variation in edaphic traits and deer population density was conflated with exclosure age, we do not present data on how age influenced community change, although we recognize that it may be an important source of variation (Collard et al. 2010; Royo et al. 2010; Hidding, Tremblay & Cote 2013; Habeck & Schultz 2015; Pendergast et al. 2016).
Results from this study and Frerker, Sabo & Waller (2014) support the contention that regional abundances of many native forb, shrub and palatable tree species may be threatened by overabundant deer. Deer impacts are not only due to the direct effects of herbivory on susceptible plants but also the alteration of a broad suite of environmental conditions. When deer increase light availability and soil resources by reducing sapling abundance, they also alter soil structure and morphology. Such indirect effects act in concert with browsing damage to alter understorey communities in complex and potentially long-lasting ways. However, we acknowledge that inferences drawn from our data are based primarily on correlation, and firmer conclusions regarding causes vs. consequences of the various deer impacts on forest understoreys await the results of future studies manipulating specific
© 2017 The Authors. Journal of Ecology © 2017 British Ecological Society, Journal of Ecology, 105, 1386-1398
1396  A. E. Sabo et al.
aspects of the environment and/or plant community. Additional research could examine how species' distributions along environmental gradients change across a continuum of deer population densities, paying particular attention to how species perform at their environmental extremes. Results from such studies will help us understand and predict how species and community assemblages are likely to respond to both the direct effects of herbivory and the complex effects that deer can have on environmental conditions.
Authors' contributions
A.S., D.W., K.F. and E.K. designed the study. A.S. ad K.F. collected and processed data. A.S. and E.K. conducted data analyses. A.S. generated the first draft of the manuscript, and all authors contributed substantially to revisions.
Acknowledgements
We thank Friends of Peninsula State Park, the UW-Madison Botany Department and the Mclntire-Stennis Program for funding our research. We appreciate that Kemp Natural Resources Station, Wisconsin Department of Natural Resources, U.S. Forest Service and Plum Creek Timber Company built, maintained and allowed us access to their exclosures. We also thank our field assistants, Andy Jandl, Marion Lee and Tyler Schappe, and two anonymous reviewers.
Data accessibility
Data used in this study deposited in the Dryad Digital Repository https://doi. org/10.5061/dryad.4nf75 (Sabo et al 2017).
References
Amatangelo, K.L., Johnson, S.E., Rogers, D.A. & Waller, D.M. (2014) Trait-environment relationships remain strong despite 50 years of trait compositional change in temperate forests. Ecology, 95, 1780-1791.
Anderson, R.C. (1994) Height of white-flowered trillium (Trillium grandiflo-rum) as an index of deer browsing intensity. Ecological Applications, 4, 104-109.
Anderson, M.J. (2001) A new method for nonparametric multivariate analysis of variance. Austral Ecology, 26, 32—46.
Anderson, M.J. (2004) PERMDISP: A FORTRAN Computer Program for Per-mutational Analysis of Multivariate Dispersions (For Any Two-Factor ANOVA Design) Using Permutation Tests. Department of Statistics, University of Auckland, Auckland, New Zealand.
Arkle, R.S., Pilliod, D.S. & Welty, J.L. (2012) Pattern and process of prescribed fires influence effectiveness at reducing wildfire severity in dry coniferous forests. Forest Ecology and Management, 276, 174-184.
Augustine, D.J. & Frelich, L.E. (1998) Effects of white-tailed deer on populations of an understory forb in fragmented deciduous forests. Conservation Biology, 12, 995-1004.
Augustine, D.J., Frelich, L.E. & Jordan, P.A. (1998) Evidence for two alternate stable states in an ungulate grazing system. Ecological Applications, 8, 1260-1269.
Baiser, B., Lockwood, J.L., La Puma, D. & Aronson, M.F. (2008) A perfect storm: two ecosystem engineers interact to degrade deciduous forests of New Jersey. Biological Invasions, 10, 785-795.
Balgooyen, CM. & Waller, D.M. (1995) The use of Clintonia borealis and other indicators to gauge impacts of white-tailed deer on plant communities in northern WI, USA. Natural Areas Journal, 15, 308-318.
Barrios-Garcia, M.N. & Ballari, S.A. (2012) Impact of wild boar (Sus scrofa) in its introduced and native range: a review. Biological Invasions, 14, 2283-2300.
Bodmer, R.E. (1990) Ungulate frugivores and the browser-grazer continuum. Oikos, 57, 319-325.
Bradshaw, L.B. & Waller, D.M. (2016) Impacts of white-tailed deer on regional patterns of forest tree recruitment. Forest Ecology and Management, 375, 1-11.
Bressette, J.W., Beck, H. & Beauchamp, V.B. (2012) Beyond the browse line: complex cascade effects mediated by white-tailed deer. Oikos, 121, 1749-1760.
Burton, J.L, Mladenoff, D.J., Clayton, M.K. & Forrester, J.A. (2011) The roles of environmental filtering and colonization in the fine-scale spatial patterning of ground-layer plant communities in north temperate deciduous forests. Journal of Ecology, 99, 764-776.
Chantigny, M.H. (2003) Dissolved and water-extractable organic matter in soils: a review on the influence of land use and management practices. Geoderma, 113, 357-380.
Chips, M.J., Magliocca, M.R., Hasson, B. & Carson, W.P. (2014) Quantifying deer and turkey leaf litter disturbances in the eastern deciduous forest: have nontrophic effects of consumers been overlooked? Canadian Journal of Forest Research, 44, 1128-1132.
Collard, A., Lapointe, L., Ouellet, J.-P., Crete, M., Lussier, A., Daigle, C. & Cote, S.D. (2010) Slow responses of understory plants of maple-dominated forests to white-tailed deer experimental exclusion. Forest Ecology and Management, 260, 649-662.
Cote, S.D., Rooney, T.P., Tremblay, J.-P., Dussault, C. & Waller, D.M. (2004) Ecological impacts of deer overabundance. Annual Review of Ecology, Evolution, and Systematics, 35, 113-147.
Dävalos, A, Nuzzo, V. & Blossey, B. (2014) Demographic responses of rare forest plants to multiple stressors: the role of deer, invasive species and nutrients. Journal of Ecology, 102, 1222-1233.
Dävalos, A., Nuzzo, V. & Blossey, B. (2015) Single and interactive effects of deer and earthworms on non-native plants. Forest Ecology and Management, 351, 28-35.
Dobson, A. & Blossey, B. (2015) Earthworm invasion, white-tailed deer and seedling establishment in deciduous forests of north-eastern North America. Journal of Ecology, 103, 153-164.
Dufresne, M., Bradley, R.L., Tremblay, J.-P., Poulin, M. & Pellerin, S. (2009) Clearcutting and deer browsing intensity interact in controlling nitrification rates in forest floor. Ecoscience, 16, 361-368.
Fletcher, J.D., McShea, W.J., Shipley, L.A. & Shumway, D. (2001) Use of common forest forbs to measure browsing pressure by white-tailed deer (Odocoileus virginianus Zimmerman) in Virginia, USA. Natural Areas Journal, 21, 172-176.
Frelich, L.E. & Lorimer, C.G. (1985) Current and predicted long-term effects of deer browsing in hemlock forests in Michigan, USA. Biological Conservation, 34, 99-120.
Frerker, K., Sabo, A. & Waller, D. (2014) Long-term regional shifts in plant community composition are largely explained by local deer impact experiments. PLoS ONE, 9, el 15843.
Godefroid, S. & Koedam, N. (2004) Interspecific variation in soil compaction sensitivity among forest floor species. Biological Conservation, 119, 207-217.
Gömez, J.M. (2005) Long-term effects of ungulates on performance, abundance, and spatial distribution of two montane herbs. Ecological Monographs, 75, 231-258.
Habeck, C.W. & Schultz, A.K. (2015) Community-level impacts of white-tailed deer on understorey plants in north American forests: a meta-analysis. AoB Plants, 7, 1-12.
Hanley, T.A. (1997) A nutritional view of understanding and complexity in the problem of diet selection by deer (Cervidae). Oikos, 79, 209-218.
Heckel, CD., Bourg, N.A., McShea, W.J. & Kalisz, S. (2010) Nonconsumptive effects of a generalist ungulate herbivore drive decline of unpalatable forest herbs. Ecology, 91, 319-326.
Hidding, B., Tremblay, J. & Cote, S.D. (2013) A large herbivore triggers alternative successional trajectories in the boreal forest. Ecology, 94, 2852-2860.
Hobbs, N.T. (1996) Modification of ecosystems by ungulates. The Journal of Wildlife Management, 60, 695-713.
Holmes, S.A. & Webster, CR. (2011) Herbivore-induced expansion of generalist species as a driver of homogenization in post-disturbance plant communities. Plant Ecology, 212, 753-768.
Horsley, S.B., Stout, S.L. & DeCalesta, D.S. (2003) White-tailed deer impact on the vegetation dynamics of a northern hardwood forest. Ecological Applications, 13, 98-118.
Hutchinson, G.E. (1959) Why are there so many kinds of animals. American Naturalist, 89, 145-159.
Hutchinson, T.F., Boerner, R.E.J., Iverson, L.R., Sutherland, S. & Sutherland, E.K. (1999) Landscape patterns of understory composition and richness across a moisture and nitrogen mineralization gradient in Ohio (U.S.A.) Quercus forests. Plant Ecology, 144, 177-189.
Jenkins, L.H., Jenkins, M.A., Webster, CR., Zollner, P.A. & Shields, J.M. (2014) Herbaceous layer response to 17 years of controlled deer hunting in forested natural areas. Biological Conservation, 175, 119-128.
© 2017 The Authors. Journal of Ecology © 2017 British Ecological Society, Journal of Ecology, 105, 1386-1398
Deer effects on forest understoreys 1397
Jensen, N.R., Webster, C.R., Witt, J.C. & Grant, J.B. (2011) Ungulate winter habitat selection as a driver of herbaceous-layer heterogeneity in northern temperate forests. Ecosphere, 2, 1-16.
Kain, M., Battaglia, L., Royo, A. & Carson, W.P. (2011) Over-browsing in Pennsylvania creates a depauperate forest dominated by an understory tree: results from a 60-year-old deer exclosure. Journal of the Torrey Botanical Society, 138, 322-326.
Kalisz, S., Spigler, R.B. & Horvitz, C.C. (2014) In a long-term experimental demography study, excluding ungulates reversed invader's explosive population growth rate and restored natives. Proceedings of the National Academy of Sciences of the USA, 111, 4501^1506.
Kirschbaum, CD. & Anacker, B.L. (2005) The utility of Trillium and Maianthemum as phyto-indicators of deer impact in northwestern Pennsylvania, USA. Forest Ecology and Management, 217, 54-66.
Knapp, A.K., Blair, J.M., Briggs, J.M., Collins, S.L., Hartnett, D.C., Johnson, L.C. & Towne, G. (1999) The keystone role of bison in North American tall-grass prairie: bison increase habitat heterogeneity and alter a broad array of plant, community, and ecosystem processes. BioScience, 49, 39-50.
Knight, T.M., Dunn, J.L., Smith, L.A., Davis, J. & Kalisz, S. (2009) Deer facilitate invasive plant success in a Pennsylvania forest understory. Natural Areas Journal, 29, 110-116.
Krueger, L.M., Peterson, C.J., Royo, A. & Carson, W.P. (2009) Evaluating relationships among tree growth rate, shade tolerance, and browse tolerance following disturbance in an eastern deciduous forest. Canadian Journal of Forest Research, 39, 2460-2469.
Kuijper, D.P.J., Cromsigt, J.P.G.M., Jgdrzejewska, B., Miscicki, S., Churski, M. & Jgdrzejewski, W. (2010) Bottom-up versus top-down control of tree regeneration in the Bialowieza Primeval Forest, Poland. Journal of Ecology, 98, 888-899.
Long, Z.T., Pendergast, T.H. IV & Carson, W.P. (2007) The impact of deer on relationships between tree growth and mortality in an old-growth beech-maple forest. Forest Ecology and Management, 252, 230-238.
Loss, S.R., Hueffmeier, R.M., Hale, CM., Host, G.E., Sjerven, G. & Frelich, L.E. (2013) Earthworm invasions in northern hardwood forests: a rapid assessment method. Natural Areas Journal, 33, 21-30.
Lundstrom, U.S., van Breemen, N. & Bain, D. (2000) The podzolization process. A review. Geoderma, 94, 91-107.
McCune, B. (2006) Non-parametric habitat models with automatic interactions. Journal of Vegetation Science, 17, 819-830.
McCune, B., Grace, J.B. & Urban, D.L. (2002) Analysis of Ecological Communities. MjM Software Design, Gleneden Beach, OR, USA.
McCune, B. & Mefford, M.J. (2011) HyperNiche v. 2.10. MjM Software, Gleneden Beach, OR, USA.
McEwan, R.W. & Muller, R.N. (2011) Dynamics, diversity, and resource gradient relationships in the herbaceous layer of an old-growth Appalachian forest. Plant Ecology, 212, 1179-1191.
McGraw, J.B. & Furedi, M.A. (2005) Deer browsing and population viability of a forest understory plant. Science, 307, 920-922.
Michigan Department of Natural Resources and Environment (2010) Michigan Deer Management Plan. Wildlife Division Report 3512, Michigan Department of Natural Resources and Environment.
Miller, S.G., Bratton, S.P. & Hadidian, J. (1992) Impacts of white-tailed deer on endangered and threatened vascular plants. Natural Areas Journal, 12, 67-74.
Mudrak, E.L., Johnson, S.E. & Waller, D.M. (2009) Forty-seven year changes in vegetation at the Apostle Islands: effects of deer on the forest understory. Natural Areas Journal, 29, 167-176.
Murray, B.D., Webster, CR. & Bump, J.K. (2013) Broadening the ecological context of ungulate-ecosystem interactions: the importance of space, seasonality, and nitrogen. Ecology, 94, 1317-1326.
N0rnberg, P., Sloth, L. & Nielsen, K.E. (1993) Rapid changes of sandy soils caused by vegetation changes. Canadian Journal of Soil Science, 73, 459-468.
Nuttle, T., Ristau, T.E. & Royo, A.A. (2014) Long-term biological legacies of herbivore density in a lands cape-scale experiment: forest understoreys reflect past deer density treatments for at least 20 years. Journal of Ecology, 102, 221-228.
Nuttle, T., Yerger, E.H, Stoleson, S.H. & Ristau, T.E. (2011) Legacy of top-down herbivore pressure ricochets back up multiple tropic levels in forest canopies over 30 years. Ecosphere, 2, 1-11.
Oksanen, J., Blanchet, F.G., Kindt, R., Legendre, P., O'Hara, R.B., Simpson, G.L., Solymos, P., Stevens, M.HH. & Wagner, H. (2013) Vegan: Community Ecology Package. R Package Version 2.0-7. Available at: https://CRAN. Rproject.org/package-vegan (accessed 1 September 2013).
Pastor, J., Dewey, B., Moen, R., Mladenoff, D.J., White, M. & Cohen, Y. (1998) Spatial patterns in the moose-forest-soil ecosystem on Isle Royale, Michigan, USA. Ecological Applications, 8, 411—424.
Pendergast, T.H., Hanlon, S.M., Long, Z.M., Royo, A. & Carson, W.P. (2016) The legacy of deer overabundance: long-term delays in herbaceous understory recovery. Canadian Journal of Forest Research, 46, 362-369.
Persson, I., Danell, K. & Bergstrom, R. (2000) Disturbance by large herbivores in boreal forests with special reference to moose. Annales Zoologici Fennici, 21, 251-263.
Promis, A., Gartner, S., Butler-Manning, D., Duran-Rangel, C, Reif, A., Cruz, G. & Hernandez, L. (2011) Comparison of four different programs for the analysis of hemispherical photographs using parameters of canopy structure and solar radiation trans mittance. Waldokologie, Landschaftsforschung und Naturshutz, 11, 19-33.
Putman, R. (1988) The Natural History of Deer. Cornell University Press, Ithaca, NY, USA.
Rice, W.R. (1989) Analyzing tables of statistical tests. Evolution, 43, 223-225.
Rooney, T.P. (2001) Deer impacts on forest ecosystems: a North American perspective. Forestry, 74, 201-208.
Rooney, T.P. & Dress, W.J. (1997) Patterns of plant diversity in overbrowsed primary and mature secondary hemlock-northern hardwood forest stands. Journal of the Torrey Botanical Society, 124, 43-51.
Rooney, T.P. & Waller, D.M. (2003) Direct and indirect effects of white-tailed deer in forest ecosystems. Forest Ecology and Management, 181, 165-176.
Royo, A.A., Stout, S.L., deCalesta, D.S. & Pierson, T.G. (2010) Restoring forest herb communities through lands cape-lev el deer herd reductions: is recovery limited by legacy effects? Biological Conservation, 143, 2425-2434.
Russell, F.L., Zippin, D.B. & Fowler, N.L. (2001) Effects of white-tailed deer (Odocoileus virginianus) on plants, plant populations and communities: a review. The American Midland Naturalist, 146, 1-26.
Sabo, A.E., Frerker, K.L., Waller, D.M. & Kruger, E.L. (2017) Data from: Deer-mediated changes in environment compound the direct impacts of herbivory on understorey plant communities. Dryad Digital Repository, https://doi.org/10.5061/dryad.4nn5
Sanborn, P., Lamontagne, L. & Hendershot, W. (2011) Podzolic soils of Canada: genesis, distribution, and classification. Canadian Journal of Soil Science, 91, 843-880.
Seagle, S.W. (2003) Can ungulates foraging in a multiple-use landscape alter forest nitrogen budgets? Oikos, 103, 230-234.
Shelton, A.L. & Inouye, R.S. (1995) Effect of browsing by deer on the growth and reproductive success of Lactuca canadensis (Asteraceae). American Midland Naturalist, 134, 332-339.
Shelton, A.L., Henning, J.A., Schultz, P. & Clay, K. (2014) Effects of abundant white-tailed deer on vegetation, animals, mycorrhizal fungi, and soils. Forest Ecology and Management, 320, 39—49.
Tahtinen, B., Murray, B.D., Webster, C.R., Tarasoff, CS. & Burton, A.J. (2014) Does ungulate foraging behavior in forest canopy gaps produce a spatial subsidy with cascading effects on vegetation? Forest Science, 60, 819-829.
Tanentzap, A.J., Kirby, K.J. & Goldberg, E. (2012) Slow responses of ecosystems to reductions in deer (Cervidae) populations and strategies for achieving recovery. Forest Ecology and Management, 264, 159-166.
Tanentzap, A.J., Bazely, D.R., Koh, S., Timciska, M., Haggith, E.G., Carelton, T.J. & Coomes, D.A. (2011) Seeing the forest for the deer: do reductions in deer-disturbance lead to forest recovery? Biological Conservation, 144, 376-382.
Tilghman, N.G. (1989) Impacts of white-tailed deer on forest regeneration in northwestern Pennsylvania. The Journal of Wildlife Management, 53, 524-532.
Tremblay, J.-P., Huot, J. & Potvin, F. (2007) Density-related effects of deer browsing on the regeneration dynamics of boreal forests. Journal of Applied Ecology, 44, 552-562.
Tripler, C.E., Canham, CD., Inouye, R.S. & Schnurr, J.L. (2005) Competitive hierarchies of temperate tree species: interactions between resource availability and white-tailed deer. Ecoscience, 12, 494-505.
USDA, NRCS (2013) Web Soil Survey. Soil Survey Staff. Available at: http:// websoilsurvey.sc.egov.usda.gov (accessed 1 February 2013)
USDA, NRCS (2016) The PLANTS Database. National Plant Data Team, Greensboro, NC, USA. Available at: http://plants.usda.gov (accessed 1 September 2016).
Waller, D.M. & Alverson, W.S. (1997) The white-tailed deer: a keystone herbivore. Wildlife Society Bulletin, 25, 217-226.
Webster, C.R., Jenkins, M.A. & Rock, J.H. (2005) Long-term response of spring flora to chronic herbivory and deer exclusion in Great Smoky Mountains National Park, USA. Biological Conservation, 125, 297-307.
© 2017 The Authors. Journal of Ecology © 2017 British Ecological Society, Journal of Ecology, 105, 1386-1398
1398  A. E. Sabo et al.
White, M.A. (2012) Long-term effects of deer browsing: composition, structure and productivity in a northeastern Minnesota old-growth forest. Forest Ecology and Management, 269, 222-228.
Wilde, S.A. & Leaf, A.L. (1955) The relationship between the degree of soil podzolization and the composition of ground cover vegetation. Ecology, 36, 19-22.
Wisconsin Department of Natural Resources (2010) Deer Abundance and Densities in Wisconsin. State of Wisconsin. Available at: http://dnr.wi.gov/topic/ hunt/maps.html (accessed 1 March 2012).
Wright, D.M., Tanentzap, A.J., Flores, O., Husheer, S.W., Duncan, R.P., Wiser, S.K. & Coomes, D.A. (2012) Impacts of culling and exclusion of browsers on vegetation recovery across New Zealand forests. Biological Conservation, 153, 64-71.
Yamada, H. & Takatsuki, S. (2015) Effects of deer grazing on vegetation and ground-dwelling insects in a larch forest in Okutama, western Tokyo. International Journal of Forestry Research, 2015, 687506. doi:10.1155/2015/ 687506
Received 1 May 2016; accepted 16 January 2017 Handling Editor: Matthew Heard
Supporting Information
Details of electronic Supporting Information are provided below.
Table SI. Means (with standard deviations) and ranges, by location and deer treatment, for the 20 biotic and abiotic environmental variables used in our environmental NMDS ordinations.
Table S2. Scores of common species along vegetation NMDS axes for saplings, seedlings and shrubs, and herbs.
Table S3. Correlation coefficients of the relationships between environmental NMDS axes and their individual environmental components.
Table S4. The abundances of several species and species groups were significantly related to thickness of the soil E horizon.
Table S5. The abundances of several species and species groups were significantly related to soil compaction.
Table S6. The frequencies of several groundlayer species and species groups were significantly related to openness.
Table S7. The frequencies of several groundlayer species and species groups were significantly related to total basal area.
Table S8. Relationships between exclosure effects on the frequency of a particular groundlayer species or group and those on individual environmental parameters.
Table S9. Results of nonparametric multiplicative regression (NPMR) for species or groups in which patterns of abundance across our ordination of environmental space differed between exclosure and deer-access treatments.
Table S10. Results of nonparametric multiplicative regression (NPMR) for common species in which abundance varied across our ordination of environmental space but did not differ between exclosure and deer-access treatments.
© 2017 The Authors. Journal of Ecology © 2017 British Ecological Society, Journal of Ecology, 105, 1386-1398