Vegetation data were collected in the Crau nature reserve, hereafter called Crau area (Figure 1). This area is characterized by a Mediterranean climate (dry summers, mild winters, with rainfall peaking in autumn, some detailed temperature and precipitation data are available for the 1997–2006 period in Wolff et al., 2015), frequent strong northwesterly winds (Saatkamp et al., 2021) and is dominated by “pseudo-steppic” vegetation (Tatin et al., 2013). Dominant species composition is characteristic of Mediterranean Thero-Brachypodietea habitat “Pseudosteppe with grasses and annuals,” a high priority habitat for conservation (Corine Biotopes Code 34.512, Bissardon et al., 1997; Natura 2000 habitat type 6,220, San Miguel, 2008; EUNIS habitat E1.311 and E1.312, Louvel et al., 2013). This pseudo-steppic vegetation located in the fossil delta of the Durance river is characterized by a high stone cover (mainly pebbles) that limits the installation and the development of dense grasslands. Plant communities are highly diverse and dominated by low or very short herbaceous vegetation (49.4% of all plant species are annuals, Devaux et al., 1983; Helm et al., 2019; Pfeiffer et al., 2019), as a result of the combined effects of climate, soil, fire, and grazing. Grazing by itinerant sheep flocks has been reported since at least 2000 BP on the basis of Roman archeological remains (Badan et al., 1995; Henry et al., 2010; Robin et al., 2018; Saatkamp et al., 2021). This extensive grazing, occurring mostly from February to June, is heterogeneous throughout the steppe in the Crau area (Buisson and Dutoit, 2006; Saatkamp et al., 2017). These different parameters, Mediterranean climate, strong northwesterly winds, high stone cover, and two millennia of extensive grazing, all participate to shape this ecosystem.

We selected four study sites locally named in French “Coucou, Négreiron, Petit Carton, and Le Merle,” along a South-West/North-East gradient following the geomorphological-pedological-climatic gradient to cover all the environmental variability in the area (climate, topography, geology, pedology, and water accessibility) and in grazing pressures (Molinier and Tallon, 1950; Tatin et al., 2013; Table 1; Figure 1). An increase in fine elements (clay, silt) relative to decrease of coarse elements (sand, pebbles) has been shown to drive the presence of the biological crust micro-patterns in the south and their absence in the north of the plain (Rieux et al., 1977). Finally, the four sites are grazed by sheep kept in enclosures or led by a shepherd at an average grazing pressure of 400.8 ± 150 grazed day sheep/ha/year (see below, unpublished data from the Agricultural Chamber of the Bouches-du-Rhône department and the Le Merle Domain). Two out of the four sites had long-term, established ungrazed areas (exclosures): one of approximately 0.2 ha since 1982 at Le Merle and another of approximately 3 ha since 2001 at Négreiron (Bourrelly et al., 1983; Pilard and Brun, 2001).

Within the four sites described above, three micro-patterns in plant communities have been primarily reported: (i) Mediterranean grasslands, locally named “coussoul” or Asphodeletum fistulosi (Thero-Brachypodietalia) following the phytosociological typology by Molinier and Tallon (1950), hereafter named as lightly grazed reference communities, (ii) biological crusts, areas dominated by lichens on shallower soils, locally named “tonsure” or Crassuletum tillaeae (Tuberarietalia guttati) following the phytosociological typology by Molinier and Tallon (1950), corresponding to present or historically high grazing pressure areas (Rieux et al., 1977), and (iii) edge communities dominated by the tussock perennial grass [Brachypodium retusum (Pers.) P.Beauv.] of this grassland, locally named “grossier” and associated to the absence of grazing or very low grazing pressure areas, often localized at the edge of grazing plots (Dureau and Bonnefon, 1998). A total of 70 relevés, i.e., a list of each vascular plant species in a community for a defined spatial area, were carried out in the three Mediterranean vegetation micro-patterns: 20 in biological crusts (5 in exclosures), 30 in reference communities (10 in exclosures) and 20 in edge communities (Table 1). This difference in the number of relevés was due to the absence of long-term exclosures at two sites (i.e., no exclosure in sites Coucou and Petit Carton) or the absence of biological crusts in the north of the plain (i.e., in sites Le Merle and Petit Carton). Each relevé delimitation was defined taking into account the floristic, structural, and ecological homogeneity necessary to set up a phytosociological survey and on previous phytosociological studies already realized in the Crau area (Molinier and Tallon, 1950; Rieux et al., 1977). A minimum distance of 10 m between relevés located in a same site was observed to avoid any influence of each type of plant community on one another, and relevés were done at least 10 m away from fences (Buisson and Dutoit, 2006; Saatkamp et al., 2017). The floristic data were gathered from April to May 2018, which is the period of the year during which most plants can be identified to the species level in South-East France. All vascular plants were recorded in each relevé, and abundance data were estimated for each species according to the six abundance classes defined by Braun-Blanquet (+, 1, 2, 3, 4, 5; Braun-Blanquet, 1932). The sampled areas of the 70 relevés ranged from 4 to 40 m² (average area of grazed biological crusts 8.5 ± 7.9 m² vs. ungrazed biological crusts 7.2 ± 2.2 m²; of grazed reference communities 14.4 ± 9.2 m² vs. ungrazed reference communities 6.6 ± 2.8 m²; of edge communities 13.4 ± 7.8 m², sampled area of each relevé are detailed in Supplementary Table A1). The variability of the sampled areas arose from differences in the area occupied by the three different micro-patterns and to their species richness, as we stopped samples when 95% of the species of the plant community was incorporated into the relevé, ratio also based on previous phytosociological surveys realized in the Crau area (Molinier and Tallon, 1950; Rieux et al., 1977). To avoid potential bias, the area was included in all subsequent analyses.

For each relevé, we defined five independent indices that captured the main local environmental parameters: (i) an index of vegetation volume, based on visual estimates of vegetation cover, bare ground, and maximum height of the herbaceous stratum; (ii) a soil fertility index, based on 12 parameters related to soil chemistry (cation exchange capacity, pH water, and pH KCI, available phosphorus, potassium oxide, magnesium oxide, calcium oxide, sodium oxide, total nitrogen, organic carbon, and carbon-to-nitrogen ratio); (iii) an amount of fine soil particles, measured as the amount of clay and fine silt; (iv) an amount of coarse soil particles, i.e., fine sand; and (v) an amount of pebbles (>5 cm), based on visual cover estimates. These indices capture known gradients at the scale of the Crau area (Details in Supplementary Material Appendix I; Supplementary Figure A1) that broadly influence vegetation characteristics. All visual estimations have been collected by the same observer to reduce subjective variability. Grazing pressure was retrieved from inquiries of pastoralists realized by the Agricultural Chamber of the Bouches-du-Rhône department and the Merle Domain for the four study sites concerning the 55 grazed relevés out of the 70 relevés (15 were located in exclosures) in 2018. We calculated the number of sheep grazing days per ha and per year as Ns * Nd/ha for the grazed season of the year 2018 (autumn 2017 and spring 2018, before the transhumance to the mountain pastures), where Ns the number of sheep in a herd, Nd the number of days of sheep presence on the grazed area, and ha the area of the site (Table 1). It is important to note that we used two grazing variables. We used first a quantitative variable, described above and named hereafter “grazing pressure.” This variable was mostly used in statistical analyses. Second, we compiled a qualitative variable, grazed vs. ungrazed in exclosures, hereafter named “treatment” and used in some graphical representations where appropriate.

We calculated the species richness and Pielou’s measure of species evenness, hereafter named evenness, which is weakly sensitive to rare taxa, not sensitive to median or dominant taxa, and independent of species richness for sites with more than 20 species (Beisel et al., 2003; Gosselin, 2006). Evenness was calculated as J = H′/ln(S) where H′ is the Shannon diversity index based on the abundance of species in a relevé and S is the total number of species in a relevé.

To explore the possible relationships between the three micro-patterns and the considered response variables, i.e., amount of clay and fine silt, amount of fine sand, soil fertility, pebble cover, vegetation volume, species richness, and evenness, we fitted linear mixed-effects models using the restricted maximum likelihood algorithm as recommended in case of unbalanced design associated to random effect (Pinheiro and Bates, 2000; Bolker et al., 2009). We included three fixed-effect predictors: micro-patterns (qualitative, three levels), grazing pressure (quantitative), and the interaction between both variables, and included the site identity as a random effect on the intercept. We included the logarithm of the sampled area as a control, fixed-effect variable for the models that included species richness, evenness, vegetation volume, and pebble cover as response variables. For the model with pebble cover as response, we did not include grazing pressure and the interaction between grazing pressure and micro-pattern as explanatory variable. Equations are detailed in Supplementary Material Appendix II. We adjusted the p-values obtained from the linear mixed-effects models (Kuznetsova et al., 2017) and from the ANOVA applying the Holm correction (Sinclair et al., 2013). We checked visually using quantile-quantile plots that the residuals were consistent with the assumption of a normal distribution. We calculated the pairwise comparison Wilcoxon tests between micro-pattern group levels according to each response variables with corrections for multiple testing (Benjamini and Hochberg, 1995; Benjamini and Yekutieli, 2001) and between grazed and exclosures micro-patterns.

Multivariate analyses were used to study the relationship between the species composition of plant communities, environmental variables, and grazing pressure. In order to analyze plant community data, Braun-Blanquet scores were transformed into percentage cover on the basis of the average percentage of cover of each class: + = 1%; 1 = 5.5%; 2 = 17.5%; 3 = 37.5%; 4 = 62.5%; 5 = 87.5%. We performed Adonis and pairwise Adonis post-hoc tests with Bonferroni adjustments to check whether the plant communities were significantly different between and within groups defined by the micro-pattern, the treatment (grazed vs. exclosure), and the interaction of both variables. Moreover, we implemented a redundancy analysis (RDA) to investigate the variability in plant community composition at the micro-local scale in the Crau area constrained by centered environmental variables and grazing pressure. Note that we considered site and area as conditional variables, to account for the spatial nesting of sampled areas in sites and for the variable area sampled per relevés, respectively. As recommended by Borcard et al. (2011), we Hellinger-transformed the abundance data to reduce the wide disparity in magnitude between species abundances prior to the analysis. We retained all environmental variables (i.e., amount of clay and fine silt, amount of fine sand, soil fertility, and pebble cover) and grazing pressure as explanatory variables according to the Variance Inflation Factor testing collinearity between explanatory variables (all VIF < 5; Borcard et al., 2011). We finally performed a global permutation test to assess the significance of the RDA and a test of the canonical axes significance.

Finally, in order to evaluate and compare direct and indirect effects of grazing versus soil variables on species richness, evenness, and vegetation volume of the three micro-patterns communities, we performed a single multigroup structural equation modeling (SEM), a probabilistic model including multiple predictors and response variables in a single causal network. It provides a means of representing complex hypotheses about causal networks, accounting for factors that are both causes and effects, taking into account interaction effects, and testing for model and data consistency (Grace, 2006; Lefcheck, 2019). We proposed an initial conceptual model containing six regression analyses based on an a priori knowledge about the relations among the variables that was both consistent with our data and which made biological sense (Supplementary Figure A2). We considered species richness, vegetation volume, and evenness as endogenous variables, i.e., responses variables, and amount of clay and fine silt, amount of fine sand, soil fertility, pebble cover, and grazing pressure as exogenous variables, i.e., predictors. Note that an endogenous variable can also predict other variables in the model and be characterized simultaneously as endogenous and exogenous variables. This is the case for vegetation volume, species richness, and amount of clay and fine silt which were both endogenous and exogenous variables. A logarithm term of the area was considered as an exogenous variable influencing species richness, vegetation volume, and evenness to account for the variable area sampled per relevé (Supplementary Figure A2). As expected, this term had a direct significant and positive effect on species richness of biological crusts (Standardized estimate = 0.72, value of p < 0.001), reference communities (Standardized estimate = 0.3, value of p < 0.001), and edge communities (Standardized estimate = 0.6, value of p < 0.001), i.e., the larger the sampled area, the higher the species richness. It has not been represented in figures as this variable is a sampling artifact. In each regression analyses of the multigroup SEM, we included an independent random effect to account for the spatial nesting of sampled areas in site. We performed a single multigroup SEM testing relationships between variables for the three micro-patterns of the Mediterranean dry grassland of the Crau area. This single multigroup SEM allows testing for significant differences in the relationships between variables across and between micro-patterns. According to the test of directed separation, no relation was missing in the conceptual model. To evaluate the model and the goodness-of-fit with the data, we used a Fisher’s C test (high p-values indicate a good fit, with values well above 0.05 considered acceptable since the null hypothesis here tested is that the data predicted by the model and the observed data are consistent; Lefcheck, 2019).

All analyses and figures were done using R version 3.2.4 (R Development Core Team, 2008), and the packages vegan v2.5.6 (Oksanen et al., 2019), nlme v3.1.142 (Pinheiro et al., 2017), piecewiseSEM v2.1.0 (Lefcheck, 2016), ggplot2 v3.3.1 (Wickham, 2009).

We found differences between the three micro-patterns in terms of species composition. All groups constituted by grazed and ungrazed communities of the three micro-patterns of the Crau area were different from each other in terms of species composition (Pairwise Adonis test, adjusted value of p < 0.01; Supplementary Tables A2, A3) except grazed reference communities and grazed edge communities (Pairwise Adonis test, F value = 2.21, adjusted value of p = 0.41; Supplementary Tables A2, A3). We found a significant effect of the three micro-patterns (Adonis test, F value = 11.33, value of p < 0.001), of the treatment (Adonis test, F-value = 6.82, value of p < 0.001) and of the interaction of the both variables micro-patterns and the treatment (Adonis test, F-value = 0.38, value of p = 0.015) on the species composition of all groups.

We found differences between the three micro-patterns in terms of soil characteristics, vegetation volume, and species composition. Considering the RDA analysis, permutation tests supported the significance of the RDA (F-value = 4.03, value of p < 0.001; Supplementary Table A4): variance in plant community composition was explained to 20.55% by the explanatory variables and to 29.35% by the conditional variables (RDA adjusted R² = 16.4; Supplementary Table A4). In addition, the first and second canonical axes were significant (Permutation test, RDA 1 F-value = 12.16, value of p < 0.001; RDA 2 F-value = 7, value of p <0.001; Supplementary Table A4) and accounted together for 79.17% of the total explained variance. The first RDA axis captured differences in composition of micro-patterns and vegetation volume: it contrasted biological crusts, which had higher level of fine sand and soil fertility, whereas reference communities and edge communities had higher levels of clay, fine silt, and vegetation volume (Figure 2). The second RDA axis separated grazed plant communities from ungrazed plant communities, i.e., exclosures.

We found evidence of differences between the three micro-patterns in terms of species richness, evenness, vegetation volume, and soil characteristics. Linear mixed models showed significantly higher levels of fine sand for biological crusts compared to the other micro-patterns, significantly higher levels of clay and fine silt, of vegetation volume for reference communities and edge communities compared to biological crusts (Table 2; Figure 3; Supplementary Tables A5, A6). Thus, biological crusts were characterized by coarse silt and fine sand, whereas reference communities and edge communities were characterized by a mixture of variables element’ size: clay and fine silt and coarse sand. Edge communities were characterized by the lower amount of fine sand in soil than for other micro-patterns (i.e., the higher amount of coarse sand (Figure 3). Finally, linear mixed models did not show significant differences among micro-patterns concerning soil fertility and pebble cover. We found significantly higher levels of clay and fine silt for ungrazed communities compared to grazed communities, higher levels of soil fertility for grazed biological crusts, lower levels of species richness, and evenness for ungrazed reference communities (Supplementary Figure A3).

We investigated the direct and indirect effects of grazing pressure and soil variables on species richness, vegetation volume, and evenness of plant communities of the three micro-patterns and the fitting parameters of the adequate multigroup SEM were correct (Fisher’s C = 9.75, value of p = 0.78, ddl = 14, n = 70; Figure 4). This analysis showed that the direct effect of grazing varied across micro-patterns (significant interaction between grazing pressure and micro-patterns against species richness, value of p < 0.05), whereas effects of other predictors were similar across micro-patterns (value of p > 0.05; Figure 4). Effect sizes revealed three distinct functioning, each one specific of a micro-pattern. Soil variables and grazing had stronger direct and indirect effects on plant communities of biological crusts than on reference communities and edge communities. Direct effects of grazing on species richness differed among each micro-pattern, being negative in biological crusts and positive in reference and edge communities. Finally, soil fertility and pebble cover did not appear to influence plant communities of any micro-patterns according to the multigroup SEM.

This work highlights the contrasted functioning across micro-patterns, as grazing had a variable impact on species richness and soil granulometry. Micro-patterns with higher vegetation cover and volume are those that tend to be less grazed (although this has not been quantified here)—biological crusts are thus likely to be where grazing is the highest, followed by the reference communities and edge communities. Grazing exhibited a negative effect on species richness in biological crusts versus a positive one on richness in reference communities and an indirect and smaller positive one on richness of edge communities. This is consistent with the theoretical impact of grazing on species richness—intense grazing decreases diversity (in biological crusts) but lower levels favor diversity, which is what we observe in other micro-patterns (Milchunas et al., 1988; Tatin et al., 2013).

Differences among micro-patterns were particularly stark between biological crusts and the two other micro-patterns, as the former showed very different soil granulometry and vegetation volume. Biological crusts appeared primarily driven by soil variables, and second by grazing. They were more homogeneous in terms of species composition, and had an upper soil containing more fine particles (fine sand, coarse silt) than the other two micro-patterns. Such differences in soil characteristics suggest that biological crust communities may have their origin in present or former grazing practices (Saatkamp et al., 2021). Indeed, soil characteristics are influenced by grazing, suggesting that trampling could be the main driver of local soil differences, which compound with already-existing larger-scale differences, such as the combined geomorphological-pedological-climatic gradient of the Crau area and the spatial heterogeneity of the vegetation (Molinier and Tallon, 1950; Adler et al., 2001; Tatin et al., 2013). Soil trampling and compaction by animals alter the soil structure by compressive deformation: particles are brought closer together by the applied charge, which reduces the total pore space and permanently expels air or water from the soil pores (Patto et al., 1978; Silva et al., 2003). These alterations of the upper part of the soil affect soil physical properties, density, water dynamics, mineral element cycle, and plant growth (Lowery and Schuler, 1991; Schlotzhauer and Price, 1999; De Neve and Hofman, 2000; Drewry et al., 2008). Biological crusts could correspond to past high-impact areas of grazing, which are known to have a long-lasting effect, and to finely textured soils, which are more compressible than others (Larson et al., 1980), especially under stronger grazing pressures. Indeed, the higher amount of fine sand in soil composition could promote biological crusts by compaction, which may impair the installation of the dominant and key species in the Crau area Brachypodium retusum. In contrast, the heterogeneity of soil particle sizes in reference and edge communities seems to promote Brachypodium retusum establishment, as observed in field. Also, the increased abundance of this perennial graminoid species has been related to lightly grazed edge communities (Porensky et al., 2017). The installation of this key species facilitated by soil grain size characteristics and fire resilience could then result in changes in the microenvironment by improving water infiltration, reducing radiation to the soil, and accumulating seeds, fine soil particles, and nutrients between plant roots (Aguiar and Sala, 1999; Vidaller et al., 2019; Davies et al., 2022). This could also reduce herbivory through a phenomenon of facilitation between species (survival of seedlings protected by the dominant species, Aguiar and Sala, 1999). Taken altogether, these processes could explain the persistence of such contrasted community types at a local scale (<10 m).

Reference and edge communities appeared much more similar to each other than to biological crusts, as they had similar species compositions and abundances, species richness, evenness, and vegetation volume, all consequences of variable levels of herbivory (Génin et al., 2021). However, they differ in terms of soil composition (amount of fine sand). Edge communities appeared directly driven by soil variables and in a minor way by grazing, thus contrasting with reference communities which were directly influenced by grazing. Edge communities are often associated with abiotic modifications such as water infiltration (Masson et al., 2015) and not only to low grazing pressure, which does not affect plant community composition unlike high grazing pressure. One notable implication in terms of conservation and restoration is that it is not enough to remove grazing to produce a transition from reference to edge community in the short term because changes must also occur in soil granulometry (i.e., high amount of coarse sands in the soils of edge community). In fact, only long-term grazing removal affects will change the species composition of reference communities into that of edge communities (Saatkamp et al., 2017). Despite being more similar, switching from one reference to edge community (or the other way around) may thus require more than short-term changes in grazing pressure.

No effect of pebble cover has been found, although it is widely thought to impact positively species richness, vegetation volume, and evenness, by constituting a root shelter allowing the installation of Brachypodium retusum (Bourrelly et al., 1983) and other species, improving microenvironmental conditions and reducing herbivory by physically protecting the species from sheep’s mouths (Dureau and Bonnefon, 1998; Aguiar and Sala, 1999). This is probably due to the differences in sampled area among micro-pattern (more pebbles in the less grazed areas are also the largest sampled areas), which is not independent from the variable size and density of pebble related to the gradient of the alluvial deposits of the Durance river (Devaux et al., 1983), the largest pebbles at the highest recovery being encountered in the north-west where we did not sample biological crusts (Le Merle site; Figure 1), while the fine elements in the superficial part of the soil are more present in the south of the Crau area (Coucou site; Figure 1). Likewise, no effect of soil fertility has been found and no difference among each micro-pattern even though the former over-grazed areas could be suspected to be characterized by accumulation of nutrients and animal waste, beneficial to plant growth, long-term fertility, and known to cause a shift in community diversity and composition (Matches, 1992; Liu et al., 2016). These persistent legacies of historical sheep concentration, as Roman ones, on soil properties have persistent effects on today’s vegetation in this area (Saatkamp et al., 2021).

This work showed that ecosystem functioning varies over very small spatial scales: the three distinct communities at the micro-local scale are maintained by different direct and indirect effects of grazing and soil characteristics. Thus, there is no one-size-fits-all strategy to influencing the trajectory of the plant communities and one of the main levers of restoration strategies, changing grazing pressure, may have variable effects depending on the community type it is applied to. Grazing has a long-term residual effect on plant communities and soils, so it is not sufficient to remove it in the short term to restore a particular Mediterranean dry grassland community (Buisson et al., 2020). Indeed, grazing-induced changes in soil characteristics and fertility via long-term persistent eutrophication (soil phosphorus, calcium, and nitrogen) and trampling have effects on historical time-scales, as there is little opportunity to uncompact soil structure (Greenwood and McKenzie, 2001; Saatkamp et al., 2021). Grazing thus acts as a one-way lever: given the relative similarity between reference and edge community, it is likely that increased grazing pressure in edge communities would make them shift toward reference communities. However, the sole removal of grazing in reference communities is unlikely to drive them into edge-type communities on short time-scales, because of the need to reverse soil characteristics too. A similar and more pronounced effect is also likely in biological crusts: driving them toward other community types would require undoing the high soil compaction, which only occurs on very long time-scales (centuries) without drastic human intervention. Indeed, effect of past sheepfolds and intense grazing are still detectable on vegetation and soils, even after two millennia (Saatkamp et al., 2021).

In this situation, there is little opportunity to reverse such effects for restoration purposes. However, authors have already shown that a species of harvesting ant, Messor barbarus, could play an important role in the restoration of plant communities in the Crau area, by displacing a large quantity of plant seeds and changing physical and chemical soil characteristics toward that of the reference communities (Cramer et al., 2008; García-Ruiz and Lana-Renault, 2011; De Almeida et al., 2020). Typical dry grassland species are known to have low capacity to colonize new habitats or persist in the soil seed bank (Buisson et al., 2006; Metzner et al., 2017). Thus, increasing pebble cover and re-introducing Messor barbarus may provide an interesting and low-cost tool to reverse compaction and disperse seeds (Bulot et al., 2014; De Almeida et al., 2020).

In this work, we presented the first overview of the multiple effects of grazing and soil characteristics on plant communities of micro-patterns at the micro-local scale. Our results should enable conservation planners to use grazing and/or ecological engineer species as effective tools to influence plant community composition and soil structure. A more complete picture may emerge by measuring additional processes influencing the distribution of micro-patterns such as water infiltration capacity in case of compacted soil, wind-based redistribution of seeds, and fine soil elements that may accumulate in tall vegetation patches, and fire (Burke et al., 1998; Aguiar and Sala, 1999; Vidaller et al., 2019). Evaluating these additional effects would bring an understanding of the governing forces that shape community assembly in the context of Mediterranean dry grasslands conservation and restoration (Török et al., 2018).