Cushion Shrubs Encroach Subhumid Rangelands: Fertility Islands in Fragile Soils Along a Grazing Gradient in Patagonia

Purpose This work assesses changes in soil and vegetation structure associated with grazing intensity (GI) in subhumid grasslands. We conducted the study in the Subandean district of Patagonia, Argentina. Non-degraded Subandean grass steppes have extremely erodible volcanic soils and are valuable grazing ranges. However, nowadays vast portions exhibit a heterogeneous cover that is mostly of cushions shrubs, with big eroded soil patches. Methods We selected four study sites along a GI gradient and one grazed-excluded site. Soils, vegetation cover and patches structure were characterised. We took soil samples beneath grass and shrub patches and their interpatches and in undisturbed spots. showed high heterogeneity associated with and Results also showed that and sand particles remobilised from bare to vegetated patches. Total nitrogen and organic matter increased in the same direction. Grass cover decreased as GI increases, while shrubs cover and total richness increased, until a collapse at the highest intensity. Relative cover of cushion shrubs and bare soil grow steadily with GI. Classication acronyms are P for perennial life cycle and A for annual; N for native status and E for Exotic; S for shrub life form, G for grass and grasslike, and F for forb. Species acronyms are: ACA.pin (Acaena pinnatida), ADE.cam (Adesmia campestris), ADE.obc (Adesmia obcordata), ANE.mul (Anemone multida), ARM.mar (Armeria maritima ssp andina), BER.emp (Berberis empetrifolia), BER.het (Berberis heterophylla), BRO.set (Bromus setifolius), BRO.sta (Bromus stamineus), CAL.sp (Calceolaria sp), CAM.den (Camissonia dentata), CAR.arg (Carex argentina), CAR.sub (Carex subantartica), CER.arv (Cerastium arvense), COL.hys (Colletia hystrix), ERO.cic (Erodium cicutarium), FES.pal (Festuca pallescens), GAL.apa (Galium aparine), GAL.ric (Galium richardianum), GAM.sel (Gamocarpha selliana), HOF.tri (Hoffmannseggia trifoliata), HOL.lan (Holcus lanatus), HOR.com (Hordeum comosum), HYP.inc (Hypochaeris incana), HYP.rad (Hypochaeris radicata), MUL.spi (Mulinum spinosum), MYO.sco (Myosotis scorpioides), NAR.chi (Nardophyllum chiliotrichioides), NAS.acu (Nassauvia


Introduction
Rangelands across the world are under increasing pressure by a growing demand for livestock products, land competition with agriculture and urbanisation, cultural and social drivers, and climate change (Godde et al. 2018). Large modi cations in soil and vegetation are expected, especially in drylands (Oliva et al. 2020). Shrub encroachment in grasslands is a global concern since this shift in density and cover of native shrubs represents losses from a livestock production standpoint (Beeskow et al. 1995 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 Paci c Ocean determines the prevalence of strong winds from the west in spring and summer (Anchorena and Cingolani 2002;Paruelo et al. 1998a).
Grass steppe is the dominant physiognomic units of the Subandean district, with Festuca pallescens (St.-Yves) Parodi predominating at higher altitudes and Pappostipa speciosa (Trin. & Rupr.) Romasch [= Stipa speciosa Trin. & Rupr. var. major (Speg.) Parodi] in lower altitudinal ranges (Del Valle et al. 1995). Glacial frontal and side moraines con gured a complex landscape, with hills alternated with lowlands (Andrada de Palomera 2002). Soils developed from volcanic ashes in a moisture dependant gradient, from Andisols to Mollisols (Irisarri et al. 1995;La Manna et al. 2020). Variations in landscape morphology and soils determine an intricate mosaic pattern where grass steppe coexists with scrub thickets and meadows developed in azonal environments (Del Valle et al. 1995). There are also incursions, throughout irregular boundaries, of species from subantarctic forests and occidental district from west and east respectively (León et al. 1998).
Volcanic ashes originated most soils in the Patagonian Andean region (Irisarri et al. 1995), with noncrystalline minerals giving these soils their distinctive properties (Dahlgren et al. 2004;McDaniel et al. 2012). In the subhumid ecotone between the forest and the steppe, precipitation and silica (Si) concentration quite restricts the formation of non-crystalline minerals (Par tt et al. 1984). These volcanic soils are highly erodible (La Manna et al. 2016) and have potential and actual erosion rates even higher than those soils in the more arid, eastern steppes (La Manna et al. 2019). It was suggested for these soils that desiccation, due to cover loss, affects non-crystalline materials, leading to soil organic matter decrease and an increase of soil erodibility (La ). Intense human pressure over the vulnerable soils of this ecotone may produce strong changes in soil and vegetation.
Non-degraded states of Subandean steppes have high value as grazing ranges, characterised by high homogeneity of perennial grass cover and little shrub cover (Bertiller and Defossé 1993;Bertiller et al. 1995; León et al. 1998; Paruelo and Golluscio 1993). Since the end of the XIX century, the entire area had intense pressure from human origin because of livestock grazing, rewood extraction and wild res (Aguiar et al. 1996;Willis 1914). Current vegetation in large portions of the area is visually dominated by invasive native cushion shrub Mulinum spinosum (Cav.) Pers. (León and Aguiar 1985;Roig 1998). This shrub covers almost all landscape positions except meadows and hilltops. Perennial grasses cover which complements the vegetation, consist of scarce, small and coned-shaped tussocks of F. pallescens, Poa ligularis Nees ex Steud and P. speciosa. Barren soil spaces are big and ubiquitous. There are little but relatively big and dense thickets formed mainly by Discaria articulata (Phil.) Miers, Colletia hystrix Clos, also with Berberis heterophylla Juss. Ex Poir. [=Berberis microphylla G. Forst.]. There are also sparse specimens of trees like Maytenus boaria Molina, Austrocedrus chilensis (D. Don) Pic. Serm. & Bizzarri and Nothofagus antartica (G. Forst.) Oerst.
We de ned our potential spectrum of study sites among sheep ranches with the typical productive con guration for this region (Irisarri et al. 1995). Since stocking rates registration is not a common practice in the region, we prioritised for our sites selection the availability of information about historical use. We selected a ranch, Estancia Rio Percy (14000 hectares, 730-1010 metres above sea level -masl-, 42°51' S, 71°23' W), which belongs to the same family since the early XX th century. This enabled us to get, through interviews with its owners and its manager, invaluable information orally shared from a generation to the next, and reconstruct the ranch's grazing history. Like many other ranches in this region of Patagonia, it followed the regional pattern of decreasing stocks (Golluscio et al. 1998;INTA 2003).
Four paddocks with different mean stocking densities were selected as study sites (PL, CA, AR, PM). We included a 5th site (AN) that was a small paddock (12 hectares) of grass steppe near Esquel city (570 masl, 42°55' S, 71°20' W) which remains ungrazed since 1961. This latter site, possibly disturbed by anthropic activities before its fencing and with probable border effects because of its size, represented the only steppe portion we could locate meeting the criteria of high total and tussocks cover and low shrub density characterising Subandean district non-degraded states ( (Oliva et al. 2020). Then, we estimated for each paddock the normalised difference vegetation index (NDVI) values for the same time lapse through Landsat TM 5 and Landsat 8 OLI satellite images using QGIS software (QGIS Development Team 2009). We calculated aboveground net primary production (ANPP) for the ve plots using the NDVI value for each pixel, through the method adjusted by Paruelo et al. (2004) for Patagonian steppes. We de ned grazing intensity (GI) by calculating the proportion of ANPP needed to satisfy herbivores dry matter requirements for each plot for each year (Milchunas and Lauenroth 1993). Finally, we assigned ve categories, ranging from less (non-grazed) to more intensity (A-E) based on the average GI for the two decades.

Field survey
We identi ed at each paddock sectors of grass-shrub and shrub-steppe. To better capture vegetation heterogeneity after more than a century of sheep livestock grazing in the area, we randomly selected two of them (one of each condition). Following the MARAS method (Oliva et al. 2020) we placed in each selected spot 3 transects of 50 m length (6 in total). We determined, by point-intercept lines, oristic composition, diversity and cover, and soil surface conditions in two transects. Patch-interpatch structure was measured by gap-intercept lines in the third transect. With transects raw data we calculated rangeland structure properties such as the absolute and relative cover of different vegetation traits, vegetation layer superposition, mean length of interpatches and richness (Oliva et al. 2020). Field work was in late spring and early summer when this western Patagonian rangeland exhibits its maximum productivity (Paruelo et al. 2004).

Laboratory analysis
We took 19 0-5 cm soil samples in each ranch's paddock. Sixteen samples were associated with vegetation transects. In the grass-shrub sector, 4 samples beneath F. Pallescens patches (PATCH_GRASS) and 4 more samples in the barren soil spaces next to each of those patches (BARE_GRASS). Within the shrub-steppe, we took 4 samples below shrub patches of M. spinosum (PATCH_SHRUB) and 4 from the barren soil interpatches associated (BARE_SHRUB). We also took 3 more 0-5 cm samples beneath tall shrubs canopies inside woody thickets (REFERENCE) as reference parameters. These spaces make part of the vegetation mosaic of the area and probably develop from the occurrence of early res or chopping in Austrocedrus forests (Anchorena and Cingolani 2002;Bran 2000). Since they present a very dense stem arrangement with many species having hardened spines which turn them into spaces di cult to reach, we assumed these spots as less livestock-disturbed. To measure the depth of soil horizons and their variations, and the micro-topographic relief of vegetation hummocks, we cut at each grazed paddock two soil trenches (one for grass and one for shrub) (Charley and West 1977). Each spanned from the centre of the hummock beneath the vegetation patch to the centre of the barren soil interpatch.
In the ungrazed paddock, we took 3 more 0-5 cm PATCH_GRASS soil samples beneath dominant tussocks of P. speciosa and 3 more BARE_GRASS samples from associated interpatches. We estimated they were su cient because of the size and vegetation homogeneity of this plot. No cushion shrubs were present in this area. For each soil sample, both in grazed and ungrazed paddocks, we took an undisturbed sample for bulk density determination.
At the laboratory, all soil samples were air-dried and sieved through 2 mm mesh to separate coarse fragments from ne soil particles (IRAM 1998). We determined particle sizes distribution through laser diffraction with previous organic matter elimination and sonication (Arriaga et al. 2006). Effective porosity was calculated by subtracting eld capacity to total porosity (Helalia 1993). The percentage of wind erodible fraction was determined according to Chepil (1953) technique. To estimate non-crystalline clays presence, we determined soil pH in sodium uoride solution (Fieldes and Perrot 1966 To perform an indirect and direct gradient analysis of vegetation data from the transects, we run a constrained or canonical correspondence analysis (CCA) (Ter Braak 1986). We assessed the relation of community structure and oristic changes with rangeland structure attributes and growing grazing intensity (Ter Braak 1987). We also performed a detrended correspondence analysis (DCA) (Hill and Gauch 1980)

Results
Grazing intensity classes Table 1 shows the grazing intensity (GI) classes according to average values of proportional ANPP required to satisfy the needs of the herbivore stocking density (SD) assigned at each plot. The assigned class for each paddock, from lower to higher GI, is expressed with letters A to E. It is interesting that for these study sites, GI values followed roughly SD increasing average levels.
Soil analysis Table 2 lists soil variables mean and standard deviation by vegetation and cover condition. The mean relative interaction index (RII) and its standard deviation by vegetation class are listed, together with Bonferroni a-posteriori mean comparison test.
Correlations of the 16 soil variables plus RII on the rst 3 dimensions of the PCA were analysed. The biplot in gure 1 displays the rst two dimensions. Dimension 1 (27% of variance) had strong positive correlations with silt content (SILT) (cos2=0.82), very coarse sand (VCS) (cos2=0.51) and clay content (CLAY) (cos2=0.50). It had strong negative correlations with medium sand (MS) (cos2=0.62) and ne sand (FS) (cos2=0.55). There were moderate positive correlations with soil pH in sodium uoride (NAF) (cos2=0.34) and relative interaction index (RII) (cos2=0.29). Dimension 1 had moderate negative correlations with erodible fraction composed by particles smaller than 0,84 mm (ER) (cos2=0.30) and effective porosity (EP) (cos2=0.23). Dimension 1 thus de ned a change in soil texture. There were losses of medium-sized particles (MS and FS) from vegetated patches to bare soil spaces, and gains of ner (CLAY and SILT) and coarser (VCS) fractions in the same direction. Dimension 1 also exhibited the contrast between different soil conditions related to shrubs, with covered and uncovered patches displayed distant from each other. Also, RII values followed the increments in ner and coarser particles.
Instead, grass-related soil samples were close to each other on this axis. Reference patches differed little with shrub and grasses patches or with grass related small bare soil spaces about soil textural variables.
Lawley-Hotelling and Hotelling tests showed signi cant differences (p < 0.05) between all groups, except grass-related ones (bare soil and patch).
Dimension 2 (17.3% of variance) was strongly correlated with total nitrogen (N) (cos2=0.79), and total organic matter (OM) (cos2=0.70). It had moderate positive correlations with electrical conductivity (EC) (cos2=0.48) and effective porosity (EP) (cos2=0.32). Negative correlations were with coarse fragment (CF) (cos2=0.24). Thus, chemical fertility-related properties (N and OM) de ned dimension 2, with an increasing fertility gradient from bare spaces to vegetated patches (and from grasses to shrubs inside this latter group). Unlike dimension 1, reference patches differed clearly from shrub and grass ones at this level. Also, there were increments in effective porosity (EP) from shrub to grass bare soil spaces, and from patches to reference spots. Dimension 2, therefore, differentiates through chemical properties those patches of current vegetation from the scrub/thicket reference spots. As in dimension 1, bare soil spaces associated with shrub patches showed the poorest soil condition. Lawley-Hotelling and Hotelling test showed signi cant differences (p < 0.05) between all groups, except grass-related ones (bare soil and patch).  The nest pool of soil particles (CLAY+SILT) had percentage values growing from vegetated spaces to bare ones. All grass-related spaces and shrub patches had similar values, without signi cant differences between them. Sand particles pool (FS+MS) showed lower values proportionally to higher values of ne particles, in an inverse gradient. The reference group took, for all particle categories, intermediate values between the most degradation-prone soil condition (BARE_SHRUB) and the others (Fig. 2a).
Mean values for chemical variables (OM and N) decreased from vegetated patches to bare spaces ones, and from grass to shrub inside each soil cover condition (bare and patch). The reference group had the highest levels for these two properties in all cases (Fig. 2b).
The trenches we opened showed the high soil heterogeneity existing in grazed areas, between soil mounds formed under canopies and barren interpatch areas. Barren spaces are depletion zones exposed directly to the biotic and abiotic environment. Super cial soil horizons were thicker beneath canopies, with a narrowing trend towards the barren space. Grassy areas had more homogeneity between patches and bare spaces than trenches opened in shrub area, where erosion created more contrast ( g. 3)
We found some major community composition changes with the indirect gradient analysis performed over CCA. Table 3 shows the best-represented species for each axis branch (negative and positive) together with its life form, cycle and status. The species gradient length for CCA axes was 3.09 standard deviations (sd) for the rst, and 2.32 sd for the second. The gradient showed an important community change, but also that some species remain common between sites. We expected a full turnover in species composition in about 4 sd, and a half-change in that composition in about 1 sd (Hill and Gauch 1980).
The rst axis showed, on its positive branch, a shrub-grass community with 57% of shrubs and 43% of grasses, all being perennial and native. On the negative side of the axis, half of the species were forbs (52%), with grasses (19%) and shrubs (29%) sharing the other half. Species at this end were 24% annual, while 33% were exotic species. Gradient length for this axis showed that if not a full turnover (sd below 4), there was a major change in composition throughout it. Forbs have a dramatic role in this change since they were absent at one edge and turned into the dominant life form at the other. Most of the forbs were perennial (71%) and native (65%). There was also a shift in the grass community, with 22% being annual and 33% exotic species. All annual grasses were exotic at this axis end.
The second axis had a gradient length of 2.32 sd. From positive to negative scores, forbs kept stable with a participation from 47 to 45%, grasses decreased from 41 to 36%, and shrubs raised from 12 to 18%. At the positive end of the gradient, annuals represented 18% of total species, while at the negative end they represented 27%. Similar to the rst axis, there was an important change in exotics from 12% to 55%. This change was best represented by forbs (25 to 80%), with a shift from native to exotic dominance in this group. The abundance of grasses decreased and exotics incremented dramatically within this group (going from 0 to 50%), being annuals half of them. All shrubs were native species.
Thus, the rst CCA axis determined a shift from a shrub-grass mixed community to a forbs-dominated one, with some losses in perennial and native species. The second axes showed decay in the participation of grasses and an increment in shrubs in the same direction, with a dramatic increase of forbs and annual grasses in the community.
Our model explained most of the variation since the rst residual axis (CA1) had a much smaller eigenvalue (0.07). Lawley-Hotelling and Hotelling test showed signi cant differences between GI categories. The extremes of grazing intensity range (A and E) were signi cantly different between them and also compared to middle intensities (B, C and D), which didn't differentiate between them (p > 0.05).
Shrub cover increased from less to more GI, with a nal decay in the most intensely-used class (E).
Superposition of vegetation layers and total richness followed a similar trend (Fig. 5a). Total grass cover followed an inverse trend compared to shrubs, decreasing with a higher GI. Bare soil and relative cover of cushions were coupled in their trends. Cushions behaved as increaser with increments in grazing intensity levels (Fig. 5b).

Discussion
Soil and fertility islands Through estimations of ANPP, a requirement to calculate GI, we could see the consequences of the grazing history in each paddock. Indeed, aerial biomass production calculated for our sites re ected fairly good the historical use contrasts we identi ed among plots with the qualitative survey made via interviews. We found that the non-grazed site produced more biomass than all grazed ones and that within them, the most heavily grazed site yielded far less than the others. Overgrazing affects ecosystem structure and functioning (López et al. 2011). ANPP has been proposed as a comprehensive indicator of system functioning through condensing information about energy ow and nutrient cycling (Gaitán et al. 2014). Though the assessment of ecosystem functionality was not an explicit aim of our work, our results concur with the observation of ecosystem misfunctioning at high grazing loads.
Our results evidenced the formation of fertile islands, with a differential expression between the cover categories evaluated. We could identify a signi cant land use effect on soil properties ( Fig. 1 and table 2). Also, it was related to indicators of degraded states in grasslands as shrub cover, large interpatches, and annual grasses (Fig. 4). The effect was linked to shrubs and their surrounding bare soil spaces. Medium sand and ne sand fractions were remobilised preferentially over other fractions. Important chemical variations occurred on organic carbon and total nitrogen, two of the most important variables for soil fertility characterization. Many rangelands across the world developed fertile islands, with changes in many soil properties and vegetation arrangement (Allington and Valone 2014; Ridol et al. 2008). We found that bare soil areas were bigger when associated with shrub patches, compared to those associated with grasses. This resulted in enlarging the bare soil vs. patch contrast and consequently in a bigger island effect detected for shrub encroached sites. The magnitude of the heterogeneity of the resources developed in cushion dominated patches is really important. This becomes clear when comparing our reference spots with the grazed areas (Table 2). Even being the long-term less disturbed conditions, we found that soil heterogeneity between the scrubs and their paddocks was lower (but not signi cant) than between cushions and their interpatches. We suppose the differential trapping ability together with the greater distance from redistributed particles sources contributed to maintaining less altered particles distribution of our reference spots compared to areas dominated by cushion shrubs and grasses. Our results showed that soil of grass and cushion shrub patches differ to a large degree from reference spots in the particles intervening in remobilisation (Fig. 2). The higher content of organic matter could also mediate this effect (Figs. 1 and 2). The organic matter reduces soil susceptibility to erosion by enhancing the stability of its aggregates, which turns into lesser particle mobilisation (Cambardella and Elliott 1993; La ). Thus, because of less particle entrapment and also to fewer losses, we assume that the soil of these spots exhibits a situation closer to the existent before intense disturbance from sheep grazing.
Our results showed preferential remobilisation of soil, with differential effects at deposition patches (sink) and source areas. Sand medium and ne particles were higher in the soil beneath canopies. In source patches (bare soil spaces), we found consequent increments of ne particles (clay and silt) and gravel. In Patagonia, these redistribution effects were found for climatic regimes ranging from arid to subhumid (Kröp et  Conversely, vertical stems of shrubs or trees from tribe Colletieae (Tortosa et al. 1996) enables sunlight to reach the basement. We propose that soil fertility under high shrubs re ects the light-enhanced degradation of the high-quality litter produced by these plants.
We didn't nd extensive indications of amorphous materials at our study sites. We expected to have these minerals, in the same line of previous assessments for soils with similar moisture regimes in this region 2008). Besides a grass cover contraction, we also found an increment in total shrub cover with increases in disturbance (grazing) intensity. The most relevant species (>1% total frequency) found in sites with high grass cover were Pappostipa speciosa (Pap.spe) and Bromus stamineus (Bro.sta). On the opposite edge, where grass cover reduced, Mulinum spinosum (MUL.spi) and Festuca pallescens (Fes.pal) were the most notorious species (Table 3 and g. 4). Increments in cushion shrubs were more associated with grazing intensity and bare soil than total shrubs, which shows their encroacher trait in disturbed sub-Andean rangelands (León and Aguiar 1985). Grazing is also a major driver of fertility island formation (Cai et al. 2020; Allington and Valone 2014). As discussed before, we found a greater fertility island effect associated with shrub patches and their paired large bare soil interpatches ( Fig. 1 and Table 2). Moreover, we also veri ed a determinant effect related to increments of shrub in community contribution and cover (Fig. 4).
We detected that larger bare interpatches accompanied the increase in total shrub cover ( g. 4). Lin et al.
(2010) proposed patch edge as a key indicator of water and sediments trapping e ciency. Larger patches like those formed by shrubs have less edge by area than smaller patches like those of grasses and forbs. Our survey results support this asseveration since grass patches were higher in sand remobilised particles, nitrogen and organic matter than those of shrubs (Fig. 2). Also, the bare interspace between shrubs was bigger (Fig. 4) and more heterogeneous with their corresponding patch than those of grasses ( Fig. 1, Table 2). These increased losses in shrub interspace areas together with soil amelioration under both grass and shrub canopies created the island effect we found.
We also found that the percentage of uncovered area incremented with growing disturbance intensity (Fig. 4). Our assessment detected structural changes, such as increments of forbs and annuals (in cover and proportional contribution to community) (Fig. 4). Literature mentions them as indicators of grazing disturbance for sub-Andean rangelands (Anchorena and Cingolani 2002; Cesa and Paruelo 2011; Paruelo and Golluscio 1993). These changes were both associated with grazing pressure increases but were not much related to each other. Community with forb increments had Poa ligularis (POA.lig), Acaena pinnati da (ACA.pin) and Carex argentina (CAR.arg) as the most relevant species (>1% total frequency). The most relevant in communities where annuals increased was Vulpia sp (VUL.sp) ( Table 3 and Fig. 4). When a patch collapse after heavy disturbance, the resources ow to bare spaces promoting recovery processes in rangelands (Kröp et al. 2013). Resources are lost if surrounding patches and interpatches cannot catch them e ciently. Following Lin et al. (2010) we assume that having less proportional patch edge, shrub patches trap resources depleted from bare soil with lesser e ciency than grasses. Thereby, these encroached ranges reached soil conditions limiting vegetation growth. In that situation, some species can take advantage of spatial heterogeneity created by disturbance.
At our study sites, annual species accompanied shrubs cover growth, an effect previously reported for moderate and heavily grazed rangelands (Noy-Meir et al. 1989). Annuals increment usually follows perennials decrease (Milchunas and Lauenroth 1993). These species prevail in bare interpatches because of advantages such as tolerance to harsh environmental conditions, higher seed production, shorter and prostrate canopy, secondary compounds and life cycles suitable to bene t from any benign microclimatic We veri ed clay increments in super cial soil of bare shrub interpatches (Fig. 2). This causes a xeri cation effect in Subandean highly eroded soils, diminishing water available for plants uptake (Paruelo and Golluscio 1993). The situation described here is comparable to the shrub-steppe of NW Chubut named state IV, with a very unlikely return to a more productive state because of irreversible soil changes (Paruelo and Golluscio 1993). In such harsh conditions, reinforced by a high content of coarse fragments (Table 2 and g. 3), we assumed as possible the constitution of niche gaps in the highly eroded big barren interspaces, where only annual grasses could complete an entire life cycle. Annuals are the most successful trait for such levels of resources depletion and spatial heterogeneity.
We found forbs gained participation in cover, being most of them perennial and native. We didn't observe the increment in exotics with grazing intensity reported for other subhumid rangelands (Lyseng et al. 2018). Many forbs share key traits with successful exotic invader plants, like grazing tolerance/avoiding.
Forbs are usually smaller than shrubs and grasses, which enables them to explore microhabitats created by disturbing activities such as herbivory, increasing plant richness (Jobbágy et al. 1996). Vegetation cover may become more fragmented because of destructive grazing and trampling (Lin et al. 2010), creating spatial heterogeneity and opening niche gaps, thus encouraging species coexistence (Gao and Carmel 2020). We believe that grazing intensity could lay behind the independent response of forbs to abiotic heterogeneity reported by Jobbágy et al. (1996). We observed vegetation layers superposition together with increasing forbs cover (Fig. 4) in our study sites, variables which are showing the mentioned niche heterogeneity. Indeed, in those sites with less bare soil and shrub cover, and smaller barren spaces, growth several life forms of perennial vegetation other than shrubs, like native forbs (Fig. 4). This situation pairs with grass and grass-shrub steppe of SW Chubut named as states II and III by Bertiller and Defossé (1993). Maintaining current grazing pressure will cause a probable transition to shrub-grass steppes of M. spinosum or Acaena spp. We consider that grazing management should focus on sites like these since they exhibit less soil heterogeneity and still sustain life forms other than shrubs and annual grasses.
Our results showed that increasing grazing intensity exerted a positive effect on total richness. This effect manifests along all grazing intensities evaluated with exception of the highest, were richness nally declines (Fig 5). While shrubs and grasses responded by replacement in species of one life form by the other, perennial forbs (and secondarily annual grasses and forbs) represented an addition to the community (Fig. 4). Golluscio and Mercau (1995) identi ed shrubs and cushions as increasers, some grasses and palatable forbs as decreasers and other grasses and forbs with maximum cover values at mid grazing intensities. Grazing induces changes in interspeci c competition, producing a shift in species abundance that otherwise won't succeed because of resistance of non-grazed grasslands to invasion (Lyseng et al. 2018). The effects of grazing on grasslands richness have been controversial, and attempts to generalise them had arrived at different conclusions (Eldridge et al

Conclusions
Our results indicates that historical livestock grazing activity caused a vegetation cover loss, which turned into organic matter decrease and a further increase of soil erodibility. The soil heterogeneity was in part created by soil resources remobilisation and magni ed by differential erosion. Indeed, it was stressed by greater losses in shrub interspace areas. Bare soil areas were larger and more eroded when associated with shrub patches. In contrast, shrub patches exhibited the biggest island effect. Thus, fertility islands were associated with indicators of degraded states in grasslands like high shrub cover, large interpatches, and annual grasses cover. This loss of resources turned into limiting conditions for vegetation growth.
Perennial grasses had a sustained cover decrease with a higher intensity of grazing. Grasses and shrubs reacted to grazing by replacement in species of one life form of species from the other. Higher grazing intensities incremented barren soil area and favoured the establishment of shrubs. Grazing provoked cushion shrub increases, even more than shrubs showing their encroacher trait at disturbed sub-Andean rangelands. Cushions increments with grazing intensity were paired with growth of bare soil area. Shrub total cover decay in the most intensely grazed class.
Increasing grazing intensity exerted a positive effect in vegetation richness of these subhumid rangelands. Increments of forbs and annuals were both associated with grazing pressure increases. Livestock grazing also favoured the increment of vegetation and soil resources heterogeneity. Moderate grazing intensities created niche gaps which favoured the coexistence of more plant species and resulted in perennial forbs species addition. Once surpassed an intensity threshold, this effect fades out. Our results con rmed early warnings about the drastic magnitude of changes in humid grasslands of Patagonia.
Regardless of the intensity, livestock grazing favoured shrub encroachment and a decrease in aerial biomass production. There was a loss in system services provision from a range management and conservation point of view. However, richness increases at moderate grazing intensities could represent an enhancement in system functionality. Bare soil increment was part of the shrub encroachment process, which resulted in challenging conditions for plant growth. However, shrub encroachment is aiding to keep within the system some resources depleted from bare spaces. Highest grazing intensities must be avoided, since causes shrub fertile island collapse and the eventual loss of these resources.     Diagrams of soil pro le corresponding to a cushion shrub patch (Mulinum spinosum) (Fig. 3a) and native perennial palatable grass patch (Festuca pallescens) (Fig. 3b). Dot plot showing mean values of selected cover and structure variables with high correlations on the rst two dimensions of CCA, against grazing intensity classes (GI) displayed in X-axis. Fig. 5a shows shrubs total cover (SHRUBS) together with total richness (RICH) and percentage of vegetation layers superposition (LAYERS). Fig. 5b shows grasses total cover (GRASSES) together with bare soil percentage (BARE) and cushion shrubs relative cover (CUSHIONS). Line bars represent standard errors

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. MANUSCRIPTWITHTABLESANDFIGURESEMBEDDED.pdf