Decreasing deer browsing pressure influenced understory vegetation dynamics over 30 years

Abstract

Key message

Thanks to the concomitant recordings of vegetation and deer browsing sampled first in 1976, then resurveyed in 2006, we show that forest plant communities shifted in response to deer population dynamics, stand management and eutrophication.

Context and aims

High deer populations alter forest understory dynamics worldwide. However, no study ever attempted to rank the importance of deer herbivory relatively to other environmental drivers. In the Arc-en-Barrois National Forest (France), we investigated whether (i) deer browsing is a critical driver of vegetation composition and dynamics, (ii) the vegetation communities recover after a decrease in deer populations.

Methods

In 2006, we resurveyed 321 plots from a network of 1027 plots where vegetation composition and browsing pressure was first assessed in 1976. We used coinertia analysis to identify the gradients in vegetation composition in 1976, when abiotic variables were also recorded. We assessed shifts in plant community composition using mean Ellenberg indicator values, analysed plot scores shifts along the axes of the coinertia analysis and correlated these shifts with changes in browsing pressure.

Results

Two major gradients determined vegetation composition in 1976: edaphic variables (nitrogen availability and soil moisture) and browsing pressure. Over the next 30 years, we noticed a strong increase in nitrophilous plant species frequency and community composition shifted towards lightly browsed characteristics, accompanying a decrease in browsing pressure. Shifts in community composition were significantly correlated with the intensity of changes in browsing pressure, showing that deer population dynamics were a determinant driver of changes in plant assemblages.

Conclusion

Our results provide evidence for a structuring effect of deer browsing on vegetation composition, once forest site variations (soil moisture and nitrogen) were accounted for. We observed an incomplete recovery of the communities 25 years after the reduction of deer densities, suggesting a delayed response to deer population reduction. Long-term monitoring of forest biodiversity should therefore include browsing pressure assessment to control for potential effects of wild ungulates.

1 Introduction

Because plant communities are relevant bio-indicators for environmental conditions and provide evidence for the consequences of global changes in the biosphere (Gilliam 2007), the mid- and long-term dynamics of herbaceous vegetation have received considerable attention for several decades. Since the early 1990s, the effect of atmospheric nitrogen deposition on changes in vegetation composition has been widely described, particularly for herbaceous forest vegetation (see review in the introduction of Verheyen et al. (2012)). In the last decade, the focus has shifted towards the effects of global warming (e.g. Bertrand et al. (2011)). But, mid- or long-term variations in herbivory pressure are another potential driver of vegetation communities that has rarely been addressed in Europe. The eradication of large predators since the early twentieth century combined with strict hunting regulation laws has led to a global increase in ungulate populations in the last few decades in Western Europe (Apollonio et al. 2010). Consequently, the pressure on forest vegetation due to high deer populations has grown (Estes et al. 2011).

The effects of high deer populations on forest ecosystems have long been documented worldwide (see (Rooney 2001) for a review). Thus, deer feeding directly affect the structure of the plants they consume, both horizontally—by reducing their cover percentage and abundance—and vertically—by limiting their height growth (Horsley et al. 2003). Moreover, through selective feeding, deer may affect species composition by avoiding unpalatable species or favouring species particularly resilient to browsing (Augustine and McNaughton 1998). Deer also enhance the dispersion of both epi- and endo-zoochorous species (Myers et al. 2004) and, conversely, may alter the reproductive process of other plants by consuming flowers (Mårell et al. 2009) or fruits, like acorns, that are preyed upon by deer (Picard and Gégout 1992). These effects of deer activity directly impact plant communities and more generally alter environmental conditions that could induce secondary changes in vegetation composition. For example, the reduction in shrub cover caused by browsing leads to an increase in light availability for ground-layer vegetation (Gill and Beardall 2001) and thus enhances tree seedling germination (Traveset et al. 2007); this may drive vegetation towards more heliophilous communities. Inversely, de la Cretaz and Kelty (1999) observed that herbivore removal of competing plants led to the accelerated growth of the hay-scented fern (Dennstaedtia punctilobula), which in turn inhibited tree seedling establishment and development through allelopathic mechanisms. Monitoring vegetation over decades is particularly relevant to detecting such indirect, long-term effects of deer in field studies (Barrett et al. 2006).

Even though deer browsing often appears to be a probable cause of changes in plant community composition, most long-term forest vegetation monitoring programmes have no precise measurements to support this hypothesis (Corney et al. 2008). Fenced exclosures, either with controlled deer densities or without any deer, are the most rigorous method to infer deer impact. Though often of small size and few in number, such experimental exclosures appear to be complementary to large-scale vegetation surveys focusing on spatial gradients in ungulate impacts (Waller et al. 2009). Watkinson et al. (2001) clearly point out the need for a local assessment of browsing pressure to help to explain a part of the variation in ground flora composition; however, such quantification is rarely done in vegetation surveys (but see Veblen et al. (1989)). Consequently, the influence of deer browsing on vegetation composition and dynamics in field studies remains hypothetical most of the time (Decocq et al. 2005; Verheyen et al. 2012).

In 1976, due to high deer populations, the northern part of the Arc-en-Barrois forest (the Arc-Châteauvillain forest) in north-eastern France was deemed suitable for studying the impacts of deer browsing on vegetation communities (Allain et al. 1978). Indeed, roe deer, a browser, is likely to impact shrub cover and composition, and red deer, an intermediate feeder that consume not only ligneous but also herbaceous plants, is likely to alter ground flora. One year later, the southern part (the Arc-Carrefour de Joinville Forest), which hosted lower deer densities and was moderately impacted by deer browsing, was inventoried to serve as a reference. On 1027 plots, both vegetation composition and browsing pressure were surveyed at that time. Thirty years later, we resampled 321 of those plots in an attempt to reveal the links between spatio-temporal variations in deer browsing pressure and vegetation changes. The first aim of this study was to identify the spatial gradients that structure vegetation communities. We expected that abiotic conditions such as water availability, temperature and mineral nutrition would be the main determinants of vegetation composition. We also focused on whether or not and how herbivory pressure shapes plant communities. Thanks to the original data available on browsing pressure and vegetation communities, we were able to carry out comparisons over time. A decrease in deer population that occurred in the Arc-Châteauvillain Forest at the end of the 1970s over a 5-year period (Boulanger et al. 2009) gave us the opportunity to study the vegetation response to a decrease in herbivory pressure under natural conditions. Making the most of this context, we analysed the intensity of the shifts in vegetation community composition and more precisely the correlations with the changes in herbivory pressure. We predicted (i) that the vegetation composition would recover following the decrease in deer population in the Arc-Châteauvillain Forest and (ii) that the magnitude of the shifts in plot vegetation composition would be correlated with the change in local browsing pressure. Finally, we discuss the ranking of the various drivers of the change in vegetation composition (climate warming, eutrophication, forest stand management, and deer browsing) that occurred in this forest over 30 years.

2 Material and methods

2.1 Study area

The Arc-en-Barrois state forest (48°0′N 5°00′E) is located in the Haute-Marne administrative region in north-eastern France. It is divided into two forests: (i) the 8797-ha “Arc-Châteauvillain” forest to the North (hereafter “North Forest”) and (ii) the 2191-ha “Arc-Carrefour de Joinville” to the South (hereafter “South Forest”). The Arc-en-Barrois state forest is situated in the ecological region of the Champagne-Ardenne limestone plateaux on a calcareous tertiary sub-stratum; altitudes span 238–426 m a.s.l. Mean annual temperature and rainfall were 9.5 °C and 860 mm/year during the 1957–1976 period, and 10.5 °C and 882 mm/year during the 1987–2006 period at the nearby Langres meteorological station (Meteo-France), indicating a +1 °C warming trend.

Since the 1970s, the previous coppice-with-standards regime has been gradually converted to even-aged high forest. Forest stands are dominated by two oaks (in 2006, Quercus robur L. and Quercus petraea L. respectively occurred in 79 and 12 % of the plots in the North Forest, 59 and 60 % in the South Forest) and beech (Fagus sylvatica L.), mixed with hornbeam (Carpinus betulus L.) and hazel (Corylus avellana L.) in the coppice layer.

Free-ranging populations of wild boar (Sus scrofa L.), red deer (Cervus elaphus L.) and roe deer (Capreolus capreolus L.) are present in both forests. Since the early 1980s, forest management has shifted from game-oriented objectives, targeting mainly red deer and wild boar populations, to wood production. This refocusing required a substantial decrease in game populations to allow forest regeneration. No reliable data was available to estimate deer population densities, and their temporal trends as frequent changes in game management policies made hunting counts unrepresentative of population levels (Morellet et al. 2007). Consequently, we focused on changes in herbivory pressure, measured both in 1976 and 2006 during browsing surveys, which provide reliable information about deer impacts on vegetation (Frerker et al. 2013). The surveys revealed that average browsing pressure decreased substantially from 40 % in 1976 to 10 % in 2006 in the North Forest (Figure S1). In the South Forest, it remained consistently low and below 10 %.

2.2 Sampling design

In 1976, a network of 1027 plots was established at the nodes of a regular grid (10-ha cells, plots spaced 333 m apart). The North Forest was sampled in 1976 from July 15th to October 15th (820 plots) and the South Forest in 1977 from September 15th to October 15th (207 plots). Hereafter, we will refer to this baseline sampling as the 1976 campaign. Plots were offset to the nearest tree; a cross was painted on the tree to identify the plot for later relocation. In the summers (July 15th to October 15th) of 2005 and 2006 (hereafter “2006”), we resampled 321 of these plots (248 in the North Forest and 73 in the South Forest) in mature stands that had not been regenerated since 1976 nor harvested during the year previous to the survey. These sites were selected to match a partial resampling conducted in 1981 (not studied here) and to increase spatial coverage by filling in the largest gaps in this 1981 sampling. One hundred and twenty-two plots out of the 321 resampled were relocated exactly (painted cross on the tree), while the others were relocated thanks to precise points recorded on the initial 1/10,000 maps and notes taken during the original surveys (Fig. 1).

Fig. 1

Map of the Arc-en-Barrois Forest (embedded its location in France), showing the two parts (North and South) and the location of the plots (small points sampled in 1976 only, large points resampled in 2006)

At each plot, all plant species were recorded and separated into four vegetation layers: herbaceous (up to 50 cm in height), low shrubs (from 50 cm to 2 m in height), high shrubs (from 2 m in height to suppressed trees) and trees (co-dominant or dominant). A cover value based on an ordinal scale with six classes (Braun-Blanquet 1932) was assigned to each species in each layer. Because the limit between the herbaceous and low shrub layers was not clearly defined in the 1976 protocol, we pooled these two layers and retained the maximum cover mark (hereafter called the “<2 m layer”) before all subsequent analyses. In 1976, the field observers had also collected general information for each plot: slope, exposure, topographic position and soil variables (humus type, texture and stoniness in successive 20-cm layers) to establish a typology of the forest plots. These latter variables are hereafter called “abiotic site variables”.

In 1976, vegetation surveys at each plot were conducted on areas varying from 400 to 1000 m², according to the minimal area definition (Braun-Blanquet 1932), without keeping any record of this area. In 2006, we standardized both sampling time (to 30 min with two operators) and sampling area (to 400 m²) in order to correct for lack of exhaustiveness in plant inventories (Archaux et al. 2006). In addition, we visually estimated the global cover percentage for each vegetation layer.

In both 1976 and 2006, a browsing score was assigned to each woody or semi-woody species (n = 56) recorded in the vegetation survey and exhibiting twigs available to deer (below 2 m in height). These marks are based on the proportion of shoots browsed and the global shape of the species (Table S1). Before carrying out the 2006 browsing survey, we conducted a calibration phase based on advice from the leading expert during the 1976 survey (J.F. Picard).

2.3 Data analysis

2.3.1 Index of local browsing pressure

For each species i and each plot j, cover and browsing marks (C _ij and B _ij , respectively) were first converted to percentages (median value of each class, supplementary material). We then calculated an average browsing pressure index per plot and sampling date as the ratio between the amount of browsed vegetation (sum of C _ij  × B _ij ) and the total amount of vegetation available (sum of C _ij ). We used this variable after square root transformation as an index of browsing pressure at the plot scale (Boulanger et al. 2009; Frerker et al. 2013).

2.3.2 Indices of stand composition, stand structure and environmental conditions

We performed a correspondence analysis (CA) to identify a smaller set of general variables describing stand structure and tree composition on the 1027 plots recorded in 1976. We grouped the species occurring in the tree and high-shrub layers into six categories, retaining the highest cover-abundance coefficient: beech, hornbeam, oaks, pioneer tree species (Betula sp., Populus tremula L. and Salix caprea L. pooled together), fruit trees (Malus sylvestris Mill., Prunus avium L., Pyrus pyraster Burgsd., Sorbus aria (L.) Crantz and Sorbus torminalis (L.) Crantz pooled together) and other species (pooled together). In this way, we avoided deleting any rare species and were able to correct for possible misidentifications (especially inconsistencies between 1976 and 2006). We interpreted the first two axes of the CA (accounting for 17.3 and 14.4 % of the total variance) as (1) the ratio between hornbeam and beech (“Stand.1”) and (2) stand maturity (opposing pioneer tree species to oaks and fruit trees, “Stand.2”). We projected the 321 plots surveyed in 2006 as supplementary individuals in the same PCA to measure the changes that occurred in stand characteristics. The coordinates of the plots positioned along the first two axes were subsequently used as stand variables.

We used a multiple correspondence analysis for abiotic site variables because some of them were discrete and non-ordered. We interpreted the first axis (11.3 % of total variance) as a gradient of water availability (“Soil.Water”), opposing superficial stony soils with low water reserves and deep soils with high water reserves. The second axis only explained 5.8 % and did not represent any obvious ecological gradient.

2.4 Community changes

2.4.1 Plot mean Ellenberg indicator values

For the two sampling dates at each plot, we calculated the unweighted mean Ellenberg indicator values (L light, T temperature, F water availability, R pH, and N Nitrogen) for the community of species occurring in the <2 m layer (Diekmann 2003). As these values were non-normally distributed, we tested for changes in plot indicator values between the two inventories using Wilcoxon rank-signed tests to account for the pairwise structure of the data that resulted from resampling. We tested the relation between these changes in indicator values and the change in browsing pressure with Spearman’s rank correlation tests.

2.4.2 Plant community composition

The floristic table of the 1976 sample contained 1027 relevés and 142 species; we removed rare species (occurring in less than 10 plots). The environmental table of the 1976 sample contained 1027 relevés and 11 variables (deer browsing index, the two general stand axes from the CA “Stand.1” and “Stand.2”, the global soil water availability Soil.Water, tree cover, high shrub cover and mean Ellenberg indicator values for light, temperature, soil moisture, pH and nitrogen). First, we applied a correspondence analysis to the floristic table. Cover percentages were square-root-transformed to limit the influence of high-cover species; rare species were made supplementary, i.e. a posteriori positioned along the axes according to the plots in which they occurred. Second, a principal component analysis (PCA) was applied to the environmental table. Third, in order to analyse the links between the vegetation community structure and environmental gradients in 1976, we used a co-inertia analysis corresponding to a joint analysis of the CA on the floristic table and the PCA on the environmental table. Co-inertia is particularly suited for the identification of structures in vegetation communities that are associated with environmental gradients, especially when these environmental variables are numerous and possibly correlated (Dray et al. 2003). As the total forest area was surveyed in 1976, we considered the vegetation communities and their gradients as a baseline for the comparison with the 2006 sampling.

Next, we projected the 321 floristic relevés from 2006 as supplementary sites in the coinertia analysis in order to obtain their scores in the multivariate space defined by the 1976 data. In this way, temporal shifts in vegetation composition could only be related to ecological variables if these shifts followed the same patterns as the spatial gradients identified through the co-inertia analysis applied to 1976 data. We analysed the changes in plot scores on the factorial axes using paired t tests. We computed the length of the vector linking the two positions of a relevé in 1976 and 2006. This provides a Euclidean measure of the dissimilarity quantifying the 30-year shift in community composition; this shift could be broken down along the factorial axes for further interpretation. We identified and ranked the drivers of the changes in vegetation composition using multiple linear regressions linking these shifts (vector length and decompositions along the axis of the co-inertia analysis) to the standardized changes in browsing pressure, Ellenberg plot indicator values, stand composition and tree and high shrub layer cover values.

All the analyses were done using R 3.0.2 and ade4 package (Dray and Dufour 2007). In the following statistical tests, the threshold level of rejection was set at 5 %, but if a meaningful trend appeared, significance at 10 % was also indicated.

3 Results

3.1 Gradients of variation in vegetation composition in 1976

Co-inertia analysis of the 1976 vegetation and environmental data tables returned two major axes of variation accounting for 84.2 % of the total variability.

The first axis (52.9 %) represented a gradient of soil water availability and acidity, opposing dry calcareous plots to plots with a better water supply (Fig. 2) and deeper soils. A gradient of stand maturity is also correlated to this edaphic axis. The species that obtained the highest scores on the first axis are (meso-) xerophilous and light-demanding species, as opposed to acid-tolerant or nitrogen- and water-demanding species that obtained the lowest scores (Fig. 3).

Fig. 2

Co-inertia analysis of the 1027 plots × 142 species sampled in 1976, showing plot positions and the correlations with environmental variables along the first and second axes. The environmental variables are Ellenberg mean indicator values for light (EIV.L), temperature (EIV.T), soil moisture (EIV.F), soil pH (EIV.R) and nitrogen (EIV.N), water availability constraint (Soil.Water), cover of the tree (Tree.cover) and high shrub (H.shrub.cover) layers, balance between Hornbeam and Beech (Stand.1), degree of maturity of the Oak stands (Stand.2) and browsing pressure (Deer.Browsing). See “Materials and Methods” for details on the computation of the environmental variables

Fig. 3

Coordinates of the species position projected along the first and second axes of the co-inertia analysis. Only the species that contributed the most to the definition of the axis are shown. Full bars positive values, open bars negative values

The second axis (accounting for 31.3 % of the total inertia) could be interpreted as a gradient of browsing pressure and high shrub cover, the two variables being negatively correlated (Spearman’s ρ = −0.31). The species associated with high browsing pressure at the bottom of this axis were mainly broadleaf forbs and grasses; on the opposite end, woody species and sedges were associated with low browsing pressure (Fig. 4, p < 0.01, F value = 5.08). The correlation between browsing pressure and the position of the plots along the second axis was high (ρ = −0.61, p < 0.001) for the whole forest (North and South pooled together). Within the North Forest, the correlation remained high (ρ = −0.49, p < 0.001) whereas in the South Forest, the correlation was not significant (ρ = 0.06, p > 0.1). Analysis of covariance of the effects of browsing pressure and forest (North or South) on the position of plots along the second axis revealed the significant effects of both browsing pressure (p < 0.001), forest (p < 0.001) and the interaction between forest and browsing pressure (p < 0.05). Tree cover and indicator values for light and for nitrogen were also significantly correlated to browsing pressure but to a lesser extent (ρ = 0.24, ρ = 0.17 and ρ = 0.26, respectively, all p < 0.001).

Fig. 4

Distribution of the coordinates of plant species along the second axis of the coinertia analysis according to the four main forage types (broadleaf herbs, grasses, sedges and woody species)

3.2 Shift in vegetation community composition between 1976 and 2006

Between 1976 and 2006, the mean Ellenberg indicator value for nitrogen (N) increased both in the North Forest (V = 2135, p < 0.001; mean difference: 0.27) and the South Forest (V = 537.5, p < 0.001; mean difference: 0.16) (Fig. 5). In the South Forest, no other changes in mean Ellenberg indicator values were significant. In the North Forest, the mean value for humidity (F) increased significantly (V = 5889, p < 0.001; mean difference: 0.08) and the mean values for light and temperature decreased significantly (mean difference: −0.06, p < 0.01 and −0.03, p < 0.001, respectively).

Fig. 5

Mean Ellenberg indicator values (error bars standard error on the mean) per year and per forest (North Forest 248 plots, South Forest 73 plots). Significance of plot paired Wilcoxon rank-signed test (*p < 0.05; **p < 0.01; ***p < 0.001). Full symbol North Forest, open symbol South Forest, triangle 1976, circle 2006

The projection of the relevés conducted in 2006 revealed very significant shifts on both the first and the second axis (Fig. 6). On the first axis, the scores of the relevés were statistically lower in 2006 than in 1976 for the plots in the North Forest (paired t test: t = −9.56, p < 0.001, mean difference = −0.27, i.e. 50 % of the standard deviation of the coordinates in 1976 on the first axis). This was also true for the South Forest (paired t test: t = −6.87, p < 0.001, mean difference = −0.27, i.e. 55 % of the standard deviation of the coordinates on the first axis). The shift was towards wetter, more shade tolerant and nitrophilous vegetation communities in both the North and South Forests. On the second axis, the scores were significantly higher in 2006 than in 1976 in the North Forest (paired t test: t = 6.77, p < 0.001, mean difference = 0.14, i.e 37 % of the standard deviation of the coordinates in 1976 on the second axis). In the South Forest, scores were not statistically different (paired t test: t = −1.24, p = 0.2, mean difference = −0.04). There was a clear shift towards less browsing pressure and more high shrub cover in the North Forest, whereas the communities in the South Forest did not significantly change.

Fig. 6

Changes in the position of the 321 plots sampled both in 1976 and 2006 along the first two axes of the co-inertia analysis with bivariate SD-ellipses and their centroids grouping plots according to sampling year and forest (North and South separately). Full symbol North Forest, open symbol South Forest, triangle 1976, circle 2006. Embeded mean scores (and confidence interval) along the two axes of the North and South Forest; arrows indicate the trajectory from 1976 to 2006

In 2006, when much lower (Figure S1), deer browsing pressure was far less influent on plant community composition. A coinertia analysis conducted on 2006 data, similarly to the one conducted in 1976, reveals also two major gradients; as in 1976, deer browsing pressure contributes mainly to the second axis but with a low contribution (3 % of explained variance in 2006 whereas is was 30 % in 1976).

3.3 Correlation between shifts in community composition and changes in browsing pressure

The correlation between the change in the Ellenberg indicator value for light and the change in browsing pressure was significant and positive in the North Forest (ρ = 0.31, p < 0.001), but the trend was less significant in the South Forest (ρ = 0.22, p < 0.1). No other significant correlations were identified between the changes in Ellenberg indicator values and the changes in browsing pressure.

Both in the North and the South Forests, once all the changes in the other environmental variables were accounted for, changes in browsing pressure still explained a significant part of the change in community composition (“vector length” Table 1). In the North Forest, changes in the cover of high shrub layer explained the major part of the vector length and changes in tree cover were also highly significant. Once all other changes in the ecological variables used were accounted for, the change in deer browsing pressure was still a significant explanatory variable for the length of the shift along the first axis in the multivariate analysis (Table 1). Estimated standardized regression coefficients indicate that the decrease in the light and soil moisture indicator values influenced the intensity of the shift along this axis the most. However, changes in browsing pressure also explain a significant part of the upward shift along the second axis. For both the first and the second axes, the models invariably returned changes in tree and high shrub cover as highly significant variables. In the South Forest, the shifts in community composition were poorly explained by changes in the selected environmental variables. However, a significant relationship remained between the intensity of the shift along the first axis and browsing pressure.

Table 1 Relationships between changes in ecological conditions and plant community shifts (length of shift along axis 1 and axis 2, Euclidean distance in the first factorial map (axis 1, 2): VectLength)

Based on a linear model linking deer browsing pressure and plot scores along the second axis in 1976 (R ² = 0.38, p < 0.001), we computed predicted plot positions according to deer browsing pressure in 2006 (this latter variable was not included in the coinertia). These predicted positions were significantly higher than the projected positions of the 2006 vegetation relevés in the North Forest (paired t test: t = 14.9, p < 0.001, mean difference = 0.22), indicating a lag in the response of vegetation communities to changes in deer browsing pressure. In the South Forest, the difference was also significant but lower (paired t test: t = −3.03, p < 0.01, mean difference = −0.064).

4 Discussion

4.1 Eutrophication but no warming trend in vegetation composition

Over 30 years, the most obvious change in plant species composition was an increase in the nitrophilous character of the ground flora, in both the North and South Forests. The increase in mean Ellenberg values for nitrogen in the communities as well as the pattern of shift in the multivariate analyses indicates a higher frequency of nitrophilous species (see also Table S2). This directional change in species composition highlights an eutrophication of the forest ecosystem. This is consistent with the effects of atmospheric nitrogen deposition observed over the last few decades in Europe (Dirnböck et al. 2014) and particularly in north-eastern France (Thimonier et al. 1992). Moreover, ungulates intervene in the nitrogen cycle in various ways. Ungulate excreta, namely, urine and dung, contain urea nitrogen directly usable by plants, and its input affects plant species composition at the plot scale (Murray et al. 2013). High deer populations occurring in the past could have had a persistent impact on soil nitrogen dynamics in the Arc-en-Barrois forest; however, the total balance of the nitrogen flux remains uncertain as it depends on herbivore type (Pastor et al. 2006) and presumably on the time scale considered. A non-exclusive hypothesis is that deer could have introduced nutrients—notably nitrogen—into the forest interior from surrounding agricultural crops (Abbas et al. 2012). Since the North Forest is surrounded by crop fields, while the South Forest is included inside a larger forested area, this possible nitrogen transfer from croplands to forest habitats by deer and the consequences on plant species composition could be tested through the analysis of deer movements.

Mean annual temperature increased by 1 °C over the 30-year time period, but plant species composition did not show any shift towards more thermophilic communities. On the contrary, mean Ellenberg indicator values for temperature decreased significantly (although not very significantly) in the North Forest, and did not significantly change in the South Forest. In lowland forests, Bertrand et al. (2011) explains this low, even undetectable, response of plant communities to climate warming as a delayed response due to the plasticity of plant species to changes in temperature and geographical constraints that limit their migration. Moreover, plant community characteristics at the plot scale are also strongly dependant on microclimate, be it determined by topography or influenced by forest management (De Frenne et al. 2013).

4.2 Changes in forest management impacted plant species communities

The shifts in community composition were strongly influenced by changes in tree and high shrub cover. We identified a very significant positive shift along the first axis of the co-inertia analysis, which was not only correlated with trophic conditions (pH, water, nitrogen and temperature indicator values) but also with forest canopy closure (light indicator value, stand maturity). The decrease in mean Ellenberg indicator values for light indicates a closing of the forest ecosystem. This could in turn explain the concomitant decrease in mean indicator value for temperature and the increase in mean soil moisture indicator value. Canopy closure is the result of both an increase in the age of high forests and an abandonment of coppicing in coppice with standards regimes; these dynamics lead to increasing high shrub (including the dominated trees in the coppice) and tree cover, decreasing light availability for ground flora and changing microclimatic conditions. Consequently, the canopy closure that appears to be an important driver of the shifts in understory community composition could be attributed to the conversion from coppice with standards to even-aged forest that only begun in the late 1970s, a trend previously observed in lowland temperate forests (Van Calster et al. 2007; Baeten et al. 2009).

4.3 High deer browsing pressure shaped the composition of vegetation communities

In 1976, the dominant gradient in the composition of vegetation communities was clearly related to soil conditions (nutrient and water availability), which strictly limit forest stand growth on the Champagne-Ardennes calcareous plateau. This underlines that ground flora remains a reliable indicator for site conditions, even in a context of high deer populations and the heavy browsing pressure associated.

Nevertheless, browsing pressure, negatively correlated to high shrub cover, was the second gradient that structured ground-layer vegetation communities. Along this gradient, woody species and sedges were associated with low browsing pressure and, conversely, grasses and broadleaf herbs with heavy browsing pressure. This pattern is characteristic of the impact of deer browsing on vegetation in temperate forest environments (Rooney 2009). In addition, several species associated with high browsing intensity are known for their toxicity (Atropa bella-donna L. and Agrimonia eupatoria L.); this unpalatability provides them with a competitive advantage within the plant community (Augustine and McNaughton 1998). The negative correlation between browsing pressure and both tree and high shrub cover suggests that deer browsing contributes to open and to maintain open forest ecosystems by limiting the recruitment of shrubs over the browse line (Tanentzap et al. 2009), thereby impacting the composition of ground flora. Indeed, species associated with high deer browsing pressure are also mainly light-demanding species.

In 2006, the correlation between community composition (position of the plots along the second axis) and browsing pressure was weaker than in 1976. It was significant for the whole forest, but was not in the North or the South Forest considered separately. This result shows that deer browsing impacts the spatial structure or shape of forest ground flora communities when browsing pressure reaches high levels. Indeed, deer browsing is likely to be a significant gradient in vegetation communities when variations are strong. At low levels of herbivory, persistent gradients might be due to long-term legacy effects of high deer browsing on understory vegetation (Nuttle et al. 2014).

4.4 Response of vegetation communities to decreasing deer browsing pressure

Between 1976 and 2006, the communities in the North Forest moved significantly upward along the second multivariate axis, indicating that, over 30 years and following the decrease in deer browsing pressure, communities shifted back towards a composition similar to the more lightly browsed communities identified in 1976, namely, in the South Forest. In the South Forest, the shift was not significant, in agreement with the observed small variation in browsing pressure. According to forest managers involved in the present study, the South Forest can be taken as a reference, particularly in 1976 when moderate deer populations allowed forest regeneration. With this in mind, the shifts in the North Forest vegetation communities along the second axis suggest that the communities are recovering from a high level of deer browsing disturbance. However, our results reveal a lag in plant community response to a reduction in deer browsing pressure. Recently published studies of forest ecosystem monitoring following the reduction of (over-)abundant deer populations generally highlight the slowness of the recovery process (Tanentzap et al. 2009; Royo et al. 2010; Wright et al. 2012). In their review, Tanentzap et al. (2012) explain this slow recovery of forest understory through (i) the long-term consequences of plant biomass consumption that reduces individual plant species’ ability to react to reduced deer densities, (ii) consumption of reproductive structures and propagules that limit plants’ ability to reproduce and colonize and (iii) the formation of recalcitrant understories that prevent from the recolonization of other species, particularly palatable species. However, changes in canopy closure (see above) are also involved. Disentangling the effects of canopy closure and variations in deer browsing pressure on vegetation dynamics would require an ad hoc sampling design.

One practical application of our results is that plant community species composition could be used as a phyto-indicator of deer browsing pressure. Previous work on phyto-indicators have mainly dealt with quantitative measurements such as height or flowering abundance for specifically selected species (see review in Waller (2014)). The concomitant use of several species makes sense, however, since each indicator species is not always present and because the averaging of measurements on several species is likely to increase the reliability of the indicator. Recently, experimental work based on controlled browsing enclosures has identified species combinations that could be used as indicators of browsing pressure and has suggested using them to monitor ecosystem recovery after reducing herbivore densities (Bachand et al. 2014). Here, we raise the issue of identifying assemblages of species that vary along a gradient of browsing pressure, not only spatially as in the previous studies to date but also temporally. Our results suggest that this is possible only if there is enough variation in the intensity of the browsing pressure. Given the delayed response of vegetation communities to the decrease in deer browsing we observed, the kinetics of plant species dynamics appears to be a critical issue for such a subset of indicator species.

5 Conclusion

Ranking the main drivers of shifts in plant community composition is a key issue for the global assessment of the changes currently occurring in the biosphere and more specifically in vegetation communities (Corney et al. 2006). Yet, the deer browsing factor has not been precisely tested in many cases because of a lack of adequate measurements.

Our 30-year monitoring study of vegetation composition has revealed several directional changes in vegetation composition. The dominant shift was an increase in nitrophilous species frequency, pointing to the ongoing eutrophication of terrestrial ecosystems in Western Europe. Secondly, the canopy closure caused by the abandonment of coppicing and stand ageing induced a shift in vegetation community composition. Finally, a change in deer browsing pressure affected existing species assemblages. A decrease in herbivore populations, particularly in the North Forest, caused a clear shift in plant communities along the gradient of browsing pressure. It is noteworthy that the temporal shifts during the phase of decreasing browsing pressure occurred all along the spatial gradient of browsing pressure from the first date of observation. Obviously, plant community assemblages are sensitive to temporal changes in deer abundance, but even 25 years after the reduction in deer populations, their recovery appears incomplete. Assessing the resilience of vegetation communities to episodes of heavy deer browsing, the kinetics involved and the ability of communities to recover requires prolonging the monitoring of currently available long-term study sites.