Compatibility of regeneration silviculture and wild ungulates in a Mediterranean pine forest: implications for tree recruitment and woody plant diversity

Abstract

Key message

Small-scale forest interventions (< 0.75 ha) promoted advanced regeneration and woody plant beta-diversity without increasing ungulate habitat use and detrimental browsing damage. Rubbing damage by ungulates was higher in the treated areas and no effect was found on woody plant alpha-diversity.

Context

Adapted silviculture is needed to promote forest persistence and plant diversity in the current context of wild ungulate overabundance.

Aims

This study examines the ungulate effects on tree recruitment and woody plant diversity after silviculture treatments (small-scale regeneration fellings on Pinus species).

Methods

We compared tree recruitment, browsing/rubbing damage, and woody plant diversity on 17 pairs of control/treated areas in an ungulate-dominated Pinus halepensis forest.

Results

Recruitment levels were significantly higher in the treated areas as compared to intact (control) plots only for large saplings and juveniles (> 130-cm high). Ungulates did not use the treated areas more often than the control plots but caused significantly greater rubbing damage in the treated areas. Silvicultural treatments did not have a significant effect on alpha woody plant diversity but did promote beta-diversity, with a 49.7% woody species turnover. We did not find any clear patterns indicating that the treated areas suffered heavier browsing damage across all woody plant species.

Conclusion

This study highlights that small-scale forest interventions (< 0.75 ha) are small enough to avoid greater habitat use and browsing damage by ungulates but sufficiently large to promote advanced regeneration (large saplings and juveniles), with the additional benefit of increasing woody plant heterogeneity and structural diversity.

1 Introduction

Understanding the effects of wild herbivores on plant composition, tree regeneration, and forest persistence is crucial to define adequate management practices, particularly in the current context of increasing rates of wild ungulates (Côté et al. 2004; Massei et al. 2015). Forest management through silvicultural practices can help promote forest regeneration and persistence in areas where wild ungulates are abundant or overabundant (Heikkilä and Härkönen 1996; Bugmann and Weisberg 2003; Mansson et al. 2010; Edenius et al. 2014). In fact, new ungulate-adapted silvicultural practices have been proposed to reduce conflicts between wild ungulates and forest managers (Kuijper et al. 2009; Apollonio et al. 2010; Edenius et al. 2014; Faison et al. 2016a). In this context, a better understanding of ungulate habitat use, their preferences for plant species, and their overall effects on woody vegetation is crucial to create “opportunity windows” for tree regeneration (Reimoser 2003; Didion et al. 2009) and to achieve sustainable forest management (Morellet et al. 2007; Perea et al. 2015).

Silvicultural treatments create heterogeneity in the forest canopy, change the environmental conditions, and modify vegetation structure and composition (Decocq et al. 2005; Harrod et al. 2009). Changes in the understorey vegetation after silvicultural practices can be strongly conditioned by ungulate pressure, with browsing the most studied process (Gill 1992; Bugmann and Weisberg 2003). Several studies have analyzed the effects of ungulate browsing on tree regeneration to assess forest management (Bugmann and Weisberg 2003; Perea et al. 2016; Faison et al. 2016b). Browsing on young trees was found to be highly dependent on the silviculture treatment (e.g., logging), the stem density in the stand, or the species palatability (Faison et al. 2016a). For instance, in boreal forests, ungulate-adapted slash treatment increased the available forage biomass (Edenius et al. 2014) and felled pines have been used as supplemental food for wild ungulates, potentially reducing browsing pressure on young trees (Mansson et al. 2010).

Previous studies have also shown that gaps opened through regeneration silviculture are attractive to ungulates since the herbaceous cover as well as the number, diversity, and growth rate of young trees is higher than under the surrounding closed-canopy areas (Kuijper et al. 2009; Royo et al. 2010). Other studies highlight that woody forage quality is lower within the forest gaps, although this could be compensated by for greater volumes of lateral branches (Edenius et al. 1993; Hartley et al. 1997). In Mediterranean environments, the presence of shrub cover plays a crucial role in facilitating tree recruitment (Pugnaire et al. 1996; Gómez-Aparicio et al. 2008) and nurse shrubs have been used as a defense against ungulate browsing (Castro et al. 2002), particularly non-palatable evergreen shrubs (Perea and Gil 2014a). Rubbing damage by ungulates on young and mature trees has been shown to be important (Nielsen et al. 1982; Ramos et al. 2006; Gerhardt et al. 2013; Charco et al. 2016) but its relationship with silvicultural practices has not been extensively studied in Mediterranean environments. Thus, we still lack sufficient understanding of the interactions between ungulates and silvicultural treatments regarding ungulate habitat use and their effect on tree recruitment and woody plant diversity, especially in highly diverse systems (e.g., Mediterranean), where conservation of vegetation is usually the primary management goal.

This study investigates the ecological sustainability and compatibility of silvicultural practices in ungulate-dominated environments. Particularly, we examine the effects of a large wild ungulate (Ammotragus lervia) on a highly diverse Mediterranean vegetation following silvicultural treatments (regeneration fellings). These silvicultural practices are aimed at improving natural regeneration of the main tree species (Pinus spp., Quercus spp., and Juniperus spp.), as well as increasing the diversity of forest structure and woody plants. Specifically, we compare long-term plant damage (herbivory and rubbing) on shrubs and tree regeneration in treated areas and their adjacent closed-canopy areas (control plots). We hypothesize that tree regeneration will be more abundant in the treated areas. However, given the expected ungulate preference for forest gaps, we predict greater ungulate habitat use in the treated areas, causing higher browsing and rubbing damage as compared to the control areas. We also expect that regeneration structure will be more heterogeneous in the treated areas, with greater woody plant diversity as compared to control plots. Finally, we assessed whether browsing pressure is sustainable across all woody species and compared the proportion of species with unsustainable damage between treated and control areas. In this way, we aim to evaluate whether regeneration fellings are ecologically sustainable and compatible with wild ungulate populations to promote plant diversity and ensure the conservation of highly diverse woody systems.

2 Materials and methods

2.1 Study area

The work was carried out in Sierra Espuña Regional Park, a 17,804-ha protected area located in southeastern Spain (Murcia province; 2°4′-2°14′ N, 37°47′–37°57′ W). The elevation varies between 300 and more than 1500 m above sea level. The climate is Mediterranean semi-arid, with long summer droughts (> 4 months). Annual rainfall is usually low, between 200 and 500 mm, with strong seasonal and annual variability. Temperatures are warm, without frosts (thermo-Mediterranean thermotype), in the lower areas, but relatively cold in winter in the upper slopes and summits (supra-Mediterranean thermotype), where continentality is high and frosts are frequent between November and February. The meso-Mediterranean thermotype is the most abundant, covering 59% of the study area (Fernández-Olalla et al. 2016). The lithological substrate is basic, usually limestone, dolomite, or marle.

The vegetation is dominated by pines (Pinus halepensis Mill., P. pinaster Ait., and P. nigra Arn.), holm oak [Quercus ilex subsp. ballota (Desf.) Samp. = Quercus rotundifolia Lamk.], kermes oak (Quercus coccifera L.), and shrub-like junipers (Juniperus phoenicea L. and J. oxycedrus L.), with extensive patches of evergreen shrubs (Rosmarinus officinalis L., Lithodora fruticosa (L.) Griseb, Thymus vulgaris L., Cistus spp.). Xerophytic perennial grasslands, mostly dominated by Macrochloa tenacissima (L.) Kunth or Helictotrichon filifolium (Lag.) Henrard, are also widespread all over the area. The upper parts are covered by low cushion-shaped scrub formations under scattered pines and junipers. The vegetation shows high levels of biodiversity, as a consequence of both the high number of vascular species and the high rate of endemisms (Alcaraz et al. 2008). It also includes some protected habitat types included in the European Natura 2000 Network, such as natural and seminatural grasslands (Habitat codes 6110, 6220, 6170, 6420), Mediterranean scrubs (1430, 4090, 5210, 5330), and forests (9340, 9540) (Alcaraz et al. 2008).

The aoudad or Barbary sheep (Ammotragus lervia) is a large caprin ungulate (males up to 160 kg of body weight; females up to 90 kg). This species is native to large areas of semi-arid and arid mountain regions of northern Africa and is threatened as a vulnerable species (Cassinello et al. 2008). Its phylogenetic position has been considered intermediate between sheep and goat, although closer to goats (Valdez and Bunch 1980). Barbary sheep were released in 1970 in the study area with the main purpose of adding game diversity (San Miguel et al. 2011). Barbary sheep are the only browsing ungulate species in the study area with current density estimations of 5.9–8.4 individuals km⁻² (Eguía et al. 2015 for the period 2012–2014). Management of this population includes supplementary feeding, water provision, and population control (hunting).

2.2 Forest management and study sites

The main aim of the Forest Management Plan for Sierra Espuña Natural Park, approved in 2001, was to develop sustainable forest management according to the limitations of the semi-arid conditions. An adaptive silvicultural framework was implemented to enhance forest persistence and make it compatible with small-scale forest harvesting, biodiversity conservation, and park visitors.

The regeneration fellings analyzed in this study were organized in a group selection system to fulfill the following aims: (1) ensure and improve natural regeneration of the main tree species (Pinus spp., Quercus spp., and Juniperus spp.), (2) diversify forest structure, and (3) increase woody plant diversity of both tree and shrub species. The first implementation plan began in autumn 2001 and continued until 2010. In the group selection system defined for regeneration achievement, management practices consisted of felling all mature pines within a stand dominated by Pinus spp. Gap size, where cutting was carried out, varied from 0.05 to 0.75 ha depending on environmental restrictions and biodiversity requirements. In parallel, coppicing with standards was applied to Quercus spp. in the felled areas, and woody residues were added to the soil through chipping (Velamazán et al. 2006).

We selected six zones with similar environmental conditions (meso-Mediterranean thermotype) where Pinus halepensis was the dominant species. A total of 17 sites were selected within the six zones. Each site comprised two plots: the regeneration gap (with silvicultural practices) and the control area (with no management). Plots of each site were separated by at least 30 m from each other to avoid edge effects. We also ensured that the control plot had the same ecological conditions (slope, aspect, topography) and similar stand density and structure as the treated plot before cutting.

2.3 Data collection

Field sampling was carried out in three consecutive years (2013–2015), in late winter and early spring, the best season for quantifying browsing damage in Mediterranean woody species (Perea et al. 2014). Each site was measured only once within this 2013–2015 period, 6–13 years after treatment. In both plots of each site, we characterized the forest structure by recording the DBH (diameter at 130-cm height) of all standing trees in a 10-m-radius circular area, taken from the center of the plot. Young plants (< 7.5 cm DBH) of pines, holm oaks and shrub-like junipers were quantified in a concentric 5-m-radius area and classified in four regeneration categories according to height and DBH: 1 = seedlings, with plant height < 30 cm; 2 = small saplings, with plant height between 30 and 130 cm; 3 = large saplings, with plant height > 130 cm and DBH < 2.5 cm; 4 = juveniles, with plant height > 130 cm and DBH of 2.5–7.5 cm.

To estimate browsing damage by Barbary sheep on each woody plant, we established a 5-m-radius circular area within the above-mentioned 10-m plot (hereafter, browsing plot). This area (78.5 m²) corresponds to the minimum area concept (Braun-Blanquet 1951; Pfeifer et al. 1996) to survey woody plant communities in our study area (Fernández-Olalla et al. 2016). Browsing preferences were studied by comparing the utilization of every woody species occurring in the plot with their availability through an adaptation of the forage ratio index (Perea et al. 2014, 2015). The availability of each species was estimated with the Braun-Blanquet abundance scale (Braun-Blanquet 1951; Guisan and Harrell 2000), using the ground cover percentage of each plant species in the plot. Total shrub cover in each plot was also estimated with the Braun-Blanquet abundance scale. Utilization for each woody species was estimated by analyzing browsing damage using a 6-rank (0–5) browsing score: 0 = no browsing evidences; 1 = very few twigs (< 10%) browsed; 2 = low browsing, with 10–30% of twigs browsed and signs of regeneration (flowering/fruiting) in unconsumed twigs; 3 = intense browsing (30–60% of the twigs) but sustainable (flowering/fruiting signs); 4 = heavy browsing (> 60% of the twigs), with clear modification of the plant shape and no signs of flowering/fruiting; 5 = maximum browsing, no or almost no browsable twigs available. More details about this 0–5 rank and the preference index can be found in Perea et al. (2015). Unsustainable browsing was defined when browsing damage was > 3 since it clearly limits regeneration success and plant growth (Perea et al. 2014; Velamazán et al. 2017).

Plots were located using an accurate GPS (GARMIN etrex ® 10, Garmin Co., Kansas, USA), digital topographic maps (1:25000), orthophotos (50-cm resolution), and Geographical Information Systems (software ArcGis® 10.0, Esri, Redlands, USA). Digital maps and orthophotos were obtained from the Regional Geographic Institute (www.murcianatural.org).

To estimate the habitat use (relative density) of Barbary sheep in each of the browsing plots (5 m circular), we used pellet counting. Only fecal pellet groups containing six or more pellets were recorded by the standing crop method following Smart et al. (2004) and Acevedo et al. (2010). Additionally, in each plot, we also recorded rubbing damage on the tree regeneration as the percentage of young plants (DBH < 7.5 cm) with rubbing damage. All values are mean ± SE unless otherwise indicated.

2.4 Data analysis

All the analyses were performed using the R programming environment (R Core Team 2016; http://www.r-project.org/). To analyze regeneration abundance, we used a Generalized Linear Mixed Model (GLMM) where number of young plants was the response variable (count data; Poisson error distribution with a log link function). Treatment (control vs. cutting), tree species, and their interaction were the predictors. Site and age after treatment were included in the random structure of the model. We used a Cumulative Link Mixed Model (CLMM) to compare the browsing effects between treated and control plots. CLMM’s allow for regression methods similar to linear models while respecting the ordered categorical nature (0–5 rank) of our observations (Greene and Hensher 2010; Perea et al. 2015). Browsing damage for each species at each sample site was the response variable (ordinal variable). We included treatment (control vs. cutting), plant species, and their interaction as predictors. Site and age after treatment were included as random effects. The model was fitted by the Laplace approximation with a Probit link function, using the “clmm” function within the “ordinal” package (Christensen, 2013) for R 3.1.0 software (www.r-project.org). We used another GLMM to analyze rubbing damage (presence vs. absence of impact; binomial error distribution and logit function) with the same predictors and random effects. Relative density of Barbary sheep (pellet counts) was analyzed with another GLMM, where number of pellet groups was the response variable (count data; Poisson error distribution) with treatment (control vs. cutting) as the only predictor in the model.

We compared structural/species diversity between control and treated plots (pairwise). Structural measures included shrub cover and stand tree densities (through basal area). We estimated the stand basal area from the DBH of all standing trees within each plot. We measured structural heterogeneity through (1) variance of plot-level DBH for adult trees (> 7.5 cm DBH), (2) Shannon diversity applied to 5 cm DHB classes for adult trees, and (3) Shannon diversity applied to tree regeneration categories (del Río et al. 2016; Fahey et al. 2015). We calculated two woody species alpha-diversity indicators (Shannon index and species richness) and a beta-diversity index (Jaccard) for presence-absence data (Koleff et al. 2003) to compare more intuitively changes in species composition across all woody plants. Generalized Linear Mixed Models (GLMM) with different transformations and family distributions were used to analyze the following response variables: Shannon applied to DBH classes (squared root transformation; family = gaussian), variance of plot-level DBH (squared root transformation; family = gaussian), Shannon applied to tree regeneration categories (family = gamma), shrub cover percentage (family = binomial), Shannon index (family = gaussian), and species richness (family = poisson).

We used the Jaccard similarity index (Jaccard 1912; Koleff et al. 2003), which has a minimum value of zero (completely different communities) and a maximum of 1 (identical communities in terms of species presence/absence) using the “vegan” package of R. Similarly, we calculated all the possible pairwise Jaccard indexes across the 17 control areas [C (17, 2) = 136 pairwise combinations] and compared them with those obtained between the treated and control areas to see if there is higher species similarity across the control areas (at different ecological conditions) than between treated and control areas at each site. To validate these comparisons, we performed Exact Permutation Tests estimated by Monte Carlo (9999 replications), using the library “perm” and the function “permTS” in the R software.

Data availability

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

3 Results

3.1 Regeneration abundance

Overall, the abundance of tree regeneration was significantly greater in felled areas than in control areas (Table 1; Fig. 1). Regeneration abundance was higher in felled areas for all regeneration categories except for category 1 (seedlings), which showed similar values (Table 2). Regeneration categories 3 and 4 (plants with height > 1.3 m) showed the greatest abundance difference between treated and control areas, with regeneration densities seven times greater in the treated areas (Table 2). Most of the tree regeneration (all categories) was comprised of Pinus halepensis young plants (42.88%), followed by Quercus ilex (25.92%) and Juniperus oxycedrus (30.08%).

Table 1 Summary of the Generalized Linear Mixed Model (GLMM) to analyze the effect of treatment (regeneration fellings) on the abundance of regeneration of the three main species (Pinus halepenis, Quercus ilex, and Juniperus oxycedrus)

Fig. 1

Boxplot representing the abundance of tree regeneration (all categories) in the felled and control areas. Each box shows the median (band in the middle of the box) and the first and third quartiles (edges). Whiskers represent the lowest and highest datum within the 1.5 interquartile range of the lower and upper quartile

Table 2 Summary of the regeneration values (mean ± SE) for each regeneration category (1–4) and for each paired site (treated vs. control plot). N total number of individuals, BD browsing damage, R rubbing percentage (%). Plant categories: 1 = plants < 30-cm high, 2 = plants 30–130-cm high, 3 = plants > 130-cm high and DBH < 2.5 cm, and 4 = plants > 130-cm high and DBH of 2.5–7.5 cm

3.2 Browsing damage on regeneration

Browsing damage on young plants was generally low in both control and felled areas (Table 2). CLMM models revealed that, overall, there were no significant differences in browsing damage between the felled and control areas (Table 3). Browsing damage on Q. ilex was significantly higher than in Pinus and Juniperus (Table 3). In addition, browsing damage on young plants was only significantly higher in the felled areas for Quercus ilex (browsing damage of 1.83 ± 0.11 for the treated areas vs. 0.92 ± 0.12 for the control plots). Pinus halepensis showed greater browsing damage in the control areas (0.21 ± 0.06 for the treated areas vs. 0.44 ± 0.21 for the control; Table 3) as compared to Juniperus oxycedrus and Quercus ilex (Table 3). Juniperus oxycedrus showed very similar values of browsing damage in the treated and control plots (0.37 ± 0.06 for the treated areas vs. 0.36 ± 0.05 for the control plots).

Table 3 Summary of the Cumulative Link Mixed Model (CLMM) to analyze the effect of treatment (regeneration felling) on the browsing damage of the three main tree species (Pinus halepenis, Quercus ilex, and Juniperus oxycedrus). Values of tree species are shown against Juniperus oxycedrus

3.3 Rubbing damage

The percentage of trees damaged by rubbing was consistently higher in felled areas than in control areas (χ²₁ = 5.92; P = 0.015; Fig. 2). Interestingly, occurrence of rubbing damage increased with regeneration size in both treated and control areas (Table 2). No rubbing damage was found at all for the smallest regeneration categories (plants < 130 cm) in control areas (Table 2). Rubbing damage was mainly focused on Pinus halepensis (98.26% of all rubbed young plants; χ²₂ = 22.91; P < 0.001).

Fig. 2

Proportion of young plants (%) with rubbing damage in the felled and control areas. Whiskers represent standard error

3.4 Relative use by Barbary sheep

We counted a mean of 4.29 ± 1.18 fecal pellet groups per plot in felled areas and 3.53 ± 1.02 in control areas. The GLMM analysis revealed that ungulate relative density did not significantly respond to the treatment (χ²₁ = 0.02; P = 0.892).

3.5 Forest structure

The mean basal area in felled areas was almost half (15.34 ± 2.78 m²/ha) than that in the control areas (29.15 ± 2.46 m²/ha). Stand structural diversity indicators also showed significant differences between the control and felled areas. Mean diameter variance was higher in the felled areas (66.02 ± 14.51 in felled areas vs. 48.02 ± 6.28 in the control areas; χ²₁ = 17.54; P < 0.001) whereas the Shannon index applied to 5 cm DBH classes was greater in control areas (0.68 ± 0.14 in felled areas vs. 1.09 ± 0.13 in control areas; χ²₁ = 74.30; P < 0.001).

Regeneration structural diversity estimated by Shannon index applied to the regeneration categories was significantly higher in the felled areas than in the control plots (1.09 ± 0.05 vs. 0.67 ± 0.10; χ²₁ = 273.68; P < 0.001; Fig. 3). Finally, total shrub cover in treated areas was, approximately, 10% greater than in control areas (57.35 ± 6.57 vs. 48.21 ± 7.51%; χ²₁ = 17.54; P = 0.004).

Fig. 3

Regeneration structural diversity estimated by Shannon index for the felled and control areas. Whiskers represent standard errors

3.6 Woody plant diversity

A total of 41 woody species were found in the 34 plots, with 26 species in at least three plots (Table 4). Neither Shannon Index (t = 0.948; P = 0.350) nor Species Richness (t = − 0.863; P = 0.394) were significantly different between the felled and control areas. The Jaccard similarity index (beta-diversity index) between treated and control plot (pairwise analysis; n = 17 pairs) was 0.497 ± 0.033. However, the Jaccard similarity index between all control areas (pairwise analysis; n = 136 pairs of control areas) was 0.754 ± 0.138, indicating significantly higher similarity in woody plant composition across control plots than between treated and control plots on the same site (P < 0.001 estimated by 9999 Monte Carlo replications).

Table 4 Comparison of the number of individuals (N), ground cover (%), browsing damage (BD), and preference index (PI) for the woody species that were present in the treated and control plots. Only species that were found in, at least, three surveys are shown. Values indicate mean ± SD

3.7 Browsing damage on woody plants

We compared browsing damage on each woody plant species and six plant species showed unsustainable browsing levels (browsing damage > 3; Table 4). Two of them (Lonicera sp. and Bupleurum fruticosum) showed unsustainable browsing damage in both areas, felled and control (Table 4). Overall, five and three species showed unsustainable browsing levels in control and treated areas, respectively (Table 4).

4 Discussion

Our results show that, overall, regeneration abundance was higher in treated areas where regeneration fellings were applied. This positive effect of silvicultural treatments on natural regeneration was particularly strong for larger young plants (> 130-cm high), which indicates that silvicultural practices mostly favored the survival and development of large saplings but not that of small saplings and seedlings (< 130-cm high), suggesting that germination and establishment also took place in the intact (control) areas. This agrees with other studies that emphasize the positive effect of releasing young plants from competition by removing adjacent trees (Rodríguez-Calcerrada et al. 2008), particularly for light-demanding tree species (Zavala et al. 2011).

The results did not support our prediction regarding ungulate habitat use. Although we found greater number of pellets in felled areas, the difference was not significant, contrary to other studies (Kuijper et al. 2009; Wattless and Stefano 2013) which state that canopy gaps are visited by ungulates more often. The non-significant difference in our study could be due to the small size of the forest gaps, as their characteristics are probably not different enough from the forest matrix to become significantly more attractive for herbivores. However, as expected, we did find greater ungulate effect (particularly rubbing) in the treated areas (Table 2). Interestingly, occurrence of rubbing damage was proportionally higher than occurrence of browsing when comparing treated and control plots. This proportionally higher rubbing effect might be related to the greater abundance of large saplings and juveniles (> 130-cm high) in the treated areas, which were the preferred tree size for rubbing. In fact, no individuals below 130-cm high were found with rubbing damage in the control plots. Trees of intermediate size, such as saplings and juveniles, have also been demonstrated to be preferred for rubbing by other wild ungulates (e.g., red deer, Cervus elaphus) in similar Mediterranean pine forests (Charco et al. 2016). These results do not contradict previous studies that found more rubbing damage on intermediate size trees since they considered all diametric classes (regeneration and mature trees) in the stand (Johansson et al. 1995; Massei and Bowyer 1999; Ramos et al. 2006; Charco et al. 2016). Our results also highlight that rubbing damage by ungulates might be more important than browsing damage for the recruitment and growth of some species (e.g., pine trees in this study). Only for Quercus species did we find significantly greater browsing damage in treated areas as compared to control plots. These findings reveal the importance of discriminating between rubbing and browsing damage and the differential ungulate preference for conifers (mostly used for rubbing) in comparison to other species (mostly used for browsing).

As expected, the mean basal area in the treated areas was reduced by the implementation of regeneration fellings although structural diversity indicators for adult trees responded differently, depending on the index used. However, the structural regeneration diversity, estimated by Shannon index across regeneration categories, was higher in treated areas as compared to the surrounding forest canopy. We suggest that this positive effect was mostly due to the silvicultural treatment per se (increase of light and reduction of competition) and not to the ungulate pressure since we found very low browsing damage on pines (the dominant tree species) in line with previous studies in this area (Fernández-Olalla et al. 2016).

Apart from the influence on regeneration structure, silvicultural treatments did not significantly affect any woody plant diversity indexes (both species richness and Shannon). However, Jaccard similarity indexes revealed an approximately 50% replacement (species turnover) for woody plants. Interestingly, species turnover was even higher between control-treated areas on the same site than between control areas on different sites, suggesting that silvicultural treatments promote beta-diversity at a local scale.

Finally, we found unsustainable browsing damage (browsing score ≥ 3) on six shrub species (Bupleurum fruticosum, Genista valentina, Coronilla juncea, Coronilla glauca, Lithodora fruticosa, and Lonicera spp.) that are mostly considered highly preferred plant species for browsing by Ammotragus lervia (Fernández-Olalla et al. 2016 and Table 4). However, we did not find any clear patterns indicating that shrubs in treated areas suffered significantly heavier browsing damage. In fact, we found a higher number of species with unsustainable browsing damage in the control plots (five species) than in treated areas (three species). This might be due to the lower shrub cover and, thus, lower food availability, in control (closed-canopy) plots as compared to treated areas but further studies on the relationship between forest practices, shrub cover, and browsing would be desirable.

4.1 Conservation and management implications

This study confirms that small-scale forest interventions are suitable practices to improve biodiversity, as they promote structural heterogeneity with positive effects on trees and understory vegetation (Torras and Saura 2008). We highlight the fact that forest gaps smaller than 0.75 ha appear to be small enough to avoid higher habitat use by ungulates but sufficiently large to increase sunlight and recruitment growth. Although small-scale regeneration silviculture does not seem to increase woody plant species richness, it does increase beta-diversity and provides structural diversity and understory heterogeneity by replacing many woody species that thrive in the adjacent closed-canopy areas. In addition, small-scale interventions contribute to increased shrub cover, which is considered essential to facilitate tree recruitment in dry environments (Pugnaire et al. 1996; Castro et al. 2004) particularly in ungulate-dominated areas (Perea and Gil 2014b). Therefore, small-scale interventions might benefit from shrub encroachment to increase favorable microsites for seedlings and, hence, enhance the protection against ungulates.

In order to reduce the possible detrimental effect of wild ungulates (rubbing and browsing), we recommend selecting the regeneration areas sufficiently far away from supply points or attractive foraging areas (water, mineral correctors, feeding areas, crops, pastures, etc.) since relative high densities of wild ungulates have been proven to decrease regeneration probability and woody plant diversity (Perea et al. 2015; Velamazán et al. 2017). The use of standing dead trees or the addition of artificial decoy posts may be an adequate alternative to reduce the rubbing impact of ungulates on saplings and juveniles, which seems to be the most deleterious effect of small-scale interventions in ungulate-dominated environments. Further studies should analyze the effect of forest gap size on ungulate habitat use as well as the consequences for woody plant diversity and tree regeneration in order to estimate the most appropriate scale for intervention in areas where wild ungulates are abundant.