When disturbances favour species adapted to stressful soils: grazing may benefit soil specialists in gypsum plant communities

Introduction

Herbivory by domestic and wild ungulates is one of the main drivers of global vegetation dynamics. Grazing mammals affect plant performance by biomass removal (Huntly, 1991) and, accordingly, plants have developed a wide array of adaptations to cope with grazing disturbance throughout evolution (Díaz et al., 2007). Grazing is usually considered a crucial biotic filter by restricting the range of trait values to those of species that survive and establish successfully (Violle et al., 2007). Grazing can also exert contrasting effects on plant community properties. In productive environments, grazing may alleviate plant-plant competition through changes in competitive hierarchies (Noy-Meir, Gutman & Kaplan, 1989; Louda, Keeler & Holt, 1990), biomass removal and trampling associated with grazing can also create spatial heterogeneity (Moret-Fernández et al., 2011), and thus allow species coexistence due to niche differentiation (Rosemond, Mulholland & Elwood, 1993). However, the selective removal of less grazing-tolerant species may also result in a reduction in species richness and diversity (Milchunas, Sala & Lauenroth, 1988), particularly when soil resources are scarce and plant productivity is low (Cingolani, Noy-Meir & Díaz, 2005).

Extreme soils, such as saline, limestone, serpentine or gypsum, have particular physical and chemical characteristics that restrict plant growth and species distribution (Rorison, 1960; Kazakou et al., 2008; Munns & Tester, 2008; Moore et al., 2014). Due to these constraints, atypical substrates are also major drivers of plant evolution (Hulshof & Spasojevic, 2020), leading to the development of specialised floras with numerous edaphic-endemics (Braun-Blanquet, 1932). Plants have developed different mechanisms to cope with the harsh conditions of extreme soils, and edaphic-endemics are usually soil specialists with substrate-specific strategies (Kruckeberg & Rabinowitz, 1985). These strategies allow them to optimise their performance and growth over other plants in their singular habitat (Cody, 1978), but may render them less competitive on non-extreme soils, which would explain why edaphic-endemics are frequent in plant communities associated with extreme soils (Rajakaruna, 2004).

The assemblage of plant communities on extreme soils has traditionally been explained in relation to plant adaptation to the limiting conditions of these special substrates (Caçador, Tibério & Cabral, 2007; Kazakou et al., 2008; Luzuriaga, González & Escudero, 2015; Luzuriaga et al., 2020). Plant communities associated with extreme soils are open shrublands or grasslands (Brady, Kruckeberg & Bradshaw, 2005; Mota et al., 2017), generally associated with large mammal herbivory and livestock grazing (Asner & Levick, 2012; Bakker et al., 2016). Consequently, the effect of grazing on plant community composition could be combined with soil restrictions in these systems. Evidence of the effect of grazing on vegetation in extreme soils is controversial. Ballesteros et al. (2013) and Pueyo et al. (2008) reported negative consequences on the abundance of two gypsum endemic species due to grazing, while other studies found that livestock grazing favoured edaphic-endemics over other plant species in extreme soils (Bonis et al., 2005; Beck et al., 2015). Most of these studies analysed the effect of grazing on individual species growing on extreme soils, but more studies at the community level are lacking.

Gypsum plant communities offer an excellent study system to evaluate the joint effect of extreme soils and grazing on plant community assemblage. The weathering of gypsum rock generates an unusual soil with high Ca and S content and low P availability that severely limits plant life due to nutrient imbalances (FAO, 1990; Casby-Horton, Herrero & Rolong, 2015). Gypsum soils occur worldwide in drylands, where extensive grazing almost always occurs (Pueyo et al., 2008; Akhani, 2015). The flora associated with gypsum is a unique endemic flora identified as an international conservation priority (Moore et al., 2014; Escudero et al., 2015), composed by edaphic-endemics, which are soil specialists (Palacio et al., 2007), and also species with wide ecological ranges (Meyer, 1980). Gypsum endemics show generally high affinity for gypsum soils (Luzuriaga, González & Escudero, 2015), also referred as gypsophily value (Mota, Sánchez-Gómez & Guirado, 2011; Musarella et al., 2018), and have foliar chemical composition characterized by high foliar S-accumulation (Duvigneaud & Denaeyer-De Smet, 1966; Salmerón-Sánchez et al., 2014; Merlo et al., 2019). Braun-Blanquet & de Bolòs (1957) suggested that plant communities rich in edaphic-endemics might be favoured by moderate grazing in gypsum soils, as grazing hardens soil conditions and impedes the weathering of gypsum, preventing the formation of a more organic soil that would favour other plant communities without gypsum endemics. However, no previous studies have evaluated the effect of livestock grazing on the assembly of species with different affinity for gypsum soils, or its relationship to the foliar composition of gypsum plants.

In this context, functional traits may play a key role in species composition on extreme soils with different intensities of herbivory. Plant functional traits that favour persistence under grazed conditions can be classified into avoidance and tolerance mechanisms (Briske & Richards, 1995). Avoidance mechanisms include traits that reduce plant accessibility and palatability, whereas tolerance traits lead to increased growth rate to compensate for biomass loss due to grazing. Leaf C and N content indicate species tolerance to grazing (Capó et al., 2021), where high N content is generally related to high growth rates (Pérez-Harguindeguy et al., 2013) and thus, to species that are able to compensate the biomass lost by herbivory (Grime, 2006). Moreover, leaf S content is a functional trait related to plant specialisation to gypsum soils (Merlo et al., 2019; Cera et al., 2021), although its ecological significance remains unknown (Palacio et al., 2007). High leaf S content could play a significant role in grazing-avoidance in gypsum environments (Palacio et al., 2014), because foliar S accumulation has been related to herbivore-deterrent compounds in Brassicales as glucosinlates (Ernst, 1990), and in some species of Acacia as crystals with S (He et al., 2015).

The aim of this study was to evaluate the extent to which the assembly of perennial plant communities under different grazing intensities and on high gypsum soils, is mediated by the affinity of species for gypsum, by species traits related to gypsum adaptation (leaf S concentration) or by traits related to plant tolerance to herbivory such as leaf C and N concentrations. If stressful conditions derived from atypical gypsum soils are the main selective force shaping plant communities on gypsum, we would expect gypsum-endemics to dominate over other species, independently of the grazing pressure. In this sense, the assembly of plant communities in gypsum environments would mainly depend on the gypsum affinity of plants (Luzuriaga, González & Escudero, 2015, Luzuriaga et al., 2020) and species with high S leaf content would be always dominant in plant assemblages. However, if grazing is the main factor driving plant assembly on gypsum plant communities, species with a rapid-growth strategy (i.e., species with high leaf N content and tolerance mechanisms, Grime, 2006) and/or with deterrent traits (avoidance mechanisms, Briske & Richards, 1995) would be dominant under high grazing pressure in relation to non-grazed conditions. This would lead to increased values of community weighted mean (CWM) for leaf N and decreased leaf N functional diversity in the most grazed plant communities. Finally, if as expected, herbivores have historically played an important role in shaping plant assembly in extreme soils of the Mediterranean region (Braun-Blanquet & de Bolòs, 1957; Montserrat-Martí & Gómez-García, 2019), edaphic-endemics should have developed both mechanisms to tolerate or avoid grazing and mechanisms to persist in restrictive soils. In this context, species with high affinity for gypsum and traits to cope with restrictions typical of gypsum soils (i.e., leaf S-accumulation), also fitted with traits to tolerate (i.e., high growth rates) or avoid herbivory would be favoured under grazing conditions, increasing CWMs values of leaf S and gypsophily value (GV) and decreasing their functional diversity in the most grazed conditions.

We also evaluated the contribution of intraspecific trait variability to cope with different levels of herbivory during the species assembly process. For this purpose, we applied the method proposed by Lepš et al. (2011) and de Bello et al. (2011) to disentangle the effects of interspecific vs intraspecific trait variability on species assembly. In this context, if species adapted to gypsum are also adapted to grazing disturbance, we would expect responses at the population level (i.e., intraspecific variability), since species would have enough intraspecific variability to cope with changes in grazing intensity. Otherwise, if species adapted to gypsum are not necessarily adapted to cope with grazing disturbance, we would expect those species that cannot withstand grazing would disappear when grazing intensity increases. This would lead to a shift in species composition (inter-specific variability) depending on their ability to cope with herbivory regardless of their soil affinity, so we would expect species assembly to be mainly due to species turnover.

Materials and Methods

Results

Variation in species composition

We found 52 perennial plant species (Table S2). Total plant cover and canopy height decreased from low to high grazing intensities (Fig. 1, Table 2). Plant composition was significantly different among grazing levels (F-ratio = 17.213, P-value = 0.001, TVE = 9.68%), as well as in response to soil features (F-ratio = 10.37, P-value = 0.001, TVE = 2.92%; Figure S2, Table S3). Specifically, Soil C (R² = 0.41) and gypsum content (R² = 0.32) explained the largest proportion of variability in plant species composition (Table S3).

Table 2:

Effect of grazing intensity on the main features of plant community properties and the gypsophily value after GLMMs.

Effect of grazing intensity on the main features of plant community properties and the gypsophily value after GLMMs. Chi-square values obtained by Wald tests based on generalised linear mixed models plus the family of error distributions and link functions assumed in the models are indicated. Id: identity link function; Log: logarithmic link function. Superscripts indicate GLMMs were run with the glmmTMB function to correct for dispersion of residuals.

	Family (link)	Df	Chis-square	Pr(>Chisq)
Vegetation cover (%)	Gamma (Log)1	2	41.889	0.001
Canopy height (cm)	Gamma (Id)1	2	76.541	0.001
Simpson Index	Binomial1	2	1.506	0.471
Gypsophily value	Gaussian	2	30.32	0.001

DOI: 10.7717/peerj.14222/table-2

Gypsum affinity and leaf traits (C, N, S-contents) at the community level

Gypsum specialists were favoured over non-specialists with increasing grazing (Fig. 2, see means and SE between locations in Figure S3). Our results showed higher CWM of gypsophily values (referred to GV) on the most heavily grazed plots than on the less grazed ones using a standard test (Table 3). However, these results may be slightly optimistic: when we applied a correction to avoid type I error (i.e., max test) the row-based test was statistically significant (P-value = 0.001), but the column-based test was not (P-value = 0.413). Plants with a rapid-growth strategy (high leaf N, low leaf C) and high leaf S were more likely to assemble in the most grazed plots, since CWM values of leaf C decreased and CWM values of leaf N and leaf S increased in higher grazing intensities (Fig. 2, Table 3). Medium grazing intensity was the most heavily filtered level on leaf C and leaf N, since plant communities showed the lowest FD values. Contrastingly, FD values of leaf S were statistically different (Table 3) and decreased with increasing grazing pressure (Fig. 2, P-adjusted = 0.06 between High and Low intensity). However, all plots along the gradient were heavily filtered, since FD values of leaf traits were generally lower than the Simpson Index (leaf C = 0.22 ± 0.01, leaf N = 0.25 ± 0.01, leaf S = 0.18 ± 0.01), which was equal along the gradient (Mean = 0.78 ± 0.01).

Table 3:

Results of generalised linear mixed models (GLMMs) for community weighted mean (CWM) and functional diversity indices.

Results of generalised linear mixed models (GLMMs) for community weighted mean (CWM) and functional diversity indices. Chi-square values obtained by Wald test based on GLMMs. The family of error distributions and link functions assumed in each model are also indicated. Id: identity link function; Log: logarithmic link function. FD: functional diversity values. Superscript indicate analysed with glmmTMB function to correct for dispersion of residuals.

	Community-weighted mean
	Family (link)	Df	F-ratio	Pr(>F)
Leaf C	Gaussian	2	18.23	0.001
Leaf N	Gaussian	2	67.34	0.001
Leaf S	Gaussian	2	31.81	0.001
	Functional diversity
	Family (link)	Df	F-ratio	Pr(>F)
Leaf C	Gaussian	2	22.06	0.001
Leaf N	Gaussian	2	31.995	0.001
Leaf S	Gamma (Id)	2	6.7121	0.035

DOI: 10.7717/peerj.14222/table-3

Overall, the process of species assembly in different grazing intensities was mainly determined by intraspecific shifts in functional traits (Fig. 3), since CWMs of leaf traits along the grazing intensity gradient were largely explained by ITV. Contrastingly, the amplitude of the range of trait values that occur in different grazing intensities was explained by species turnover, since shifts in FDs of studied leaf traits were explained mainly by interspecific trait variability (Fig. 3). Further, the contribution of interspecific and ITV on CWM and FD values were positively correlated for leaf N and negatively for leaf S. Leaf C showed a positive correlation for FD and negative for CWM values (Fig. 3).

Discussion

Our results provide evidence that livestock grazing may benefit soil specialist species during community assembly. Specifically, the relative abundance of gypsum specialists and of the leaf S content of species increased in species assemblages under medium and high grazing intensities. Although some studies found that certain soil specialists are more vulnerable to herbivory than their non-specialist relatives (Dechamps et al., 2008; Kay et al., 2011; Strauss & Boyd, 2011), our results highlight that under medium grazing pressure species with higher affinity for gypsum soils were favoured. Our results showed a clear trend of increasing relative abundance of gypsum specialists along the grazing gradient, although some caution should be applied to the statistical interpretation of this result. The sampling design should include more replicates of grazing level per site rather than using pseudoreplicates (i.e., plots), which could result in low β-diversity within sites, and increase type I error rates in the CWM approach (Zelený, 2018a). Further, the trend of increasing relative abundance of gypsum specialists was different between the three locations studied, although the highest values of GV were always reported in the medium and high grazing intensity site. Despite these concerns, our study is a remarkable result that aligns with previous studies of plant communities growing on extreme soils, such as white sands of tropical forests (Fine, Mesones & Coley, 2004), serpentine grasslands (Beck et al., 2015) and saline soils (Bonis et al., 2005). All these studies found that the occurrence of edaphic specialists was dependent on herbivory, most likely because herbivores modified competitive plant-plant interactions (Louda, Keeler & Holt, 1990; Grover & Holt, 1998) and may benefit the less competitive soil specialist species.

Our study showed that species prone to accumulate S in their leaves were favoured in medium and highly grazed conditions (i.e., larger CWM values for the leaf-S trait). Leaf S-accumulation is usually related to gypsum specialist species (Merlo et al., 2019), but the specific ecological role of S-accumulation in gypsum plants remains unknown (Palacio et al., 2007). It has been described that a high translocation of the excess element in soil to plant tissues (i.e., S in gypsum soils) could be a strategy to optimise plant growth in extreme soils by avoiding interference of that specific element with plant metabolism (Kabata-Pendias, 2010; Tran et al., 2020). Other studies proposed that high leaf S content could be related to a herbivore-deterrent strategy to avoid biomass loss in nutrient-limited habitats (Ernst, 1990; Palacio et al., 2014), as proposed for some species of Brassicales and Acacia (He et al., 2014; Tuominem et al., 2019). These results point at a selection of increased foliar S accumulation in plants growing on gypsum as a mechanism to deter herbivores (Boyd, 2007; Hoerger, Fones & Preston, 2013). These results are compatible with a potential role of herbivory as a selective force underlying the evolution of edaphic specialists (Fine et al., 2006; Lau et al., 2008), promoting the selection of foliar S accumulation in gypsum specialists. Therefore, it could be interpreted that foliar S accumulation could be an adaptive strategy of gypsum specialists to obtain better performance in gypsum environments disturbed by herbivores, where nutritional and water stresses are also present.

Some authors suggested that when resource availability is scarce, the costs of losing plant tissues due to herbivory are high, and plants that invest in chemical defences should be selected (Coley, Bryant & Chapin, 1985), mainly because the ability to compensate biomass losses is dependent on resource availability (Strauss et al., 1999). Thus, tolerance to herbivory is expected to be low in edaphically stressful substrates. Nevertheless, our results showed that species with higher leaf-N contents were selected in highly grazed compared to low grazed conditions (i.e., greater CWM of leaf-N and low leaf-C in highly grazed sites), suggesting that species with comparatively higher growth rates under the low resource availability of gypsum soils may have been favoured under grazing. These species could also be soil specialists. In the case of gypsum environments, soil specialists are expected to perform better in their atypical substrate than other soils (Cera et al., 2021). Also, they are more likely to assemble in gypsum soils than other soils (Luzuriaga, González & Escudero, 2015), as we observed high relative abundance of gypsum specialists along the gradient, regardless of grazing pressure. High foliar N concentrations could also be related to the accumulation of N-rich unpalatable compounds, more aligned with a grazing-avoidance strategy (Tuominem et al., 2019). However, such compounds are also frequently C-rich (Tuominem et al., 2019), which would have resulted in higher C concentrations, contrary to our observation of decreased foliar C in high grazing plots.

The effect of grazing on gypsum plant communities varied with grazing intensity. Low and medium grazed plots showed similar plant cover, but medium grazed plots displayed lower FD values (narrower range of trait values) than low and high grazed plots for leaf C and leaf N. Both leaf traits are linked to plant growth strategy (Grime et al., 1997; Pérez-Harguindeguy et al., 2016), indicating a narrow range of growth rate values of the species that assemble in medium grazed plots. These results can be explained by the effect of different grazing intensities on plant competiveness in the stressful conditions of gypsum soils. The higher FD in low grazed plots can be explained by the hypothesis of limiting similarity (Abrams, 1983). In these plots, there is high competition for resources due to the nutrient scarcity of gypsum soils (Boukhris & Lossaint, 1970). Species with different N requirements can coexist, because N acquisition niches do not overlap (Montesinos-Navarro et al., 2017), leading to higher FD values of traits related to growth strategy than medium grazed plots. Contrastingly, the top-down effect of sheep in high grazed plots seems to be due to disturbance associated to grazing that may reduce the biomass of dominant species (Noy-Meir, Gutman & Kaplan, 1989), alter the harsh physical crust typical of gypsum soils (Moret-Fernández et al., 2011) and eventually create new gaps for colonisation (Rosemond, Mulholland & Elwood, 1993). These processes may allow species with contrasting growth rates to coexist in heavily grazed conditions.

The relative contribution of species turnover and intraspecific trait variability on the species assembly process at different grazing intensities (Lepš et al., 2011) has been poorly analysed. Our results showed that intraspecific variability was the main component of species assembly for CWM values, while shifts in FD resulted from species turnover. The design of this study did not allow checking whether intraspecific variation was due to plastic responses of plants or heritable differences between genotypes (Bolnick et al., 2011). Furthermore, the design of our study limited the exploration of interspecific trait variation, as we analysed leaf traits in only a subset of 14 species growing along the gradient. The high relevance of ITV indicates that the species that finally conformed gypsum plant communities had enough functional variability among individuals to cope with different grazing intensities. Considering that livestock has been a usual anthropogenic activity since the Neolithic in the Iberian Peninsula (Balaguer et al., 2014), our results are compatible with the hypothesis that grazing has acted as an evolutionary driver over time, promoting a regional species pool fitted with successful strategies and enough functional variability among individuals to cope with herbivory. Consequently, in our study system, grazing did not act strictly as a biotic filter selecting certain species and jeopardizing others, but it acted in a subtler way, modifying the range of values of each plant trait selected during plant community conformation. This is a remarkable result that aligns with previous studies on the effect of herbivory on plant community assembly in environments with a long evolutionary history of grazing, such as Tibetan alpine meadows (Niu et al., 2016), and Inner Mongolia grasslands (Zheng et al., 2015).

Conclusions

To sum up, our results are compatible with the notion that grazing has likely been a powerful evolutionary driver in the conformation of plant assemblages on gypsum soils. Our results seem to indicate that herbivores may promote plant edaphic specialists in gypsum soils. Extensive livestock grazing should be considered as key tool to promote plant communities in gypsum soils where endemics persist, and should be evaluate with long-term studies, monitoring the intensity of livestock to avoid over-grazing, and linked other factors, such as climatic, that could affect the assemblage of communities.

Supplemental Information