Mountain sheep grazing systems provide multiple ecological, socio-economic, and food quality benefits

Pastoral systems face increasing pressure from competing global markets, food sector industrialization, and new policies such as Europe’s post-2020 Common Agriculture Policy. This pressure threatens the use of extensive sheep-grazing systems in mountain areas of low productivity but high natural value. Using information gathered at a long-term research setting in a mountainous area of the Basque Country (northern Spain), we assessed the multiple benefits of extensive dairy sheep grazing systems from multiple perspectives using indicators pertaining to ecological, socio-economic, and food quality domains. In this way, we address the benefits that would be lost if sheep grazing abandonment persists in mountain regions. Our results show that the benefits of extensive dairy sheep grazing in the research area include the production of healthy and high-quality foods and multiple ecological benefits including biodiversity conservation. Extensive dairy sheep grazing also contributes to rural development by generating employment and income in marginal, low-productivity lands that can support few economic alternatives. In particular, we found that sheep farmers who produce high-value products, such as cheese, have enhanced their economic profitability and are less dependent on public subsidies. However, careful attention to sustainable practices, support for new generations of farmers, and streamlined supply chains are required. These would contribute to ensure socio-economic benefits for farmers, avoid the ecological costs associated with grazing abandonment, and enhance ecosystem services for the whole society.


Introduction
Rangelands make up 30-40% of the Earth's land surface (Sayre et al. 2013) and are intertwined with extensive livestock systems. These systems depend on grazing of marginal lands, providing multiple services, income, and employment for 1 to 2 billion people around the world (Herrero et al. 2013). Small ruminants alone provide vital resources to 300 million smallholders worldwide using mainly rangelands unsuitable for crops (FAO 2018).
Grassland-based systems include a variety of different management approaches in diverse habitats. This heterogeneity of contexts makes the management of grasslands a key driver of land use changes with multiple, often contradictory, consequences. For example, in some parts of the world, grazing is a main cause of deforestation, while in other areas, overgrazing of grasslands leads to land degradation (Herrero et al. 2013). Elsewhere, grazing abandonment threatens the existence of species-rich semi-natural grasslands (MacDonald et al. 2000). In general, grassland-based systems in temperate mountainous areas are used in marginal locations of low economic profitability but high natural and cultural value (Lomba et al. 2017;Navarro & López-Bao 2018).
Extensive grazing is threatened by intensification and abandonment, especially in low-productivity zones such as mountainous areas, which are characterized by ecological, economic, and political marginality (Sayre et al. 2013). In Europe, where traditional farming systems are in decline, grazing abandonment has been observed in various mountain regions (MacDonald et al. 2000). Factors leading to grazing abandonment include meager incomes, tedious working conditions, competing land uses, integration into global markets, and changes in Europe's Common Agriculture Policy; these factors have been discussed elsewhere in detail (e.g., MacDonald et al. 2000). The consequences of grazing abandonment are also diverse, although most studies have focused on associated costs and dynamics from a single perspective. Recent efforts have integrated multiple fields of knowledge into the study of extensive grazing systems (e.g., Bernués et al. 2011;Ripoll-Bosch et al., 2012;Benoit et al. 2019). Our aim is to address the full range of socio-ecological costs associated with mountain sheep grazing abandonment by assessing extensive grazing systems benefits from multiple perspectives that integrate all available knowledge.
With this aim, this study applies multiple ecological, socioeconomic, and food quality indicators to a mountain sheep grazing system at a long-term research setting in the Aralar Natural Park (Basque Country, northern Spain). This area, which is representative of Atlantic European temperate marginal lands, provides a unique laboratory for the study of mountain pastoral systems. In addition to a long tradition of sheepherding that goes back to Neolithic times, the area has a strong and dynamic shepherd's community, includes species and habitats of high conservation value, and produces highly valued sheep products (Andonegi et al., 2021).

Study area
Aralar Natural Park is an 11,000-hectare Special Area of Conservation (in the European Natura 2000 network) located in the Atlantic Basque Country of northern Spain. The area has an oceanic climate, with a mean annual temperature of 7°C and annual precipitation of about 1,330 mm. Although calcareous substrates predominate, the oceanic climate and abundant precipitation throughout the growing season have led to the development of acidic soils. The soils contain an average of 6,534 ppm Kjeldahl total nitrogen (SD 607.3;min. 4,500;max. 9,700) and 8.3 ppm available phosphorus (SD 2.7 ppm; min. 3 ppm; max. 23 ppm) and have a mean soil water content of 46.3 % (SD 4.7 %; min. 27 %; max. 65%) at a 10cm depth and a mean pH of 4.8 (SD 0.3;min. 4.2;max. 7.4) (Odriozola et al. 2017a). The vegetation in the park comprises a mosaic of gorse-heather shrublands and semi-natural grasslands.
Livestock herds are managed extensively in transhumance, using lowland farms in winter and upland grasslands from May to November. The number of livestock in Aralar Natural Park has remained almost stable for the last 25 years (17,000 Latxa ewes; average weight of a female Latxa ewe 55 kg); however, as in other Spanish regions, there has been a shift towards fewer sheep farm holdings with higher livestock loads. Most transhumance flocks range between 200 and 500 lactating ewes. Latxa sheep produce one lamb per year after a gestation period of 5 months. During the 140-day lactation period, a ewe produces an average of about 1.3 kg of milk per day. The milk is used primarily for making cheese: the Idiazabal Protected Denomination of Origin (PDO) is recognized worldwide as a high-quality cheese, but most of the production is sold in local markets and/or directly from farms.

Methodology: multiple indicators, data collection, and analysis
The application of multiple indicators in agri-food systems in general and pasture-based livestock systems in particular has a long trajectory (e.g., Bernués et al. 2011;Ripoll-Bosch et al., 2012;Botreau et al. 2014). These indicators can be used for assessments (e.g., Singh et al. 2012) or as inputs for multicriteria modeling approaches to choose, sort, or rank alternative options (see Gésan-Guiziou et al. (2020) for a recent review on the diversity and potentiality of multi-criteria decision analysis for agri-food research). In this study, we adopt the first approach and use multiple indicators pertaining to the ecological, socio-economic, and food quality domains to assess the benefits of a mountain sheep grazing systems from multiple perspectives. In this way, we address, with input from multiple disciplines, the benefits that would be lost if mountain sheep grazing abandonment takes place, a trend that has been observed in many European regions (e.g., MacDonald et al. 2000). To cover the wide range of benefits associated with these systems in a balanced way, we assess an equal number of indicators for each of the dimensions. Unlike modeling approaches, our assessment uses actual data collected over a long period from a research setting. The incorporation of multiple disciplines in the data collected at the research setting has occurred progressively. While the ecological assessment is based on more than 10 years of data, food quality and socio-economic data were incorporated into the analysis more recently. Data collection methods, including sample size and sampling period for each indicator, are provided in the following subsections.
The analytical strategy of this study, for identifying and evaluating the benefits of extensive grazing systems (see Fig. 2), considers the existence of emerging properties of complex systems at multiple scales, which provide nonequivalent information to understand the performance of the analyzed systems (Giampietro 2004). That is, when analyzing the performance of sheep dairy systems at the level of the whole (Level n in Fig. 2), we would analyze the total production of dairy products and consumption of resources. By moving down one level of analysis (Level n-1 in Fig. 2), we differentiate the relationship between a grazing system and the ecological context (e.g., by comparing grazing systems with non-grazing systems, as we have done, or with partial-grazing systems). To determine the nutritional benefits of dairy products associated with different grazing systems in this research area, we analyze the effects of mountain versus valley grazing during the grazing season (Level n-2 in Fig. 2). Finally, to understand the socio-economic effects of different production systems on mountain herders, we examine the two main typologies of mountain grazing production systems in the area: cheese and milk producers (Level n-3 in Fig. 2). Our study focuses solely on grazing systems through an analysis of Levels n-1, n-2, and n-3. A broader Level n analysis, which would assess the contribution of both grazing and non-grazing systems to the whole sheep dairy system, requires aggregated socio-economic and biophysical data related to the entire dairy system and is beyond the scope of this article.

Ecological indicators
Over a 10-year period, we measured four essential ecological indicators that are sensitive to grazing in Jasiono-Danthonietum grasslands. These indicators respond (positively or negatively) to grazing management or can be affected significantly by grazing disturbance (Odriozola et al. 2014(Odriozola et al. , 2017a.  & Species richness (S) is the number of species in a given surface unit. The effects of grazing on the biodiversity and species composition of plant communities depend on productivity, evolutionary history, and grazing pressure.
According to the generalized model, productive grasslands with a long evolutionary history of grazing have experienced divergent selection for species traits that increase resistance to grazing and the capacity to compete for light, resulting in a species pool of both short and tall species. Under such conditions, grazing intensity strongly affects both species composition and community structure, as rapid competitive exclusion of short species for light would occur with a hypothetical cessation of grazing. Grazing mediates spatial heterogeneity in grasslands by modulating plant inter-and intra-specific interactions; the patterns generated by species interactions are directly related to diversity (Deléglise et al. 2011). To determine the contribution of grazing to biodiversity, we assessed plant diversity by measuring species richness. This is one of the most frequent and recommended quantitative measures of diversity (Wang et al. 2019). The analysis was developed at a 0.25 m 2 -quadrat level, the recommended sampling size in grassland studies. & Soil microbial metabolic quotient (qCO 2 ) is the soil microbial CO 2 emissions (CO 2 flux, expressed in g CO 2 m −2 h −1 ) per microbial biomass (a measure of the mass of bacteria and fungi living in soil organic matter, expressed in mg microbial C kg −1 soil). This soil ecophysiological trait indicates the ecological well-being of the soil ecosystems and is used to monitor ecosystem health. As microbial biomass responds quickly to soil management changes, the soil microbial metabolic quotient provides information about the efficiency of microbial metabolism. Lower values reflect a more efficient microbial metabolism typical of mature, stable ecosystems, while higher values may indicate more immature or stressed communities (Aldezabal et al. 2015;Epelde et al. 2017). Since this quotient demonstrates the efficiency of microbial communities in utilizing C substrate in the soil, it also assesses indirectly the soil potential for C sequestration. & Soil compaction (resistance to penetration, Mpa) reflects the increase in bulk density or decrease in porosity of soil due to externally or internally applied loads. This indicator provides a useful measure of the soil's resistance to root growth (Correa et al. 2019). Soil compaction, which may be caused by trampling under high grazing intensity, can adversely affect nearly all physical, chemical, and biological properties and functions of soil. For example, compacted soils may hinder root penetration and consequently the development of rhizosphere and edaphic invertebrates, decreasing soil health. This indicator, in combination with the soil microbial metabolic quotient, provides a proxy to assess soil health. & Carbon/nitrogen ratio (CN ratio ), a quotient (nondimensional ratio) of leaf carbon content and leaf nitrogen content per leaf dry matter, is used to estimate forage nutritional value (a proxy of forage quality). Leaf nitrogen content is positively correlated to forage nutritional value and livestock productivity, while leaf carbon content is related to fiber content and is negatively correlated with livestock productivity. This means that the lower the carbon/nitrogen ratio, the higher the forage digestibility. Livestock grazing enhances the forage quality at moderately high stocking rates (i.e., 3.2 livestock units ha −1 d −1 : 13% beef cattle, 52% dairy sheep, and 35% horses; see Odriozola et al. 2014) by increasing white clover abundance in grasslands (Aldezabal et al. 2019). As a consequence, the production of enteric methane per unit intake diminishes, especially in ruminant livestock (Waghorn et al. 2002, Lee et al. 2004).
We estimated the ecological indicators described above using data collected in a manipulative experiment conducted in the study area. For this study, permanently fenced plots (50 m × 50 m each) were installed in May 2005 at four experimental sites to ensure grazing exclusion. Next to each of these exclusion plots, there was a grazed plot, where sheep, cattle, and horses were allowed to graze continuously from May to November (Odriozola et al. 2014) (see Fig. 3).
To estimate species richness, we located 100 sampling points systematically in both exclusion (E) and grazing (G) plots at three experimental sites (sites 1, 2, and 3). We took samples during the growing season of 2014, 9 years following livestock exclusion. Sampling points were spaced 2 m apart using the triangulation method, as a lag of 2 m proved suitable for pattern analysis in a pilot study of these grassland plots (unpublished data). At each sampling point, we measured flora composition and structure using two types of abundance metrics in quadrats (overlaid on the sampling points) of 0.5 m × 0.5 m: (i) species frequencies in 49 sub quadrats of 0.07 m × 0.07 m each and (ii) visual estimation of species percentage cover (Odriozola et al. 2017b). From these data, we calculated species richness for each experimental plot.
We measured soil health and ecological indicators (soil microbial metabolic quotient and soil compaction) in E and G plots at two experimental sites (sites 1 and 4) in 2010. After cutting the grass to ground level, we collected a total of 48 random soil samples (12 samples per site and plot) from a depth of 0-10 cm using a 3-cm diameter core. Each soil sample comprises a pool of 10 cores (subsamples). All laboratory analyses were carried out in duplicate within the following 2 months. We determined microbial biomass C, following the original method of Vance et al. (1987), and soil CO 2 emissions (CO 2 flux), using a PP-Systems® EGM-4 IRGA linked to a cylindrical soil respiration chamber SRC-1 (in triplicate around each sampling point, to give an average value). The soil microbial metabolic quotient (qCO 2 ) was inferred directly by dividing the CO 2 flux by the microbial biomass C (Aldezabal et al. 2015). To determine soil compaction (resistance to penetration, measured with Rimik CP40® cone penetrometer), we carried out five in situ cone index measurements around each sampling point and obtained an average value.
Additionally, we defined structural plant species as those which account for at least the 80% of the community cover and measured leaf carbon content (C), leaf nitrogen content (N), and CN ratio in seven randomly selected (7 in G and 7 in E plots), fully grown, undamaged individuals from each species (common bent Agrostis capillaris L., chewing fescue Festuca nigrescens subsp. microphylla [St. -Yves] Markgr.-Dannenb., heath bedstraw Galium saxatile L., field wood-rush Luzula campestris [L.] DC., and white clover Trifolium repens L.), following the procedure described in Cornelissen et al. (2003). After removing the petiole or rachis, we prepared powdered leaf samples and used these to determine the C and N contents using combustion elemental analysis (AOAC 1990). We obtained the CN ratio by dividing the C content by the N content per leaf dry mass. Finally, we calculated the community weighted mean, obtaining indicator estimates per site × treatment combination (for detailed sampling protocol, see Odriozola et al. (2014)).
We used the statistical package SPSS (Version 25.0, IBM Corp., Armonk, New York) for our statistical analyses. To investigate the effect of treatment (exclusion versus grazing) on ecological indicators (site was included as random factor), we used the generalized linear mixed model. For species richness, we used data available from Odriozola et al. (2017b). The raw data we used for calculating soil microbial metabolic quotient (qCO2) and soil compaction are presented in Supplementary Material (Table S1 and S2). For leaf nitrogen content (LNC), leaf carbon content (LCC), and LCC to LNC ratio (CN ratio ), we used raw data available from Appendix C (trait.csv) of Aldezabal et al. (2019). Statistical significance was declared at P ≤ 0.05.

Food quality indicators
To assess the contribution of an extensive mountain dairy sheep farming system in terms of healthy and high-quality food production, we focused on Idiazabal PDO cheese, the reference product of the area and main destination of the milk produced by shepherds in lowland and highland farms.
Cheeses made with milk from extensively managed flocks provide a higher content of healthy fatty acids (e.g., vaccenic, rumenic, and α-linolenic acids) as compared to flocks managed indoors and/or under part-time grazing that are supplemented with concentrate (Valdivielso et al. 2015). Moreover, as flock levels of concentrate and oil supplementation increase, the content of non-healthy trans-10 octadecenoic acid (trans10-18:1) in ruminantderived products increases (Aldai et al. 2013). This indicator assesses the cheese fat nutritional quality in relation to the forage intake of animals. The grass-related fatty acid accumulation (non-dimensional ratio) is calculated as the content ratio between certain unsaturated fatty acids in cheese, as follows: The fat healthiness indicator, a nondimensional ratio that assesses the healthiness of the cheese, is adapted from the atherogenicity index defined by Ulbrich and Southgate (1991). As previously reported, sheep milk fat from flocks managed by mountain extensive grazing shows lower atherogenicity values than that from flocks managed by indoor feeding and/or part-time grazing (Valdivielso et al. 2015(Valdivielso et al. , 2016b. Likewise, the ratio between polyunsaturated and saturated fatty acids is higher in milk from extensive mountain grazing flocks compared to that from valley grassland flocks (Bravo-Lamas et al. 2018). The fat healthiness indicator is calculated as the content ratio between unsaturated and saturated fatty acids in cheese, as follows: ) assesses the presence in cheese of antioxidant compounds beneficial for human health that come directly from the animal diet. Previous studies report significantly higher levels of tocopherols (vitamin E) and terpenoids in cheese from mountain-grazed flocks compared to those from flocks managed indoors and/or through part-time grazing, although the retinoid content was similar for each (Valdivielso et al. 2015(Valdivielso et al. , 2017. Antioxidant capacity is calculated as the content ratio between fat-soluble antioxidants in cheese, as follows: The sensory typing indicator (nondimensional ratio) assesses the typical cheese flavor in terms of sensory quality. The typical sensory characteristics of cheeses are strongly linked to the type of forage and to specific botanical species. Some studies have found significant differences in the sensory profile of the mountain cheese compared to valley cheese. Overall, Idiazabal PDO cheese made on mountain farms has more intense flavor and odor than valley cheese (Valdivielso et al. 2016a;Amores et al. 2021). Sensory typing is calculated as the ratio between specific sensory attributes in cheese, as follows: ST ¼ acid flavor þ 4 x rennet odor ½ þ 4 x piquant flavor ½ ½ = milky odor þ toasty odor ½ We based the food quality indicators on chemical and sensory data collected at 22 small rural dairies in the study area during the Idiazabal PDO cheese making seasons of 2009-2015. During 2009 and 2010, 10 farms located in valley grasslands were monitored monthly from February (indoor feeding) to July (extensive grazing) (Virto et al. 2012). In 2011, 6 flocks were monitored in February (indoor feeding), April (part-time valley grazing), and June (extensive mountain grazing), after having been moved from valley to mountain farms in May (Valdivielso et al. 2015). In 2015, a comparative study during May and June was carried out with 3 flocks located in valley farms and 3 other flocks that were moved to mountain farms. Both valley and mountain flocks were managed under extensive grazing, with valley flocks grazed on semi-natural pastures (Bravo-Lamas et al. 2018). The Jasiono-Danthonietum grassland dominates the grazing areas of the mountain flocks (Valdivielso et al. 2016b). In all studies, replicate raw bulk milk samples and cheeses ripened for 3 to 6 months were sampled from each dairy. Fatty acids, tocopherols, retinoids, terpenoids, and sensory attributes were determined in the samples according to previously described methods (Valdivielso et al. 2016a(Valdivielso et al. , 2016b. More details on analytical data from the different samplings used for estimation of food quality indicators are presented in the Supplementary Material (Tables S3-S6).
We performed a statistical analysis using the statistical package SPSS and used one-way ANOVA to investigate the effect of the livestock system (indoor feeding in valley, parttime grazing in valley, extensive grazing in valley, and extensive grazing in mountain) on food quality indicators. We included the livestock system as a fixed effect in the model and used the Tukey test for pairwise comparisons. Numerical values of all dependent variables were log transformed for statistical data treatment, and we checked normality and homoscedasticity of the data within each group according to the different levels of the main factors. Statistical significance was declared at P ≤ 0.05.

Socio-economic indicators
We measured four socio-economic indicators using data collected from 20 Aralar farms during 2016. The representativeness of this year was checked against data available for 7 farms since 2012.
& Net margin. Livestock grazing generates economic revenues in marginal lands of low productivity (Sayre et al. 2013). To assess such benefits, we estimated the net margin (in euros) generated by cheese, milk, meat, and livestock production at each local farm. We calculated net margin for a farm as follows: income (including revenues from sales plus other incomes such as subsidies)variable costs (feed, herd replacement, commercialization, veterinary)fixed costs (land rental, salaries, maintenance). & Employment. Small ruminants provide vital resources and employment opportunities in remote rural landscapes with limited alternatives (FAO 2018). For this study, we measured the employment generated directly by sheep grazing in annual work units (AWU). & Economic profitability. The long-term viability of mountain grazing systems depends on its ability to provide return to labor and total family income. To assess economic viability, we divided the net margin of each farm by the number of AWUs it uses. & Dependency on premiums. Another important socioeconomic feature of agricultural sector activities is their dependency on public subsidies. This is of particular relevance in European mountain grazing systems, where Common Agricultural Policy support income schemes are often employed to compensate activities that have multiple socio-ecological benefits but low productivity.
To assess such dependency in the study area, we estimated the percentage of subsidy over the net margin and compared it with other agricultural holdings in the same region.
Two types of dairy producers were identified in the study area: those who transform milk into artisanal Idiazabal PDO cheese and sell it directly to consumers (cheese makers) and those who sell milk to cheese factories (milk sellers). To capture the diversity of shepherds operating in the area, we surveyed 40% of them, 8 cheese makers and 12 milk sellers. Of these, 7 owned 150 to 300 sheep, 11 owned 300 to 500 sheep, and 2 owned more than 500 sheep. On average, cheese makers use approximately 55.9 ha of mainly communal land and produce an annual average of 34,611 L of milk, while milk sellers use 57.3 ha of land and produce 37,375 L of milk annually per farm.
We used the statistical package SPSS for our statistical analysis. To investigate the effect of farm typology (cheese maker versus milk seller) on socio-economic indicators, we used the General Linear Model of Analysis of Variance (GLM-ANOVA) as follows: where Y ijk = dependent variables, μ = intercept, D i = farm typology as fixed effect, L j = livestock (flock size) and M k = total milk production as covariates, and ε ijk = residual random effects. Numerical values of all dependent variables were log transformed for statistical data treatment. We checked normality and homoscedasticity of the data within each group according to the different levels of the main factors. Statistical significance was declared at P ≤ 0.05.

Ecological assessment
All ecological indicators measured at our experimental sites were affected significantly (P ≤ 0.05) by grazing abandonment (Table 1). Species richness was higher in the grazing plots, with an average of seven more species than in the exclusion plots. This increase in species diversity was due to the presence of more forbs and legumes in the grazing plots. This was especially true for Trifolium repens, which benefits from reduced competition for light in grazed areas.
The soil microbial metabolic quotient was around 0.27 units lower in the grazing plots than in the exclusion plots. Higher basal soil respiration (CO 2 flux: 1.55 mg CO 2 m −2 h −1 in exclusion plots vs 1.32 mg CO 2 m −2 h −1 in grazing plots) and lower microbial carbon biomass (1315.0 mg microbial C kg −1 soil in exclusion plots vs. 1451.6 mg microbial C kg −1 soil in grazing plots) were responsible for the higher soil microbial metabolic quotient observed in exclusion plots.
Soil compaction was higher in grazing plots, although noticeable variation was found among sites, probably due to differences in stocking rates or grazing intensities. In our study, soils under grazing regimes exceeded 3 MPa. The high grazing pressure in these grasslands clearly contributes to soil compaction, especially at site 4, where horses and cattle are present in higher densities because of the lower altitude. After 5 years, grazing exclusion allowed soil to recover to compaction values close to 1-1.5 MPa.
Finally, the carbon/nitrogen ratio was lower in grazing plots because of the higher nitrogen content of plants under a grazing regime (there were no significant differences in carbon content).
All indicators, except soil compaction, benefit from grazing (Table 1). Similar trends for the Jasiono-Danthonietum community are reported by Epelde et al. (2017), who observed that species richness decreased from 16 to 14 after short-term grazing exclusion. Likewise, Odriozola et al. (2017a) detected a negative linear relationship between species richness and abundance of competitive species, suggesting that strong aboveground competitors outcompeted other species (mainly forbs) when herbivores were excluded. This means that one important benefit of livestock grazing is the maintenance of plant diversity, which, in turn, enhances the forage quality of Jasiono-Danthonietum grasslands. Regarding forage quality, similar trends were observed by Aldezabal et al. (2019), who noted that an increase in carbon/nitrogen ratio after grazing abandonment was related to a reduction in white clover. Odriozola et al. (2017a) demonstrated that strong aboveground competitors (mainly tall grasses and long-spreading stoloniferous species) outcompeted small species such as white clover via a light-competition mechanism (by creating large intra-specific aggregations when herbivores were excluded). Sheep grazing exerts an equalizing effect on competing relationships among species because sheep preferentially consume grasses (more than 60% of their diet) and nitrogen-rich white clover to a lesser extent (Valdivielso et al. 2016b). Table 1, medium-to long-term grazing abandonment decreased soil compaction at the 0-10 cm soil depth, improving conditions for root development. This is consistent with previous studies (Epelde et al. 2017). Decreased soil compaction in exclusion plots could promote water infiltration capacity, gas circulation through the soil, and the accumulation of plant necromass due to the absence of herbivores. However, grazing abandonment also lowers the amounts of available nutrients (e.g., nitrogen, phosphorus) from animal excreta and urine, which changes the quantity and quality of root mass. Changes in organic matter dynamics and nutrient supply are likely to have a profound influence on microbial community structure as well as on plant physiological responses to grazing. Differences in soil microbial metabolic values for grazed and excluded plots suggest that soil microbial community functions are sensitive to the impacts of livestock grazing exclusion. According to Aldezabal et al. (2015), soil microbial metabolic quotients are lower in grazed plots due to reduced CO 2 emissions and increased microbial biomass. This effect could be related to higher soil temperatures observed in grazed plots (Odriozola et al. 2014), which also positively affect microbial biomass. In the absence of herbivores, a rise in microbial oxidative activity implies an undesirable increase in CO 2 emission into the atmosphere from microbial respiration.

As shown in
At this point, it is worth noting that a complete climate impact assessment would consider other sources of GHG fluxes (e.g., sheep enteric fermentation). For instance, in the study area, Batalla et al. (2015) found that the C-footprint of sheep milk in Latxa breed farms grazing in mountain and valley areas was greater (2.8-15 kg CO 2 -e/kg fat-and protein-corrected milk (FPCM)) than that from dairy sheep farms under confined management (2.3 kg CO 2 -e/kg fatand protein-corrected milk (FPCM)). Nevertheless, if soil organic carbon (SOC) changes are considered, the C-footprint of milk produced by Latxa breed farming systems can be significantly reduced, offseting GHG emissions, on average, about 55% (Batalla et al. 2015), of the total C-footprint. Last, but not least, in the light of pledges calling for afforestation of pastoral lands (e.g., mountain areas) as a means of mitigating climate change (e.g., Harwatt et al. 2020;Hayek et al. 2021), it is important to note that the resultant changes in C stocks or Table 1 Mean values and standard deviations for ecological indicators. a,b Means with different superscripts in the same row indicate statistically significant differences (P ≤ 0.05) between treatments (exclusion, grazing).

Indicators
Grazing treatment Exclusion treatment Species richness (number of species) 32.33 a ± 2.33 25.00 b ± 2.33 Soil microbial metabolic quotient (mg CO 2 m −2 h −1 /mg microbial C kg −1 soil) 0.93 a ± 0.28 1.20 b ± 0.26 Soil compaction (Mpa) 3.31 a ± 0.86 1.28 b ± 0.21 CN ratio 24.51 a ± 0.35 29.06 b ± 0.88 albedo are highly context-dependent (Bond et al. 2019), potentially counterproductive (Nuñez et al. 2021), and subject to large uncertainties as landscapes will be colonized by wild ruminants and other methane emitters (Manzano and White 2019). In terms of biodiversity, our results corroborate that livestock grazing exerts important and direct benefits for plant species richness, which also benefits the diversity of other aboveground and belowground organisms (e.g., insects, worms, microorganisms) (see also Wang et al. 2019). In this regard, enhancing shepherds' skills on guided grazing (Meuret and Provenza 2015) can help ensure these benefits across the entire landscape. Active livestock management plans, such as the new payments for ecosystem services schemes being promoted to increase grazing on abandoned (under-grazed) areas in the Navarran side of the park, support the conservation of biodiversity-rich Jasiono-Danthonietum grasslands. Table 2 shows mean values of food quality indicators, measured in milk and cheese samples from commercial Latxa flocks managed under different livestock systems on valley and mountain farms. Feeding management (indoor, part-time grazing, extensive grazing) and/or farm location (valley or mountain) influenced all indicators significantly (P ≤ 0.05).

Food quality assessment
As expected, milk and cheese from commercial flocks fed indoors on valley farms yielded the lowest values for grassrelated fatty acid accumulation and fat healthiness indices. These indices increased in flocks managed with extensive grazing systems, indicating a healthier fatty acid profile for milk and cheese from the latter. Furthermore, mean values for the grassrelated fatty acid accumulation index were significantly (P ≤ 0.05) higher for milk and cheese produced from extensive mountain grazing as compared to extensive valley grazing. In both types of extensive grazing, sheep were allowed to graze freely outside for the entire day, except during milking (twice a day, 12 h apart) (Bravo-Lamas et al. 2018).
The antioxidant capacity index is related to the presence of fat-soluble antioxidant compounds, which are beneficial for human health. Several studies have reported that plant biodiversity and botanical composition of mountain grasslands contribute to a high-quality animal diet (Falchero et al. 2010;Revello Chion et al. 2010;Valdivielso et al. 2016b). We found significantly (P ≤ 0.05) higher AC values in milk and cheese from extensively grazed flocks as compared to those from flocks managed under part-time grazing and indoor feeding ( Table 2). These results confirm the beneficial effect of extensive grazing management on the nutritional quality of milk and cheese, particularly when produced from mountain farms. The higher nutritional quality and healthiness of mountainproduced dairy foods are related to increased content of healthy mono-and poly-unsaturated fatty acids and fatsoluble antioxidants and decreased specific unhealthy trans and saturated fatty acids (Aldai et al. 2013;Valdivielso et al. 2015;Bravo-Lamas et al. 2018). Therefore, the abundance of specific botanical species in pastures and their fatty acid, tocopherol, retinoid, and terpenoid contents strongly influence the sheep milk and cheese composition (Valdivielso et al. 2016b).
With respect to sensory typing, the typical odor and flavor of Idiazabal PDO cheese have been described and defined. Our study is the first, however, to report scientific information on the specific sensory properties of Idiazabal PDO mountain cheeses. Table 2 compares the mean sensory typing values for ripened cheeses made on valley and mountain farms during the extensive grazing period. Cheeses produced on mountain farms show significantly (P ≤ 0.05) higher values than those produced on valley farms. The fact that Idiazabal PDO cheeses from extensive mountain grazing can be sensory differentiated from those made on valley farms may be an important factor for increasing consumers' willingness to pay more for mountain cheeses. Compared to valley farms, the incorporation of automatic milking machines, small cooling tanks, automatic or semi-automatic vats, and small ripening chambers with controlled conditions (temperature and humidity) in mountain farms could also contribute to improved Table 2 Mean values and standard deviations for food quality indicators measured in milk and cheese samples from commercial Latxa flocks managed under different livestock systems on valley and mountain farms. The number of farms (n) sampled in each livestock system was different depending on the year. a,b,c,d Means with different superscripts in the same row indicate statistically significant differences (P ≤ 0.05) among different livestock systems. 1 Index not calculated due to lack of tocopherol and retinoid data in milk and cheese samples. 2 Index not calculated due to lack of sensory attribute data in cheese samples.

Indicators
Indoor 18.6 ± 2.5 b ---1 46.9 ± 5.4 a Sensory typing (ST, non-dimensional units) ---2 ---2 7.8 ± 5.3 b 12.1 ± 7.6 a production of high-quality products that are better valued in the market. These technological improvements will undoubtedly improve the cheese making process in mountain facilities through increased control of the manufacturing process, better working conditions, and time savings for the farmer although they will also imply an additional economic cost for the shepherds. Finally, considering the intensification trend that has been observed in some European regions, it would be worthwhile to assess the benefits of extensive mountain sheep grazing in comparison to more intensive types of lowland grazing systems than those addressed in this study. A comparison could combine multiple indicators, such as the ones considered in this study as well as other indicators defined in collaboration with relevant stakeholders. Multi-criteria aggregation algorithms could be used to rank and compare different intensification levels of farming systems. In addition, further research is needed to assess the impacts of different degrees of abandonment in order to identify critical socio-ecological thresholds and tipping points at which the multiple benefits noted in this study might be lost irreversibly.

Socio-economic assessment
In socio-economic terms, the net margin and employment of the farms we analyzed show that sheep grazing systems have the potential to generate economic benefits in marginal mountain regions with harsh working conditions and few economic alternatives. In this regard, we also identify significant differences (P ≤ 0.05) among cheese makers and milk sellers for the selected indicators (see Table 3). Cheese makers show a superior net margin, due to the high added value of the cheese. This is reflected in a higher income from sales that compensate higher variable and fixed cost of cheese makers in comparison to milk sellers. The net margin of cheese makers is 2.5 times higher than of milk sellers.
On average, each farm employs 1.49 shepherds (AWU), which is close to the averages for Basque (1.7 AWU) and Spanish (1.6 AWU) small ruminant farms (RECAN (Red Contable Agraria Nacional), 2015). Because of their higher workload, Aralar cheese makers employ 1.8 AWU per farm, while milk sellers employ only 1.3 AWU (Table 2). These numbers are similar to those reported in the Basque Country for cheese makers (1.9 AWU) and milk sellers (1.5 AWU). On average, cheese makers require 216 ewes to generate 1 AWU, while milk sellers require 257. The 50 farms registered in the area, both full and part-time grazing, account for in total 99.5 AWU. These results show that, although mountain dairy sheep farming systems may have little impact on national economies, they can create jobs and increase income in remote mountainous areas such as Aralar. This study considered only the direct creation of jobs, but there are multiplier effects in local income and employment derived from both backward linkages with local suppliers and forward linkages stemming from successive expenditures created by increased sales in the region. In terms of employment, these linkages imply that for every job created in dairy sheep farming, indirect and induced jobs are created in other sectors. A multiplier analysis, based on Social Accounting Matrices, found that backward linkages are greater than forward ones for the dairy sector in Spain (Philippidis et al. 2014). The analysis estimated an employment multiplier of 13 in Spain's dairy sector. The additional jobs induced by farming vary at local and national levels but are usually low compared with other activity sectors (tourism, industry) (Hostiou et al. 2020). However, farms that sell directly in local markets, as those of Aralar, tend to have higher forward multipliers and consequently generate more income and employment in the region (Ekanem et al. 2016;Malagon-Zaldua et al. 2018).
The average economic profitability of cheese makers is 1.7 times higher than that of milk sellers and is higher than that of Basque (28,938 €/AWU) and Spanish (28,768 €/AWU) small ruminant farms (RECAN (Red Contable Agraria Nacional), 2015 ). The dependency of cheese makers on subsidies is in line with the average for Basque (30%) and Spanish (34%) . The dependency of milk sellers is much higher, almost double that of cheese makers, due to lower income from sales. The topography, climate, and remoteness of mountainous areas detract from the economic competitiveness of extensive grazing systems unless these are linked to food products with high added value (Madelrieux et al. 2018), such as Idiazabal PDO label cheeses. In the study area, the net margin per AWU of cheese makers is almost double that of milk sellers, and, the per liter price for milk used to manufacture cheese sold directly by farms is almost double the price of milk sold to dairy factories. This is in line with other studies (e.g., Renting et al. 2003) that address farmers' benefits associated with short supply chains and direct commercialization. Fresh milk sales direct to consumers could also increase farm incomes, but at present, this type of sales has no market niche. Despite its advantages, cheese making is not a silver bullet for all shepherds. Cheese makers assume higher variable costs (purchased feed), higher fixed costs (investment in equipment), and higher dedication in terms of working hours for manufacturing and marketing cheese. Note that the additional income is often not enough to hire extra personnel and/or to provide attractive salaries for qualified people. These requirements may discourage new generations from livestock farming . As in other Spanish regions, there are few young farmers in Aralar. The average age of Aralar farmers is 51 years, and only 11% of the shepherds aged 60 years and older have ensured generational replacement.
The high dependency on subsidies, especially in the case of milk sellers, makes the economic viability of mountain sheep grazing systems vulnerable to external changes. Without subsidies, the number of flocks might be reduced. Even so, the dependency of sheep grazing systems on subsidies is much lower than that for horse (mares) or dairy/meat production systems, whose subsidies exceed in both cases their net income. Considering social preferences for multifunctional agricultural landscapes and the multiple benefits of mountain sheep grazing systems documented in this and other studies (e.g., Dumont et al. 2019), subsidies may be justified. At the European level, payments for ecosystem services are gaining popularity as a way to recognize good livestock management practices that support high natural value systems benefiting the entire society. To promote socioecological benefits, however, subsidies should ensure farmers engage in sustainable practices (Navarro and López-Bao 2019). This requires active management planning, the spatial distribution of grazing livestock, or restoring abandoned grazing lands to avoid under-grazed areas.

Conclusion
This study confirms that extensive mountain dairy sheep grazing systems provide multiple ecological, socio-economic, and nutritional benefits. Interestingly, in contrast to the trade-offs between ecological and economic dimensions found in other studies (e.g., Ripoll-Bosch et al. 2012), we found significant synergies among ecological, socio-economic, and nutritional benefits. These results demonstrate that it is possible to enhance the economic profitability of farms based on sustainable extensive grazing practices that provide socio-ecological benefits for the whole society (e.g., high-quality food and biodiversity conservation). Indeed, the economic profitability is most remarkable for on-farm cheese makers, who produce a high-quality food product (Idiazabal PDO cheese) that depends on the maintenance and grazing of ecologically rich mountain pastures (see also Valdivielso et al. 2016b).
Even so, European extensive mountain sheep grazing faces significant challenges: a high risk of economic failure, an aging workforce, low social regard, competing land use, and integration into global markets (MacDonald et al. 2000). In addition to subsidies that will help sustain the socio-ecological benefits of extensive sheep grazing systems in mountain areas, measures are needed to attract new generations of sheep farmers (e.g., improved quality of rural life and working conditions), enhance short supply chains, and improve the marketing and recognition of value-added mountain products (e.g., through specific labeling).
Code availability Not applicable.

Declarations
Ethics approval Not applicable.
Consent to participate Not applicable.
Consent for publication Not applicable.

Conflict of interests
The authors declare no competing interests.
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://creativecommons.org/licenses/by/4.0/ .