Plant Diversity in Sardinian Mountain Rangelands: Analysis of Its Relationships with Grazing, Land Management, and Pastoral Value

In this study, we analyzed the effects of grazing on native and endemic plant diversity, as well as its relationship with pastoral value along a gradient of abiotic and biotic factors and types of land management in a mountainous area of central-eastern Sardinia, Italy. Plant diversity was estimated by conducting a floristic survey within plots. In total, 231 plant species were recorded in 63 plots distributed within the study area, and this total number included 20 endemic species. Species richness was mainly affected by the type of management, soil attributes, altitude, and bioclimate. Pastoral value was strongly affected by nutrient availability and bioclimate. Our results suggest that the cover of endemic species increases with altitude. Finally, in Sardinian rangelands, a negative effect of grazing pressure on endemic species was observed.


Introduction
Rangelands cover approximately one-third of the Earth's land area and at least one billion people depend on these lands to live, both directly through livestock production and indirectly from other resources and ecosystem services [1,2]. They are traditionally characterized by highly heterogeneous natural vegetation communities, which often have a high conservation value, frequently grow in harsh environments [3], are associated with wildlife and/or domestic grazing, and are managed by ecological/traditional rather than agronomic methods [4]. One of the most common forms of use of rangelands with domestic livestock is continuous grazing [5]. Nevertheless, rotational grazing (i.e., moving livestock among paddocks) has become common and has provoked discussion among rangeland ecologists and managers regarding the best management practices [6,7]. Consequently, the management and preservation of rangelands, which are important ecosystems for the maintenance of biodiversity [8], must be supported by a deep knowledge of natural resources and by an accurate level of grazing [9].
Rangelands in the Mediterranean area and on Sardinia are predominantly man-made grasslands, oak-wood grasslands, and oak-dominated woodlands [10], of which the stability and productivity largely depend on the amount of species diversity [11]. Importantly, there is increasing evidence that maintaining low or sustainable levels of grazing livestock is essential to preserve plant species and habitat types linked to pastures [11]. Furthermore, particularly in Sardinia, the presence of wild animals (e.g., mouflons and wild boars), overgrazing, soil-ploughing in natural or semi-natural land, mechanical bush removal, and high frequency of fires can lead to degradation of natural vegetation, particularly where soils are shallow, topography is mountainous, and the climate is thermo-Mediterranean [12]. In addition, the intensive exploitation of rangelands under high stocking rates promotes an increase in unpalatable species (e.g., poisonous, toxic, and spiny species), and a decrease in forage biomass and value, which are common indicators of land degradation [12,13].
The Mediterranean Basin is a renowned global biodiversity area with high plant diversity, sheltering 20% of the world's total floristic richness [14]. Thus, understanding the drivers that modify plant diversity, particularly in rangelands, and its relationships with grazing and land management is a fundamental question in applied ecology and an important factor for conservation purposes, management practices, and the adoption of suitable strategies for the preservation of plant diversity. The aim of the present study was to evaluate changes in plant diversity (native and endemic species richness) and in pastoral value (based on species pabularity) due to grazing pressure, abiotic and biotic factors, and land management types in a mountain area on a Mediterranean island. More precisely, the present research addresses the following questions: (1) How does plant diversity, plant community composition, and pastoral value respond to the effect of grazing? (2) Which abiotic, biotic, and management variables drive the highest differences in plant community composition?

Study Area
The research was carried out on Sardinia, Italy, the second largest island in the Mediterranean basin (24,100 km 2 ), in the mountainous area of the Ogliastra region, which is located in central-eastern Sardinia (Figure 1), at an altitude ranging between 290 and 1350 m a.s.l. The climate is typically Mediterranean, and is characterized by two main seasons, a hotdry and a cold-humid season [15]. The vegetation is composed of meso-Mediterranean silicicolous shrublands, Sardinian holm-oak forests and supra-Mediterranean holm-oak forests, Mediterranean riparian elm forests, Southern and Sicilian Italian Quercus pubescens woods, Mediterranean xeric grasslands, evergreen oak matorral and Juniper matorral, and conifer plantations, according to the Map of the Habitats of Sardinia [16].

Vegetation Surveys
Vegetation surveys were carried out in 63 plots, located according to a random sampling design, with field surveys from May 2013 to July 2015, when the floristic richness is higher and it is possible to observe higher plant diversity. The surveys were scheduled taking into account the altitudinal gradient, starting from the plots at the lower altitudes. The rangelands in the study area are grazed throughout the year, but more frequently from the beginning of spring until the end of autumn (i.e., for about 200-250 days of grazing per year). Each plot, sized 1 × 1 m, was georeferenced with a portable GPS device. The Ogliastra region was selected as a suitable study area due to its rich flora and due to its being predominantly characterized by extensively managed rangelands, which contain agro-forestry land used for grazing various animals (e.g., cattle, donkeys, goats, horses, pigs, and sheep) throughout the year. The investigated rangelands are partly included in territories managed by the Sardinian Forest Agency (FoReSTAS). These rangelands are also partly included within Special Areas of Conservation (SAC) and Special Protection Areas (SPA) according to the Council Directive no. 92/43/EEC of 21 May 1992 on the conservation of natural habitats and of wild fauna and flora.

Vegetation Surveys
Vegetation surveys were carried out in 63 plots, located according to a random sampling design, with field surveys from May 2013 to July 2015, when the floristic richness is higher and it is possible to observe higher plant diversity. The surveys were scheduled taking into account the altitudinal gradient, starting from the plots at the lower altitudes. The rangelands in the study area are grazed throughout the year, but more frequently from the beginning of spring until the end of autumn (i.e., for about 200-250 days of grazing per year). Each plot, sized 1 × 1 m, was georeferenced with a portable GPS device. Plant species presence and abundance were visually estimated. Plant abundance was recorded according to the Braun-Blanquet scale (r, +, 1-5).
Plant specimens were collected only if was not possible to identify the species in the field, then pressed, dried, and identified in the laboratory according to Arrigoni's Flora [15]. Endemic species were identified following Arrigoni's Flora [15]. This was also the reference used for the diverse biogeographical categories of the endemic species observed.

Estimation of Response Variables Plant Diversity and Pastoral Value
Species richness (native and endemic), and Pielou's evenness (J) were used to express plant diversity. Species richness was the number of plant species observed in each vegetation plot. The function 'diversity' finds the most used diversity indices and the R package 'vegan' was used to calculate the diversity indices, estimated in each plot [17].
The pastoral value (PV) for each plot was estimated using the methodology described by Vacca et al. [18]. The abundance value assigned to each plant species according to the Braun-Blanquet scale was converted into an average cover value per class, as proposed by [19], r = 0.1; + = 2.5; 1 = 7.5; 2 = 17.5; 3 = 37.5; 4 = 62.5; 5 = 87.5, to calculate the specific coverage coefficient (CSR). Therefore, PV was calculated as: PV = 0.2 × (ΣCSPi × SI), where CSPi is the presence specific contribution of a single species and SI is the specific index (from 0-species of no forage interest, to 5-species excellent for quality, palatability, and productivity) [20][21][22]. CSPi = (CRSI/CRStot) × 100, where CSRi is the specific coverage coefficient of eaten species and CSRtot is the specific coverage coefficient of the community. PV ranges from 0 to 100 and estimates the forage potentiality of a pasture area [21].
Finally, we considered for each plot the following functional plant traits: life form, life span, presence of thorns, and pollination type.

Estimation of Explanatory Variables Abiotic and Biotic Predictors
A set of predictor variables in GIS format was retrieved from the Sardinian regional geo-database (Sardegna Geoportale, http://www.sardegnageoportale.it, accessed on April 2020). Topographic variables-namely, altitude, aspect, and slope-were determined for each plot from a digital terrain model (DTM, raster 10 × 10 m resolution), using QGIS 3.16.3 software (QGIS Development team, 2018, Open Source Geospatial Foundation). The type of geological substrate was selectedfrom the following: intrusive rocks (granites, granodiorites, and leucogranites), metamorphic rocks (schists, arenaceous, and shales), effusive rocks (basalts), and acid-effusive rocks (andesites and rhyolites) [23]. The bioclimate map of Sardinia was used to identify the main bioclimatic categories found in the 63 plots of the study area (Supplementary Materials, Table S1) [24]. The Special Areas of Conservation (SAC) vector map (https://www.regione.sardegna.it, accessed on April 2020) was used to identify the plots that were inside SACs. The visual assessment of the presence and proximity of trees and shrubs for each plot was carried out and classified as presence and absence. The effects of past fires on the plant community were visually evaluated for each plot and classified as absent or very low, low, medium, or high.
Topsoil samples were collected at a depth of approximately 30 cm for each plot after removing litter, air-dried, and the analyses were conducted on the <2 mm soil fraction after sieving. Soil chemical characteristics were measured in the laboratory for pH, organic carbon (OC) using the Walkley-Black method, nitrogen (N) using the Kjeldahl procedures, and phosphorous (P) as P 2 O 5 using the Olsen method [25].
Grazing pressure was visually assessed (VLU) for each plot and classified as absent or very low, low, medium, or high. In addition to this visual assessment, the spatial variability of livestock units (LU) was evaluated through 49 pasture visual random observation sites (GPS located) within the study area. In each of the 49 sites, all grazing animals within a radius of about 1 km were recorded. Stocking rate by pasture was expressed in livestock units (LU), using the conversion values reported in Regione Autonoma della Sardegna-Programma di Sviluppo Rurale 2007-2013 for the different types of livestock (e.g., cattle value as 1 LU, pigs 0.3 LU, sheep and goats 0.15 LU). Spatialized livestock unit pressure (SLU), used hereafter as a proxy for grazing pressure, was performed between interpolated maps by livestock unit (LU) and ordinary kriging (contains the predicted values with the coordinate covariates) with the R packages 'gstat' [26] and 'kriging' [27]. Spatial interpolation used the data from the 49 sampling (observation) points to predict values of the observation variable at the 63 plots where data was not collected. We checked the relationship between the visual assessment and the spatial interpolation (ordinary kriging) and retrieved a non-significant association (p-value = 0.0688). Therefore, we decided to consider only the spatial value (SLU), which is a quantitative variable, in all of the following statistical analyses involving grazing pressure.

Statistical Analysis
Plant community composition at each plot was analyzed using canonical correspondence analysis (CCA) of the vegetation matrix (231 species × 63 plots). The species scores were weighted average plot scores. CCA is a gradient analysis that shows the relationship amongst vegetation patterns under the influence of the abiotic, biotic, and management predictors described above (i.e., analyzing the ecological gradient of vegetation).
Analysis of variance (ANOVA) was performed to test the relationship between the species richness, evenness diversity index, pastoral value, and the number of endemic species, calculated in the 63 sampled plots of the study area.
We used generalized linear models (GLMs) to assess the effects of abiotic, biotic, and management predictors (explanatory variables) on species richness, pastoral value, and the number of endemic species (response variables). Response and explanatory variables are listed in Table 1 and explained in more detail in Table S2 (see Supplementary Materials). A stepwise variable selection procedure was adopted. We used an information-theoretic model selection approach based on Akaike's information criterion (AIC) to identify the final model among all of the subsets of reduced models. Models considered best had a delta AIC higher than two compared to the next-ranked model.
In addition, we also used GLM to evaluate the effect of the SLU (i.e., a proxy for grazing pressure, considered as an explanatory variable) on functional traits (e.g., life forms, life span, and the presence of thorns) (all considered as response variables) within the study area.
The ANOVA results showed that the number of endemic species per plot was positively correlated with species richness. In addition, we found a negative correlation between endemic species and evenness index (Supplementary Materials, Table S4). Additionally, the GLM procedure highlighted that the percentage cover of endemic species was significantly positively affected by altitude and negatively by SLU and the number of endemic species was positively influenced by both altitude and SLU (p-values > 0.05).
Through the GLM, we found that species richness was positively affected by the abiotic predictors of altitude, soil attributes (OC and N), and iso-bioclimatic type 20. Conversely, soil P 2 O 5 and iso-bioclimatic type 31 resulted in a negatively correlated species richness. Interestingly, species richness showed a highly significant positive correlation with SLU.

Ordination Analysis
The overall results of the CCA ordination plot are shown in Figure 2. A correlation among plant species and the explanatory variables was observed, with 48.73% of the variance explained by the first two canonical axes of CCA. A plot of the first two axes shows two major trends. Axis 1 (eigenvalue = 0.262) was mainly dominated, with a positive correlation, by abiotic factors such as N and the iso-bioclimatic type 37, and negatively correlated with iso-bioclimatic type 28. Soil N is the variable with a long arrow, indicating that it was more strongly correlated with the ordination axis 1 and more closely related to the pattern of community variation than those with short arrows, whereas axis 2 (eigenvalue = 0.224) was distinguished, with a negative correlation, by abiotic factors such as P 2 O 5 and biotic factors such as the proximity of plots to trees. In addition, the plot shows that the plant communities' composition had a significant negative relationship with grazing pressure (Table 2, Figure 2

Pastoral Value
The pastoral value evaluated in the 63 plots of the study area ranged between 18.05 and 46.27 (mean = 29.43 ± 6.23). The results of the GLM are presented in Table 3. The most significant abiotic predictors that influenced PV were soil N and iso-bioclimate "lower mesomediterranean-lower subhumid-euoceanic weak" (isobioclim type 20, Table S1). Pastoral value was not correlated with the altitudinal gradient. The ANOVA results showed that the pastoral value was not significantly correlated with the number of endemic species per plot (p-value = 0.83) (Supplementary Materials, Table S4), or with the cover of endemic species per plot (p-value = 0.902). In fact, the endemic species recorded in the plots had a lower specific index (used to assess the PV) compared to non-endemic species, but often a lower percentage of cover (Supplementary Materials, Figure S2, Table S4). Table 3. Stepwise generalized linear model (GLM). The table shows the GLM analysis of pastoral value, species richness, and the cover of endemic species per plot (response variables) using abiotic, biotic, and management predictors within the 63 sampled plots of the study area. Significance codes: *** 0.001; ** 0.01; * 0.05. Iso-bioclimatic type 20: lower mesomediterranean, lower subhumid, euoceanic weak; Iso-bioclimatic 31: upper mesomediterranean, lower humid, semicontinental weak; N: nitrogen; OC: organic carbon; P 2 O 5 : phosphoric anhydride; SLU: spatialized livestock unit pressure.

Management Type
Of the total number of 63 plots, 38 were under FoReSTAS management and 25 plots were on common land, i.e., they were grazed without any management plan. The GLM procedure highlighted that the land managed by FoReSTAS had a significant negative influence on species richness (Table 3), whereas it did not have a significant effect on the cover of endemic species and pastoral value, or on plant community composition (Figure 2).

Effect of Grazing on Plant Traits
The grazing values (SLU) ranged between 0.11 and 115.9 (mean = 35.32 ± 23.63). We found, through the GLM, a significant positive effect of SLU on the percentage cover of hemicryptophytes and biannual species per plot (Table 4).

Discussion
The range of abiotic, biotic, and management conditions explored in this study allowed us to evaluate the relative importance of these factors in affecting plant diversity (native and endemic species richness) and pastoral value in a Mediterranean mountain rangeland, and reinforced previous studies that showed the effect of grazing pressure on plant community diversity in other parts of Sardinia.

Plant Diversity, Endemic Species, and Pastoral Value
In the present study, higher species richness was evident in plots at higher altitudes, with a high content of organic carbon and high nutrient availability (soil N). Soil attributes can have a positive relationship with plant diversity due to the creation of optimal growth conditions, and species richness is prone to change due to organic matter, providing a suitable place for the activity of soil microorganisms [30].
Mountains areas in the Mediterranean are characterized by the presence of endemic species [31][32][33], including narrow endemics. As already described for Sardinia [34,35] and other Mediterranean islands [36], and according to our results, altitude was the factor that mainly affected the cover of endemic species per plot, suggesting that the increase in endemic species is highly dependent on the increase in altitude.
We found a clear predominance of therophytes in all of the investigated plots, which is in accordance with previous studies in the Mediterranean region [37]. On the contrary, although the authors of a study [38] in a Mediterranean dwarf shrubland found a positive correlation between pH and species richness concerning annuals, geophytes, and chamaephytes, we found no significant correlation among species richness and pH.
Ordination analysis (CCA) showed that the plant community composition was mainly positively or negatively affected by abiotic (P 2 O 5 , N, and iso-bioclimatic type) and biotic factors (grazing pressure and proximity of plots to trees). Some of these variables have a prevailing effect on plant community structure on the island of Sardinia, according to Drissen et al. [39]. Nevertheless, the results of the CCA indicated that the abiotic and biotic factors considered in the study do not explain all of the floristic differences in the study area. Furthermore, our results highlighted that the pastoral value was influenced by different abiotic factors, e.g., nutrient availability (soil N) and the lower mesomediterranean-lower subhumid-euoceanic weak iso-bioclimatic type.
Management type-in our case land managed by FoReSTAS, and therefore subject to some level of regulation of the grazing activities-had a significant importance in conditioning species richness.

Effect of Grazing
Grazing is known to be a key factor for ecosystem dynamics in Mediterranean island landscapes [40]. Accordingly, in our study and that reported by Noy-Meir and Oron [41], grazing pressure promotes the abundance of geophytes and hemicryptophytes. The abundance of these two lifeforms is related to the production of bulbs, rhizomes, and buds protected in different ways at the ground level, thus acting as an effective barrier against grazing [42]. In addition, recurrent fires can promote these plant traits. However, we did not find any significant relationship in the present study, which can be related to the fact that land managed by FoReSTAS is also actively protected by the risk of wildfires.
In general, moderate grazing pressure optimizes pasture productivity, increases plant diversity [43], and might lead some species to become more dominant and others to become more relegated. On the contrary, overgrazing threatens plant diversity, especially on islands [44]. On the island of Sardinia, a negative effect of high grazing pressure on endemic species was demonstrated by Pisanu et al. [45], which is consistent with our results. However, in a study carried out by Camarda et al. [42], the authors found that 137 (45.4%) endemic species are not affected by grazing, whereas 116 (38.4%) endemic species are promoted by grazing. Grazing reduces the abundance of herbaceous palatable species, whereas it can promote the presence and abundance of species with defensive traits (e.g., spiny), and this may also promote a number of the endemic species, such as Astragalus genargenteus and Genista morisii. Therefore, the effects of grazing pressure can be twofold, both promoting the local richness of endemic species and reducing their cover.

Conclusions
Our results highlighted that altitude, soil attributes, and grazing are the main drivers of the plant diversity in the investigated Sardinian mountain rangelands, although the type of management had a notable influence, and pastoral value was highly affected by soil N. Hemicryptophytes were significantly affected by grazing, which must be considered an important component in the maintenance of the landscapes. The effects of altitude and grazing on endemic species should also be taken into account in rangeland management plans in Mediterranean mountain areas.
Supplementary Materials: The following are available online at https://www.mdpi.com/2673-4 133/2/1/9/s1, Figure S1: The pie charts show the number of species distributed according to the lifeforms, lifespan, presence of thorns, pollination type, and the number of endemic species in the 63 sampled plots of the study area. Lifeforms-Ch: Chamaephytes, G: geophytes, H: hemicryptophytes, He: helophytes, NP: nano-phanerophytes, P: phanerophytes, T: therophytes. Lifespan-A: annual, B: biennial, P: perennial. Presence of thorns-Y: yes, N: no. Pollination-A: anemophily, E: entomophily. Endemic species-Y: yes, N: no. Figure S2: Correlation between the presence of endemic species and the specific index for the evaluation of forage pastoral value (from 0-species of no forage interest, to 5-species excellent for quality, palatability, and productivity, according to Delpech, 1960). Table S1: The main iso-bioclimatic types, found in the 63 plots of the study area, extracted from the bioclimate map of Sardinia. The code defines the categories in the original vector map of the bioclimates in Sardinia (Italy). Table S2: Predictors used in the generalized linear models (GLMs), with their detailed descriptions. Table S3: Floristic list and plant traits. Table S4: Analysis of variance (ANOVA) was performed to test the effect of the number of endemic species (explanatory variable) on species richness, evenness diversity index, and pastoral value, calculated in the 63 sampled plots of the study area.