This is the accepted version of the article: Urbina Barreto, Ifigenia; Grau Fernández, Oriol; Sardans i Galobart, Jordi; [et al.]. «Encroachment of shrubs into subalpine grasslands in the Pyrenees changes the plant-soil stoichiometry spectrum». Plant and Soil, 2020. DOI 10.1007/s11104-019-04420-3 This version is avaible at https://ddd.uab.cat/record/218420 under the terms of the license 1 Encroachment of shrubs into subalpine grasslands in the Pyrenees 2 changes the plant-soil stoichiometry spectrum 3 Ifigenia Urbina 1,2, Oriol Grau 1,2,3, Jordi Sardans 1,2, Josep M. Ninot 4 and Josep Peñuelas 1,2 1 4 CSIC, Global Ecology Unit, CREAF-CSIC-UAB, Cerdanyola del Vallès, Catalonia, Spain 2 5 CREAF, E08193 Cerdanyola del Vallès, Catalonia, Spain 3 6 Cirad, UMR EcoFoG (AgroParisTech, CNRS, Inra, Univ Antilles, Univ Guyane), Campus Agronomique, 7 97310 Kourou, French Guiana 4 8 Biodiversity Research Institute (IRBio) and Department of Evolutionary Biology, Ecology and 9 Environmental Sciences, University of Barcelona, Barcelona, Catalonia, Spain 10 11 Correspondence: 12 1. Ifigenia Urbina B - CSIC, Global Ecology Unit, CREAF-CSIC-UAB- Bellaterra, Spain. 13 ORCID ID: 0000-0002-8541-950X 14 Email: ifigeniaurbinab@gmail.com 15 16 17 Acknowledgments 18 We thank all the people who were involved in the field and laboratory work and Jose V. Roses-Díaz for his 19 comments that helped to improve this manuscript. We thank two anonymous reviewers for their contribution 20 to improve this manuscript. This research was supported by the Spanish Government grant [CGL2016- 21 79835], Catalan Government grant SGR 2017-1005, European Research Council Synergy grant ERC-SyG- 22 2013-610028, IMBALANCE-P, and the Spanish National Parks agency (OAPN) through the ARBALMONT 23 project [634S/2012]. 1 24 Keywords: shrub encroachment, plant-soil stoichiometry, nutrient stocks, subalpine grassland succession, 25 plant strategy. 26 ABSTRACT 27 Aims 28 Shrub encroachment has been reported over a large proportion of the subalpine grasslands across Europe and 29 is expected to have an important impact on the biogeochemical cycle of these ecosystems. We investigated 30 the stoichiometric changes in the plant-soil system along the succession (e.g. increase in encroachment from 31 unencroached grassland to mature shrubland) at two contrasting sites in the Pyrenees. 32 Methods 33 We analyzed the chemical composition (C, N,15N, P, K, Ca, Mg and Fe) in the soil and in the aboveground 34 plant compartments (leaves, leaf-litter and stems) of the main herbaceous species and shrubs at three 35 contrasting stages of the succession: unencroached grassland, young shrubland and mature shrubland. 36 Results 37 The plant-soil stoichiometry spectrum differed between the successional stages. Shrub encroachment 38 generally increased the concentration of C and Ca and the C:N ratio and often reduced to concentrations of N, 39 P and K in the leaves and leaf-litter, while several soil nutrient concentrations (N, P, K Ca and Mg) decreased. 40 The stocks of C, N, P, Ca, and Mg in the total aboveground biomass increased with encroachment. 41 Conclusions 42 Shrub encroachment favored the dominance of long-lived species with low concentrations of N and P in the 43 plant-soil compartments, high C:nutrient ratios in the aboveground biomass and increase the uptake of N 44 through ericoid or ectomycorrhizal fungi. We highlight the role of shrubs in the sequestration of C and 45 nutrients through the allocation to the aboveground biomass. The changes in plant-soil elemental composition 46 and stocks suggest a slowdown of the biogeochemical cycles in the subalpine mountain areas where shrub 47 encroachment occurred. 48 2 49 INTRODUCTION 50 The encroachment of shrub species into grasslands causes important changes in many grass-dominated 51 ecosystems at landscape and regional scales from low to high latitudes (Van Auken 2009, Myers-Smith et al. 52 2011, Naito and Cairns 2011, Komac et al. 2013, Formica et al. 2014). The expansion of shrubs is mostly 53 caused by changes in climatic conditions and land use (Eldridge et al. 2011). For example, the ongoing 54 increase in mean annual temperature and the thinning of permafrost at high latitudes are promoting the 55 expansion of shrubs across the Arctic and subarctic tundra (Tape et al. 2006; Hallinger et al. 2010; Myers- 56 Smith et al. 2011). Shrub expansion has also been observed in most mountain ranges and massifs across 57 Europe at lower latitudes, where humans have used the subalpine and alpine grasslands for traditional 58 activities such as extensive livestock herding (Roura-Pascual et al. 2005; Anthelme et al. 2007; Hallinger et 59 al. 2010; Targetti et al. 2010; Ferré et al. 2013) for many centuries or even millennia (Gassiot and Jiménez 60 2006; Pèlachs et al. 2007; Gassiot et al. 2016). The progressive abandonment of these practices in recent 61 decades has favored the encroachment of shrubs into subalpine and alpine grasslands (MacDonald et al. 2000; 62 Dullinger et al. 2003; Komac et al. 2013; Ameztegui et al. 2016). Climatic warming is also expected to 63 promote the expansion of woody species in subalpine and alpine regions in the Pyrenees (Grau et al. 2013; 64 Peñuelas et al. 2016; Angulo et al. 2019) and across Europe in the coming decades (Sanz-Elorza et al. 2003; 65 Körner and Paulsen 2004; Wookey et al. 2009; IPCC 2013). The abandonment of subalpine and alpine 66 grasslands is very apparent in the Pyrenees (northeastern Iberian Peninsula), where traditional extensive 67 livestock grazing and the frequency of intentional fires to create or maintain pasture have substantially 68 declined in the last century, mainly due to socio-economic changes in this area (Lasanta et al. 2000; 69 MacDonald et al. 2000; Serrano et al. 2000; Roura-Pascual et al. 2005; Jiménez and Pujol 2010; Barrio et al. 70 2013). 71 Shrub expansion is expected to cause several changes in ecosystem functioning, but few studies have 72 focused on the impacts of expansion on the functioning of subalpine or alpine grasslands in the Pyrenees 73 (Vitousek 1984; Montané et al. 2007, 2010; Barrio et al. 2013; Catalan et al. 2017; Grau et al. 2019) and 74 across Europe (MacDonald et al. 2000; Wookey et al. 2009). Some studies have reported that shrub 75 encroachment has increased the carbon (C) and nitrogen (N) concentrations in the soil and reduced soil pH 76 (Knapp et al. 2008; Eldridge et al. 2011), but such changes may strongly depend on shrub traits and 3 77 ecosystem features. Shrub and grass species have contrasting growth forms and differ in many functional 78 traits of their adaptive and reproductive capacities and their strategies of resource acquisition and allocation 79 (Chapin and Körner 1994). Different co-occurring shrub species may also have contrasting traits and 80 strategies (Illa et al. 2017), so the identity of shrubs is also crucial to our understanding of changes in 81 ecosystem functioning when encroachment occurs (Grau et al. 2019). Succession from grassland to shrubland 82 is thus expected to cause important changes in the distribution, cycling and stoichiometry of chemical 83 elements in the plant-soil system. Furthermore, the increase in woody biomass during the succession adds 84 complexity to the persistence and cycling of C and nutrients in the ecosystem. Wood is highly persistent and 85 stores C and nutrients such as N, phosphorus (P), potassium (K), calcium (Ca) and magnesium (Mg) (Chave 86 et al. 2009; Sardans and Peñuelas 2015). The nutrient stocks stored in wood, however, depend mainly on 87 wood density and may vary greatly among species (Sardans and Peñuelas 2013; 2015). 88 Most of the research on the impacts of shrub encroachment into grasslands has focused on soil C and 89 N balances (Hibbard et al. 2001; Jackson et al. 2002; Hood et al. 2003; Throop and Archer 2008; Van-Auken 90 2009), including studies conducted in the Pyrenees (Garcia-Pausas et al. 2007; Montané et al. 2007; Garcia- 91 Pausas 2010). However, changes in the chemical composition of both plants and soils, including key elements 92 such as P, K, Ca, Mg and iron (Fe), along the succession from grassland to shrubland, have not yet been 93 investigated in detail. We studied two common and contrasting landscapes where shrub encroachment 94 typically occurs in the Pyrenees, one facing north (hereafter North face) and one facing south (hereafter South 95 face). Each study site represents an independent case study because the grassland and shrubland species differ 96 between the two sites. In each of the two sites, we established a replicated sampling design in three 97 contrasting stages along the succession: a) the unencroached grassland (initial stage of succession, where 98 shrubs are not yet present), b) the mixed, young shrubland (intermediate stage of succession, where grassland 99 has small patches of shrubs) and c) the mature shrubland (advanced stage of succession, where the dominant 100 shrub species form large, monospecific patches). The main aim was to analyze the changes in the chemical 101 composition (concentrations and stocks) of the aboveground-biomass compartments (leaves and leaf-litter in 102 herbaceous species; leaves, leaf-litter and stems in shrubs) and the soil along the succession from grassland to 103 shrubland in the Central Pyrenees. 4 104 A previous study conducted in the same study sites in the Pyrenees highlighted that the 105 concentrations of some soil nutrients (e.g. N, P or K) generally declined along the succession from grassland 106 to shrubland (Grau et al. 2019), although there were differences between sites and shrub species. It was 107 hypothesized that shrubs may store nutrients in the biomass to take control of nutrients of the ecosystem to 108 avoid nutrient limitation and promote their further expansion along the succession, as hypothesized for 109 forested ecosystems with nutrient limitation in the soil (Grau et al. 2017). In this current study, we explore if 110 this hypothesis can be supported and investigate whether the decline in the concentrations of soil nutrients 111 along the succession was coupled with changes in concentrations and/or stocks of nutrients in the 112 aboveground biomass. Moreover, shrub expansion into grassland is expected to shape the plant-soil 113 stoichiometry spectrum and the biogeochemical cycle through changes in the abiotic and biotic conditions 114 along the succession from grassland to shrubland. We thus hypothesized that each successional stage would 115 show a contrasting plant-soil stoichiometric spectrum, resulting from changes in vegetation structure and 116 allocation strategies in each stage, and that C:nutrient ratios in the aboveground biomass and soil would 117 increase over the succession due to the lignification of plant compartments. 118 MATERIAL AND METHODS 119 Study area, species and sampling design 120 The two study sites were located in the Central Pyrenees (NE Iberian peninsula), one site on the North face 121 (Bargadèra, Val d‟Aran, 42°39′53.9″N, 00°50′07.8″E, with Cambro-Ordovician schists, sandstones and 122 quartzites) and one on the South face (Fogueruix, Pallars Sobirà, 42°35′32.9″N, 01°04′28.8″E, with a 123 Devonian slaty limestone), both on the periphery of the „Aigüestortes i Estany de Sant Maurici‟ National Park 124 (Figure SF1). The two sites are separated by only 22 km but are characterized by contrasting macro- and 125 microclimates; the site in Val d‟Aran has a strong Atlantic influence, whereas the site in Pallars Sobirà is 126 more continental. The mean annual number of days with precipitation or fog and the relative air humidity are 127 therefore higher at North face than South face (Catalan Meteorological Service, 128 http://meteo.cat/wpweb/climatologia, accessed on March 2018). The North face is a smooth north-facing 129 slope at an altitude of 1800-1900 m a.s.l. Since 2015, annual precipitation at the nearest meteorological 130 station to North face (Bagergue, 1400 m a.s.l.) has ranged from 800 to 900 mm, mean annual temperature has 5 131 ranged from 8 to 9 °C and minimum and maximum absolute temperatures were -14 and 30 °C, respectively. 132 The South face site is a steep south-facing slope at an altitude of 2000-2100 m a.s.l. Since 2015, annual 133 precipitation at the nearest meteorological station to South face (Planes de Son, 1540 m a.s.l.) has ranged from 134 600 to 800 mm, mean annual temperature has ranged from 8 to 10 °C and minimum and maximum absolute 135 temperatures were -15 and 30 °C, respectively. Both meteorological stations are at altitudes (ca. 400 m) lower 136 than the study sites, so precipitation is expected to be higher and temperatures are expected to be lower 137 (average and range) than the reported values at both sites. Both sites are in areas that have been intensively 138 used for livestock grazing for many centuries, but socio-economic changes in the Pyrenees since the 1950s 139 reduced the density of livestock, and some summer pastures have been abandoned (MacDonald et al. 2000; 140 Jiménez and Pujol 2010), which have favored the expansion of shrubs into these former grasslands. The two 141 sites differ greatly in their abiotic conditions (macro and microclimate, bedrock lithology, topography) and in 142 the functional characteristics of their vegetation (biotic conditions); the dominant shrubs are: 143 Calluna vulgaris (L.) Hull, Rhododendron ferrugineum L. and Vaccinium myrtillus L. at North face and 144 Arctostaphylos uva-ursi (L.) Spreng, Juniperus communis L. and Juniperus sabina L. at South face. 145 The grassland at North face is classified as a mesophilic, dense subalpine grassland „Nardion‟ (Galvánek and 146 Janák 2008), dominated mostly by grasses (Festuca eskia Ramond ex DC. and Festuca nigrescens Lam.) with 147 scarcer forbs (Trifolium alpinum L., and Cerastium arvense L.). The succession is mostly driven by the dwarf 148 shrubs Calluna vulgaris, Rhododendron ferrugineum and Vaccinium myrtillus. The grassland at South face is 149 classified as a xerophilic, open montane grassland „Xerobromion‟ (Carreras et al. 1983; EAE 2019), 150 dominated by grasses (Festuca ovina L. and Festuca gautieri (Hack.) K.Richt) that co-occur with small forbs 151 (such as Hieracium pilosella L., Achillea millefolium L. and Potentilla neumanniana Rchb.). The succession 152 from grassland to shrubland is mostly driven by the dwarf shrubs Arctostaphylos uva-ursi, 153 Juniperus communis and Juniperus sabina L. and to a much lesser extent by 154 Helianthemum nummularium (L.) Mill. or Thymus pulegioides L. 155 Soils at the North face site are built from Cambro-Ordovician schists, sandstone and quartzite, and, given the 156 moderate slope of the area, they are moderately deep, rich in organic matter, and acidic, corresponding in 157 general to Humudepts. Those at the South face site result mainly from Silurian pelite and also from Devonian 158 lime-rich slate. Being this site steeper, soils are shallower and irregular, moderately acidic and humic, broadly 6 159 corresponding to lithic Humudepts (Boixadera et al. 2014). Moreover, soils at the South face experience 160 summer drought episodes resulting into apparent drying off of the shallow soil level, which does not occur at 161 the North face. The soil pH in the North face tend be more acidic (pH around 4.7-5) than in South face (pH 162 5.7-6.3) (Grau et al. 2019). The percentage of soil organic matter measured in the same study sites were 163 higher at North face (herbaceous species: 21.2%, Calluna vulgaris: 37.8%, Rhododendron ferrugineum: 164 40.12% and Vaccinium myrtillus 40.6%) than South face (herbaceous species: 16.5%, Arctostaphylos uva- 165 ursi: 16.2%, Juniperus communis:13.2% and Juniperus sabina: 13.4%) (Grau et al. 2019). 166 We selected three stages at each site along the succession from grassland to shrubland: a) unencroached 167 grassland, which is still regularly grazed by domestic animals (Oriol Grau, pers. observation), is dominated by 168 herbaceous species and has no shrubs, b) young shrubland, composed of a mosaic of herbaceous species and 169 shrub patches of ca. 1-2 m2 and c) mature shrubland, where the dominant shrub species form large 170 monospecific patches (at least 10 × 10 m). The mixed shrub patches in the young shrubland covered 25-30 % 171 of the ground, whereas the shrubs in the in the mature shrubland covered 90-100% of the ground. 172 We selected the most common shrubs at each site (C. vulgaris, R. ferrugineum and V. myrtillus at North face 173 and A. uva-ursi, J. communis and J. sabina at South face; see suplementary Table ST1a, b for an ecological 174 and functional characterization of the shrub species studied). The herbaceous species and each of the shrubs 175 occurring along the succession in each site will hereafter be referred to as „vegetation types‟. 176 The replication of the successional stages was done in each site separately; the two sites were not treated as 177 replicates in the analyses but as two independent study cases. Sampling plots were established to reproduce 178 the extant structure of the vegetation in each successional stage. We established four replicate 2 × 2 m plots 179 (separated by a minimum of 10 m) in the unencroached grassland, and in each of three mature shrublands at 180 each site (Fig. 1 and Grau et al. 2019 for further details). In the young shrubland, though, the plots for each 181 vegetation type (grass, shrub 1, shrub 2, shrub 3) were grouped because all vegetation types co-occurred; in 182 this intermediate stage, we established four groups of four plots. In total we sampled 32 plots per site along 183 the succession from unencroached grassland to mature shrubland (see Fig. 1). The distance between the 184 successional stages or among the mature shrublands of each shrub species was >100 m and the distance 185 between the four replicates inside each successional stage was 10 m. The plots in the unencroached grassland 7 186 were placed in areas that represented the mixture of grass species that co-occurred in this successional stage, 187 whereas the plots in the mature shrubland were placed in large patches dominated by each of the shrub 188 species. 189 We collected plant and soil samples towards the end of the growing season (September 2015) for the analysis 190 of their elemental compositions. The litter layer was removed prior to the soil sampling; root samples could 191 not be collected so this belowground compartment could not be accounted for this study. We collected a total 192 of 144 shrub samples: 2 sites × 2 successional stages where shrubs were present (young shrubland and mature 193 shrubland) × 3 shrub species × 4 replicates × 3 plant compartments (leaves, leaf-litter -dead leaves still 194 attached to the plant- and stems). We collected a total of 32 samples of a mixture of herbaceous species in the 195 unencroached grassland and the young shrubland: 2 sites × 2 successional stages where herbaceous species 196 were present (unencroached grassland and young shrubland) × 4 replicates × 2 plant compartments (leaves 197 and leaf-litter). We also collected soil samples to a depth of 10 cm next to each plant sample with a 5-cm 198 diameter soil corer, 48 samples for shrubs (2 sites × 2 successional stages where shrubs were present × 3 199 shrubs species × 4 replicates) and 16 samples for herbaceous species (2 sites × 2 successional stages where 200 herbaceous species were present × 4 replicates), making a total of 64 soil samples. The concentrations of C, 201 N, P, K, Mg, Ca and Fe (see details of the chemical analyses in the next section) and 15N were measured in 202 leaves, leaf-litter, lignified stems and soils. The samples were transported to the laboratory in paper 203 envelopes, the soil samples were sieved (2 mm) and all samples were oven-dried at 60 °C for 48 h. They were 204 then ground with a ball mill (Retsch, model MM400. RestchGmbH. Haan, Germany), weighed with an 205 AB204 Mettler Toledo (Mettler Toledo, Barcelona, Spain) and analyzed in the chemistry laboratories at 206 Servei d‟Anàlisi Química, Autonomous University of Barcelona (Barcelona, Spain), where the percentages of 207 P, K, Mg, Ca and Fe in dry weight were determined, and at the University of California Davis Isotope Facility 208 (Davis, USA), where the isotopic compositions (δ15N) and percentages of C and N in dry weight were 209 determined. 210 We revisited the same sampling areas in September 2016 and collected more samples to characterize the 211 mean biomass of each plant compartment and the bulk density of the soil, which we used to estimate the 212 stocks of C and nutrients in the aboveground compartments (g m-2) of the vegetation types and the soil. We 213 assumed that the aboveground biomass of the herbaceous species was similar between the two years, because 8 214 their aboveground structures are entirely or nearly newly produced each year and the climatic conditions 215 (mean annual temperature and precipitation) did not vary much (Catalan Meteorological Service, 216 http://meteo.cat/wpweb/climatologia, accessed on March 2018). The changes in biomass of shrubs are very 217 limited in this ecosystem (Ninot et al. 2010a), so we do not expect large significant variations between two 218 consecutive years with similar climatic conditions. We collected the aboveground vegetation within 25 × 25 219 cm quadrats in the young shrubland and 50 × 50 cm quadrats in the unencroached grassland and mature 220 shrubland. The biomass was then transported in plastic bags to the laboratory and sorted manually into leaves 221 and woody stems (only for shrubs). We collected a total of 48 foliar samples and 48 stem samples from 222 shrubs (2 sites × 2 successional stages where shrubs were present × 3 shrub species × 4 replicates) and 16 223 foliar samples from herbaceous species (2 sites × 2 successional stages where herbaceous species were 224 present × 4 replicates). The plant material was sorted and oven-dried at 60 °C to a constant weight, and the 225 dry weight was measured. We collected soil samples with a bulk-density corer (9 cm diameter to a depth of 226 10 cm) to estimate the bulk density for each vegetation type. The quantity of soil (g cm-2) was calculated as 227 the product of the bulk density multiplied by the core depth. We also calculated the annual leaf-litter 228 production based on estimates of foliar persistence conducted in nearby locations (Ninot et al. unpublished 229 data). Leaf-litter productivity was used to estimate the annual input of nutrients from the leaf-litter for each 230 vegetation type. 231 Chemical analyses 232 For leaves, leaf-litter and stems, 4.5 mg of dry sample were weighed and encapsulated in tin capsules. For 233 soils, 8.6 mg of each sample were used. The samples were then analyzed for C and N with an Elementar 15 234 Cube system (Elementar Analyzen system GmbH, Hanau, Germany). N isotope was analyzed by an 235 Elementar Vario EL Cube or Micro Cube elemental analyzer connected to a PDZ Europa 20-20 isotope-ratio 236 mass spectrometer (Sercon Ltd., Cheshire, UK). The plant samples were combusted at 1000°C in a reactor 237 packed with chromium oxide and silvered copper oxide and soils were combusted at 1080 °C in a reactor 238 packed with copper oxide and tungsten (VI) oxide. The oxides were then removed in a reduction reactor 239 (reduced copper at 650 °C). The samples were interspersed during the analysis with several replicates of at 240 least two laboratory standards. These standards, which were selected for their compositional similarity to the 241 samples, had been previously calibrated against National Institute of Standards and Technology (NIST) 9 242 Standard Reference Materials (IAEA-N1, IAEA-N2, IAEA-N3, USGS-40 and USGS-41). Preliminary 243 isotopic ratios were measured relative to the reference gases analyzed for each sample. These preliminary 244 ratios were refined by correcting for the entire batch based on the known ratios of the laboratory standards. 245 The long-term standard deviation was 0.3‰ for 15N. 246 The concentrations of P, Mg, K, Ca and Fe in the leaves, leaf-litter, stems and soil were determined by 247 inductively coupled plasma mass spectrometry (ICP-MS) after digestion. For leaves, leaf-litter and stems, 248 0.25g of dry material were diluted in 5 ml of concentrated HNO3 and digested in a MARSXpress microwave 249 system (CEM GmbH, Kamp-Lintfort, Germany). The solution generated was analyzed by ICP-MS to 250 determine the elemental concentrations. The accuracy of the biomass digestions and analytical procedures 251 was assessed using certified biomass NIST 1573a (tomato leaf, NIST, Gaitherburg, USA) standards and 252 regularly analyzing blank solutions with no sample (5 mL of HNO3 and 2 mL of H2O2 with no sample). Dry 253 soil subsamples of 0.1 g were dissolved in an acidic mixture of HNO3, HCl and HF and digested as described 254 for the plant samples. The solutions were diluted with 1% HNO3 (v/v) before injection into the spectrometer. 255 Blank solutions (5 mL of HNO3 and 2 mL of H2O2 with no sample) were regularly analyzed. Total stocks of 256 the elements (g m-2) for each plant compartment and the soil were calculated as the biomass (leaves, leaf-litter 257 or stems) and soil weights multiplied by the concentration of each element analyzed. The C:N, C:P and N:P 258 ratios were calculated on a mass basis. 259 Statistical analyses 260 We conducted principal component analyses (PCAs) of the elements to visualize the overall differences in the 261 chemical composition of the plant and soil samples collected for each vegetation type (shrub and herbaceous 262 species) and successional stage for each of the two sites. This allowed us to reduce the number of variables of 263 leaf, leaf-litter and soil chemistry and integrate the information in a more comprehensive way on a two- 264 dimensional plot. We standardized the variables and performed a PCA for the elements of all compartments 265 (leaves, leaf-litter and soil) for each site separately to assess the differences in the chemical compositions 266 among the successional stages and vegetation types within each site. If two variables were highly positively 267 correlated (Pearson‟s product moment correlation coefficient > 0.6), we excluded one of the two in the PCA 268 to avoid the overfitting of variables. We also conducted PCAs of each plant compartment (leaves, leaf-litter, 10 269 and stems) and soil using data from both study sites to identify differences between vegetation types and sites. 270 To identify links between the functional traits of the shrubs and the foliar elements concentrations (C, N, P, K, 271 Mg, Ca and Fe), we conducted PCAs of the foliar elemental concentrations and some key functional 272 measurements, as seed weight, specific leaf area (SLA), wood density, biomass per unit area, height and 273 annual leaf-litter production. All PCAs were conducted with the „ggbiplot‟ (Wickham 2009) and 274 „FactoMineR‟ (Le et al. 2008) packages in R. The significance of the differences between vegetation types 275 and successional stages in each PCA were tested by permutational multivariate analyses of variance 276 (PERMANOVAs) with the R „vegan‟ package (Oksanen et al. 2017). β dispersion homogeneity tests were 277 conducted prior to PERMANOVA analyses. 278 We analysed the differences in elemental concentrations (dry weight (dw) of nutrient/dw of sample × 100) 279 and stocks (g m-2) in the plant and soil and between vegetation types (herbaceous species and shrub 1, shrub 280 2, and shrub 3) and successional stages at each site. Although in our sampling design we maximised the 281 independence of samples as much as possible, the extant structure of vegetation along the succession made it 282 not possible to find a completely independent, random distribution of combinations of successional stage × 283 vegetation type in the field sites. We therefore firstly checked if there were any underlying patterns of spatial 284 aggregation in the elemental concentration and stocks data with autocorrelation semivariograms „nlme‟ 285 package (Pinheiro et al. 2017) that are used for measuring the degree of spatial dependence between 286 observations as a function of distance. As in Grau et al. (2019) the initial visual interpretation of the 287 semivariograms indicated that in most cases the data were not autocorrelated despite the plots were not 288 completely randomized. Yet, we preferred to statistically check if autocorrelation should be accounted for in 289 the models. To do so we built two types of models: 1) Generalised least square (GLS) models that did not 290 account for spatial autocorrelation, 2) GLS models that included the „corSpatial‟ R function to account for 291 potential spatial autocorrelation. For this second type of GLS models, the coordinates (longitude and latitude) 292 of each plot were specified and we fitted five different models using a different autocorrelation structure each 293 time (e.g. exponential, linear, gaussian, rational quadratic or spherical). After running all the models (one 294 without and five with autocorrelation structure) for the element concentrations and stocks in each soil or plant 295 compartment, we compared and selected the best model using scores of the Akaike Information Criterion 296 (AIC). We used maximised log-likelihood in all cases and checked the normal distribution of the data before 11 297 run the models. The unencroached grassland was the initial stage of the succession from grassland to 298 shrubland, and we considered it as the reference level in all models. The two study sites were analyzed as two 299 independent study cases, so „site‟ was not included in the models as a factor. All analyses were performed 300 with R 3.2.4 (RCoreTeam 2017). 301 RESULTS 302 Significant differences in the elemental concentrations and ratios in the plant compartments (leaves and leaf- 303 litter and stems (only for shrubs) and soil for all the vegetation types along the encroachment succession in 304 both study sites are summarized in Table 1 and 2, the means  sd values in Table ST2. Significant differences 305 for element stocks are summarized in ST3, ST4, ST5 and ST6. 306 Patterns of plant-soil stoichiometry along the succession at North face 307 The plant-soil stoichiometry spectrum at North face varied along the two axes of the PCA (Fig. 2a). The first 308 PC axis explained 32.6% of the variance mostly integrated the changes along the succession, whereas the 309 second axis explained 23% and mostly showed the variability within each successional stage, particularly for 310 shrublands. Such variation between successional stages and also vegetation types (PERMANOVA tests 311 p<0.001 in both cases). Because the β-dispersion homogeneity test showed significant differences between 312 groups (successional stages or vegetation types), the significance of these PERMANOVA tests may thus 313 result from differences in the position of the centroids as well as from differences in the dispersion of the data 314 around the centroids. 315 Leaves and leaf-litter in patches of herbaceous species in the young shrubland had a higher C:N ratio and 316 lower [N], [K] and [Fe] than the herbaceous species in the unencroached grassland; the soil [P] was lower in 317 patches with herbaceous species in the young shrubland than in the unencroached grassland (Table 1a). 318 Generally, shrub leaves in the young and mature shrubland had higher [C], C:N ratio and lower [N], [K], [Fe] 319 and N:P than herbaceous species in the unencroached grassland (Table 1a,b); leaf-litter [C] was higher in 320 shrubs than in the unencroached grassland. Foliar [Ca] and [Mg] in the mature shrubland were higher than in 321 the unencroached grassland. Regarding the soil, all vegetation types in young shrubland had higher soil C:N 322 ratio and lower [P] and [Mg] than in unencroached grassland. (Tables 1a). Soil N was also lower in the young 12 323 shrubland of R. ferrugineum and V. myrtillus and soil K, Ca and Fe were lower in the mature shrubland of C. 324 vulgaris compared to the unencroached grassland (Table 1b). Foliar 15N and [N] were lower for all the 325 vegetation types in the young and the mature shrubland than herbaceous species in unencroached grassland 326 (Fig. 3). 327 We found stoichiometric differences within shrub species between the young and mature shrublands (Table 328 2). In the mature shrubland, C. vulgaris had higher [K] and C:N ratio and lower [Mg] in the leaf-litter, higher 329 C:N ratio and lower [N] and [K] in the stems. R. ferrugineum had higher foliar N:P, higher [N], [K] and lower 330 C:N ratio in leaf-litter, higher stem [N] and N:P and lower C:N and higher soil [Mg], N:P and C:P ratio in the 331 mature than the young shrubland. V. myrtillus had higher foliar [C] and [Fe], lower leaf-litter N:P and C:P 332 ratios and higher soil [Mg] in the mature than the young shrubland (Table 2). 333 The stocks of nutrients were generally much larger (generally 100-fold in all successional stages) in the soil 334 than in the biomass. The stocks of N (9.4 g m-2), P (0.71 g m-2), K (8.5 g m-2) and Fe (0.06 g m-2) in leaves of 335 herbaceous species in the unencroached grassland were significantly larger than in the leaves of all vegetation 336 types in the young and mature shrublands (Fig. 4 and Table ST3 for significant differences). Nevertheless, 337 shrub stems generally stored larger amounts of C and nutrients compared to leaves, so that the stocks of C and 338 nutrients in the aboveground biomass in the young and mature shrubland were generally much larger than in 339 the grassland (except for K). C. vulgaris and R. ferrugineum in the mature shrubland had higher stocks of soil 340 C (8953.5 g m-2, 7388.7 g m-2, respectively) and soil N (507.5 g m-2, only C. vulgaris) than the unencroached 341 grassland (C: 3645.7g m-2, N: 294.7g m-2). The stocks of soil P and K stocks in C. vulgaris in the mature 342 shrubland tend to increase but not significantly. Herbaceous species in unencroached grassland showed the 343 highest annual leaf-litter production (509.4 g m-2 year-1) (Fig. SF2, see Table ST3 for significant differences); 344 also, the annual leaf-litter production in V. myrtillus was higher in the mature (99.8 g m-2 year-1) than the 345 young shrubland (53.8 g m-2 year-1) (Fig. SF2, see Table ST4 for significant differences), implying higher 346 inputs of C (49.5 g m-2 year-1), N (1.1 g m-2 year-1), P (0.12 g m-2year-1), K (0.68 g m-2 year-1), Mg (0.2 g m-2 347 year-1), Ca (0.7 g m-2 year-1) and Fe (0.007g m-2 year-1) in the mature shrubland (Tables ST2 and ST4 for 348 significant differences). Leaf-litter production did not differ significantly for the other vegetation types. 349 Patterns of plant-soil stoichiometry along the succession at South face 13 350 The plant-soil stoichiometry spectrum at South face varied greatly along the two axes of the PCA based on 351 the chemical composition of leaves, leaf-litter and soil (Fig. 2b). The first PC axis explained 30.7% of the 352 variance and mostly integrated the differences along the succession, whereas the second PC axis explained 353 20.1%, which highlights the variability within the successional stages and vegetation types. We found 354 significant differences between successional stages (PERMANOVA p<0.001) and vegetation types 355 (PERMANOVA p<0.004) (Fig. 2b). In general, herbaceous species in the young shrubland had lower [N] and 356 [K] in leaves, higher [C], [Fe], N:P and C:N ratio and lower [N] and [P] in the leaf-litter and higher [Fe] and 357 lower [N], [Ca] and [Mg] in soil than the unencroached grassland (Tables 1a). Shrub leaves in the young 358 shrubland generally had higher [C] and C:N ratios and lower [N], [K] and N:P ratio than herbaceous species 359 in unencroached grassland, and shrub leaf-litter generally had higher [C] and lower [N] than herbaceous 360 species in unencroached grassland. The soil in young shrubland had lower [N], [Ca], [Mg] and N:P ratio than 361 soil under herbaceous species in unencroached grassland. Leaves in the mature shrubland often had higher 362 [Ca] and [Mg] and lower [N], [K] and N:P ratio, and leaf-litter with higher [C] and lower [Fe] than the 363 unencroached grassland. Soil [N], [Ca], [Mg] were always lower in the young shrubland than the 364 unencroached grassland (Tables 1); soil [P] was also lower in the mature shrubland of A. uva-ursi and J. 365 sabina. Foliar 15N and [N] were also lower for all vegetation types in the young and the mature shrubland 366 than herbaceous species in unencroached grassland (Fig. 3), as observed at North face. 367 We also found some stoichiometric differences within shrub species between the young and mature stages 368 (Table 2). A. uva-ursi had lower foliar [P] and higher C:P ratio and higher leaf-litter N:P ratio in the mature 369 than the young shrubland, J. communis had lower [Ca] and C:N ratio in the leaf-litter and lower [Ca] in the 370 stems in the mature than the young shrubland and J. sabina had higher foliar [C], leaf-litter [C], [N] and [P] 371 and lower [Ca], C:N and C:P ratios and higher stem [Mg] and lower C:N ratio in the mature than the young 372 shrubland (Table 2). Soil stoichiometry did not differ significantly between the young and mature shrublands. 373 The stocks of nutrients were generally much larger in the soil (100-fold in unencroached grassland and young 374 shrubland and 10-fold for mature shrubland) than in the biomass. The stocks of C, N, P, Mg and Ca in the 375 leaves of J. communis and J. sabina in the young and mature shrubland were higher than in the leaves of the 376 herbaceous species in the unencroached grassland (C: 191.1 g m-2, N: 8.2 g m-2, P: 0.56 g m-2, K: 6.7 g m-2, 377 Mg: 0.45g m-2, Ca: 1.7g m-2), and the leaves of A. uva-ursi had lower K stocks (Fig. 5 and Table ST5 for 14 378 significant differences). As observed at North face, shrub stems generally stored larger amounts of C and 379 nutrients compared to leaves, so that the overall stocks of C and nutrients in the aboveground biomass were 380 generally much larger than in the grassland. The soil for most shrubs had smaller stocks of Mg and Ca in the 381 young and mature shrublands than in the unencroached grassland (Mg 77.1 g m-2 and Ca: 431.1 g m-2). The 382 soil for A. uva-ursi (348.19 g m-2) and J. communis (315.76 g m-2) in the mature shrubland had lower N stocks 383 than in young shrubland. Soil P stocks tended to decrease along the succession, but not significantly (Fig. 5 384 and Table ST6 for significant differences). Annual leaf-litter production by J. sabina was higher in the mature 385 (938.9 g m-2 yr-1) than the young shrubland (566.2 g m-2 year-1) (Fig. SF2 and Table ST6 for significant 386 differences). Leaf-litter production did not differ significantly for the other vegetation types. 387 General patterns in the plant-soil stoichiometry spectrum at North and South face 388 The PCA of the functional traits and foliar elemental composition of the shrub species from both study sites 389 (Fig. SF3) showed that the shrub species at North face differed from those at South face (PERMANOVA 390 p<0.01) along the PC2 (31.2% of variance explained). Moreover, differences between shrub species within 391 each site were also discriminated (PERMANOVA p<0.01) along the PC1 (39.1% of variance explained). 392 The PCA based on the foliar stoichiometry including both sites showed significant differences between shrubs 393 and herbaceous species (PERMANOVA p<0.5), mostly along the first axis (39,4% of variance explained, 394 SF4a). The chemical composition of leaf-litter showed a similar pattern (PERMANOVA p<0.001) but 395 varying mostly along the second axis (24% of variance explained, Fig. SF4b). The chemical composition of 396 the soil, however, differed more between the two sites (PERMANOVA p<0.001) (41.4% of variance 397 explained along the first axis, SF4c) than between herbaceous species and shrubs. In this case, though, 398 because the β-dispersion homogeneity test in Fig. SF4c was significant, the significance of the 399 PERMANOVA test may result not only from differences in the position of the centroids but also from 400 differences in the dispersion of the data around the centroids. Finally, the chemical composition of the stems 401 differed mostly between sites (PERMANOVA p<0.001) (39.3% of the variance explained along the first axis, 402 Fig. SF4d). 403 404 DISCUSSION 15 405 The plant-soil stoichiometry spectrum differed greatly along the succession from unencroached 406 grassland to mature shrubland (Fig. 2), as initially hypothesized. The shift in plant-soil stoichiometry 407 spectrum was more apparent and consistent between the unencroached grassland and the young shrubland 408 than between the young and mature shrubland (Tables 1 and 2), suggesting that changes in growth forms 409 (from herbaceous-dominated to shrub-dominated) play a greater role in shaping the plant-soil stoichiometry 410 spectrum than differences in successional stage as such. In fact, the plant-soil stoichiometry spectrum is 411 expected to be primarily shaped by changes in abiotic and/or biotic factors, such as changes in vegetation 412 along the succession. This is evidenced by the shift between herbaceous species and shrubs, from high foliar 413 [N] and [K] for herbaceous species to high foliar [C], [Ca], [Mg], C:N and C:P ratios and low [N], [P] or [K] 414 for shrubs (Table 1, Table ST2); or from low to high [C], C:N and C:P ratios in litter, or from high to low [N], 415 [P], [Mg] and [Ca] in soil. Such changes in stoichiometry may have major consequences on the functioning of 416 ecosystems (Eldridge et al. 2011). For example, the high [C] and C:N and C:P ratios in the shrub leaves and 417 leaf-litter promote the formation of recalcitrant organic matter accumulation with slow decomposition rates in 418 the top soil of shrubland compared to the grassland (Ninot et al. 2010b; Garcia-Pausas et al. 2017). 419 The changes in foliar chemical spectrum enabled us to detect shifts in the biogeochemical niche 420 (Urbina et al. 2017; Peñuelas et al. 2019) (Fig. SF4). We found that shrubs and herbaceous species differed 421 significantly in their biogeochemical niche, as expected from these two contrasting growth forms. For 422 example, herbaceous species in the unencroached grassland had higher foliar nutrient concentrations (N, K), 423 characteristic of plants with fast growth rates; whereas shrubs had higher C concentration and C:N and C:P 424 ratios in leaves, characteristic of plants with slow growth rates (Ågren 2004; Sardans et al. 2012; 425 Zechmeister-Boltenstern et al. 2015). 426 The changes in stoichiometry in the plant aboveground compartments along the succession were 427 coupled with a decrease of the soil concentrations of some essential nutrients needed for plant development. 428 Grau et al. (2019) suggested that this decrease in nutrient concentrations in the soil was possibly due to the 429 decline in nutrient inputs from excrements of domestic herbivores and to the transfer of nutrients from the soil 430 to the biomass of the shrubs (Horton et al. 2009). The allocation of nutrients to the biomass of shrubs could be 431 a mechanism of nutrient accumulation by which shrubs control nutrients in an ecosystem through a positive 432 feed-back (Chapin et al. 1997) to outcompete herbaceous species and expand into the grassland (Grau et al. 16 433 2019). In our study we found that the concentrations of N, P and K often decreased in the aboveground 434 biomass along the succession in both study sites (Table 1); at the same time, the concentrations of soil N, P or 435 K often decreased in the mature shrubland of several shrubs. The stocks of N and P in the aboveground 436 biomass (leaves and stems), though, increased for several shrubs along the succession, particularly in R. 437 ferrugineum, J. communis and J. sabina (Figs. 3 and 4). This indicates that the dilution of these nutrients in 438 the aboveground biomass of shrubs was outweighed by their increase in aboveground biomass. This 439 evidences that shrubs may act as net reservoirs of essential nutrients in the biomass even if the concentration 440 in the biomass and soil decrease along the succession. The amount of C and nutrient stocks in aboveground 441 biomass differed notably amongst shrub species, which implied differences in the final stocks in the mature 442 shrubland within each site. Shrub identity was therefore crucial to understand the variability in concentrations 443 and stocks of nutrients in the vegetation and in the soil within each site. However, we also found differences 444 between the North and South face. Shrubs were generally smaller and less sclerophyllous at North face (Fig. 445 SF3), where solar radiation and evapotranspiration are lower, the slope is smoother, and the soil is deeper than 446 at South face (Ninot et al. 2010a). The increases in the C and nutrient stocks (g m-2) in the aboveground plant 447 compartments along the succession from grassland to shrubland were less evident at North face than at South 448 face. In contrast, shrubs developed more woody tissue, were taller and had higher wood density and more 449 sclerophyllous leaves at South face. This is probably because shrubs are adapted to harsher conditions, with 450 more solar radiation, higher evapotranspiration, less precipitation, a steeper slope and shallower soil with 451 more rocky outcrops than at North face. The differences in shrub traits therefore, seem to also play a key role 452 in the C and nutrient cycles with the advance of the encroachment, as previously described by other authors 453 (Eldridge et al. 2011). 454 We also found that there was a marginally significant decrease (p-value 0.08) of P stocks in the 455 young shrubland at North face, but not at South face, when vegetation and topsoil stocks were summed 456 together. This specific pattern could be related to differences in the way P is mobilized across the plant-soil 457 compartment in each site after the abandonment of the grassland, leading to putative losses of P from the 458 ecosystem. Moreover, A. uva-ursi at the South face site showed lower storage of N, P and K in the plant-soil 459 compartments in the mature shrubland stage, which also suggests a possible loss of nutrients from the system 17 460 with the advance of the encroachment dominated by this species. These results though should be explored in 461 more detail to fully understand the underlying mechanism of this loss of nutrients from the ecosystem. 462 Finally, herbaceous species in young shrubland and shrubs in young and mature shrubland showed 463 lower foliar δ15N and [N] than herbaceous species in unencroached grassland (Fig. 3). This was often coupled 464 with a decrease in soil [N] in the young and/or mature shrubland for several shrub species (Table 1). Ruiz- 465 Navarro et al. (2016) argued that low soil N and low foliar δ15N were good indicators of a decline in soil 466 fertility in Mediterranean ecosystems, with major implications on the biogeochemical cycling rates and 467 nutrient availabilities. In the same way, our results indicate that N may become limiting in these subalpine 468 ecosystems, particularly once the shrub cover is high (Angulo et al. 2019). Moreover, lower δ15N generally 469 indicates more N uptake by ectomycorrhizal or ericoid mycorrhizal fungi compared to arbuscular mycorrhizal 470 plants such as herbaceous species (Michelsen et al. 1998), with more recycled N leading to lower N losses 471 from the ecosystem (Garten 1993; Robinson 2001; Craine et al. 2009; Anadon-Rosell et al. 2016). All shrubs 472 at North face are ericaceous and have ericoid mycorrhizae; in the South face, A. uva-ursi hosts both ericoid 473 and ectomycorrhizal fungi, whereas J. communis and J. sabina host arbuscular mycorrhizal fungi. The lower 474 δ15N in young and mature shrubland possibly indicates that more N had been taken up by ectomycorrhizal or 475 ericoid fungi (Angulo et al. 2019). The N uptake by ectomycorrhizal or ericoid mycorrhizal fungi possibly 476 promotes the uptake of organic N (e.g. from litter) by shrubs (Akhmetzhanova et al. 2012) as inorganic N 477 becomes more limiting in the soil. The fact that the arbuscular mycorrhizal species occurring in the young or 478 mature shrubland (herbaceous species at North and South face and J. communis and J. sabina at South face) 479 also show low δ15N values compared to the unencroached grassland possibly indicates that these species may 480 have taken up part of their N in organic form derived from litter of surrounding ericoid or ectomycorrhizal 481 shrubs (e.g. from patches of A. uva-ursi), but this should be further investigated. These results suggest that the 482 advance of the succession possibly changed the mechanisms by which plants uptake N to overcome N 483 limitation. 484 CONCLUSIONS 485 Shrub expansion has a major influence on the plant-soil stoichiometry spectrum along the succession from 486 grassland to shrubland in the Pyrenees. Grassland ecosystem are dominated by species with faster turnover of 18 487 nutrients between the plant and soil compartments, high concentrations of N, P and K in the plant 488 aboveground biomass and topsoil, but limited capacity to store biomass. The expansion of shrubs, though, 489 favored the dominance of long-lived species, with a more conservative strategy, high C: nutrient ratios and 490 low concentrations of N, P and K in the aboveground biomass and low nutrient concentrations in soils. The 491 total stocks of C and nutrients in the aboveground biomass of the shrubland were nevertheless high because 492 the biomass of the mature shrubland was very high compared to the grassland. We thus highlight the role of 493 shrubs in the sequestration of C and nutrients, through the allocation to the aboveground biomass. Moreover, 494 shrub encroachment altered the strategy with which N was acquired, possibly through an increased uptake of 495 N through ectomycorrhizal or ericoid fungi. Our results suggested that shrubs play an important role in the C 496 and nutrient sequestration in the aboveground biomass (through the allocation into the plant compartments) 497 along the succession, where the woody tissues play a main role as storage compartment. The changes in plant- 498 soil elemental composition and stocks suggest a slowdown of the biogeochemical cycles in the subalpine 499 mountain areas where shrub encroachment occurred, but experimental manipulations are still needed to better 500 understand the mechanisms involved. 501 502 503 504 505 506 507 508 509 510 19 511 512 513 514 515 516 517 518 519 520 Figure captions 521 Fig. 1 Experimental design. The dotted lines represent the sampling area in each successional stage: (a) 522 Unencroached grassland, (b) young shrubland and (c) mature shrubland. The squares inside each successional 523 stage represent the replicate plots (approximately 2 × 2 m). In the young shrubland we established four groups 524 of four plots; the herbaceous species, shrub 1, shrub 2 and shrub 3 were grouped because all vegetation types 525 co-occurred. The distance between the successional stages or among the mature shrublands of each shrub 526 species was >100 m and the distance between the four replicates inside each successional stage was 10 m. 527 Fig. 2 Principal component analysis (PCA) based on the concentrations of chemical elements in leaves 528 (lf), leaf-litter (ltt) and soil (0-10 cm) for the successional stages and vegetation types at (a) North face and (b) 529 South face. Only variables that were poorly correlated with each other (Pearson‟s coefficients <0.6) were 530 included in the analyses. Different colors indicate the successional stages and different shapes indicate the 531 vegetation types. The ellipses represent the dispersion around the centroid for each successional stage with a 532 normal probability of 0.85. 533 Fig. 3 Polynomial regression between foliar N and 15N for herbaceous species and shrubs along the 534 succession from unencroached grassland to mature shrubland in both study sites. Different shapes show the 20 535 different vegetation types (herbaceous species or shrub) and colors show the successional stage 536 (unencroached grassland, young shrubland or mature shrubland) in each study site (NF: North face, SF: South 537 face). 538 Fig. 4 Stocks of chemical elements (C, N, P, K, Ca, Mg and Fe) in leaves, stems and soil (0-10 cm) (g m-2) 539 for each successional stage and vegetation type at the North face. The left column in each figure 540 represents the unencroached grassland, the columns in the middle represent the vegetation types in the young 541 shrubland, and the columns on the right represent the vegetation types in the mature shrubland. Significant 542 differences are summarized in Table ST3. Leaf-litter data are not included because the data available for leaf- 543 litter could only be used to estimate the productivity (g m-2 yr-1), not the stocks. See Fig. SF2 and Table ST2 544 for the leaf-litter productivity values. 545 Fig. 5 Stocks of chemical elements (C, N, P, K, Ca, Mg and Fe) in leaves, stems and soil (0-10 cm) (g m-2) 546 for each successional stage and vegetation type at the South face. The left column in each figure 547 represents the unencroached grassland, the columns in the middle represent the vegetation types in the young 548 shrubland, and the columns on the right represent the vegetation types in the mature shrubland. Significant 549 differences are summarized in Table ST5. Leaf-litter data are not included because the data available for leaf- 550 litter could only be used to estimate the productivity (g m-2 y-1), not the stocks. See Fig. SF2 and Table ST2 551 for the leaf-litter productivity values. 552 21 553 Fig. 1 Experimental design. The dotted lines represent the sampling area in each successional stage: (a) Unencroached grassland, (b) young shrubland and (c) 554 mature shrubland. The squares inside each successional stage represent the replicate plots (approximately 2 × 2 m). In the young shrubland we established four 555 groups of four plots; the herbaceous species, shrub 1, shrub 2 and shrub 3 were grouped because all vegetation types co-occurred. The distance between the 556 successional stages or among the mature shrublands of each shrub species was >100 m and the distance between the four replicates inside each successional stage 557 was 10 m. 558 559 22 560 Fig. 2 Principal component analysis (PCA) based on the concentrations of chemical elements in leaves (lf), leaf-litter (ltt) and soil (0-10 cm) for the 561 successional stages and vegetation types at (a) North face and (b) South face. Only variables that were poorly correlated with each other (Pearson‟s coefficients 562 <0.6) are included in the analyses. Different colors indicate the successional stages and different shapes indicate the vegetation types. The ellipses represent the 563 dispersion around the centroids for each successional stage with a normal probability of 0.85. 564 565 566 23 567 Fig. 3 Polynomial regression between foliar N (%) and foliar delta 15N (15N) for herbaceous species and 568 shrubs along the succession from unencroached grassland to mature shrubland in both study sites. 569 Different shapes show the different vegetation types (herbaceous species or shrub) and colors show the 570 successional stage (unencroached grassland, young shrubland or mature shrubland) in each study site (NF: 571 North face, SF: South face) 572 573 574 575 576 577 578 579 580 24 581 Fig. 4 Stocks of chemical elements (C, N, P, K, Ca, Mg and Fe) in leaves, stems and soil (0-10 cm) (g m-2) 582 for each successional stage and vegetation type at the North face. The left column in each figure 583 represents the unencroached grassland, the columns in the middle represent the vegetation types in the young 584 shrubland, and the columns on the right represent the vegetation types in the mature shrubland. Significances 585 are summarized in table ST3. Leaf-litter data are not included because the data available for leaf-litter could 586 only be used to estimate the productivity (g m-2 y-1), not the stocks. See Fig. SF2 and Table ST2 for the leaf- 587 litter productivity values. 25 588 Fig. 5 Stocks of chemical elements (C, N, P, K, Ca, Mg and Fe) in leaves, stems and soil (0-10 cm) (g m-2) 589 for each successional stage and vegetation type at the South face. The left column in each figure 590 represents the unencroached grassland, the columns in the middle represent the vegetation types in the young 591 shrubland, and the columns on the right represent the vegetation types in the mature shrubland. Significances 592 are summarized in table ST5. Leaf-litter data are not included because the data available for leaf-litter could 593 only be used to estimate the productivity (g m-2 y-1), not the stocks. See Fig. SF2 and Table ST2 for the leaf- 594 litter productivity values. 26 595 Table 1. Significant changes in the element concentrations for the different vegetation types in (a) the young shrubland and (b) mature shrubland 596 compared to the herbaceous species in the unencroached grassland. Changes (increases (  ) or decreases (  )) in the elemental concentrations (C, N, P, K, 597 Ca, Mg and Fe) (dw/dw) and in the C:N, C:P and N:P ratios were identified by generalized least square models. The n.s. abbreviation indicates that the results 598 were not 599 Increase or decrease in the elemental concentration and ratios between the young shrubland and statisticall (a) unencroached grassland stage Leaves Leaf-litter Soil (0-10 cm) 600 y Chemical changes (dw/dw)       Ca 601 significan herbaceous species C, N, K, Fe N, K, Mg, Fe P C:N C:N C:N P, Mg 602 C, Ca, Mg Ca t (P > Calluna vulgaris N, P, K C, N, Ca, Mg K, Fe N:P, C:N, C:P, N:P, C:N, C:P C, Ca, Mg C, Ca, Mg Ca, K North Rhododendron ferrugineum N, K, Fe N, K, Fe N, P, Mg603 0.05). C:N C:N, C:P C:N Face C, Ca, Mg K, Ca Vaccinium myrtillus N, K, Fe C, Ca, K, Mg Fe N, P, Mg C:N 604 Means  C, Fe Fe herbaceous species n.s N, K N, P N, Ca, Mg N:P, C:N 605 sd are C N, K C, K, Mg N, Fe Fe N, K, Ca, Mg Arctostaphylos uva-ursi C:N, C:P N:P C:N N:P C:N N:P South C, Ca N, K C, Ca N, Ca, Mg,606 Fe summariz Juniperus communis N, P, Fe n.s Face C:N, C:P N:P C:P N:P Juniperus sabina C, Ca, Mg N, K C, K, Ca, Mg n.s Fe N, Ca, Mg607 ed in C:N N:P N:P (b) 608 Increase or decrease in the elemental concentration and ratios between the mature shrubland and Table unencroached grassland stage K, Ca, Fe609 C, Ca, Mg ST2 and Calluna vulgaris N, P, K, Fe C, Ca, Mg n.s N:P, C:N, C:P C:N, C:P C, Ca, Mg C, Ca, Mg C, Ca North Rhododendron ferrugineum C:N, C:P N, K, Fe C:N, C:P N, P, K, Fe N:P, C:N, C:P Mg 610 significan Face C, Ca, Mg N, K, Fe Ca Vaccinium myrtillus C, P, K, Ca, Mg Fe n.s C:N N:P C:N 611 ce values C, Ca, Mg N, P, K C, Mg N, P, Fe Arctostaphylos uva-ursi N, P, K, Ca, Mg N:P, C:N, C:P C:N, C:P N:P C:N C, Ca N, K C, Ca N, Ca, Mg 612 are shown Juniperus communis N, P, Fe n.s South C:N N:P C:N, C:P N:P Face Juniperus sabina C, Ca, Mg N, K C, P, K, Mg, Ca Fe n.s N, P, Ca, Mg 613 in Tables C:N N:P N:P 27 614 ST3 and ST5. 615 616 617 618 619 620 621 622 623 624 625 626 Table 2. Significant changes in the element concentrations for the vegetation types in the mature shrubland compare to young shrubland. Changes 627 (increases or decreases) in the elemental concentrations (C, N, P, K, Ca, Mg and Fe) (dw/dw) and in the C:N, C:P and N:P ratios were identified by generalized 28 628 least square models. The n.s. abbreviation indicates that the results were not statistically significant (P > 0.05). Means  sd are summarized in Table ST2 and 629 significance values are shown in Tables ST4 and ST6. 630 631 Changes in the elemental composition between the mature and young shrubland stage 632 Leaves Leaf-litter Stems Soil (0-10 cm) Chemical changes        633 Calluna vulgaris K n.s n.s n.s Mg N, K n.s C:N C:N North n.s634 Rhododendron ferrugineum N Mg n.s N, K C:N face N:P C:N N:P N:P, C:P Vaccinium myrtillus C, Fe n.s n.s n.s n.s Mg n.s N:P, C:P 635 Arctostaphylos uva-ursi P n.s n.s n.s n.s n.s C:P N:P South Juniperus communis n.s Ca n.s Face n.s n.s Ca n.s636 C:N Juniperus sabina C Ca n.s C, N, P Mg n.s n.s C:N, C:P C:N 637 638 639 640 29 641 Reference 642 Ågren GI (2004) The C : N : P stoichiometry of autotrophs - theory and observations. Ecol Lett 7:185–191. 643 doi: https://doi.org/10.1111/j.1461-0248.2004.00567.x 644 Akhmetzhanova AA, Soudzilovskaia NA, Onipchenko VG, et al (2012) A rediscovered treasure: mycorrhizal 645 intensity database for 3000 vascular plant species across the former Soviet Union. Ecology 93:689–690. 646 doi: 10.1890/11-1749.1 647 Ameztegui A, Coll L, Brotons L, Ninot JM (2016) Land-use legacies rather than climate change are driving 648 the recent upward shift of the mountain tree line in the Pyrenees. Glob Ecol Biogeogr 25:263–273. doi: 649 https://doi.org/10.1111/geb.12407 650 Anadon-Rosell A, Palacio S, Nogués S, Ninot JM (2016) Vaccinium myrtillus stands show similar structure 651 and functioning under different scenarios of coexistence at the Pyrenean treeline. Plant Ecol 217:1115– 652 1128. doi: 10.1007/s11258-016-0637-2 653 Angulo MA, Ninot JM, Peñuelas J, et al (2019) Tree Sapling Responses to 10 Years of Experimental 654 Manipulation of Temperature, Nutrient Availability, and Shrub Cover at the Pyrenean Treeline. Front 655 Plant Sci 9:1–13. doi: 10.3389/fpls.2018.01871 656 Anthelme F, Villaret J-C, Brun J-J (2007) Shrub encroachment in the Alps gives rise to the convergence of 657 sub-alpine communities on a regional scale. J Veg Sci 18:355–362. doi: 10.1111/j.1654- 658 1103.2007.tb02547.x 659 Barrio IC, Bueno CG, Nagy L, et al (2013) Alpine Ecology in the Iberian Peninsula: What Do We Know, and 660 What Do We Need to Learn? Mt Res Dev 33:437–442. doi: https://doi.org/10.1659/MRD-JOURNAL- 661 D-13-00052.1 662 Boixadera J, Ascaso E, Herrero C, Jaquet J (2014) Els sòls de Catalunya Val d‟Aran. In: ICGC-DAAM- 663 Conselh Generau d‟Aran 664 Carreras J, Carrillo A, Font X, et al (1983) Els prats de l‟alianca Xerobromion als Pirineus Catalans. Collect 665 Bot 14:151–209 30 666 Catalan J, Ninot JM, Aniz MM (2017) High Mountain Conservation in a Changing World, 1st edn. Cham 667 Chapin FS, Chapin FS, Walker BH, et al (1997) Biotic control over the functioning of ecosystems. Science 668 (80- ) 277:500–504. doi: 10.1126/science.277.5325.500 669 Chapin FS, Körner C (1994) Arctic and alpine biodiversity: patterns, causes and ecosystem consequences. 670 Tree 9:45–47. doi: 10.1016/0169-5347(94)90266-6 671 Chave J, Coomes D, Jansen S, et al (2009) Towards a worldwide wood economics spectrum. Ecol Lett 672 12:351–366. doi: 10.1111/j.1461-0248.2009.01285.x 673 Craine JM, Elmore AJ, Aidar MPM, et al (2009) Global patterns of foliar nitrogen isotopes and their 674 relationships with climate, mycorrhizal fungi, foliar nutrient concentrations, and nitrogen availability. 675 New Phytol 183:980–992. doi: 10.1111/j.1469-8137.2009.02917.x 676 Dullinger S, Dirnbock T, Greimler J, et al (2003) A resampling approach for evaluating effects of pasture 677 abandonment on subalpine plant species diversity. J Veg Sci 14:243–252. doi: 10.1111/j.1654- 678 1103.2003.tb02149.x 679 EAE (2019) Middle Europena Xerobromion grasslands. https://eunis.eea.europa.eu/habitats/3643 680 Eldridge DJ, Bowker MA, Maestre FT, et al (2011) Impacts of shrub encroachment on ecosystem structure 681 and functioning: Towards a global synthesis. Ecol Lett 14:709–722. doi: 10.1111/j.1461- 682 0248.2011.01630.x 683 Ferré A, Ninot JM, Carreras J, et al (2013) From vegetation mapping to a prediction for landscape evolution 684 in the Catalan Pyrénées. In: Colloques Phytosociologiques 685 Formica A, Farrer EC, Ashton IW, Suding KN (2014) Shrub Expansion Over the Past 62 Years in Rocky 686 Mountain Alpine Tundra : Possible Causes and Consequences. Arctic, Antarct Alp Res 46:616–631 687 Galvánek D, Janák M (2008) Management of Natura 2000 habitats.Species-rich Nardus grasslands. 688 Garcia-Pausas J (2010) Factors determinants de l ‟ acumulació de carboni edàfic als prats de muntanya del 689 Pirineu. Universitat de Lleida 31 690 Garcia-Pausas J, Casals P, Camarero L, et al (2007) Soil organic carbon storage in mountain grasslands of the 691 Pyrenees: effects of climate and topography. Biogeochemistry 82:279–289. doi: 10.1007/s10533-007- 692 9071-9 693 Garcia-Pausas J, Romanyà J, Montané F, et al (2017) Are Soil Carbon Stocks in Mountain Grasslands 694 Compromised by Land-Use Changes? In: Catalan J, Ninot JM, Aniz MM (eds) High Mountain 695 Conservation in a Changing World, 1st edn. Cham, pp 207–230 696 Garten CT (1993) Variation in foliar 15N abundance and the availability of soil Nitrogen on Walker Branch 697 watershed. Ecology 74:2098–2113 698 Gassiot E, Jiménez J (2006) El poblament prefeudal de l‟alta muntanya dels Pirineus occidentals catalans 699 (Pallars Sobirà i Alta Ribagorça). Trib. d‟ Arqueol. 89–122 700 Gassiot E, Rodríguez-Antón D, Burjachs F, et al (2016) Poblamiento, explotación y entorno natural de los 701 estadios alpinos y subalpinos del Pirineo central durante la primera mitad del Holoceno. Cuaternario y 702 Geomorfol 30:119–130. doi: 10.17735/cyg.v30i1-2.41464 703 Grau O, Ninot JM, Cornelissen JHC, Callaghan T V. (2013) Similar tree seedling responses to shrubs and to 704 simulated environmental changes at Pyrenean and subarctic treelines. Plant Ecol Divers 6:329–342. doi: 705 10.1080/17550874.2013.810311 706 Grau O, Peñuelas J, Ferry B, et al (2017) Nutrient-cycling mechanisms other than the direct absorption from 707 soil may control forest structure and dynamics in poor Amazonian soils. Sci Rep. doi: 708 10.1038/srep45017 709 Grau O, Saravesi K, Ninot JM, et al (2019) Encroachment of shrubs into subalpine grasslands in the Pyrenees 710 modifies the structure of soil fungal communities and soil properties. FEMS Microbiol Ecol 1–16. doi: 711 10.1093/femsec/fiz028 712 Hallinger M, Manthey M, Wilmking M (2010) Establishing a missing link: warm summers and winter snow 713 cover promote shrub expansion into alpine tundra in Scandinavia. New Phytol 186:890–899. doi: 714 10.1111/j.1469-8137.2010.03223.x 32 715 Hibbard KA, Archer S, Schimel DS, Valentine DW (2001) Biogeochemical changes accompanying woody 716 plant encroachment in a subtropical savanna. Ecology 82:1999–2011. doi: 10.1890/0012- 717 9658(2001)082[1999:BCAWPE]2.0.CO;2 718 Hood EW, Williams MW, Caine N (2003) Landscape controls on organic and inorganic nitrogen leaching 719 across an alpine/subalpine ecotone, Green Lakes Valley, Colorado Front Range. Ecosystems 6:31–45. 720 doi: 10.1007/s10021-002-0175-8 721 Horton JL, Clinton BD, Walker JF, et al (2009) Variation in Soil and Forest Floor Characteristics Along 722 Gradients of Ericaceous, Evergreen Shrub Cover in the Southern Appalachians. Castanea 74:340–352. 723 doi: 10.2179/08-016R3.1 724 Illa E, Ninot JM, Anadon-Rosell A, Oliva F (2017) The role of abiotic and biotic factors in functional 725 structure and processes of alpine shrub communities. Folia Geobot 52:199–215. doi: 10.1007/s12224- 726 017-9296-x 727 IPCC (2013) Summary for Policymakers. Clim Chang 2013 Phys Sci Basis Contrib Work Gr I to Fifth Assess 728 Rep Intergov Panel Clim Chang 33. doi: 10.1017/CBO9781107415324 729 Jackson R, Banner J, Jobbágy E (2002) Ecosystem carbon loss with woody plant invasion of grasslands. 730 Nature 277:623–627. doi: 10.1038/nature00952. 731 Jiménez MB, Pujol AFT (2010) Els canvis socioambientals al Pirineu català. La Vall Fosca com a escenari 732 representatiu de les mutacions en les economies tradicionals de muntanya. Doc. d‟Analisi Geogr. 733 56:557–572 734 Knapp AK, Briggs JM, Collins SL, et al (2008) Shrub encroachment in North American grasslands: Shifts in 735 growth form dominance rapidly alters control of ecosystem carbon inputs. Glob Chang Biol 14:615– 736 623. doi: 10.1111/j.1365-2486.2007.01512.x 737 Komac B, Kefi S, Nuche P, et al (2013) Modeling shrub encroachment in subalpine grasslands under different 738 environmental and management scenarios. J Environ Manage 121:160–169. doi: 739 10.1016/j.jenvman.2013.01.038 33 740 Körner C, Paulsen J (2004) A World-Wide Study of High Altitude Treeline Temperatures. J Biogeogr 741 31:713–732 742 Lasanta T, Vicente-Serrano SM, Cuadrat-Prats JM (2000) Marginación productiva y recuperación de la 743 cubierta vegetal en el Pirineo: Un caso de estudio en el valle de Borau. Boletín la A.G.E. 29:5–28 744 Le S, Josse J, Husson F (2008) FactoMineR: An R Package for Multivariate Analysis. J. Stat. Softw. 25 745 MacDonald D, Crabtree J., Wiesinger G, et al (2000) Agricultural abandonment in mountain areas of Europe: 746 Environmental consequences and policy response. J Environ Manage 59:47–69. doi: 747 10.1006/jema.1999.0335 748 Michelsen A, Quarmby C, Sleep D, Jonasson S (1998) Vascular plant 15N natural abundance in heath and 749 forest tundra ecosystems is closely correlated with presence and type of mycorrhizal fungi in roots. 750 Oecologia 115:406–418. doi: 10.1007/s004420050535 751 Montané F, Romanyà J, Rovira P, Casals P (2010) Aboveground litter quality changes may drive soil organic 752 carbon increase after shrub encroachment into mountain grasslands. Plant Soil 337:151–165. doi: 753 10.1007/s11104-010-0512-1 754 Montané F, Rovira P, Casals P (2007) Shrub encroachment into mesic mountain grasslands in the Iberian 755 peninsula: Effects of plant quality and temperature on soil C and N stocks. Global Biogeochem Cycles 756 21:1–10. doi: 10.1029/2006GB002853 757 Myers-Smith IH, Forbes BC, Wilmking M, et al (2011) Shrub expansion in tundra ecosystems: dynamics, 758 impacts and research priorities. Environ Res Lett 6:045509. doi: 10.1088/1748-9326/6/4/045509 759 Naito AT, Cairns DM (2011) Patterns and processes of global shrub expansion. Prog Phys Geogr 35:423–442. 760 doi: 10.1177/0309133311403538 761 Ninot JM, Anadon-Rosell A, Caminal M, et al (2010a) Los arbustos de alta montaña y los cambios de paisaje 762 en el parque natural de Aigües tortes i Estany de Sant Maurici. Proy Investig en Parques Nac 2012-2015 763 209–228 764 Ninot JM, Batllori E, Carreras J, et al (2010b) Reforestación natural en el dominio de Pinus uncinata del 34 765 parque nacional de AigüesTortes i Estany de Sant Maurici. Proy Investig en Parques Nac 2007-2010 766 139–158 767 Oksanen J, Blanchet FG, Friendly M, et al (2017) vegan: Community Ecology Package. R package version 768 2.4-5. 769 Pèlachs A, Soriano JM, Nadal J, Esteban A (2007) Holocene environmental history and human impact in the 770 Pyrenees. Contrib to Sci 3:421–429. doi: 10.2436/20.7010.01.19 771 Peñuelas J, Fernández-Martínez M, Ciais P, et al (2019) The bio-elements, the elementome and the 772 “biogeochemical niche.” Ecology. doi: 10.1002/ecy.2652 773 Peñuelas J, Sardanas J, Filella I, et al (2016) Tercer informe sobre el canvi climàtic a Catalunya 774 Pinheiro J, Bates D, DebRoy S, et al (2017) nlme: Linear and Nonlinear Mixed Effects Models. R package 775 version 3.1-131 776 RCoreTeam (2017) R: A language and environment for statistical computing. R Foundation for Statistical 777 Computing, Vienna, Austria. 778 Robinson D (2001) δ15N as an integrator of the nitrogen cycle. TRENDS Eol Evol 16:153–162 779 Roura-Pascual AN, Pons P, Etienne M, et al (2005) Transformation of a Rural Landscape in the Eastern 780 Pyrenees Between 1953 and 2000. Montain Res develoment 25:252–261 781 Ruiz-Navarro A, Barber GG, Albaladejo J, et al (2016) Plant d 15 N reflects the high landscape-scale 782 heterogeneity of soil fertility and vegetation productivity in a Mediterranean semiarid ecosystem. New 783 Phytol 1030–1043 784 Sanz-Elorza M, Dana ED, González A, Sobrino E (2003) Changes in the high-mountain vegetation of the 785 central Iberian Peninsula as a probable sign of global warming. Ann Bot 92:273–280. doi: 786 10.1093/aob/mcg130 787 Sardans J, Peñuelas J (2015) Trees increase their P: N ratio with size. Glob Ecol Biogeogr 24:147–156. doi: 788 10.1111/geb.12231 35 789 Sardans J, Peñuelas J (2013) Tree growth changes with climate and forest type are associated with relative 790 allocation of nutrients, especially phosphorus, to leaves and wood. Glob Ecol Biogeogr 22:494–507. 791 doi: 10.1111/geb.12015 792 Sardans J, Rivas-Ubach A, Peñuelas J (2012) The C:N:P stoichiometry of organisms and ecosystems in a 793 changing world: A review and perspectives. Perspect Plant Ecol Evol Syst 14:33–47. doi: 794 10.1016/j.ppees.2011.08.002 795 Serrano V, Lasanta T, Cuadrat J (2000) Transformaciones en el paisaje del Pirineo como consecuencia del 796 abandono de las actividades económicas tradicionales. Pirineos 155:111–133. doi: 797 10.3989/egeogr.2001.i245.267 798 Tape K, Sturm M, Racine C (2006) The evidence for shrub expansion in Northern Alaska and the Pan-Arctic. 799 Glob Chang Biol 12:686–702. doi: 10.1111/j.1365-2486.2006.01128.x 800 Targetti S, Staglianò N, Messeri A, Argenti G (2010) A state-and-transition approach to alpine grasslands 801 under abandonment. IForest 3:44–51. doi: 10.3832/ifor0525-003 802 Throop HL, Archer SR (2008) Shrub (Prosopis velutina) encroachment in a semidesert grassland: Spatial- 803 temporal changes in soil organic carbon and nitrogen pools. Glob Chang Biol 14:2420–2431. doi: 804 10.1111/j.1365-2486.2008.01650.x 805 Urbina I, Sardans J, Grau O, et al (2017) Plant community composition affects the species biogeochemical 806 niche. Ecosphere 8:1–14. doi: https://doi.org/10.1002/ecs2.1801 807 Van-Auken OW (2009) Causes and consequences of woody plant encroachment into western North American 808 grasslands. J Environ Manage 90:2931–2942. doi: 10.1016/j.jenvman.2009.04.023 809 Vitousek PM (1984) Litterfall, Nutrient Cycling, and Nutrient Limitation in Tropical Forests. Ecology 810 65:285–298. doi: 10.2307/1939481 811 Wickham H (2009) ggplot2: Elegant Graphics for Data Analysis. ggplot2 115–137 812 Wookey PA, Aerts R, Bardgett RD, et al (2009) Ecosystem feedbacks and cascade processes: Understanding 813 their role in the responses of Arctic and alpine ecosystems to environmental change. Glob Chang Biol 36 814 15:1153–1172. doi: 10.1111/j.1365-2486.2008.01801.x 815 Zechmeister-Boltenstern S, Keiblinger KM, Mooshammer M, et al (2015) The application of ecological 816 stoichiometry to plant-microbial-soil organic matter transformations. Ecol Monogr 85:133–155. doi: 817 10.1890/14-0777.1 818 37 819 Supplementary material 820 Figure SF1. Map of the study sites. Location of (a) the study area in the Central Pyrenees (red square) in southwestern Europe and (b) the two study sites (red 821 dots). 822 823 824 825 826 827 828 829 830 38 831 Figure SF2. Leaf-litter productivity (g m2 y-1) for each vegetation type in the two study sites: (a) North face and (b) South face in the unencroached grassland 832 (left panel in each figure), young shrubland (middle panel in each figure) and mature shrubland (right panel in each figure). Leaf-litter productivity was calculated 833 using the foliar persistence analysis from Ninot et al. (unpublished data). The upper and lower line of the boxes indicates the second and third quartile, 834 respectively; the black line inside de box indicates the median; the bars represent the first and fourth quartiles and dots represent the outliers. Note the different 835 vertical scales for North face and South face. 836 837 838 39 839 Figure SF3. Principal component analysis of the functional traits and foliar elemental composition of 840 the shrub species in the mature shrubland including both study sites. The ellipses represent the dispersion 841 around the centroids for each shrub species with a normal probability of 0.85. 842 843 844 845 846 847 848 40 849 Figure SF4. Principal component analysis (PCA) based on the elemental concentrations for (a) leaves, 850 (b) leaf-litter, (c) soil and (d) stems for all vegetation types and successional stages. Only variables that 851 were poorly correlated with each other (Pearson‟s coefficients <0.6) were included in the analyses. Different 852 colors indicate the successional stages and different shapes indicate the vegetation types. The ellipses in (a) 853 and (b) represent the dispersion around the centroids for herbaceous species and shrubs, and in (c) and (d) for 854 the different sites, with a normal probability of 0.85. 855 856 857 41