Moss cover and browsing scores as sustainability indicators of mountain ungulate populations in Mediterranean environments

Here, we aimed to define ecological indicators of environmental change for monitoring the effect of a reintroduced species, the Iberian ibex (Capra pyrenaica), over the vegetation of a natural protected area with contrasting habitats. This species was reintroduced 30 years ago in the Sierra de Guadarrama National Park, Spain, and its population has since grown exponentially, reaching high densities in some areas. We tested the suitability of two complementary indicators: browsing damage on woody species and the percentage of moss cover loss on rocky outcrops. For this purpose, we used a mathematical approach to zone the natural protected area according to the historical presence of the species and established five different Iberian ibex pressure classes (Classes I–V). Our results showed a direct link between Iberian ibex pressure and vegetation status. Model predictions suggested 10% moss loss in the areas with low pressure (Class I) but 64% loss in the areas with high pressure (Class III), mostly due to ungulate trampling on rocks. Similarly, browsing damage increased with increasing pressure classes for different woody species. Low palatable plant species such as Pinus spp. showed little change in browsing damage for increasing ungulate pressure whereas more palatable plants (e.g., Erica arborea) showed significantly greater browsing damage variation across pressure classes (i.e., high sensitive to herbivory pressure). Both, moss cover on rocks and widely distributed palatable woody plants proved useful ecological indicators for monitoring ungulate populations in mountainous areas with a diverse array of habitats (e.g., rocky vs. woody). This methodology may help select the most suitable indicators for each type of ecosystem or dominant habitat in ungulate-dominated ecosystems.


Introduction
Over the last few decades, populations of wild ungulates have experienced an unprecedented population increase across Europe to the point of becoming locally overabundant (sensu Caughley 1981) in many regions (Apollonio et al. 2010;Carpio et al. 2021). This situation has been mainly driven by anthropogenic global changes such as land use changes (i.e., rural abandonment, access to crop fields, extensive livestock decrease), reduction of natural predators and artificial management practices (i.e., supplementary feeding, fences, translocations; see Carpio et al. 2021 for a review). Furthermore, some natural protected areas or NPAs in Europe where ungulates have been reintroduced to restore ecological processes are now experiencing overabundance situations (Lovari et al. 2014) as many NPAs lack ungulates natural predators or any other population control measures such as culling or commercial hunting (van Beeck Calkoen et al. 2020).
At moderate density levels, reintroduced wild ungulates can help achieve the conservation objectives of NPAs as they contribute to landscape heterogeneity (i.e., creating and maintaining open grasslands; Rooney and Waller 2003), seed dispersal (Malo and Suárez 1995;Gill and Beardall 2001;Perea et al. 2013), soil nutrient cycle and wildfire control (Hobbs 1996), among others. By modifying plant species abundance and distribution , ungulates in turn benefit animal diversity (see Côté et al. 2004 for a review) and, for these reasons, they are considered keystone species (sensu Paine 1995). However, in NPAs with a long history of human landscape modification and low densities of natural predators, reintroduced ungulates can become overabundant and cause a cascade effect on lower trophic levels, exerting a negative, and often irreversible, impact on different components of ecosystems (Garrott et al. 1993;Apollonio et al. 2010). They may threaten the conservation of other native animal and plant species (Morellet et al. 2001;Lovari et al. 2014), create soil composition changes, erosion or even have a detrimental effect on various fitness components (Gaillard et al. 2000;Bonenfant et al. 2009) and overall health status (Gortázar et al. 2007;Putman et al. 2011). This is particularly worrying if we consider NPAs conservation objectives such as the protection of biodiversity and natural processes, which may not be compatible with a single species exerting a disproportionate negative impact on the system (van Beeck Calkoen et al. 2020). Therefore, in NPAs, it is critical to identify tools and indicators for monitoring and diagnosing the status of the environment after a reintroduction, to be able to establish, as soon as possible, priority areas for action.
The impacts of overabundant wild ungulates on native flora have been reported in multiple scientific studies both globally (Gill and Beardall 2001;Rooney and Waller 2003;Côté et al. 2004) and in Mediterranean environments (Soriguer 1983;Gómez 2005;Perea et al. , 2015Miranda et al. 2015;Fernández-Olalla et al. 2016;Velamazán et al. 2017). These studies indicate that, at high densities, native ungulates can cause negative impacts on the floristic composition and diversity, modifying the natural dynamics of the vegetation (Johnson and Cushman 2006;Seddon et al. 2007;Armstrong and Seddon 2008;Lovari et al. 2014;Perea et al. 2015). Furthermore, the effects of overabundance will be found first on plants, before any density-dependence effect is detected on the ungulate population (Noy-Meir 1975). Therefore, changes on vegetation can be used as indicators of ecological change, helping to monitor, through time, the impacts of large herbivores on the environment (Morellet et al. 2007).
Among plants, woody species are one of the most affected groups as they grow through apical meristems, on the contrary to most grass-like plants with basal meristems, and thus, they are less adapted to browsing (Kirby 2001). The low tolerance of woody species to herbivory causes a sequence of changes on the plant and on the overall ecosystem (Mysterud 2006), starting by changes in plant morphology, and eventually causing plant death and regeneration failure (Rooney and Waller 2003;Perea and Gil 2014), which, in extreme cases, produce biotic homogenization and diversity loss (Rooney 2009;.
Other groups such as bryophytes are also interesting in the context of global change as significant changes in bryophyte distribution can be expected due to climate warming (He et al. 2016). Bryophytes are important for ecosystem functioning, as they constitute the primary edaphic sustenance, controlling nitrogen and carbon cycles which are required for more evolved soils (Elbert et al. 2012). Furthermore, moss cover on rocky outcrops help diminish water runoff, stabilizing soils and trapping sediments and water (Evans 1996;Belnap 2006). In this context of global change, their distribution may also be affected by overabundant large herbivore populations (Rooney 2001;Moore and Crawley 2015;Morales-Molino et al. 2019). Although for many North American and European ungulate species it was reported little moss consumption (low palatability) and even a positive relationship between wild ungulates density and moss cover (Rooney 2001;Marozas et al. 2009;Chollet et al. 2013), there are some herbivores living at high latitudes or altitudes for which bryophytes may constitute a significant component of their diet (i.e., reindeer Rangifer tarandus and Soay sheep Ovis aries; Herbert and Prins 1982). Furthermore, bryophytes have proven to be sensitive not only to consumption but also to trampling (Chollet et al. 2013), which may be particularly true in environments with high percentage of rock cover and high slopes (Evans 1996). In Mediterranean environments, most studies on ungulatebryophyte interactions have been conducted in areas dominated by livestock (Concostrina-Zubiri et al. 2014 and little is known for wild ungulate populations thriving on Mediterranean mountains dominated by rock outcrops. In addition, assessing the impact of reintroduced wild ungulates on flora biodiversity and ecosystem functioning, or even their direct causal link to these changes, is rather difficult. This is because the negative consequences of their increasing abundance through time produce slow cumulative effects (Hui 2006) that cannot be properly evaluated unless the state of the system has been measured before the reintroduction. Unfortunately, this information is often not available. Another approach would involve the study of similar areas with different historic ungulate pressure (Fehér et al. 2014;Vild et al. 2017). This approach would be really effective if these areas could be located as close as possible (i.e., inside the same protected area), allowing the zonation of the area according to the historic level of ungulate pressure. In this sense, zoning NPAs may help evaluate the impact of ungulates on the environment by monitoring the historical presence of the species based on two main criteria: (1) the abundance of individuals on each area and (2) the time since the species was first recorded in a specific area. These zoning methods are frequently used in other disciplines such as environmental risk assessment (i.e., impact of energy and transportation infrastructures on wildlife; Tucker 1996; Goodale and Stenhouse 2016;Wolpert 2012), ecosystem services valuation (Warren 2011;Yokomizo, and Takimoto 2014) or to optimize resource investment in conservation (Benayas and de la Montaña 2003). However, these zoning methods have seldom been used to study the impact of wild populations over the environment (but see Pfab et al. 2018;Agüera and Byrne 2018;Akman et al. 2018).
On this study, we used a zoning method to assess the effect of a reintroduced ungulate, the Iberian ibex (Capra pyrenaica, Schinz, 1838) over the vegetation in the Sierra de Guadarrama National Park, central Spain, characterized by a marked elevational range. The species was reintroduced 30 years ago in the National Park (~ 67 individuals) and, since then, populations have reached high densities in some areas (Refoyo et al. 2015), 1 3 although the pressure has been very heterogeneous within the National Park limits (i.e., there has been a high spatio-temporal variation in the presence and pressure of the species). Therefore, we used current and historical information on the presence of the species in the National Park to zone the park according to ungulate pressure. Afterwards, and based on the ungulate pressure classification, we tested the suitability, as ecological indicators, of two different group of plant species (woody and moss species) that potentially have low tolerance to the pressure exerted by ungulates (browsing, trampling or rubbing).
For this purpose, we will: (1) zone the national park according to the current and historical presence of the Iberian ibex, its abundance and seasonality of habitat use to allow the establishment of pressure categories within the National Park; (2) analyze and compare the effect of Iberian ibex on woody vegetation and moss cover on rocks outcrops (3) define ecological indicators of population sustainability based on the relationship between ungulate demography and both groups of plants.

Study area and species
The study area is located within the boundaries of the "Sierra de Guadarrama" National Park, SGNP, (33,960 ha, central Spain), characterized by a marked elevational gradient, ranging from 1000 to 2400 m a.s.l. Climate is continental sub-Mediterranean, with high temperature variation between seasons but also along its elevational range (i.e., in the lower areas, hot and dry summers are typical Supra-Mediterranean while summers are milder and shorter in the higher areas, Oro-temperate). Average minimum temperature were − 0.4 °C and 17.0 °C (January and July, respectively) and mean annual precipitation was 1223 mm (Puerto de Navacerrada, 1894 m a.s.l.  (Lara et al. 2005), accounting for 12.6% of all the plant species in the region (Morales 2003). Geologically, SGNP offers a wide landscape diversity. It is formed predominantly by gneiss and granite substrates (Pedraza et al. 2014), which define the distinctive geomorphological features. The landscape comprises fresh outcrops, areas dominated by individual boulders partially buried in sandy deposits and areas composed by loose materials from the alteration of the geological substrate.
Iberian Ibex also known as Iberian wild goat (Capra pyrenaica) is a wild ungulate native to the Iberian Peninsula, with four known subspecies, of which only two are extant (C. p. victoriae and C. p. hispanica; Acevedo and Cassinello 2009). Nowadays, both are widely distributed along the Iberian Peninsula thanks to past reintroduction efforts (Refoyo et al. 1 3 2015; Herrero et al. 2021). In the SGNP, 67 individuals of C. p. victoriae were reintroduced in the 90 s from nearby areas. Since then, the population has grown exponentially (Refoyo 2012;Refoyo et al. 2015) as most predators were absent (only golden eagles and recently a few wolves) and no management controlled the ibex population in the SGNP (except some captures for translocations). Population growth has been monitored since the reintroduction using the distance sampling methodology (Buckland et al. 1993). Thus, the founding population of 67 individuals has increased from 359 individuals (year 2000), up to 5.403 in 2017, reaching up to 52 ind/km 2 in some areas (Refoyo et al. 2019). These densities are unprecedentedly high if we compare with other Iberian or Alpine ibex populations (0.4-15 ind/km 2 ; Acevedo and Cassinello 2009) or European populations of Capra ibex (0.2-16 ind/km 2 ; Carnevali et al. 2009;Dupré et al. 2001). Until year 2007, the species barely dispersed concentrating in the reintroduction area (approximately 4000 ha; Refoyo et al. 2015;Refoyo et al. 2016). In that area, the pressure on the environment has been constant and has produced sustained damage on woody plant as regeneration is highly compromised (see Perea et al. 2015) and significant soil loss and erosion (García-Rodríguez 2018). Afterwards, the species started to colonize progressively other areas of the SGNP, which allowed us to differentiate among areas according to the historic pressure of the species.

Zoning of SGNP according to Iberian ibex historical presence
Firstly, we stablished a square grid of 100 × 100 m cell covering the full extent of the National Park (33,960 units). For each cell of the grid, we calculated the following pressure indicators (Table 1) To estimate the value of each pressure indicator for each cell and survey, we used georeferenced data from all Iberian ibex observations on each survey. Further details on the Iberian ibex survey technique (years 2000-2017) can be found in Refoyo et al. (2015Refoyo et al. ( , 2019. Unfortunately, georeferenced data only provided an instantaneous picture of the quantity and specific location of individuals at the time of the survey instead of a grid map covering the full extent of the National Park. However, using these discrete data, we were able to calculate habitat use of the species for each survey, creating a continuous (gridded) map from survey points. We based our inference on Tobler Law, which assumed that cells that are close to one another are more alike than those that are farther apart. This methodology has been previously used to assess habitat use of species (Goovaerts 2010;González-Jiménez 2010). Thus, a geostatistical interpolation method (Inverse Distance Weighting; IDW) was performed based on the assumption that each value of a measured point or cell (i.e., real location of individuals or groups) has a local influence on its surrounding area that diminishes with distance. This interpolation gave greater weights to cells closest to the real observation location, diminishing the weights as a function of distance. For this purpose, we used the function "IDW" from Spatial Analyst Tools from ArcMap 9.3.
Finally, we calculated the overall pressure value (PV) of each cell by using an algorithm that related differently each of the five pressure indicators (Refoyo et al. 2019). The formula is as follows: Table 1 Description and values for each of the historic pressure parameters a Note that on alternative years population surveys were performed during different seasons. Spring surveys (years 2000, 2003, 2007, 2017), autumnal surveys (2005,2010) and summer surveys ( where for each cell, PV is the pressure value; N: the average number of individuals observed in a cell during all surveys; NC: the number of surveys where individuals appeared; PT: the population trend; S: seasonal presence; P: temporal proximity to the last record. We acknowledge the necessary subjectivity within the formula following Benayas and de la Montaña (2003).
Finally, based on its overall pressure value (PV), we classified each cell into five categories or classes of pressure (Table 2) using a dynamic table (Refoyo et al. 2019; Table S1).

Plant data collection
To assess ibex damage on the vegetation, a total of 95 survey plots were stablished across 8 different areas of the national park (Fig. 1). This survey plots were distributed proportionally to the total area of each class of pressure but considering a minimum number of 5 survey plots per class to allow statistical power (Table S2).
The plot area sampled to assess damage on woody plants and mosses was established according to the minimal area concept (Braun-Blanquet 1951;Pfeifer et al. 1996). In this case, a 10 m-radius circle (314 m 2 ), following previous studies on sub-Mediterranean woody vegetation (Perea et al. 2015). Within the entire plot, we recorded the number of pellet groups of Iberian ibex as a proxy for herbivore abundance. We also recorded GPS coordinates of each plot, slope, aspect, rock cover and elevation of each plot using a highaccuracy GPS (GARMIN GPSmap 62S), maps (1:25,000), orthophotos (50 cm resolution) and Geographical Information Systems (ArcGis 10.0).

Damage on the moss ground cover
The effect of the Iberian ibex over the bryophyte community was assessed by surveying the rock outcrops located inside each plot. Rocks are of great interest as the Iberian ibex heavily use them. Ungulates mostly cause damage to mosses through the mechanical effects of trampling (Kirby 2001;Pellerin et al. 2006) but also through consumption (Herbert and Prins 1982) where there is little availability of other food sources (i.e., browse or grasses). Thus, the possible damage produced by Iberian ibex was assessed by analyzing the surface of the rock layer. If a specific rock had patches where the moss had totally or partially disappeared, it showed a lighter colour that clearly marked the area and shape of the eliminated moss (see Supplementary material, Figure S1). Then, the percentage of the moss cover lacking was determined visually and with the use of images assigning a continuous value between 0 and 1. Value = 0 indicated that there was no moss loss in the rock surface and value = 1 indicated that 100% of the moss over the rock had been removed (García-Rodríguez et al. 2021). When the survey plot had little rocky outcrops, we used the nearest rock outcrops in the area to complete the plot area (300 m 2 approx.). However, some plots had no rock outcrops nearby and thus, they could not be assessed.
In addition, in order to attribute the moss damage to Iberian ibex, we met the following three criteria: (a) moss cover was present on nearby rocks with the same aspect but inaccessible for Iberian ibex, (b) browsing damage was found on woody species in the plot and (c) pellet groups were found in the area. Furthermore, for values below or equal to 20% of moss loss, we assumed that the changes could be attributed to other causes as we observed some degree of moss loss in some plots where Iberian ibex was not present. Finally, the disappearance of mosses could also be attributed to the influence of long drought periods. However, it has been observed that rocky outcrops retain all moss cover even in  1953,2228] south-facing areas of the SGNP below 900 m a.s.l. (drier and with less precipitation) where Iberian ibex is absent.

Herbivory damage on woody plants
The effect of Capra pyrenaica on the woody vegetation was studied through browsing evidences. For each circular plot, we estimated browsing intensity for each woody species, using a 6-points score (0-5) of an ordinal scale, following Perea et al. (2015): 0 = No browsing evidence; 1 = Light browsing, with < 10% twigs browsed; 2 = Low browsing, with 10-30% of twigs browsed; 3 = Intense browsing (30-60% of the twigs); 4 = Heavy browsing (> 60% of the twigs), with clear modification of plant shape; 5 = Maximum browsing (> 90% of the twigs browsed). However, in some plots, we also found presence of other ungulates that could have inflicted the observed damages to the plants (i.e., cow, horse and roe deer) and thus, we recorded the number of dungs or pellet groups found in the plot from each species. By far, cattle dung was the most often sign found in the plots (43 plots out of 96), after Iberian ibex pellets. Furthermore, both species were often found sharing the space. On the contrary, horse manure and roe deer pellets were only found in two and ten plots, respectively, with only one plot for each species co-occurring with Iberian ibex pellets.
In addition, in order to select a specific woody species as an indicator plant, the following conditions must be met (Perea et al. 2015;Fernández-Olalla et al. 2016): (1) a species that is widely distributed along the SGNP, (2) a species that is also vulnerable or sensitive to the pressure of the ungulates (i.e., medium or high preference by ungulates).

Plant data analysis
For the statistical analyses, we performed different models using the percentage of moss cover lost and the browsing damage as response variables. For the analysis, we used R statistical software (Version 3.5.2; R Core Team 2018).

Moss cover analysis
Firstly, we performed a linear model using the percentage of moss cover lost as response variable and the interaction between pressure class (factor) and altitude (continuous) as predictors. We included altitude as a covariable in the full model because pressure classes presented an altitudinal distribution in the SGNP (see Table 2). Comparing moss loss between areas of same altitude but different herbivore pressure allowed us to overcome this problem. We also included slope, aspect and percentage of rocks in the plot to control for possible differences in habitat use of the Iberian ibex. All continuous variables were mean centered to allow meaningful interpretation of the model intercept. We did not include cattle pressure in the analysis as cattle rarely climb rock outcrops. Similarly, horse and roe deer pressure were excluded from the analysis, as signs of their presence appeared seldom in the plots.

Browsing damage analysis
To select woody plant species as possible indicators of herbivore pressure, we selected those plant species that presented a widespread distribution across all the surveyed plots and thus, were abundant throughout the national park. Firstly, we selected species that were present in at least 25% of the surveyed plots (at least 24 out of the total of 95 plots). Secondly, from those species, we analyzed their altitudinal range and selected those that presented an altitudinal range greater than 614 m (50% of the difference in altitude between the lowest and the highest plot; [1172, 2400 m]).
Here, we performed a cumulative linear model (CLM) with browsing damage as response variable (factor). These models are usually performed with categorical response variables (i.e., 0 to 5 rank; Perea et al. 2015;Velamazán et al. 2017). We consider levels of browsing above 3 as unsustainable as the natural regeneration and thus species persistence in the system is highly compromised (i.e., woody plants presenting more than 60% of twigs browsed, clear shaping of the original plant form and no flowering and fruiting; Perea et al. 2015). Response variables were woody plant species (factor), pressure class (continuous) and their interaction, species (factor) and cattle pressure (continuous) and their interaction and altitude (continuous). Similarly, to the previous moss model, altitude was included to allow comparing browsing damage between areas of the same altitude but different herbivore pressure. Cattle pressure was included in the model to control for the possible damage inflicted on plants in the places where it co-habits with Iberian ibex. We did not include horse pressure because it was found in two plots, both with much higher cattle pressure. Similarly, due to its small size, we do not expect roe deer to have a high impact on the browsing damage results at such low density, especially in the plots where it shared the space with cattle. Finally, we performed model selection following the same abovementioned method.
To perform both linear models, we used the lm and clm function from the package stats and ordinal, respectively. To perform the model selection, we ranked all possible models derived from the full model based on the weighting provided by the Akaike information criterion (AIC; Burnham and Anderson 2002). For this purpose, we used the function dredge from the package MuMIn (Barton 2009). Finally, we selected the top model (that with lowest ∆AIC).
Finally, among the candidate indicator plants, we selected the species that showed highest slope in the model and highest overall browsing damage.

Moss cover analysis
During the survey, we found that, in the plots located above 2.000 m a.s.l., moss cover was minimal or null, even in the rocky areas inaccessible for the Iberian ibex (36 out of 95 plots). Therefore, effects on moss cover above this height could not be assessed. Unfortunately, areas with highest Iberian ibex pressure (Classes IV and V) occurred almost exclusively above 2000 m a.s.l. approximately, most frequently between 2100 to 2200 m a.s.l (25th-75th percentile; See Table 2). Thus, Iberian ibex damage to moss could not be assessed for Classes IV and V as sample size was very low (N = 3 and N = 1, respectively). In addition, there were four plots below 2000 m a.s.l. either with no rocky outcrops in the nearby areas or no moss found on the rocks and thus, those plots were not included in the moss analysis. For the pressure classes II and III, we could attribute most of the moss loss to the presence of Iberian ibex as plots characteristics met the three criteria [i.e., (a) presence of moss in inaccessible areas, (b) browsing damage in plants and (c) presence of pellet groups]. However, in some plots within the pressure class I, we observed some degree of moss loss (≤ 20%) that could not be attributed to the Iberian ibex as we found no presence of pellet groups in the area (criteria c).
After ranking all possible models derived from the full model, the top model was that containing Pressure class and Altitude but not the interaction between them. In addition, percentage of rock cover and aspect (both included as controlling variables) were selected in the model (LL = 26.919, AIC = -39.8, w = 0.210). Thus, the top model suggested that there was a significant difference in the percentage of moss damaged across all pressure classes (Table 3; Classes III vs. II t = 3.696, p < 0.005), while controlling by altitude, rock cover and aspect. Indeed, the lowest percentage of moss cover loss was found in low-pressure areas (Class I) while it was highest in the areas classified as Class III. Therefore, for an average altitude of 1700 m a.s.l., 50% of rock cover and 50°N aspect the model predicted a 10% moss loss for Class I, 43% for Class II and 65% for Class III (Fig. 2). In addition, we found that for the same pressure class, the moss cover loss was lower with increasing altitude (t = − 2.969, p < 0.005). Finally, we found a significant positive relationship between the percentage of rock cover and moss loss (t = 4.999, p < 0.005) but no significant relationship between aspect and moss loss (t = − 1.392, p = 0.171).

Browsing damage analysis
Across the 95 surveyed plots, we found 53 different woody species (see species list in Table S1). However, we only found four species that met the criteria of having a widespread distribution (i.e., being present in 25% of the plot surveyed and presenting an altitudinal range higher than half the range between the lowest and the highest plot surveyed). Thus, the selected species were Pinus sylvestris, Juniperus communis subsp. nana, Cytisus oromediterraneus and Erica arborea. Thus, the full model was performed using these four species.  . 2 Model predictions for the percentage of moss cover lost on rock outcrops according to the Iberian ibex pressure class. Covariables were mean centered: Altitude = 1700 m a.s.l., rock cover = 55% and aspect = 175 N. Error bars indicate model confidence intervals (95% CI). Points represent empirical data collected on moss cover lost for each Pressure Class. Significant differences in moss cover loss were found among the three pressure classes

3
After ranking all possible models, the top model included all the predictors and interactions present in the full model: pressure class, species (and their interaction), species and cattle pressure (and their interaction) and altitude (LL = − 150.664, AIC = 335.3, w = 0.907).
Results revealed that Erica arborea was the woody species with highest range of variation in browsing damage between low and high pressure classes. Cumulative linear model predicted a browsing damage equal to one for Class I and II, but two and three for Class III and IV, respectively (Fig. 3). On the other hand, Pinus sylvestris showed the lowest values of browsing damage and the smallest variation among classes, being cero for Classes I and II and one for Class III.
Results showed contrasting browsing sensitivity for each species, being Erica arborea the most heavily browsed at increasing pressure classes (showing medium-high browsing damage for pressure Class III), followed by Cytisus oromediterraneus and Juniperus communis (medium-low damage) and finally, the least palatable Pinus sylvestris (very low damage; Fig. 3). Therefore, the high browsing damage variation of Erica arborea along the pressure classes and its higher overall palatability make it a good candidate as a plant indicator of Iberian ibex pressure.
Finally, cattle pressure also influenced browsing damage, being its damage also significantly higher on Erica arborea than on the other plant species (Table 4). Fig. 3 Model predictions of the increasing browsing damage for each species while increasing pressure class. Model predictions controlled by altitude (mean centered at 1800 m a.s.l.) and cattle presence (cattle dung = 0). In order to be able to see all values for each Pressure Class and Browsing damage combination, we used the jitter option in the ggplot2 package as a way of handling overplotting caused by discreteness of the data. Note that we only predicted browsing damage results for pressure classes where each species was present

Discussion
This study highlights the potential usefulness of combining different ecological indicators to monitor the long-term effect of Iberian ibex in a Mediterranean natural protected area. Here, we propose the use of two contrasting plant indicators: the browsing damage on palatable woody plants for forests and shrublands and the percentage of moss cover loss on rock-dominated landscapes. Using these indicators, we have clearly shown the relationship between Iberian ibex pressure and their effect on the vegetation, derived from both, the direct consumption of species with low tolerance to herbivore damage (moss and woody species) and the trampling effect (mostly on moss cover). We believe both indicators, browsing damage and moss cover loss, will allow monitoring and comparing in space and time the direct changes on vegetation produced by wild ungulate species.
Results on moss cover loss indicate a progressive reduction of moss cover on the surface of the rocks below 1.900 m a.s.l. This loss, which is increasingly higher in the areas of greater Iberian ibex pressure, is, thus, related to the overgrazing and trampling of this wild ungulate. The SGNP has a wide range of moss in their catalogue, with about 300 taxa described (Blanco Castro and Acón Remacha 1984;Vicente and Ron 1989;Fuertes et al. 1995;Lara et al. 2005;Luna and Estébanez 2008) which protection should be guaranteed. Moss cover plays a key role in moderating essential ecosystem functions such as water infiltration, soil structure and stability, and nutrient cycling (Belnap 2006;Bowker et al. 2013;Delgado-Baquerizo et al. 2013). Similar to the results obtained at high latitudes or altitudes in temperate environments (Herbert and Prins 1982), mosses in mountainous Mediterranean environments were sensitive to the consumption and trampling of Iberian Ibex. This seems to be particularly true in environments with high percentage of rock cover (Evans 1996) and thus, the assessment of moss cover loss may be particularly useful in rocky environments where few woody plants grow. Results on woody species suggest that browsing damage is increasingly higher with increasing Iberian ibex pressure classes. Similarly, a previous study in SGNP reported woody species diversity loss in some of the areas with high historical herbivory pressure, where several woody species showed no sign of regeneration and unsustainable levels of browsing (Perea et al. 2015). Here, to select a useful woody plant indicator, we need to find an indicator species that show large differences in browsing damage across the ungulate pressure levels (i.e., high sensitivity to increasing herbivory pressure), which ultimately depends on the species palatability (Perea et al. , 2020. Furthermore, it is interesting to select species that are representative of the variability that exists in the diversity and cover of woody species. For this purpose, we have selected Erica arborea, a species that could serve as a good indicator of herbivore pressure as it is not only widely distributed across the whole protected area but also sensitive to changes in herbivory due to its relatively high preference by ungulates in the study area (i.e., medium-high palatability).
Below 1900 m a.s.l., the relationship is clear between Iberian ibex pressure and both vegetation groups and, what is more, both indicators complement each other (i.e., browsing damage in forest and shrublands and mosses in rocky environments). However, our study showed some limitations for the highest altitudes as the abundance of moss and palatable woody species is low. Indeed, the woody species above 1900 m a.s.l. are of lower palatability such as Cytisus oromediterraneus, which showed low variation on browsing damage across different pressure classes. Indeed, previous studies confirmed that moss cover and diversity is highest in mesic habitats (Hylander and Dynesius 2006). Aridity also had strong negative effects on moss height and cover (Mallen-Cooper et al. 2018).
Above 2000 m a.s.l in continental mountains such as those in SGNP, extreme environmental conditions become more frequent, increasing evapotranspiration, radiation and overall aridity, where, basically, only psicroxerophyllous (dry cryophilic) grasslands can grow (Rivas-Martinez et al. 1999). Furthermore, the high ungulate pressure experimented in these areas (many of them with a pressure Class IV or even V) might benefit these types of pastures as these plants have compensatory growth upon defoliation and productivity is maintained under high grazing levels (McNaughton 1979;Perevolotsky and Seligman 1998;Kirby 2001;San Miguel 2001). Thus, the conventional definitions of overgrazing may be inapplicable to these native plant-herbivore systems (McNaughton 1979) although even pasture-dependent systems might be damaged under very high densities (i.e., soil erosion or compaction). This may indicate that the effect of wild ungulates on native flora can vary according to the altitude and habitat, which requires a careful assessment of the damages.
Indeed, Iberian ibex browses more intensively at lower altitudes and during the winter, when there is very little quality food available due to snow cover or low temperatures (between November and March; Perea et al. 2015). It is mainly during this period of nutritional constraint when evergreen shrubs or even mosses are used as an alternative source of medium to low quality food (Evans 1996;Perea et al. 2015). In fact, the main differences between domestic and wild ungulates in alpine and subalpine environments are that the first are not present in the national park during the period of food shortage (they are present from late spring through summer when there is plenty of high-quality grass). In addition, domestic ungulates movements can be controlled by enclosures or moved more easily toward desired areas while wild ungulate move freely and need to meet all their nutritional requirements within the National Park (San Miguel et al. 1999). Therefore, the damages to native flora produced by Iberian ibex may be more difficult to control.
Nevertheless, further studies should investigate the importance of lichens cover and other organisms, which could be useful indicators in harsher environments such as cold 1 3 and elevated areas (i.e., above 2000 m a.s.l; Evans 1996;Marozas et al. 2009). In addition, the effect on the soil, estimated, for example, as percentage of eroded (compacted) cover, could also be a good indicator in certain areas, such as alpine pastures where woody species or mosses are scarce (García-Rodríguez 2019). Soil is key in alpine and sub-alpine environments, as it functions as the support for the rest of plants and organisms, and it takes hundreds of years to form, being very sensitive to erosive processes (Isselin-Nondedeu and Bedecarrats 2007; Konz et al. 2012). Further studies should focus on the selection of indicators of ecological changes at high altitudes, typically above timberline, where dry cryophilic grasslands are predominant.
Finally, we demonstrated that our zonation of the area according to the historic level of ungulate pressure is a useful tool for the management of a protected area as we were able to easily link specific ecological indicators to current and historic grazing pressure. Without this type of zoning method, the assessment of each area would have been much difficult as there has been a high spatio-temporal variation in the presence and pressure of the ungulate species in the national park. The specific application of this method to assess different biological groups and natural structures will allow the establishment of management criteria to better preserve the natural values of the NPAs.
As a conclusion, we want to highlight the importance of monitoring the populations of reintroduced ungulates as the consequences of their sharp population increase in humanaltered landscapes may pose an ecological risk for other species (Lovari et al. 2014). Here, we provided some applied conservation tools for the long-term monitoring of this species regarding their potential damage to vegetation. Firstly, we have confirmed the usefulness of highly herbivore-sensitive woody plants to monitor changes in herbivore pressure and ungulate populations in environments dominated by woody species (forests, scrublands, etc.,.) but discarded this type of indicators for higher alpine pastures or in highly rocky areas with little presence of higher plants. Secondly, moss cover loss has proven especially suitable in areas dominated by rock outcrops where few woody plants grow. Finally, we developed and used a spatial zoning method for the validation of different ecological indicators in order to assess of wild ungulate pressure on different ecosystems.
We believe that the long-term monitoring of both groups (vascular and non-vascular plants) will allow a finer and more direct assessment of the Iberian ibex pressure in the different areas of the National Park. These groups of plants can be used as early warning indicators to detect priority areas for action and thus, develop timely biodiversity conservation plans and population control measures to tackle Iberian ibex overabundance problems (see Martínez-Jauregui and Soliño 2021 for a discussion on population control methods and the conflicts between different social stakeholders). The use of complementary indicators with a zoning approach will help managers identify the priority areas to reduce/increase the herbivory pressure, with a particular focus on population control measures on the areas showing unsustainable levels of herbivory pressure (e.g., browsing degree > 3). Likewise, these indicators could be used, along with others, as a tool to monitor population control outcomes (e.g., through culling, translocation, etc.) and thus, as an adaptive management strategy for the future conservation of the national park.
de Guadarrama National Park and, especially to the current director, Pablo Sanjuanbenito, for facilitating fieldwork in the National Park and providing relevant information for this study. This work has been performed thanks to the funding received from Madrid Regional Government (Consejería de Medio Ambiente) and Spanish Ministry of Science through the INCREMENTO coordinated project (RTI2018-094202-BC21 and RTI2018-094202-A-C22). The project GLOBALFOR (M190020074RPGC) supported by UPM and Comunidad de Madrid (Convenio Plurianual for young researchers) also provided financial support. MP acknowledges the support of "Margarita Salas postdoctoral Fellowship" from the Ministerio de Universidades-UPM funded by the European Union NextGenerationEU.
Author contributions RP, MG and PR contributed equally to the study conception and design. Material preparation and data collection were performed by all authors. Zoning and pressure classes were performed by PR and statistical analysis were conducted by MP. The first draft of the manuscript was written by MP and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

Data availability
The authors confirm that the data associated with this manuscript and R-Script used for the analyses are available in supplementary material section.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.