Vegetation trends over eleven years on mountain summits in NW Argentina

Abstract As global climate change leads to warmer and dryer conditions in the central Andes, alpine plant communities are forced to upward displacements following their climatic niche. Species range shifts are predicted to have major impacts on alpine communities by reshuffling species composition and abundances. Using a standardized protocol, we surveyed alpine plant communities in permanent plots on four high Andean summits in NW Argentina, which range from 4,040 to 4,740 m a.s.l. After a baseline survey in 2006–2008, we resurvey the same plots in 2012, and again in 2017. We found a significant decrease in plant cover, species richness, and diversity across the elevation gradient in the three censuses and a strong decrease in soil temperature along the elevation gradient. We found a high plant community turnover (37%–49%) among censuses, differentiating according to summits and aspects; major changes of community turnover were observed in the lowest summit (49%) and on the northern (47%) and western (46%) aspects. Temporal patterns in community changes were represented by increases in plant cover in the highest summit, in species richness in the lower summit, and in diversity (Shannon index) in the four summits, over time, together with increase in small herbs and non‐tussock grasses. We suggest that the observed trend in plant community dynamics responds to short‐term temperature and precipitation variability, which is influenced by El Niño Southern Oscillation (ENSO), and due to time lags in plant community response, it may take much longer than one decade for the observed trends to become stables and statistically significant. Our study provides an important foundation for documenting more profound changes in these subtropical alpine plant communities as global climate change continues.


| INTRODUC TI ON
High mountain ecosystems are sensitive to climate change, mainly because their organisms are governed by low-temperature conditions (Halloy, 1989;Körner, 2012;Pauli et al., 2012). High mountain ecosystems are often less affected by direct anthropogenic land use impacts and biotic factors, such as competition, which decrease with altitude (Llambi, Law, & Hodge, 2004;Mark et al., 2015), as environmental stress increases (Anthelme, Meneses, Valero, Pozo, & Dangles, 2017;Halloy, 1989). Thus, high mountain ecosystems can be considered "natural experiments" to study the impact of climate change on vegetation, and useful for global scale comparisons . GLORIA (Global Observation Research Initiative in Alpine Environments) is an international network of long-term monitoring sites that have been developed to assess the response of vegetation to climate change and compare it across high mountains all over the world Pauli et al., 2012). Within GLORIA, the GLORIA-Andes Network links sites in Andean countries (Cuesta et al., 2012(Cuesta et al., , 2017. Climate change is one of the main factors modifying high mountain vegetation and may result in species migration, adaptation, or extinction in the coming decades (Hoegh-Guldberg et al., 2008;Parmesan, 2006;Pauli, Gottfried, Reiter, Klettner, & Grabherr, 2007). A reconstruction based on tree rings from the high Andes of NW Argentina showed a consistent aridity trend for the last decades with an increase in the frequency of drought events (Morales, Carilla, Grau, & Villalba, 2015). The observed trend of increasing aridity along with projected temperature increases (Urrutia & Vuille, 2009;Vuille et al., 2017) will make high mountain summits of NW Argentina too dry and warm for alpine plant specialist (cryophilic species), and therefore, they may undergo a range contraction and/ or become locally extinct. Evidence for upward displacement of species and increase in community richness at high altitudes has been recorded in different mountain ranges of Europe Erschbamer, Kiebacher, Mallam, & Unterluggauer, 2009;Grytnes et al., 2014;Pauli et al., 2012;Wipf, Stöckli, Herz, & Rixen, 2013), Africa (Hemp, 2009), and South America (Moret, Aráuz, Gobbi, & Barragán, 2016;Morueta-Holme et al., 2015;Seimon et al., 2007Seimon et al., , 2017, and also differences in species growth among summits aspects (Sklenář, Kučerová, Macková, & Romoleroux, 2016). Further, population decline and losses in genetic diversity due to warming driven range reduction have been observed in narrow-range tropical alpine species (Chala et al., 2016). In addition, shifts in species dominance and plant cover may occur, which is the most detectable change in short-to medium-term studies (Halloy, 2002;Seimon et al., 2009;Yager, Resnikowski, & Halloy, 2008).
Long-term trends are often masked by short-term dynamics and variability (Perez et al., 2010). Historic records and multidecadal time series are needed to understand the natural range of variability (Landres, Morgan, & Swanson, 1999) and assess recent trends in ecosystems and biodiversity changes. Continued observations of vegetation dynamics (changes in species composition and cover) in the high Andes encompass short periods of monitoring, most of them from the beginning of the 21st century (Cuesta et al., 2012(Cuesta et al., , 2017. Few historic records and environmental proxies provide a wider perspective of the natural variability of the ecosystem versus changes caused by human activities (Halloy, 2002;Morales et al., 2015;Morueta-Holme et al., 2015).

Ground truthing through repeated surveys in permanent plots
is essential to answer the following questions: How is the distribution of alpine vegetation affected by climate change across an elevation gradient? Will climate change cause species to shift their range upward in elevation? Is it possible to detect changes in alpine plant community's composition after 10 years of continuous observation?
The aim of this study was to analyze alpine plant species composition along summits and aspects, and to assess plant community dynamics over a period of 9-11 years (from 2006/08 to 2017) on four summits, ranging from 4,040 to 4,740 m a.s.l., in Cumbres Calchaquíes, Tucumán, Argentina, in order to establish the foundations for understanding the long-term trends from vegetation patterns of shortterm variability. To explain the observed vegetation patterns and dynamics, we related vegetation variables and plant life forms with soil temperature. We hypothesized that (1) temperature regime, controlled by elevation and aspect is the main factor determining the prevailing vegetation and the distribution of plant species, and (2) even short-term thermal changes (over several years) lead to a shift in species composition, and changes in plant cover and diversity.

| Study area
A GLORIA target region was established in Parque Provincial Cumbres Calchaquíes, Tucumán, NW Argentina (26°40' S 65°44' W; Lomáscolo et al., 2014) in 2006-08. Geologically, this mountain range belongs to the Pampean system and is characterized by a plateau at 4,300 m a.s.l., with dozens of lakes of glacial origin, referred to as Huaca Huasi lakes. The dominant bedrock is Precambrian metamorphic rocks, from solid crag to quaternary moraine accumulation (Halloy, 1985a). The area is within the High Andean ecoregion (Cabrera, 1976), Calchaquí subdistrict with high levels of plant diversity and restricted distribution ranges (Halloy, 1985b). There are three main vegetation communities dominated by (a) Festuca ortophylla (Poaceae, Iro grasslands); (b) ground-hugging plants such as Adesmia crassicaulis (Fabaceae), Tetraglochin inerme (Rosaceae), and Pycnophyllum convexum (Caryophyllaceae), a community known as cryptofruticetum (up to 20 species per 1m 2 ; Figure 1); and (c) typical wetland vegetation communities of cushion bogs (Halloy & Laurent, 1988;Halloy, 1985a). This area presents some advantages for long-term studies, such as vegetation records since 1976 (Halloy, 2002), lake level records going back to the 1980s (Casagranda, 2010;Halloy, 2002), discontinuous soil temperatures since the 1970´s (Halloy & Laurent, 1988;Halloy, 1983Halloy, , 1985a) and a detailed physical characterization of the area. The main land uses are recreation (hiking) and cattle grazing at the lowest summit. Herbivory by native camelids (Lama guanicoe) and big rodents (Lagidium viscacia) is present in the whole area, although the limited accessibility and isolation has favored ecosystem conservation.
Direct measurement of annual precipitation at 4,250 m a.s.l. are available for 1977-1978, recording c. 385 mm on average, mainly as snow, concentrated in the summer months (December to March), and mean annual air temperature was estimated as 1.5ºC (Halloy, 1985a). A recently installed meteorological station at 4,200 m a.s.l., recorded an annual precipitation of 333 mm and an average temperature of 2.9°C in 2016.

| Sampling design, data recording, and preparation
Summit site setup and vegetation surveys followed GLORIA sampling methodology (Halloy, Ibáñez, & Yager, 2011;Pauli et al, 2015). Four summits were selected at different altitudes, with the same bedrock and under the same regional climatic conditions. Summits were coded from lower to higher altitude as: 1-ALZ, 2-HUA, 3-SIN, We recorded vascular plant species and their relative percentage cover, lichen, and bryophyte cover (total), and substrate cover: solid rock (fixed in the ground), gravel/scree (debris material), bare soil (including clay and sand), organic matter, and feces (from camelids and rodents in all summits, and cattle in the lower summit).
Vascular plants were classified into life forms, following Ramsay and Oxley (1997) adapted for the GLORIA-Andes network, as shrubs F I G U R E 1 Plant community known as cryptofruticetum, (a) composed by flattened vegetation, with great diversity of tiny colorful flowers in summer. (b) Adesmia crassicaulis (in the hand) one of the species of the cryptofruticetum community observed in the ground, in Cumbres Calchaquíes, Tucumán, Argentina F I G U R E 2 Study area. GLORIA summits indicated with a red point in a Landsat image. Tucuman, Argentina map (left corner). Sampling design scheme showing 1 m 2 plot, the position of data loggers and summit lower limit (right corner, from Pauli et al., 2015). Linear regression between mean soil temperature and elevation is superposed with the image (represent the main hypothesis). The main spatial and temporal trends in vegetation variables are represented by graphs in the right; symbols indicate increasing trends over time (including subshrubs), cushion plants (lax and hard cushions), tussock grasses, non-tussock grasses (cespitose and tufted grasses), erect herbs, prostrate herbs (including trailing and creeping herbs), rosettes (including basal and acaulescent), and ferns.
We estimated species composition and richness at 1) summit scale based on semiquantitative criteria for visual estimation, by subdividing the summits in sections and walking them thoroughly, at the baseline survey, and by using flexible points and areas method (PAF; Halloy et al., 2011), for the resurveys. PAF method is useful for quantifying species cover in relation to a total area with minimum sampling effort. However, in this study, we used PAF in order to record the complete species list of each summit. The lower limit of the summit was established at 10 m vertically down from the highest summit point (HSP). We measured species richness and plant cover at 2) 1 m × 1 m plots, with four plots for each aspect (north, east, south, and west, resulting in 16 plots for each summit). Each 1 m 2 plot was subdivided into 100 subplots of 0.01 m 2 , equivalent to 1% cover with 100 crosshair points.
We estimated species diversity using the Shannon-Weaver index

| Soil temperature records
We measured soil temperature with data loggers (tidbit, Onset) at four points (north, east, south, and west) on each summit, located at 5 m vertically down from the HSP at 10 cm below the soil surface, with a recording frequency of every two hours (i.e., 12 records per day

| Data analysis
In order to assess temporal trends in temperature series data, we used linear regression, analyzing each data logger record for 2009-2015 period (seven data loggers; function lm in R; Chambers, 1992).
We explored the relationship between the seven temperature variables with elevation using linear regression (function lm in R) and aspects using two-way ANOVA (function aov in R. Chambers, Freeny, & Heiberger, 1992).
We analyzed trends over three vegetation metrics ( We used linear mixed models and generalized linear mixed models (LMM and GLMM; Bolker et al., 2009) to analyze variability in plant cover, species richness, and diversity index among aspects within summits, among summits, and among censuses, with aspect and census as fixed effects and summit as a random effect. For model construction, we considered Gaussian distribution for plant cover (square root transformed data) and Shannon index, using the package lme4 in R (Bates, Kliegl, Vasishth, & Baayen, 2015), and followed by lmerTest (Kuznetsova, Brockhoff, & Christensen, 2017), to obtain p values for fixed effects, and Poisson distribution for species richness using the package GLMM (Knudson, 2016).
To assess compositional similarity among summits, we used nonmetric multidimensional scaling (NMDS) to visualize differences between each combination of summit, aspect, and year. We created a matrix based on the total percent cover of each species in the four 1 m 2 plots within each aspect/summit combination in each year (95 species in four aspects × four summits × three surveys = 48 points).
We calculated the Bray-Curtis distance (Legendre & Legendre, 1998) between these summit-aspect-year totals after taking the square root of the raw percent cover data, which gives more weighting to low-abundance species (McCune & Grace, 2002). For the NMDS and the Bray-Curtis distance, we used PC-ORD 5.0 (McCune & Mefford, 1999). We indicated trajectories in species composition over the three surveys with "successional vectors" represented by arrows.
For Poaceae species, we used the genus level taxa to avoid the noise from uncertainties in taxonomical identification. Therefore, the 31 species of Poaceae were lumped into eight genera for the ordination analysis (Supporting Information Appendix S1). We used a two-dimensional configuration with a final stress of 16 which indicates a quite satisfactory agreement between graph configuration and the similarity matrix (Legendre & Legendre, 1998;McCune & Grace, 2002), and it was significantly different from chance (Montecarlo: 250 runs with randomized matrix, p = 0.002). We report the per- gives p < 0.0005 and 0.002, respectively, as significant (Legendre & Legendre, 1998).
We calculated the relative changes in vegetation variables in 9/11 years as the difference between the second resurvey and the baseline, divided by the baseline value. We used Kruskal Wallis nonparametric tests (Sokal & Rohlf, 1995) to test for differences in the relative changes of vegetation variables among summits and Friedman non -parametric test for dependent samples (Siegel, 1956) to test for differences in percentage cover of life forms categories in the 9/11 years between summits, and between census years across summits.

| Changes in species number and composition
Across botanical families (Supporting Information Appendix S1). The most represented families were Asteraceae with 22% and Poaceae with 20% of the total species pool. Poaceae was also the dominant family in terms of plant cover. We recorded 125 species in the baseline survey, 131 species in resurvey 1 and 114 in resurvey 2. In the 9/11 years of analysis, 14 species were newly recorded, and 13 species were no longer recorded. The number of new species arriving from the baseline to resurvey 2 decreased with altitude (Table 1) and Rubiaceae) were not recorded in resurvey 2, while one family (Cyperaceae) appeared in resurvey 2 for the first time (Table 1; Supporting Information Appendix S1, Data S1).

| Spatial pattern of vegetation, substrates, and temperature
Plant cover, species richness, and Shannon diversity significantly increased with mean soil temperature in the three censuses ( Figure 3a-c), as well as, with minimum average soil temperature.
The three vegetation metrics showed significant differences among aspects within summits, among summits, and surveys: plant cover, species richness, and diversity were highest in the lower summits, with 28% of plant cover, 12 species per m 2 on average (with a maximum of 23 species in 1-ALZ) and diversity index of 2.4, averaged across the three censuses. The highest plant cover tended to occur on eastern aspects, whereas lower species richness and diversity index tended to occur on the southern aspect, the colder side of the summits (Supporting Information Table S2). Census (time of survey) had a significant positive effect on diversity (in the four summits) and richness (in the lower summit; Supporting Information Table S3 for LMM and GLMM models).
Relating vegetation with soil temperature across the three cen-

| Patterns in species richness, plant cover, and composition over 11 years
Plots tended to segregate along two NMDS dimensions based on vascular plant composition for different summits and aspects in three censuses (Figure 4). The two axes represented 78% of the variance in species composition (axis 1, 43% and axis 2, 35%). Summits were clearly separated by soil temperature, mainly mean and minimum (and thus elevation); lower summits (associated with positive axis 1 and negative axis 2) presented higher plant cover and richness in species and life forms, and were more homogenous within aspects (lower coefficient of variation comparing scores between aspects in both axis), than the highest summits, associated with positive axis 2 (3-SIN) and negative axis 1 (4-ISA), where rock cover predominates ( Figure 4). Of the 95 species, 18 species were positively correlated with axis 1, and four with axis 2, four species were negatively correlated with axis 1, and 11 species with axis 2 (Supporting Information Table S4). The main secondary variable driven the ordination were soil temperatures (mean, min and max), which were associated with total plant cover (positive correlation with axis 1 and negative with axis 2) and thus elevation, in the opposite extreme of the ordination, coinciding with higher rock cover (Supporting Information Table S4;  Table 1).
TA B L E 2 Temperature data (°C) obtained from data loggers at 10 cm below the soil surface for the 2009-2012 period. Summit temperature is the monthly average temperature at the four aspects and aspect temperature is the monthly average across the four summits (mean + standard deviation). Daily thermal amplitude is the mean difference between maximum and minimum temperature each day. Seasonal thermal amplitude is the mean difference between growing season and non-growing season temperature. Note that two data loggers (1ALZ north and 3-SIN east) failed, and therefore no data are presented  Figure S4). Non-tussock grasses increased significantly in all summits, mainly in northern aspect, while erects and rosettes increased in the lowest summit, where major changes occurred (with increases of 6.9% of life form cover). In this summit, the increase was mainly due to extant and new species of rosettes (e.g., Gomphrena meyeniana, Hypocheris eremophylla, Geranium sessiliflorum), erect herbs (Olsynium junceum and Calceolaria glacialis), and non-tussock grasses (e.g., Poa kurtzii and Nasella rupestris). Other observed changes were not significant; 2-HUA, tussock grasses tended to increase (by 2.5%), in 3-SIN, shrubs (from the eastern aspect, particularly) tended to decrease (3.6% of decrease), and in 4-ISA cushion plants tended to increase by 2% (e.g., P. convexum; Supporting Information Figure   S4).

| D ISCUSS I ON
We found that temperature, particularly mean and minimum temperatures, and topographic features such as aspect, which modify the thermal conditions, were the primary environmental determinants of species richness, composition, and plant cover at the plot suggesting some upward movements, probably in response to short-term climate variability, instead of long-term trends, both important for modeling vegetation communities . In this study, vegetation variables were analyzed at a small spatial scale (1m 2 plot) which captured 78% (109 species) of the total richness recorded in the four summits (Table 1) the plot smaller scale is more reliable for temporal analysis in cover terms, because it is more accurate and conservative in detecting changes, and reflects the high turnover rate of the plant community, while, the summit broader scale is better focused on the summit species pool for interpreting changes in species richness (Pauli et al., 2015;Zimmer et al, 2014).
The environmental gradient is also reflected in aspect and can be comparable to the elevation gradient in terms of temperature properties determining vegetation community (Sklenář et al., 2016). The highest temperature and diurnal thermal amplitude in the northern aspects, associated with high species richness, diversity and high community turnover over time, may promote the coexistence of species with wider temperature range requirements, mainly species adapted to warmer diurnal temperature, or the effect of degree days on phenology and growth (Oberbauer et al., 2013). Some of these plants, particularly species in the highest summits, may be frost-tolerant species that resist high thermal amplitudes (Sierra-Almeida, Cavieres, & Bravo, 2009;Squeo et al., 1996), such as rosettes (e.g., Draba sp), cushion plants (e.g., Pycnophyllum convexum), and shrubs (e.g., Tetraglochin cristatum and 2). The western aspects, also presented low minimum and mean temperature, but higher species richness than the southern aspects, which may reflect the influence of higher cloud cover in the afternoon and wind, which comes predominantly from the west, particularly in winter, with lots of easterlies in summer, particularly on the three highest summits (Halloy pers. com, Korner 2012). These results support the hypothesis that temperature, differing with elevation and aspect determine vegetation composition and distribution through differences in insolation period.
These differences are reflected in several variables that affect vegetation communities, such as light intensity, soil temperature and moisture, and length of the growing season (Kutiel & Lavee, 1999;Maren, Karki, Prajapati, Yadav, & Shrestha, 2015).
Temperature was the main environmental filter separating species in ordination space, and it was associated with plant and rock cover, in different ways (and with elevation; Figure 4). In general, there was no clear temporal trajectory of vegetation in plots, but in most cases plots of second resurvey tended to assimilate to the baseline in the ordination space. The net trajectory between baseline and second resurvey (not shown), indicates minor changes, mainly in lower summits. This may highlight that in these mountains and at this temporal scale, vegetation responds to short-term environmental variability, reflected in temperature and precipitation (Supporting Information Figure S1), and probably the distance index  (Graae et al., 2018). 3-SIN east, similar than 4-ISA east, presented higher plant cover than other aspects of the summit, e.g. shrubs, cushion plants and also higher organic matter cover. The highest summit is characterized by few typical plant species, such as Aschersionodoxa cachensis, Arenaria rivularis, Valeriana pycnantha. In other mountains, these species are associated to crioturbed soils and rocky areas (Cano et al., 2011), thus, being highly vulnerable to global warming, so they can be useful as indicator species. Our study, as many others (aforementioned), focused its efforts in trying to understand vegetation community dynamics by inferring their climate-driven distributions, to predict future communities trajectories in pursuit of its conservation, however, prediction may be a real challenge linked to novel climates (more vulnerable in tropics and subtropics) that may promote "noanalog communities", that is, nonexistent communities under current conditions (Williams & Jackson, 2007).
Our results suggest that vegetation dynamics (mainly a general increase in diversity and promotion of non-tussock grasses and other small herbs) may be responding to short-term variability in temperature and humidity combined with local factors (from different sources) in each summit. The low plant cover at the highest summit seemed to respond quite quickly to temperature changes F I G U R E 5 Box plot of relative change in (a) plant cover, (b) species richness, and (c) Shannon diversity index in the four summits in 11 years (Resurvey 2-baseline/baseline). **p < 0.001, ***p < 0.0001 (doubling plant cover). At the lowest summit, an increase in the three vegetation metrics is observed, as in other mountains in the world Venn, Pickering, & Green, 2014), with abrupt changes in species richness in the first 5 years, probably responding to the short-term variability in temperature and humidity combined with anthropogenic variables, such as grazing or tourism (Steinbauer et al., 2018;Yager et al., 2008). Some studies evidenced that decreasing in domestic cattle, driven by changes in land use (promoted by socioeconomic factors; e.g., mining and tourism) may favor native herbivores populations (Barros, Monz, & Pickering, 2015;Izquierdo et al., In press;Steinbauer et al., 2018). This grazing transition may affect vegetation communities in the entire elevation gradient, a hypothesis to be tested in the future. At the second summit, 2-HUA, vegetation changes (decrease in species richness and increase in plant cover) may be driven by the combination of climate variability and local biological factors such as competition, for example; tussock grasses may be limiting plant colonization, by competing for local resources (Halloy, 1985a). The particular hypsography of these mountains, with a plateau at middle elevation, may also influence the linear relationship expected between elevation and species richness, observed at summit scale (Elsen & Tingley, 2015). In 3-SIN, changes in vegetation seems to respond to the climate variability combined with local abiotic effects, such as sun exposition and wind, having an effect over predominant plant cover (some shrubs located in the east aspect dried up), but after the 10 years of analysis a recovery of the three vegetation variables is observed (Table 1; Supporting Information Figure S1). Diversity index increases significantly over time across the four summits, reflecting the short-term temperature and precipitation variability, as well (lower diversity at the baseline coincided with dry and warm period; Supporting Information Figure S1), and as it encompasses both variables, species cover and richness, could be a good indicator of vegetation community dynamics.
These inter-annual variability cycles may help explain temperatures, precipitation, lake levels, and hence vegetation variations over that period, observed in other Andean systems (Carilla et al., 2013).
There are worries that plant community composition and perhaps community functioning will be negatively affected by global change in the future. Our study showed that over an 11-year period, total plant cover, species richness and diversity were positively associated with temperature. Some signals highlighted in this study suggest that temperature is changing at different temporal scales (e.g., Increase of 0.5°C in minimum temperature between 2009 and 2015) and vegetation dynamics may be responding to these changes. Although our data may be not enough to support the observed relationship, we suggest that trends in vegetation change respond to short-term in temperature and precipitation variability, which is strongly influenced by ENSO, as in the central Andes (Christie et al., 2009;Morales et al., 2015), thus, it may take much longer than one decade for trends to became significant, due to the natural time lags in plant community responses to shifts in rainfall and temperature. This hypothesis should be tested with longer term monitoring studies; however, our data provide an important baseline with which to document changes in these sensitive alpine plant communities as global climate continues to change.

| CON CLUS ION
Soil temperature and topography (aspect and elevation) were the primary environmental determinants of vegetation dynamics at these spatial and temporal scales, with soil temperature decreasing significantly with elevation and showing a positive trend over the six years of records. Greater plant cover, species richness, and diversity tended to occur at the lower summits and in the warmer northern and eastern aspects, while the lower species richness and diversity occurred in the highest summit and in the colder southern aspects.
After 9/11 years, we observed that diversity index increased in the four summits, species richness increased in the lower summit and plant cover increased in the highest summit, with a general increase in non-tussock grasses (at plot scale). Our results suggest that vegetation dynamics may be associated with short-term variability of soil temperature and precipitation, influenced by ENSO at this elevation. Tucumán partially covered researchers salaries and is the main research facilty. Dra. Cecilia Blundo assisted with statistical analysis.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R S CO NTR I B UTI O N
JC, AG, SC, and SH conceived of and designed the study, they also collected data and conducted initial data analyses. SC identified species, JC, FC, and AM contributed to statistical analysis. JC wrote the first draft of the manuscript, all authors contributed substantially to revisions.