2.1. Study area

Experiments were conducted in two semi-natural grasslands, Alotza (43°0'10.6" N, 2°5'22" W; 1,223 m elevation) and Uzkuiti (43°0′50″ N, 2°4′3″ W; 1,300 m elevation), located in the Aralar Natural Park, which is an 11,000 ha protected area located in the Atlantic Basque Country, Spain. The climate is oceanic with a mean annual temperature of 7°C and a mean annual precipitation of 1,330 mm. Aralar grasslands are used by mixed livestock in a moderate/high stocking rate (3.2 livestock units ha^‒1 d^‒1; 13% beef cattle, 52% dairy sheep, and 35% horses; Odriozola et al., 2014). Livestock herds are managed based on a short-distance transhumance system using lowland farms in winter/spring (December to April) and, more extensively, the upland grasslands in summer/fall (May to November). The park comprises productive native grasslands (mean aboveground net primary production of 1.97 ± 0.39 t dry mass ha^‒1 yr^‒1; unpublished data) with a calcareous substrate (Gibbons & Moreno, 2002). The grasslands (Jasiono laevis-Danthonietum decumbentis Loidi; Loidi, 1982) have a high conservation value (code 6230, 92/43/EEC; European Commission, 2013) and are associated with the production of a local cheese, Idiazabal, from the milk of the Latxa sheep breed (Batalla, 2015).

2.2. Experimental design

In May 2005, we installed a permanent 50 m × 50 m fence plot at each site (Alotza, Uzkuiti) to establish grazer exclusion plots (E plots). Attached to each E plot, we delineated a 50 m × 50 m plot (G plots), in which sheep, cattle, and horses were allowed to graze continuously during the vegetative period (May-November). Both sites are located in a flat terrain and have similar environmental characteristics and grazing pressures (Table 1). Species composition was not significantly different between E and G plots prior to fence establishment in May 2005 (Odriozola et al., 2017). Soil from each plot was analyzed in July 2010 to measure pH, Kjeldahl total N (mg L^‒1), available phosphorus (P; mg L^‒1), available potassium (K; mg L^‒1), calcium (Ca; mg L^‒1), and magnesium (Mg; mg L^‒1). Soil water content was measured at 15 cm soil depth and at 2 h intervals from June 2010 to October 2011, using Em50 data loggers connected to an ECH20 sensor system (Decagon Devices Inc., Pullman, WA, USA).

2.3. Plant sampling and traits

In June 2010, the floristic composition in each plot was determined in 15 randomly placed 1 × 1 m² quadrats. Species abundance was measured as canopy cover using the following index: presence (<1%) and a scale from 1 (1-10%) to 10 (91-100%). Plant species nomenclature followed that of the standard floras (Table Al).

In this study, we evaluated the traits of the five most abundant (structural) 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.), which covered 88% of the total area in Alotza G, 89% of the total area in Alotza E, 82% of the total area in Uzkuiti G, and 84% of the total area in Uzkuiti E. For the characterization of species in G plots, four exclusion cages of 1 m² were placed at least 20 m away from each other in May 2012 prior to the grazing period to collect plants with no signs of defoliation by grazers or of trampling. In July 2012, the cages were removed, and individual plants were collected from each plot for evaluation.

The following traits were used to assess the ecological response to grazing (i.e., response traits): plant height (H, cm) and specific leaf area (SLA, m² kg^‒1), since they each represent two key dimensions of plant function (Díaz et al., 2007, 2016). H is related to successional status and response to grazing (Díaz et al., 2007), and SLA is related to relative growth rate and resource utilization efficiency (Westoby et al., 2002). The traits used to estimate forage nutritional value (i.e., effect traits) were leaf N content (LNC, mg g^‒1) and carbon (C) to N ratio (CN). LNC is positively correlated to forage nutritional value and livestock productivity, and CN is related to fiber content and is negatively correlated with livestock productivity (Lee et al., 2017). SLA, LNC, and CN were measured in seven randomly selected, fully grown, undamaged individuals from each species, whereas H was measured in 25 individuals from each species using a metallic ruler (Cornelissen et al., 2003). SLA was defined as the one-sided area of a fresh leaf divided by its oven-dry mass. Hydrated and expanded leaves were scanned using an Epson Expression 10000 XL scanner (Epson, Suwa, Japan), and the area of each image (containing all the leaves of an individual) was measured using Win Folia 2008 (Regent Instruments, Ontario, Canada). The mass of the same leaves was measured using a precision balance after oven drying at 60°C for at least 72 h (Cornelissen et al., 2003). After removing the petiole or rachis, powdered leaf samples were prepared and used to determine the C and N contents using combustion elemental analysis (AOAC, 1990). CN was obtained by dividing the C content by the N content per leaf dry mass.

2.4. Data analysis

Distance-based redundancy analysis (dbRDA) was applied to test the null hypothesis of no effect of grazing abandonment on the abundances of the five structural species (Legendre & Legendre, 2012). The response matrix was created by Hellinger-transformation of the five-species composition matrix and computing Euclidean distances (Legendre & Gallagher, 2001). Analysis of variance (ANOVA) was fitted to each of the five structural species and four traits (H, SLA, LNC, and CN) to test the null hypothesis of no effect of grazer exclusion in the ITV of measured traits. Tukey's post-hoc test was used to separate the means. In all analyses, site (Alotza & Uzkuiti), treatment (G and E), and their interactions were used as explanatory variables.

The relative contribution of species composition and ITV to community response under grazer exclusion conditions was determined as described by Lepš et al. (2011). CWMs calculated using fixed trait averages (i.e., fixed mean trait value of each species for all individuals measured in the four plots) were only affected by changes in community floristic composition (species turnover), whereas CWMs calculated using specific trait averages for each field plot (i.e., specific mean trait value of each species in each plot) were affected by both changes in composition and ITV (i.e., ITV effect = specific average - fixed average). Linear models were used to test the null hypothesis of no effect of grazer exclusion on each community variable. Spearman correlation was performed to test the relationship between fixed and specific CWMs of response (H and SLA) and effect (LNC and CN) traits.

Variations in each trait were partitioned into fixed effects (species turnover), ITV effects, and their covariation by using the sum of squares (SS) decomposition of ANOVAs or linear models. Total variation in specific averages is the total variation in community trait averages (i.e., variation in fixed averages and ITV effects). Thus, if SS_specific = SS_fixed + SS_ITV, then changes in composition and ITV operate independently; if SS_specific > SS_fixed + SS_ITV, both effects are positively correlated; and if SS_specific < SS_fixed + SS_ITV, both effects are negatively correlated (Lepš et al., 2011). The R code and data are provided in Supplementary Data B and C, respectively.

3.1. Effect of grazing abandonment on structural species composition

The permutation test of dbRDA revealed significant effects of grazing abandonment on the abundance of the five structural species (Table A2). T. repens was clearly associated with G plots, whereas G. saxatile was associated with E plots (Figure 1). A. capillaris tended to be associated with E plots, whereas L. campestris and F. microphylla were not affected by the plot type (Figure 1).

3.2. Effect of grazing abandonment on ITV

The five structural species had high ITV in all plots regarding the studied traits (Tables 2, A3). In E plots, H and SLA showed high ITV since the individuals of the same species were considerably taller and had lower SLA (Table 2). In G plots, species also tended to have higher LNC and lower CN; however, the results were not consistent. For instance, G. saxatile had low LNC in the G plot of Uzkuiti, but high LNC in the G plot of Alotza (Table 2). Moreover, the LNC and CN of F. microphylla and T. repens were not affected by the plot type (Table A3c, 3d); in all plots, F. microphylla had low LNC and high CN, whereas T. repens had the highest LNC and the lowest CN (Table 2).

3.3. Relative contribution of species turnover and ITV to the response to grazing abandonment

Weak and inconsistent effects of grazing abandonment were observed when CWMs based on fixed trait averages were considered (Figure 2). H and CN tended to be higher, whereas LNC tended to be lower in E plots (Figure 2A, C, and D). However, SLA was reduced in the E plot of Alotza but increased in the E plot of Uzkuiti (Figure 2B). However, strong and consistent effects of grazing abandonment were observed when specific CWMs were considered (Figure 2). H and CN were markedly higher, whereas SLA and LNC were markedly lower in E plots (Figure 2). Additionally, strong correlations were observed between response (H and SLA) and effect (LNC and CN) traits using both fixed and specific CWMs (Table 3). Plots with high H had low SLA and LNC but high CN; accordingly, plots with low H had high SLA and LNC but low CN (Table 3). Species turnover and ITV selected similar H values after grazing abandonment, i.e., turnover and ITV effects were positively correlated; however, the contribution of turnover was almost negligible (Figure 3A). The difference in SLA between grazing regimes was mainly due to intraspecific variability; however, the interaction fraction (Site × Treatment bar in Figure 3B) showed negative covariance between turnover and ITV. This reflects the fact that while individuals of the same species presented lower SLA in the E plot of Uzkuiti, species with higher fixed SLA were more abundant in the same plot (Figure 2B). This inconsistency between turnover and ITV was not observed in Alotza (Figure 2B). For LNC and CN, turnover and ITV effects were important and consistent, i.e., both effects were positively correlated, and selected similar trait values (i.e., low LNC and high CN) after grazing abandonment (Figures 3C-3D).

4.1. Community response to grazing abandonment and effects on forage quality

The present study showed that the plant community had low aggregated values of H and high aggregated values of SLA under grazing conditions; however, this trend was reversed after grazing abandonment. These results were in agreement with those reported in previous studies (e.g., Cingolani et al., 2005; Cruz et al., 2010; Evju et al., 2009; Moreno García et al., 2014; Peco et al., 2017), suggesting that tolerance to grazing is the main strategy against defoliation by herbivores in productive grasslands with a long history of grazing (Cingolani et al., 2005; Díaz et al., 2007). Response traits H and SLA were strongly associated with the effect traits LNC and CN, revealing that the shift in trait composition following the change in grassland strategy after grazing abandonment (from tolerance to grazing, to aboveground competition) explains the decline in forage quality reported in these grassland sites (Odriozola et al., 2014). Theoretically, grazing tolerance requires high SLA, which increases shoot regrowth capacity (Westoby, 1999), and high LNC, which improves leaf quality and selectivity by herbivores (Cingolani et al., 2005). However, light competition requires increased H to capture the light and a conservative resource utilization (i.e., low SLA) for a higher investment in structural components (i.e., high CN) (Cingolani et al., 2005).

The production of forage of high quality under grazing conditions has been also reported in other sub-humid and humid grasslands where tolerance is the plant community strategy against herbivores (Cingolani et al., 2005; Cruz et al., 2010; Deléglise et al., 2015). This has also been observed in the semi-arid grasslands of South Africa (Moreno García et al., 2014), although intense grazing areas coincided with water points in this case. By contrast, in the mesic grasslands of France, where tolerance is combined with avoidance, some species showed reduced SLA and LNC with cutting (Pontes et al., 2010, 2007). Our results were in agreement with those of humid grasslands since the mean annual precipitation in Iberian Atlantic grasslands is 1,330 mm. Fertilization favors a resource acquisition strategy and may also improve forage quality (Pontes et al., 2010, 2007). Since defoliation by herbivores is always coupled with fertilization by urine and feces (Lezama & Paruelo, 2016), community response is probably a combination of tolerance to defoliation and response to fertilization, but these two effects cannot be disentangled in our study.

4.2. Species turnover and ITV effects on the response to grazing abandonment

In agreement with the results of previous studies, the traits measured in the present study showed high ITV (Albert et al., 2010; Cornelissen et al., 2003; Lepš et al., 2011; Westoby et al., 2002). Pontes et al. (2007) observed higher between-species variation than within-species variation. We also found that species could be distinguished based on the mean values of traits; however, the exclusion of ITV resulted in a clear underestimation of community response to grazing abandonment for all traits. Similarly, Lee et al. (2017) claimed that intraspecific physiological response is an important driver of forage quality worldwide.

We observed a mismatch between species turnover and ITV effects on the response to grazing abandonment for H and SLA; the effects of species turnover were very weak, whereas those of ITV were very strong. The lack of clear association of the tall grasses F. microphylla and A. capillaris with E plots along with the abundance of the short dicot G. saxatile might explain the weak species turnover effect. Five years after grazing abandonment, opposite effects of species turnover and ITV were observed for SLA in Uzkuiti, as indicated by the negative covariance between species turnover and ITV (Figure 3B). This response was probably a result of the high abundance of G saxatile (very high SLA) in the E plot. However, long-term changes in floristic composition were observed after grazing abandonment (nine years in both sites): F. microphylla and A. capillaris became the most abundant species in the E plots (Odriozola et al., 2017). Our results are in agreement with those reported by Lepš et al. (2011), who suggested that ITV determines the resistance and stability of grasslands to short-term changes, while species turnover reflects processes associated with the resilience of grasslands to long-term changes.

The contributions of species turnover and ITV were more equilibrated for LNC and CN. T. repens, which showed very high LNC and very low CN (Table 2), was the most sensitive species to grazing abandonment and was clearly more abundant on G plots. This species is highly digestible and nutrient-rich, and is widely cultivated for animal nutrition (Parsons et al., 2011). Our results demonstrated that T. repens is fundamental for maintaining the high nutritional value of Iberian Atlantic grasslands.

4.3. Management implications

Previous studies have reported reductions in forage quantity along with improvement in quality under grazing conditions (Cruz et al., 2010; Deléglise et al., 2015; Moreno García et al., 2014). Accordingly, an increase in the annual aboveground net primary productivity has been observed in Iberian Atlantic grasslands after grazing abandonment (unpublished data). Cruz et al. (2010) claimed that the reduction in forage availability under grazing conditions is incompatible with the high animal productivity, whereas Moreno García et al. (2014) argued that the improved forage quality might compensate the reduction in quantity. It is widely accepted that both forage digestibility and nutritive value have direct effects on animal performance (Herrero et al., 2015; Lee et al., 2017). Additionally, although short-term grazing abandonment positively affects grassland productivity, long-term grazing abandonment has negative effects because litter accumulation inhibits plant growth (Patton et al., 2007).

Our study revealed that ITV responded quickly to grazing abandonment and determined the short-term changes in forage quality. Strong species turnover effects may appear after long-term grazing abandonment when tall grasses with tissue of low quality outcompete high forage quality species (Odriozola et al., 2017). The use of mountain grasslands by domestic livestock has declined slowly in most of the Pyrenean-Basque-Cantabrian grassland systems (Ruiz et al., 2009), and significant management changes, including shorter duration in upland grasslands, lower stocking rate, and less shepherd control, have been identified (Aldezabal et al., 2015). Our study suggested that an active management system of livestock herds is necessary to prevent damage to grasslands associated with long periods of grazing abandonment. Short-term ITV-based changes in forage nutritional value are easily reversible; however, recovery times are much longer when changes are caused by species turnover effects, and high tissue quality species may be reduced in abundance or eliminated. This study also identified T. repens as a key indicator of forage quality and successional status in Iberian Atlantic grasslands.