Cushion shrubs encroach subhumid rangelands and form fertility islands along a grazing gradient in Patagonia

We assessed changes in soils and vegetation associated with different plant life forms and increasing grazing intensity (GI) in a subhumid grassland of Patagonia (Argentina), with vast portions of these grass steppes in degraded states. We reconstructed the historic GI gradient of the study sites and characterised the vegetation cover. We sampled the soil beneath patches of grass and shrub and their bare interspaces, as well as the height of the hummocks and the thickness of superficial horizons. We used multivariate analysis, inferential tests, simple regressions and the relative interaction index to measure the effect of historical grazing in these rangelands. The soil in interspaces between shrubs was the most degraded, with particle and fertility losses. While the soil of vegetation patches did not differ in any fraction, the soil beneath shrub patches was far more fertile. The soil of the sites with cushion shrubs developed the typical spatial heterogeneity of the fertility island effect, and their fertility decreased with increasing GI. With increasing GI, the relative cover of cushion shrubs grew and the total grass cover decreased, while the percentage of bare soil increased. The increasing grazing intensity favoured the transition of this subhumid grassland to shrubland. Grazing as an exogenous factor triggers processes of vegetation change and soil degradation, which lead to the encroachment by the cushion shrub Mulinum spinosum and fertility island development. This self-reinforced degradation process, well documented in arid and semiarid environments, also occurs in subhumid rangelands of the forest-steppe ecotone.


Introduction
Rangelands across the world are under increasing pressure. The drivers of change in grazing systems are the increasing demand for livestock products, land competition with agriculture and urbanisation, cultural and social drivers and climate change (Godde et al. 2018). Drylands comprise 41% of the earth´s emerged surface and are subjected to large modifications in soil and vegetation (MEA 2005). The increase in density, cover and biomass of native shrubby plants in grasslands (i.e., shrub encroachment sensu Van Auken 2009) is a global concern since it represents losses from a livestock production standpoint (Beeskow et al. 1995;D'Odorico et al. 2012). However, some studies highlight the benefits of the shrub encroachment process in ecosystem services, despite a decrease in forage provision (Eldridge and Soliveres 2014;Maestre et al. 2016). Most of the studies about shrub encroachment are focused on arid and semi-arid rangelands (Asner et al. 2004;Eldridge et al. 2012;Li et al. 2016), despite the importance of dry-subhumid rangelands in contributing to ecological goods and services (Lyseng et al. 2018).
High grazing intensities favour the establishment of shrubs and the replacement of a relatively homogeneous grass cover with shrub patches bordered by bare soil (D'Odorico et al. 2012;Gaitán et al. 2018). Grazing promotes floristic changes in rangelands, with a reduction in the cover of perennial grasses and replacement of mesophytic species with xerophytes (Cesa and Paruelo 2011;Eldridge et al. 2011). A high level of grazing disturbance frequently promotes increases in unwanted plants, such as shrubs, forbs, exotics and annual species (Anchorena and Cingolani 2002;Lyseng et al. 2018). Overall, grazing effects on vegetation vary with grazing livestock type (Tóth et al. 2018), climate, soil and use history (Díaz et al. 2007). Effects of grazing on plant diversity are contradictory, with studies reporting decreases (Herrero-Jáuregui and Oesterheld 2018; Paruelo et al. 2004), increases (Lyseng et al. 2018;Milchunas and Lauenroth 1993) and no changes (Howard et al. 2012;Maestre et al. 2016). Locally, although the problem of changes in vegetation in the wettest grasslands has required urgent discussion since the end of the twentieth century (Bertiller and Bisigato 1998), few studies have addressed this issue (Anchorena and Cingolani 2002;Lezama et al. 2014).
Grazing activities establish positive feedbacks with system heterogenisation, contributing to the development of fertility islands (Allington and Valone 2014;Ridolfi et al. 2008). The fertile island effect comprises the accumulation of fertile soil particles under plant canopies (Schlesinger et al. 1990). The soil beneath perennial vegetation becomes significantly more fertile than that in the surroundings (D'Odorico et al. 2013) and the spatial heterogeneity in soil resource distribution increases (Ding and Eldridge 2021;Mihoč et al. 2016). Because of the fertility island effect-also called fertility island, resource island, and island of fertility (Reynolds et al. 1999)the biotic processes become progressively confined to the soil beneath shrub canopies, resulting in augmented local fertility (Ridolfi et al. 2008). Nutrient cycling becomes restricted to sink spaces beneath canopies, and source interpatches become poorer (Reynolds et al. 1999;Videla et al. 2008). The effect is climate-dependent (Asner et al. 2004) and speciesspecific (Burke et al. 1998), but the research on fertility islands has been focused mostly on arid and semiarid environments and tall woody species (Allington and Valone 2014;Ochoa-Hueso et al. 2018;Ratajczak et al. 2012). Recent studies have emphasised cushion plants as island-driver agents (Zhang et al. 2011;Zhao and An 2021). Cushion-forming plants are chamaephytes or hemicryptophytes with singlestemmed growth form and hemispherical, subhemispherical, or flat-shaped canopy (Aubert et al. 2014). They are widespread in harsh environments and interact strongly with other species (Cavieres and Badano 2009). Their ground-contacting canopy has a positive influence on soil nutrient stocks (Chen et al. 2015;Zhao and An 2021) and high windblown-sediment trapping ability (Li et al. 2007;Zhang et al. 2011).
The formation of fertility islands is linked to accelerating aeolian and water soil erosion (Chartier and Rostagno 2006). This soil sediment redistribution reinforces the contrast between source interpatches of bare soil and sink patches of perennial vegetation (Parizek et al. 2002). This redistribution of particles leads to a decrease in plant-suitability and productivity of the source interpatches (Lal 2001;Wang et al.

3
Vol.: (0123456789) 2018). Size-differential particle redistribution magnifies the soil heterogeneity effect created by erosion (Navas et al. 2017;Romero Ovalle et al. 2021). Erosion effects are also soil type-dependent, with volcanic soils being an extreme example of the vulnerability to these processes (La . Soils in the Andean region of western Argentinian Patagonia originated mainly from volcanic ash (Irisarri et al. 1995). Because of restricted precipitation in the subhumid ecotone between the forest and the steppe, these soils are highly erodible and highly susceptible to losses of the vegetation cover (La .
Rangelands in the forest-steppe ecotone of western Argentinian Patagonia have experienced long-term problems resulting from management practices. Surveys and reports from the beginning of the twentieth century describe the high human and grazing pressure over these rangelands (Willis 1914). Extensive areas were burnt to clear woody vegetation and promote grass growth. Subsequent overgrazing resulted in species replacement and shrub encroachment (Rothkugel 1916). These rangelands have not been as extensively studied as arid and semiarid rangelands of Patagonia (Casalini et al. 2019;Castillo et al. 2020;Gonzalez et al. 2015;Oliva et al. 2020). Nevertheless, the available information about these subhumid systems describes an ongoing degradation process (Bertiller and Bisigato 1998;del Valle et al. 1998;Farias 2003;. Non-degraded states of Subandean steppes have high value as grazing ranges, characterised by the high homogeneity of perennial grass cover and little shrub cover (Paruelo et al. 2008). However, nowadays, vast portions of the ecotone are in a degraded state, showing high spatial heterogeneity with large denuded soil patches (Bertiller et al. 1995;Vogel and La Manna 2018). The most conspicuous components of perennial vegetation are cushion shrubs, with tussocks being less ubiquitous. Within this highly disturbed landscape, there are also native tall shrub scrubs, a fire-modified remnant of originally patchy woodlands (Anchorena and Cingolani 2002). Progressive and accelerating land degradation processes created a continuous loss of functional and structural attributes (La León and Aguiar 1985). This history of degradation may cause irreversible damage to the capacity of these rangelands to provide ecosystem services, including forage productivity (Rusch et al. 2017). It is necessary to deepen the knowledge about these valuable grasslands (León et al. 1998) and make efforts to comprehend how the degradation processes work in subhumid systems and how to stop and eventually reverse them.
This work aimed i) to identify changes in soil properties associated with different plant life forms and soil covers and with different historical grazing intensity, and ii) to identify changes in the vegetation composition and structure along a gradient of historical grazing intensity. We hypothesised that i) shrubs produce changes in the soil beneath them that lead to the development of fertility islands, ii) grazing activity strengthens the fertile island effect, and iii) livestock activity produces changes in vegetation cover and promotes shrub encroachment and cover of nonforage species. We predicted that the soil beneath shrub patches would exhibit greater fertility and more sediment accumulation than all other evaluated soils. Fertility under shrubs was expected to deviate even more strongly from surrounding interspaces than fertility under grasses. This heterogeneity in soil fertility was predicted to increase in response to higher grazing intensity. We also predicted that higher historical grazing intensity would be reflected in greater shrub, annual, exotic and forb cover and higher bare soil percentage, while it would reduce grass cover.

Study area
The study area is located in the Percy River basin, near the city of Esquel (Chubut Province, Argentina). The area belongs to the Subandean district of the Patagonian phytogeographic province. This district spans over 24,000 km 2 from 38º43´ S to 43º9´ S (Anchorena and Cingolani 2002). Climate is coldtemperate, with an annual precipitation range from 300 to 770 mm, concentrated in the colder months (April to September) and diminishing from west to east in a strong gradient determined by the orographic shadow of the Andes mountain range (Cesa and Paruelo 2011). Following this gradient, the moisture regime ranges from subhumid to semi-arid (Verón and Paruelo 2010). The semipermanent anticyclone of the Pacific Ocean determines the prevalence of strong winds from the west in spring and summer . Glacial frontal and side moraines configured a complex landscape, with hills alternated with lowlands (Andrada de Palomera 2002), determining an intricate vegetation pattern (del Valle et al. 1995). The Subandean district establishes, on its western subhumid boundary, a transition (approximately at 71°12´ W) with the Subantarctic forests, with incursions of forest species (León et al. 1998), configuring the trans-Andean forest-steppe ecotone of northern Patagonia (Kitzberger 2012). The soils, originating from volcanic ash (Irisarri et al. 1995), are highly erodible (La Manna et al. 2019). The loss of vegetation cover produces their desiccation and leads to a decrease in soil organic matter and an increase in soil erodibility (la . We defined our potential spectrum of study locations among the sheep ranches located in the ecotone with the typical layout of the region: extensive sheep grazing for wool production in a few large and heterogeneous paddocks with a seasonal continuous grazing regime (Borrelli and Oliva 2001;Peri et al. 2013). Since the registration of stocking rates is not a common practice in the region, we prioritised the availability of information about historical use. We selected the ranch Estancia Rio Percy (14,000 hectares, 730-1010 m above sea level, 42°51' S, 71°23' W, 650 mm mean annual precipitation), which has belonged to the same family since the early twentieth century. This enabled us, through interviews with its owners and manager, to access information that is not usually available at the ranch level. We were able to reconstruct the ranch's long-term grazing history in general terms and gained detailed knowledge about grazing loads from 1999 to 2016. We found out that the ranch followed the regional pattern of decreasing stocks, like many other ranches in this region of Patagonia (Golluscio et al. 1998;INTA 2003).
There, we identified the plant communities following Anchorena and Cingolani (2002), interspersed in the complex arrangement characterising the forest-steppe ecotone of Patagonia (Kitzberger 2012). We were able to recognise the communities "Stipa-Mulinum" (hereafter cushion shrub steppe), "F. pallescens" and "Stipa with Festuca" (hereafter grass steppe), and "shrub barren" (hereafter scrub) ( Table 1). All grass sites had grazing disturbance indicators such as cone-shaped tussocks, presence of M. spinosum, exotic plants, rocks on the surface. The scrub sites, with thorny woody species and non-preferred grasses (Anchorena and Cingolani 2002), were, in our opinion, the shrub community with probably less livestock-disturbed soils. Within each community, we searched for sites with slope < 9%, west aspect, and altitude ca. 800 m.a.s.l. to focus only on the average conditions of the area (Irisarri et al. 1995) and minimise the sources of variation. Following these decision rules, we finally delimited 12 sites within five paddocks distributed as follows: four cushion shrub steppe sites; five grass steppe sites (F. pallescens or Stipa with Festuca), and three scrub sites. The only grass steppe site that we were able to locate with no recent grazing was a small fenced paddock near Esquel (570 m.a.s.l., 42°55' S, 71°20' W, 635 mm mean annual precipitation). We decided to include it despite its lower elevation because it has remained ungrazed since 1961 and still maintains a high cover of tussocks. We calculated the Trabucco & Zomer Aridity Index (2019) Because of the absence of information on grazing distribution, we assumed it was randomly (homogeneous) distributed (Adler et al. 2001). Next, we calculated the normalised difference vegetation index (NDVI) values for each site for the same period with Landsat TM 5 and Landsat 8 OLI satellite images (30 m pixel) using QGIS software (QGIS Development Team 2009). With the NDVI values of the pixels comprised within each site, we estimated the annual aboveground net primary production (ANPP) in kilograms by hectare using the empirical relationship based on field data proposed by Paruelo et al. (Paruelo et al. 2004) for the Patagonian steppe as follows: Finally, we estimated the grazing intensity (GI) for each site by calculating the proportion of its ANPP (1) ANPP = 7951.5 * NDVI − 1088.4 needed to satisfy the DSE requirements (Milchunas and Lauenroth 1993).

Field survey
We followed the Environmental Monitoring of Arid and Semiarid Regions method (MARAS, for its name in Spanish) (Oliva et al. 2020). It is a monitoring system adapted to Patagonian conditions based on the Western Australian Rangeland Monitoring System (WARMS) (Watson et al. 2007). In each steppe site (4 cushion shrub and 5 grass sites), we placed 2 pointintercept and 1 gap-intercept transects (50 m in length each). In the point-intercept lines (250 equidistant points by transect), we determined floristic composition, total richness and soil cover. Patch-interpatch structure and soil surface conditions were measured with the gap-intercept lines (the third transect at each site). Using the point-intercept transects raw data, we calculated rangeland structure. Total vegetation cover was calculated (in percentage) as the subtraction from 250 non-vegetated point intercepts. The absolute cover of the plant life forms (shrub, grass, forb) was calculated by dividing by 250 the number of points intercepted by the life form of interest. The relative cover of different vegetation traits (cushion shrubs and annual plants) was calculated by dividing the number of points intercepted by each trait by the total number of points intercepted by any vegetation. We obtained the percentage of foliage superposition (strata) by dividing the number of double (or higher) plant interceptions at a single transect point by the number of total interceptions by vegetation. Total richness was the number of vascular plant species present at each transect. Finally, from the gapintercept transect, we calculated the mean bare interpatch length as the sum of centimetres of barren soil, divided by the total number of barren interpatches along the transect. Fieldwork was conducted in late spring and early summer of 2016 when this western Patagonian rangeland exhibits its maximum productivity (Paruelo et al. 2004).

Soil sampling and laboratory analysis
Based on previous literature about the redistribution of sediments (Funk et al. 2012;Li et al. 2009 we took samples of the top 5 cm of soil for analysis. Each sample was complemented with an undisturbed sample for bulk density determination. Each steppe site had 8 samples associated with it. In the grass sites, there were 4 samples beneath F. Pallescens or P. speciosa v. major patches (GRASS_ PATCH) and 4 samples in the barren soil spaces next to each of those patches (GRASS_BARE). Within each of the cushion shrub sites, we took 4 samples below patches of M. spinosum (SHRUB_ PATCH) and 4 from the associated barren soil interpatches (SHRUB_BARE). In order to collect samples of soil from ungrazed shrub situations (SHRUB_PATCH with GI = 0), we took 3 samples (0-5 cm) beneath tall shrubs (D. articulata, B. heterophylla and S. patagonicus) at each scrub site. To measure the variation in thickness of superficial soil horizons and the micro-topographic relief of vegetation hummocks, we cut a soil trench (Charley and West 1977) at each grazed site (GI > 0). Each soil trench spanned from the centre of the hummock beneath the vegetation patch (grass or cushion shrub) to the centre of the barren soil interpatch.

Fertility island effect quantification
We used the Relative Interaction Index (sensu Armas et al., 2004), which has recently been applied for assessing the fertile island effect through the heterogeneity in soil properties (Cai et al. 2020;Ding and Eldridge 2021;Zhao and An 2021). The numerical values for this index are between -1 and + 1 and symmetrical around zero (equal absolute value means equal intensity in the effect). Its use allowed us to quantify the difference in soil properties between a vegetation patch and bare soil (and also the inverse relation, between bare soil and a vegetation patch) and to compare this difference between different situations. Since soil properties do not change in the same direction responding to an equal force of a driver (some increase and some decrease), we considered only the absolute value of the index for assessing the magnitude of changes. Our purpose was to have a unique numerical expression of the fertility island effect through one condensed measure of the heterogeneity of 16 soil properties, for every soil sample. To do this, we calculated a one-by-sample average Relative Interaction Index (RII) as follows: where p 1-16 denotes each of the soil properties we measured and X a is the value of p measured beneath a given vegetation patch, and X b is the value taken by p on the bare soil interpatch paired with it. To calculate the RII of scrub site samples, since we did not collect bare soil samples from them, we estimated the term X b from the average value of each soil property p of all the bare soil samples obtained from the other sites in the same paddock.

Statistical analysis
We used a principal component analysis (PCA) (Hotelling 1933;Rao 1964) to make a descriptive analysis of soil variables. We included all soil samples in the PCA and identified them by plant life form (GRASS/SHRUB) and cover condition (PATCH/BARE), without pooling them by site. ( Following O'Connor (1988), we standardised all variables to unit variance. To compute the analysis, we used the package FactoMineR (Lê et al. 2008) for R software (R Core Team 2021). We plotted PCA results with the R package factoextra (Kassambara and Mundt 2020).
To test statistical differences of soil variables among categorical groups of samples, we used analysis of variance (ANOVA) and Tukey a posteriori test. Kruskal Wallis non-parametric test was utilised when the assumptions of ANOVA were not fulfilled. To account for the sampling design of our assessment, with multiple soil samples per site, we averaged them at site-level prior to all inferential tests. We proceeded in the same way to test differences in soil horizon thickness and mound height. We used InfoStat software for the analysis and for drawing the dot plots (di Rienzo et al. 2011).
The vegetation data and their variations with grazing intensity were analysed with distance-based redundancy analysis (db-RDA) (Anderson and Willis 2003;Legendre and Anderson 1999;). The db-RDA allows one to employ a user-defined matrix of dissimilarities, rather than relying on a pre-assumed model of relationships between the species and the predictive variables (Paliy and Shankar 2016). We employed the semi-metric Bray-Curtis measure to define our matrix since it is appropriate for use with ecological data of abundances (Mcardle et al. 2014). For computing the db-RDA, we used the R package Vegan (Oksanen et al. 2020).
We used regressions to analyse the response patterns of soil variables and vegetation structure along a gradient of increasing GI. We followed the advice of Kreyling et al. (Kreyling et al. 2018) about using gradient designs as the most efficient way to diminish unexplained variance in dealing with ecological responses to continuous environmental drivers like physical disturbance. To account for the sampling design of our assessment, with multiple soil samples and vegetation transects per site, we averaged them at site-level prior the regression analyses. To perform the regressions, we used the tidyverse (Wickham et al. 2019) and caret (Kuhn 2021) R packages and we plotted them with ggpubr R package (Kassambara 2020).

Grazing intensity
The estimated grazing intensities (GI) ranged from 0 to 2.89. Table 2 shows the GI, together with the estimated ANPP and the past stocking densities provided by the ranch owners and manager, which we used to calculate the GI.

Soil analysis
We grouped soil variables according to plant life form and soil cover condition for the PCA. Their acronyms, with their mean value and standard deviation, Table 2 Steppe sites assessed, ordered from lower to higher GI. The rank of ANPP is for 1999-2016. The stocking density (in DSE) is the mean value for the same period Acronyms are G for grass sites and CS for cushion shrub sites. Variables acronyms are ANPP for aboveground net primary productivity, DSE for dry sheep equivalent, and GI for grazing intensity are shown in Table 3. The biplot in Fig. 1 displays the first two PCA-Axes. There was a change in soil texture along axis 1, with increases in erodible particles from bare soil spaces to vegetated patches and a decrease in fine (clay and silt) fractions in the same direction. There was also a change in fertility-related properties (organic matter and nitrogen), showing a gradient along axis 2 of increasing fertility from bare spaces to vegetated patches (and from grasses to shrubs inside the latter group). The shrub patches were clearly differentiated from the grassy ones regarding these chemical properties. Shrub patches were enriched in organic matter and nitrogen with respect to their interpatches, while the grasses were not. The analysis exhibited a great contrast between different soil conditions related to shrubs, with covered patches and uncovered bare spaces being distant from each other, while grass-related samples were in proximity in the PCA plot. The soil heterogeneity measured through RII was higher for shrub-related soils. The detailed results of the PCA are shown in Online Resource 1. The inferential analysis of the above-mentioned variables (Fig. 2) showed that the interpatches between shrubs were enriched in fine fraction < 50 µm (clay + silt), while they had the highest losses of erodible particles (Fig. 2a), differing (p ≤ 0.05) from perennial grass patches in fine particle content and from the interpatches between grasses for the erodible fraction. Beneath patches, the erodible particle content was similar for shrubs and grasses. The fertility properties (Fig. 2b) yielded their lowest values for interpatches in shrub sites and their highest values beneath the shrub patches. We tested the variations in RII between cushion shrub and grass sites (Fig. 2c) and found the higher values of heterogeneity in soil properties for cushions (p ≤ 0.05). The complete results of inferential tests are shown in Online Resource 2. The response to the GI gradient of the soil fertility and heterogeneity (OM, N and RII), which we Table 3 Mean values and standard deviation (SD) of soil variables, grouped by vegetation type and cover condition Variables acronyms are: C_N (carbon/nitrogen ratio); CF (coarse fragments > 2000 μm); CLAY (clay content 0,1-2 μm); CS (coarse sand 500-1000 μm); EC (electrical conductivity); EP (effective porosity); ER (erodible fraction < 840 μm); FS (fine sand 100-250 μm); MS (medium sand 250-500 μm); N (total nitrogen); NAF (soil: NaF pH); OM (total organic matter); PH (soil pH in water); RII (Relative Interaction Index); SILT (silt content 2-50 μm); VCS (very coarse sand 1000-2000 μm); VFS (very fine sand 50-100 μm) analysed with linear regressions, (Fig. 3) showed that OM (Fig. 3a) and N (Fig. 3b) decreased significantly (p ≤ 0.05) with increasing GI for shrub-related samples. The regressions for RII (Fig. 3c) were not significant, but shrub sites showed higher RII than grasses along the entire GI gradient and fewer variations than grasses with increasing GI. Regression models are shown in Online Resource 3. The soil trenches (Fig. 4) showed that mounds formed beneath cushion shrubs were not  Figure 2a shows the percentual participation in the soil of erodible particles (ER %) and the pool of particles smaller than 50 µm composed of clay and silt (C + S %). Figure 2b shows the percentage of total Nitrogen (N %) and total Organic Matter (OM %). Figure 2c shows the Relative Interaction Index (RII). For Figs. 2a and 2b, the first letter of the categories in the x-axis shows plant life form (grass = G; shrub = S) and the second shows soil cover condition (bare soil = B; soil beneath vegetation patch = P). For Fig. 2c, the categories in the x-axis are grass (G) and cushion shrub (CS). Soil variable acronyms are shown in Table 3 significantly higher than mounds formed beneath grasses. Cushion shrub patches were larger (p ≤ 0.05) and accumulated more litter than grasses, but it was not a significant difference. Vegetation patches, be they cushion or grass, accreted far more sediments than bare spaces. The mean thickness of horizon A in grass sites was higher, but not significantly, than in cushion shrubs. While grass patches had the thicker A horizon, cushion shrub bare soil spaces had the thinner one. The complete results of the inferential tests are shown in Online Resource 2.

Vegetation analysis
The db-RDA analysis (Fig. 5) showed the growth in the relative cover of cushion shrubs as the main response to increases in grazing intensity (GI). The cushion shrub M. spinosum was the most relevant species related to these conditions, followed by Acaena splendens as a secondary species. The higher the GI levels, the lower the total covers of grasses and forbs, but these covers also showed differences between them. In the most grassy sites, the non-palatable Pappostipa speciosa var. major was the most relevant species, followed by the palatable Bromus stamineus, both being perennial. The higher forb cover was related to increases in total richness and in the superposition of vegetation layers. In this condition, the palatable grass Poa ligularis was the most relevant species. There was also a less relevant presence of non-palatable forbs and palatable grasses or grasslikes. The complete list of the species, with their acronyms, life form, cycle, status and palatability are listed, together with the complete db-RDA results, in Online Resource 4. The regression analysis (Fig. 6) showed three different kinds of responses of the variables to GI increments. The uncovered area (Fig. 6c) and bare interpatch size (Fig. 6h) increased considerably in a linear fashion (p ≤ 0.05), in line with the increase in cushion shrub relative cover (Fig. 6d), while grass cover (Fig. 6g) decreased in the same linear manner. The total cover of forbs (Fig. 6a) and the relative cover of annual grasses (Fig. 6b) exhibited similar behaviour. Their values increased until mid-intensities and experimented decay onwards. Their curves were similar to those of total richness of vascular plants (Fig. 6f) and the percentage of foliage superposition . The X-axis shows the distance (DIST) from the plant hummock centre in cm (Festuca pallescens and Mulinum spinosum, respectively). Y-axis shows in the positive scores the relative elevation (Height) in cm with respect to the lowest record in the trench. In the negative scores, the Y-axis shows the depth in cm (Depth) of the superficial layer of soil from the surface; Litter, Accretion and A horizon (Hz A); measured from the soil surface. Points show the mean values of all sites and bars show standard errors (Fig. 6e). Regression models are shown in Online Resource 3.

Discussion
To our knowledge, this is the first mention of M. spinosum forming fertility islands in subhumid rangelands, and even the first for any cushion shrub encroaching and forming fertility islands in subhumid and low altitude rangelands. The theoretical framework of shrub encroachment was first developed for arid and semi-arid systems (Asner et al. 2004;Schlesinger et al. 1990) and later extended to all drylands, including dry-subhumid regions (Briggs et al. 2005;Eldridge et al. 2011;Maestre et al. 2016). Despite this, most studies are still comprised between the original climatic boundaries (Howard et al. 2012;Li et al. 2016), which is also true for the region under assessment (Aguiar et al. 1996;Cesa and Paruelo 2011;Romero Ovalle et al. 2021). In addition, research on fertility island effect has been focused mostly on arid and semiarid environments too (Charley and West 1975; Ratajczak et al. 2012;Schlesinger et al. 1990). It has been also expanded to subhumid environments (Burke et al. 1998;Ding and Eldridge 2021) but like shrub encroachment, most studies are still conducted in dryer systems (Allington and Valone 2014;Ochoa-Hueso et al. 2018). In particular, cushion shrubs have been widely assessed, in terms of their role in the fertility island effect, in mountain regions well above 1000 m.a.s.l (Chen et al. 2015;Erfanzadeh et al. 2020;Zhao and An 2021) because of their great importance in high-altitude environments (Aubert et al. 2014). However, these studies have been carried out only for arid and semiarid regimes (Mihoč et al., 2016;Reid et al., 2010;Yang et al., 2017 but see also Liu and Guo, 2012).
The redistribution of soil and its spatial heterogenisation play a key role in the formation of fertility islands (Cheng et al. 2004). Increasing grazing causes a decrease in vegetation cover, particularly when sheep is the herbivore species (Tóth et al. 2018), which in turn produces the desiccation of soil superficial layers (Cavaleri and Lawren 2010). This affects the non-crystalline clays and reduces the potential of the volcanic soils to stabilise carbon (Hernández et al. 2012), leading to a decrease in OM and N contents (La ) and aggregate stability (Cambardella and Elliott 1993). The subsequent disaggregation resulted in the loss of the erodible particles from the topsoil, since light sands of volcanic soils are wind-erosion-prone (McDaniel et al. 2012;Nanzyo and Kanno 2018). This way the process of cover reduction initiated by grazing and trampling accelerated soil erosion and enlarged bare soil area, with subsequent decreases in plant cover. These positive feedbacks between the grazing-driven vegetation change and erosion reinforce the degradation and lead to desertification of the rangeland Webb et al. 2014). The particles eroded from interpaches are trapped by plant canopies with a success that varies between species (Zhang et al. 2011). We expected to find far more erodible material beneath cushion shrub patches than beneath grasses (Cai et al. 2020;Ding and Eldridge 2021), since the interpatches by the former were highly eroded. We hypothesise that our results did not confirm this expectation because the larger area of cushion patches is negatively affecting the proportional interface between them and the surrounding bare soil (Wu et al. 2000), making the source-sink ratio insufficient to prevent the leakage of these resources from the system (Kröpfl et al. 2011).
The spatial heterogeneity of the soil, in terms of its remobilization, was not due to a greater accumulation of particles beneath shrub canopies but rather due to an accentuated loss of erodible fraction from the shrub interspaces. This was not the case for soil fertility, because shrub interspaces were also degraded, but the soil beneath their canopies was more fertile than that beneath the grasses. This enrichment of soil beneath canopies is a fundamental biological process in fertility island formation (Allington and Valone 2014;Burke et al. 1998;Zhao and An 2021). We think that shrub patches were enriched in organic matter and nitrogen with respect to their interpatches because of the quality of litter inputs and their decomposition rate (Carrera and Bertiller 2010). M. spinosum accumulated more litter beneath its canopy than grasses, like other shrub species of various traits (Daryanto et al. 2012;Graff and Aguiar 2017;Liu and Guo 2012). Moreover, most of the shrub species in our study are deciduous and therefore produce litter of relatively high quality that mineralises fast (Cornwell et al. 2008;Satti et al. 2003). Compared to evergreens, they have a lower C:N ratio (Araujo and Austin 2015;Bertiller et al. 2005) and fewer secondary compounds in leaves (Carrera and Bertiller 2010;Saraví Cisneros et al. 2013). We suppose that these litter traits, with a positive effect on decomposition (de Paz et al. 2017), produced the fertility uprise and greatly influenced the heterogeneity in the soil properties that we assessed. This spatial heterogeneity of soil resources, produced also by other cushion shrub species (Chen et al. 2015), is one of the most reliable indicators of desertification (Cheng et al. 2004).
The increasing grazing intensity favoured the transition of the assessed subhumid grassland to shrubland and its encroaching by cushions shrubs, as in dryer environments (Komac et al. 2013;Kröpfl et al.  2011; León and Aguiar 1985;Ridolfi et al. 2008). It is interesting that M. spinsoum, cited as encroacher and fertile island driver for semiarid ranges (Gaitán 2009), was not even mentioned in the first botanical description of the F. pallescens steppe (León and Aguiar 1985). Shrub encroachment results from a combination of exogenic and endogenic factors (D'Odorico et al. 2012), but the driving force is high levels of herbivory (an exogenic factor), which reduces the competition from grasses (Milchunas and Lauenroth 1993;Ravi et al. 2010;van Auken 2009). The competitive advantage of grasses (Borrelli and Cibils 2005;Sala and Maestre 2014) is that they maintain their cover and prevent shrub establishment and growth even if the resource availability is suitable to sustain them (Aguiar et al. 1992;Stavi et al. 2009). The assessed reduction in grass cover with increasing grazing intensity is a well-documented effect for arid to subhumid rangelands (Eldridge et al. 2011). Changes in their cover are much more pronounced with long-term, heavy grazing (Asner et al. 2004). When the grazing intensity is enough to attenuate the competitive advantage, shrubs experience a physiological (affecting individual plants) and demographic release (affecting overall cover) with a greater effect in systems with higher precipitation than in semiarid regions (Ratajczak et al. 2012). The endogenic factors enhanced the encroachment with positive feedback processes including vegetation-microclimate feedback (Romero Ovalle et al. 2021) and soil erosion-vegetation feedback at patch scale (Cao et al. 2019), which can cause a transition to even more degraded states and make the process irreversible (Bertiller and Defossé 1993;van Auken 2000).
The physical properties of the soil in the interpatches changed due to the loss of erodible fraction and the increase in fine particles, which hinders the establishment of perennial plants (Romero Ovalle et al. 2021;Rostagno and Degorgue 2011). The low permeability of clay-enriched soils increases the potential for runoff (Rostagno and Degorgue 2011) and diminishes the water available for plant uptake (Paruelo and Golluscio 1993). These conditions were exploited by short-lived perennial forbs and annual grasses, both indicators of grazing disturbance (Anchorena and Cingolani 2002;Cesa and Paruelo 2011). In subhumid environments, the moderate removal of vegetation produces an enhanced intensity of plant interactions and a competitive release (Lezama et al. 2014;Milchunas and Lauenroth 1993;Segre et al. 2016). As predicted by the Milchunas-Sala-Lauenroth model (Milchunas et al. 1988), the richness in this subhumid area peaked at moderate grazing intensities because it encouraged species coexistence. Perennial native forbs took advantage of these situations, where the resistance of grasslands to shifts in species abundance and composition was reduced by grazing (Lyseng et al. 2018). In the places where the intensity was higher, shortly perennial forbs and annual species increased their cover, exploiting their tolerance to harsh conditions or their short life cycles to avoid them (Asner et al. 2004;Vetter 2005).

Conclusions
Increasing grazing intensity favours the transition of this subhumid grassland to shrubland. The changes that grazing produces in the vegetation and the soil lead to fertility island formation and encroachment by the cushion shrub Mulinum spinosum. These rangelands show high spatial heterogeneity, with large bare soil interspaces and nutrient cycling restricted to spaces beneath canopies. The feedback processes of vegetation change and soil degradation diminish the value of these steppes as grazing range, with negative ecological and economic implications. This study shows that a self-reinforced degradation process, well documented in arid and semiarid environments, also occurs in subhumid rangelands of the forest-steppe ecotone.