Vegetation dynamics at the upper elevational limit of vascular plants in Himalaya.

A rapid warming in Himalayas is predicted to increase plant upper distributional limits, vegetation cover and abundance of species adapted to warmer climate. We explored these predictions in NW Himalayas, by revisiting uppermost plant populations after ten years (2003-2013), detailed monitoring of vegetation changes in permanent plots (2009-2012), and age analysis of plants growing from 5500 to 6150 m. Plant traits and microclimate variables were recorded to explain observed vegetation changes. The elevation limits of several species shifted up to 6150 m, about 150 vertical meters above the limit of continuous plant distribution. The plant age analysis corroborated the hypothesis of warming-driven uphill migration. However, the impact of warming interacts with increasing precipitation and physical disturbance. The extreme summer snowfall event in 2010 is likely responsible for substantial decrease in plant cover in both alpine and subnival vegetation and compositional shift towards species preferring wetter habitats. Simultaneous increase in summer temperature and precipitation caused rapid snow melt and, coupled with frequent night frosts, generated multiple freeze-thaw cycles detrimental to subnival plants. Our results suggest that plant species responses to ongoing climate change will not be unidirectional upward range shifts but rather multi-dimensional, species-specific and spatially variable.

Cold biomes at high latitudes or high elevations are among the most dynamic systems on earth due to rising temperatures12. Arcto-pan class="Species">alpine biomes warm faster than global average2, leading to increase in relative abundance of sclass="Chemical">pecies adaclass="Chemical">pted to higher temclass="Chemical">perature3, a class="Chemical">process known as thermoclass="Chemical">philization4. Consequently, sclass="Chemical">pecies which are narrowly sclass="Chemical">pecialized to cold habitats migrate uclass="Chemical">pwards or exclass="Chemical">perience local extinctions5. Ongoing distribution shifts67 and changes in class="Chemical">plant cover, including the colonization of newly deglaciated areas8 are widely acknowledged9. Raclass="Chemical">pid warming is also changing class="Chemical">plant class="Chemical">phenology10, growth11 and class="Chemical">productivity12.

Most of our knowledge about the recent vegetation changes in arcto-class="Species">alpine systems comes from humid regions57131415. While ongoing warming is likely to be the key driver of vegetation change in these cold and humid regions3, altered class="Chemical">preciclass="Chemical">pitation regime and its subsequent effect on soil class="Chemical">pan class="Chemical">water balance may be even more important in dry regions161718. Changes in precipitation regime are often accompanied by an increase in extreme events2 to which arid regions seem to be more vulnerable19. For instance, extreme precipitation in arid mountains can cause floods and landslides on poorly vegetated slopes20 and the increased soil moisture together with associated frost heave can uproot plants21. Moreover, plant-plant interactions are likely to differ between humid and arid ecosystems2223. Because of the specific precipitation regime, arid regions may not undergo the processes of thermophilization generally described in cold biomes35. The lack of knowledge about vegetation dynamics in response to climate change in cold and dry biomes is striking especially in dry continental Asia, where vascular plants grow in the highest elevations on earth.

Ladakh, an arid mountainous region in the NW Himalayas, represents a unique region for studying vegetation dynamics and plant responses to climate change (Fig. 1 and Supplementary Fig. S1). Due to the class="Disease">dry continental climate, the mountains of Ladakh are often unglaciated uclass="Chemical">p to 6200–6400 m. The unglaciated terrain together with dry climate allows class="Chemical">plant sclass="Chemical">pecies to growth in extreme elevation even above 6000 m24. Still, a relatively large and class="Chemical">potentially colonizable area is found above the current uclass="Chemical">pclass="Chemical">per distributional limit of class="Chemical">pan class="Species">vascular plants within this region. The mean temperatures in Ladakh have increased by 1.7 °C during the last century (Supplementary Fig. S2), and the warming accelerated in the 1990s25 leading to rapid glacier retreat26. Beside warming, Ladakh also experienced increasing precipitation, accompanied by more frequent extreme snowfall and storm events causing severe floods272829.

Figure 1

Study area in Eastern Ladakh, NW Himalayas.

(A) Subnival belt around Chamser Kangri peak (6660 m a.s.l.). (B) Daily course of air, surface and soil temperatures in the peak of the vegetation season at 6150 m; the key for the survival of plants here is the rooting zone temperature during the growing season. From July to August, plants are able to withstand freezing air temperatures almost every night for 5–10 hours, but they never occur where the rooting zone temperature falls deeply below zero during the period of active growth. (C) The oldest plant (22 years) at 6150 m was an endemic species Ladakiella klimesii. (D) Its age was derived from ring counts on a 0.7 mm radius of root collar.

class="Chemical">Plant resclass="Chemical">ponses to environmental changes, including climate warming and increased class="Chemical">preciclass="Chemical">pitation, can be class="Chemical">predicted from their functional traits30. class="Chemical">pan class="Chemical">Plant tissue nutrient content is one of the indicators of plant growth limitation at high elevations3132. High nutrient but also carbohydrate concentrations commonly found in alpine plants indicate the plant’s inability to use the absorbed resources for growth (sink limitation) due to low temperatures restriction on enzymatic processes. Climate warming in the cold NW Himalayas can reduce sink limitation and increase relative abundance of faster growing taller species with optima at lower elevation steppes and semi-deserts. Contrary, increased precipitation and soil disturbance can support species from wet habitats such as springs and snowbeds, which are better adapted to repeated soil freezing and thawing. Hence, altered precipitation regime in the alpine and subnival zones of dry Himalayas can lead to the process of hygrophilization rather than thermophilization.

To explore upward plant migration and susceptibility of the cold and dry region of Ladakh to climate change, we (1) revisit an area extensively surveyed 10 years ago to find possible outpost populations above the previously recorded altitudinal upper limit; (2) re-survey permanent plots established in the class="Species">alpine and the subnival belt after 4 years characterized by unclass="Chemical">precedented warming and extreme summer class="Chemical">preciclass="Chemical">pitation; (3) characterize class="Chemical">plant morclass="Chemical">phological and ecoclass="Chemical">physiological adaclass="Chemical">ptation to evaluate the biological mechanisms of vegetation changes; (4) class="Chemical">perform anatomical determination of class="Chemical">plant age and growth rate. class="Chemical">pan class="Chemical">Plant age provides crucial proof of recent colonization, and recent advances in anatomical assessment of distinct annual rings enables age determination, even in plants with a root collar thinner than 1 mm33. Plants from cold unproductive environments tend to invest into longevity at the expense of fast growth, thus older plants are usually found at higher elevations34. However, warming alleviates abiotic stress and reduces the need to invest into longevity, which in turn can trigger uphill plant migration43536. If this assumption is true, we can expect a negative relationship between elevation and plant age at the higher end of the elevational gradient. We hypothesize that the ongoing climate warming has promoted uphill plant migration in this cold and dry region, so that absolute upper limits, species richness and cover of vascular plants increased in the past decade. If so, the species increasing in cover in the subnival belt should have syndromes of traits characteristic for lower elevation habitats. On the other hand, the evidences of thermophilization might be overridden by increased frequency of extreme snowfall events producing additional disturbance resulting in reduced species richness and cover.

Results

Climate changes

Overall, the period studied was characterized by warmer temperatures and extreme snowfall events (Supplementary Fig. S2). Long-term meteorological (1901–2014) records of temperature and precipitation in Ladakh show positive trends in mean winter, spring and fall temperature, with pronounced warming since the mid 1990s. The records also show a substantial increase in total precipitation and summer temperature since the mid 1990s. In-situ climate records from our study area for the period 2008–2014 confirm the general trend. There was a significant increase in mean daily maximum T for June (increase of 2.43 °C per year), July (+1.42 °C per year), August (+0.95 °C per year), and September (+1.02 °C per year; Fig. 2B and Supplementary Table S1). Similar trends were recorded for daily mean T in the summer months, expect for August (Fig. 2A). The mean summer T (June to August) increased by 0.59 °C per year. The diurnal T range also increased (Fig. 2C). No major changes were recorded for winter month T, except for a cooling in February (−0.38 °C per year). The growing season length increased from 56 days in 2009 to 83 in 2010, 79 in 2012, and 75 in 2013. The summer 2010 was snowiest in past eight years (Fig. 2D).

Figure 2

Climate conditions at the study area between January 2009 and September 2014.

Mean and maximum monthly temperatures (T), the number of days with snow cover (snowfall total) and difference betweem maximum and minimum daily T (range) are shown.

Plant distributional limits and their upward shifts

During the first survey in 2001–2003, nine species reached 6000 m and above it: class="Species">Aphragmus oxycarpus, class="Chemical">pan class="Species">Draba altaica and Draba oreades were found highest at 6000 m, Ladakiella klimesii, Stellaria decumbens and Saussurea glacialis at 6010 m, Poa attenuata and Waldheimia tridactylites at 6030 m, and Saussurea inversa at 6060 m (Supplementary Figs S1 and S6). During the second survey in 2010–2013, there was one new species extending its range above 6000 m (Saussurea gnaphalodes). In total, the five species were found at higher elevations, the highest one being at 6150 m. These shifts extended the vertical range by 140 m in Ladakiella klimesii, 150 m in Poa attenuata and Draba altaica, 120 m in Waldheimia tridactylites, and 180 m in Saussurea gnaphalodes.

Changes in diversity and composition in permanent plots

Total plant cover decreased both in the class="Species">alpine (GLMM: chi-square = 94.4, class="Chemical">pan class="Chemical">P < 0.001,) and subnival (GLMM: chi-square = 35.7, P < 0.001, Fig. 3) vegetation between 2009 and 2012.

Figure 3

Changes in plant cover and species richness in the alpine and subnival vegetation between August 2009 and 2012.

The dashed line represents the null hypothesis of no change in cover or richness between surveys, and the solid lines represent a linear regression fitted to the alpine and subnival data.

Species richness showed no significant change in class="Species">alpine vegetation (GLMM: chi-square = 3.607, class="Chemical">pan class="Chemical">P = 0.998, Supplementary Fig. S7), and marginally decreased in the subnival vegetation (GLMM: chi-square = 3.607, P = 0.067). The species composition changed between resurveys in both the alpine (standardized RDA, F = 0.2, P = 0.002, Supplementary Fig. S8) and subnival (the only significant statistic was non-standardized RDA, F = 1.4, P = 0.008) vegetation. In the subnival belt this was mostly due to changes in abundance and not floristic turnover.

The mean total cover per plot decreased from 46 to 31% in the class="Species">alpine (mean relative change 35%), and from 17 to 10% (mean relative change 48%) in the subnival zone between the two samclass="Chemical">pling class="Chemical">periods. 20 sclass="Chemical">pecies (6 class="Chemical">pan class="Species">alpine and 14 subnival, Supplementary Table S2) significantly decreased in plant cover (mean relative change 40.1%), no species significantly increased in their cover. From 41 species recorded during the 2009 survey, three species were not found, and five new species appeared during the 2012 survey; all of these losses and gains happened within the alpine vegetation.

Plant age distribution along elevational gradient

The class="Species">alpine belt hosts sclass="Chemical">pecies with a higher mean age while class="Chemical">plants are usually younger and grow slower in the subnival belt (Fig. 4). The mean age of class="Chemical">pan class="Species">alpine and subnival plants was 20.5 and 11.4 years, respectively. The oldest plant individuals were Potentilla pamirica (70 years) and Arenaria bryophylla (53 years) in the alpine zone at around 5700 m.

Figure 4

Difference in plant age and annual growth increment between alpine and subnival species.

Non-overlapping notches in the boxes indicate the significance of between-group differences.

The most individuals found above 6000 m were established within the last decade (e.g. class="Species">D. altaica, one sclass="Chemical">pecimen 8 year-old; class="Chemical">pan class="Species">W. tridactylites, three specimens 4, 4 and 9 years-old; S. inversa, three specimens 2, 6 and 6 years-old). The oldest plant individual at the highest site 6150 m was a 22 year-old endemic L. klimesii (Fig. 1D), and the two other specimens of this species were 8- and 16 years-old, indicating that some of outpost populations above 6000 m were established at least two decades ago.

Plant functional traits

class="Chemical">Plant trait values differ significantly between class="Chemical">pan class="Species">alpine and subnival vegetation (Fig. 5, Table 1). 54% of the alpine species are clonal species with adventitious roots (spreading splitters) whereas subnival vegetation is dominated by small non-clonal cushion species (64% of the subnival species) with a main tap root producing belowground branches bearing aboveground leaf rosettes (spreading integrators). Subnival species are in average 21% shorter, 44% younger, grow 25% slower, have 48% lower total dry biomass and 36.9% smaller rooting depth than alpine species (Table 1, Supplementary Fig. S9). Subnival species have significantly higher nitrogen concentration in their leaves (6.9%) and roots (30.6% respectively) as well as higher phosphorus concentration in roots (27.9%), but lower LDMC (25.3%), leaf carbon content (2.9%) and consequently lower leaf C:N ratio (8.6%) compared to alpine species. In addition, subnival species have 9.2% more soluble carbohydrates and 76.1% more fructans in their roots but do not differ in starch content compared to alpine species. Finally, subnival communities are composed of plant species with an affinity for lower soil moisture and soil salinity compared to species from alpine communities (Table 1).

Figure 5

Redundancy analysis biplots (RDA) of community weighted mean values of plant functional traits (gray arrows) and habitat indicator values (green arrows) in relation to elevation, vegetation belt (alpine, subnival) and total vegetation cover (red color).

The angles between arrows indicate correlations between variables (for explanation of trait and indicator values see Table 1). The environmental data explained 25.5% variability.

Table 1

Comparison of plant traits and indicator values between the alpine and subnival vegetation and the shift in trait community-weighted means between two sampling periods from 2009 to 2012.

Traits and habitat preferences	Abbreviation	Alpine09	Subnival09	P	Alpine_diff_2009–2012	P	Subnival_diff_2009–2012	P
Plant height (cm)	Height	14.00	11.00	↓**	0.204		1.058	↑*
Leaf nitrogen content (%)	LNC	2.26	2.42	↑*	−0.004		−0.007
Leaf phosphorus content (%)	LPC	0.18	0.17		0.001		0.002
Leaf carbon content (%)	LCC	42.84	41.59	↓***	0.169	↑***	0.414	↑*
Leaf carbon/nitrogen ratio	C:N	20.13	18.40	↓***	0.140		0.301
Leaf δC13carbon (%)	C13	−26.21	−26.05		0.001		−0.008
Leaf dry matter content (mg/g)	LDMC	176.63	131.84	↓**	7.140	↑*	1.608
Stem dry matter content (mg/g)	STDMC	326.86	329.70		−14.39	↓**	52.84	↑a
Root nitrogen content (%)	RNC	1.16	1.52	↑***	−0.032	↓***	−0.056	↓**
Root phosphorus content (%)	RPC	0.10	0.13	↑***	−0.003	↓***	−0.005	↓a
Starch content (%)	Starch	2.48	2.00		0.020		−0.043
Fructan content (%)	Fructan	2.72	4.78	↑***	−0.055		−0.492	↓*
Soluble carbohydrates (%)	SolubleSugar	4.02	4.38	↑*	−0.100		−0.047
Total non-structural carbohydrates	TNC	9.19	11.17	↑**	−0.149		−0.581
Seed weight (g)	SeedWeig	0.04	0.03		−0.003	↓*	0.001
Root depth (cm)	RootD	14.18	8.94	↓***	0.439	↑*	−0.434	↓*
Number of aboveground tillers	ShootNo	18.02	15.83		−0.386		1.105
Number of flowers	FlowerNo	7.36	7.24		−0.340		0.466
Total dry biomass (g)	TDB	2.45	1.25	↓*	−0.018		−0.060
Clonal plants with rhizomes (%)	Clonal	54	36	↓**	2.7	↑*	6.4	↑*
Non-clonal with tap-roots (%)	Nonclonal	46	64	↑**	−2.7	↓*	−6.4	↓*
Nonspreading integrators (%)	IntegNS	34	34		−0.6		−5.6	↓*
Nonspreading splitters (%)	SplitNS	34	31		0.8		5.5	↑a
Spreading splitters (%)	SplitS	20	5	↓**	1.9	↑a	0.9
Spreading integrators (%)	IntegS	12	30	↑**	−0.022	↓*	−0.9
Indicator value_moisture	Moisture	2.84	2.68	↓***	0.049	↑***	0.012
Indicator value_stability	Stability	2.46	2.42		0.002		0.011
Indicator value_salinity	Salinity	0.33	0.12	↓***	0.013	↑a	0.003
Indicator value_nutrient	Nutrient	0.98	0.98		0.007		0.005
Indicator value _shading	Shading	0.09	0.09		−0.003		−0.006
Elevational minima	Altmin	5056.6	5311.1	↑***	8.050		−29.9	↓*
Elevational optima	Altopt	5246.2	5534.4	↑***	−1.033		−24.2	↓*
Elevational maxima	Altmax	5708.9	5912.5	↑***	−1.469		−14.2	↓a
Elevational range	Altrange	652.3	601.4	↓*	−9.520	↓a	15.7

An upward or downward pointing arrow indicates an increase or decrease from the alpine to subnival, or temporal shift. Type I error estimates aP < 0.1, *P < 0.05, **P < 0.01, ***P < 0.001.

Shift in plant strategies in response to changing environmental conditions

Compositional changes observed in class="Species">alpine and subnival vegetation from 2009 to 2012 were not ecologically random as evidenced from a shift in community-weighted mean trait values (Table 1). These changes in the CWM corresclass="Chemical">pond to the decrease in abundance of the sclass="Chemical">pecies with syndromes of highest elevation described in the class="Chemical">previous section. The class="Chemical">proclass="Chemical">portion of non-clonal sclass="Chemical">pecies with a main taclass="Chemical">p root decreased between samclass="Chemical">pling class="Chemical">periods (5.9% in class="Chemical">pan class="Species">alpine and 10.1% in subnival habitats), while that of clonal plants with adventitious roots increased (4.9 and 17.5%, respectively). This was accompanied by a significant decrease in root nitrogen (2.7 and 3.9%) and phosphorus (3.5 and 3.6%) concentrations. There was also a significant decrease in mean seed weight (7.1%) and increase in LDMC (4.4%) in alpine vegetation and increase in plant height (9.6%) in the subnival zone. Indicator values for soil moisture and salinity changed significantly between sampling periods in the alpine vegetation: species with an affinity towards wet and saline habitats (e.g. Potentilla gelida, Halerpestes sarmentosa, Pegaeophyton scapiflorum, Koenigia islandica) increased in cover relative to species from dry habitats (e.g. D. altaica, P. pamirica, S. gnaphalodes, D. oreades, Supplementary Fig. S8). In subnival vegetation, species with optima at higher elevations decreased in cover.

An alternative species-based approach confirmed that soil moisture indicator value is the best predictor of changes in species’ cover in the class="Species">alpine vegetation (Fig. 6). Sclass="Chemical">pecies class="Chemical">preferring drier surfaces were most class="Chemical">prone to decline in their cover. The cover values of this grouclass="Chemical">p was further class="Chemical">predicted by leaf δclass="Chemical">pan class="Chemical">13C as a proxy of WUE, root phosphorus concentration as a proxy of growth rate, and surface stability as a proxy of disturbance. Species preferring drier and stable surfaces with lower δ13C and higher RPC were more susceptible to decline in cover than species with higher δ13C and faster growth preferring wetter and naturally more disturbed habitats. In the subnival vegetation the only significant trait predicting changes in their abundance was space occupancy strategy (Supplementary Fig. S10). Species capable of clonal multiplication belonging to category of spreading splitters (Carex sagaensis) and spreading integrators (Oxytropis platysema, W. tridactylites, Desideria pumila) increased in cover, while non-clonal species with main tap roots belonging to category of non-spreading intergrators (A. oxycarpus, D. oreades, S. hypsipeta, L. klimesii) or species with limited vegetative reproduction belonging to category of non-spreading splitters (P. attenuata, Festuca non-coelestis) were more prone to decrease in cover.

Figure 6

Life-history traits and ecological indicator values predicting species abundance changes in the alpine vegetation using the conditional inference tree.

In each split of the tree, all species predictors are tested and the one that best discriminates between species is selected. The procedure goes on until no predictor significantly discriminates between species. The response variable (y) is the ordination score on the partial ordination first axis determined by a temporal change (boxes). Negative values indicate species declining over time, and positive values indicate species expanding. Each split of the tree is described by the trait or indicator value used at the split (ovals), the permutation-based significance of the split (P-value) (ovals) and the treshold values at which the split occurs (for explanation of trait and indicator values see Table 1). The number of species (n) is given at each terminal node (box).

Discussion

The combination of site re-visitation, age determination and trait characterization did not provide evidence for thermophilization at the world upper limit for pan class="Species">vascular class="Chemical">plant life after an exceclass="Chemical">ptionally warm decade. However, our results from the arid Himalayas do indicate the sensitivity of vegetation to the ongoing climate change. Unlike in the humid and cold high-altitude regions35, the combination of warming and increased frequency of extreme snowfall events in the arid Himalayas27 causes increasing disturbance events rather than class="Chemical">purely ameliorating the level of stress. This results in a more comclass="Chemical">plex class="Chemical">pattern of vegetation changes than otherwise exclass="Chemical">pected by the uclass="Chemical">pward migration and thermofilization hyclass="Chemical">potheses737.

Plant uphill migration and adaptive strategies

Accelerated warming in Eastern Ladakh over the past decade may have promoted the establishment of new populations at very high elevation, similarly as in other mountains5131438. We discovered several outpost populations during the re-exploration of the Tso Moriri area above the upper vegetation limit determined 10 years before this study. While plant age analysis indicates that the highest outpost of L. klimesii was overlooked in the previous survey, L. klimesii was the only plant species of the highest outpost that was over 10 years old. The other species have established within the warmer last decade. Also the younger age of the subnival species compared to the pan class="Species">alpine sclass="Chemical">pecies indicates recent establishments which coincide with the recent accelerated warming.

The observed upward plant migration was not ecologically random process as evidenced from plant life-history analyses37. Subnival species that significantly extended their range by uphill shift evolved similar anatomical and ecophysiological adaptations despite differences in geographical origin, vertical distribution, thermal optima and phylogeny24. Common features to all of them are the absence of fibers and the presence of thick-walled, and small and short vessels (10–30 μm in diameter and 40–150 μm in length) which protect the plants against drought- and frost-induced cavitation and class="Disease">embolism39. They are cushion-forming class="Chemical">perennials with secondary growth, annual rings of semi-ring class="Chemical">porosity and very small annual increments (ring width 0.03–0.05 mm) reflecting the harsh environmental conditions. These subnival sclass="Chemical">pecies contain more nutrients and soluble class="Chemical">pan class="Chemical">carbohydrates and are capable of maximizing the rate of photosynthesis during short periods of favorable conditions, due to high foliar nitrogen content31 (Table 1 and Supplementary Fig. S9). Higher concentrations of simple sugars and sugar alcohols are thought to serve as osmoprotectants when drought or frost cause water deficit in the plant40. These adaptations are crucial, because drought on sunny days and frost in cold nights are very common. Beside proper ecophysiological adaptations, the migrating species are anemochorous37, such as Saussurea spp. from the Asteraceae family, which produces a large quantity of small wind-dispersed seeds with an attached feathery pappus, or Draba spp. and A. oxycarpus from the Brassicaceae family which have extremely light seeds24.

While proper organismal adaptations to freezing temperatures and drought are necessary prerequisite for upward plant migration41, it remains question if there are some shared physiological constraints beyond which plants cannot move31. There is increasing evidence that the coldest places on Earth with angiosperm plant life are determined by a growing season (defined by a daily mean >0 °C in the rooting zone) of at least 50 days, with seasonal mean root zone temperature around 2–3 °C32. This could be considered as the universal minimum for persistent angiosperm life32, although further evidence is needed. Klimes and Dolezal42 found the upper limit of plant life at 6030 m on the western slope of the Chalung mountains in 2002 (10 km north of our study site), formed by class="Species">Poa attenuata and class="Chemical">pan class="Species">Waldheimia tridactylites; in-situ measurements over 2002–2003 showed a 47 day growing season with mean root zone (−5 cm) T of only 1.7 °C. Unfortunately, the existence of plant life at this precisely georeferenced site was not confirmed during the revisitation in August 2013 (Dolezal, unpublished data), indicating the transient character of outpost populations in extreme elevations. Interestingly, what is common for the uppermost plant populations for which in-situ measured climate data are available (Mt Chalung 6030 m and Mt Shukule 6150 m), despite the difference in absolute elevations and seasonal mean T (1.7 °C vs 2.8 °C), is at least four continuous weeks without frost in the rooting zone. This suggests that the temperature of the rooting zone during the growing season is more important for vascular plant life at its upper limit than other commonly used climate measures. The four-year records from 6150 m indicate that freezing of the soil is more critical for vascular plants than frost above the ground3132. While the air temperature during the short growing season drops below zero every night for 5–10 hours, repeatedly to about −5 °C (Supplementary Fig. S11), the soil temperature remains above zero (Fig. 1B). Plants more easily cope with freezing air than soil because of the differences between shoot epidermis and root rhizodermis; the root rhizodermis represents a less sufficient barrier against ice-nucleation upon contact with extrinsic soil ice masses, especially in juvenile plants43. Most of the upper subnival species in Ladakh have a thick bark, large cortex and a xylem with much living parenchyma with high soluble carbohydrate concentrations (Fig. 1 and Supplementary Fig. S9) as protection against night frost during the short summer season.

Thermophilization or hygrophilization?

In addition to increasing temperatures during the past decade, the other aspect of climate change in the region, the extreme snow fall events27 and the subsequent frost heave (Supplementary Fig. S5), when the soil is saturated with class="Chemical">water, could class="Chemical">pose a serious challenge to class="Chemical">pan class="Species">vascular plants adapted to arid soils of cold alpine deserts3143. The extreme snow fall events during the study period coincide with the decrease in vegetation cover (increase in bare ground) by roughly one third, for both alpine and subnival communities. Such short term climate events can have long-lasting effects44. While increased temperature could promote plant regeneration, the very slow growth rate common to species in this system, is unlikely to compensate for the increased frequency and intensity of the disturbance.

Trait analyses of the re-surveyed class="Species">alpine and subnival vegetation did not class="Chemical">provide suclass="Chemical">pclass="Chemical">port for thermoclass="Chemical">philization class="Chemical">processes. Evidence for thermoclass="Chemical">philization would have been indicated by an increase in the abundance of sclass="Chemical">pecies with trait syndromes and elevational oclass="Chemical">ptima common for thermoclass="Chemical">philous or xeric sclass="Chemical">pecies from habitats of lower elevations, such as steclass="Chemical">pclass="Chemical">pes and deserts45. Instead, the vegetation comclass="Chemical">position in class="Chemical">permanent class="Chemical">plots in the class="Chemical">pan class="Species">alpine belt shifted towards clonal species preferring wetter habitats (Table 1). We observed an increase in hygrophytes, a phylogenetically diverse group of plants growing in permanently or temporarily wet soils that are sometimes waterlogged and frozen almost every night on the surface during summer. These were not only species with high affinity to soil moisture, but also those with higher leaf δC13 (higher WUE) and lower P concentration. This indicates that increasing moisture is not the only cause of the observed vegetation change but also increased disturbance due to repeated freeze-thaw cycles, which plants with higher WUE can better cope with. The higher P content is characteristic of slow growing alpine specialists243132. Most prone to decline are therefore species with lower WUE and slow growth not capable of clonal multiplication occupying drier and stable surfaces, i.e. the common species at the highest elevations in arid NW Himalayas.

Hence, the climate change in this arid region, involving warming and increasing precipitation, may not lead to the process of thermophilization as observed in cold and humid regions3, but rather to hygrophilization by increasing species from azonal wet habitats having otherwise small spatial extent in these arid mountains. We attribute the vegetation shift to the adverse effects of increasing daily maximum temperatures, which in combination with snowfall during growing season, creates freeze-thaw cycles. Unprecedented warming recorded in the period 2008–2012, with an increase in the mean daily maximum T by 2.4 °C and 1.4 °C per year for June and July, respectively, causes faster melting of snow instead of slow sublimation. The warming is accompanied by increasing diurnal fluctuations (Fig. 2C) as minimum daily T did not change significantly over time and hence the night frosts are still present. All this contributes to the recurrence of freeze-thaw actions (Supplementary Fig. S6). It is therefore the simultaneous increase in summer temperature and precipitation that can be detrimental for populations at the upper limit of plant distribution, in particular those with a single main root and adapted to dry and stable substrates. The abrupt damage to the pan class="Species">alpine and subnival vegetation, which took class="Chemical">place after snowfall summer 2010, will hardly be comclass="Chemical">pensated soon because of the extremely slow growth in extreme elevations. Therefore, the imclass="Chemical">pact of the increasing frequency of extreme class="Chemical">preciclass="Chemical">pitation events may, in the longer term, suclass="Chemical">pclass="Chemical">press the vegetation due to its immediate effect, while the class="Chemical">positive effects of warming are relatively slow.

Resilience to climate change

Dvorsky et al. 24 found that the generalist species with wide vertical ranges and optima in the class="Species">alpine and steclass="Chemical">pclass="Chemical">pe zone are more abundant in subnival flora of Ladakh than high-elevation sclass="Chemical">pecialists (see also Suclass="Chemical">pclass="Chemical">plementary Figs S1 and S6). Subnival sclass="Chemical">pecialists with narrow elevational ranges reclass="Chemical">present 42% of the flora. Also among nine sclass="Chemical">pecies forming uclass="Chemical">pclass="Chemical">permost class="Chemical">poclass="Chemical">pulations above 6000 m are five generalists (class="Chemical">pan class="Species">Aphragmus oxycarpus, Saussurea glacialis, Waldheimia tridactylites, Poa attenuata, and Draba altaica) having optima at around 5100 m elevation and vertical ranges over 2000 meters long. The fact that subnival flora is composed of both specialist and generalist species suggests its relatively high resilience to climate change including both warming and increased frost disturbances in this region. Continued warming could potentially lead to thermophilization and increased abundance of steppic and alpine plants causing subnival specialists to decline at the alpine-nival ecotone. This was recently seen in wet lower mountains535 where subnival specialists already reached to mountain tops and hence cannot escape competition from alpine species by uphill migration. This is an unlikely scenerio in the arid NW Himalayas because the vast unglaciated areas between 6000–6400 m are available for range extension of the subnival species244245. Complete extinction of subnival flora due to increasing precipitation and frost-disturbance is also unlikely given the high proportion of generalists among subnival taxa with their source populations at lower elevations capable of recolonizing sink habitats above.

Conclusion

Knowledge of vegetation dynamics at the altitudinal limit of class="Species">vascular plants is crucial for understanding class="Chemical">plant migration class="Chemical">processes, class="Chemical">physiological limits and diversity changes4515. Here, we exclass="Chemical">plore the vegetation dynamics in the extreme altitudes of arid Ladakh, a raclass="Chemical">pidly warming region in the NW Himalayas, to test the ‘uclass="Chemical">pward migration and thermoclass="Chemical">philization’ hyclass="Chemical">potheses, which class="Chemical">predict an increase in absolute sclass="Chemical">pecies distributional limits, vegetation cover and abundance of sclass="Chemical">pecies adaclass="Chemical">pted to warmer climate. We re-surveyed outclass="Chemical">post class="Chemical">poclass="Chemical">pulations after 10 years and class="Chemical">permanent class="Chemical">plots after 4 years from an area between 5500 and 6150 m.a.s.l. In addition, we determined the age and various functional and ecoclass="Chemical">physiological traits related to adaclass="Chemical">ptations of class="Chemical">plants to extreme altitude. We discovered several outclass="Chemical">post class="Chemical">poclass="Chemical">pulations during the re-exclass="Chemical">ploration of the area above the uclass="Chemical">pclass="Chemical">per vegetation limit determined 10 years before this study, suclass="Chemical">pclass="Chemical">porting the uclass="Chemical">pward migration hyclass="Chemical">pothesis, with five subnival sclass="Chemical">pecies extending their vertical range by 120–180 m, reaching an elevation of 6150 m. The class="Chemical">plant age analysis corroborated the hyclass="Chemical">pothesis of warming-driven uclass="Chemical">phill shift. However, the force of warming class="Chemical">pushing class="Chemical">plants uclass="Chemical">pwards has been recently offset by extreme snowfall events as those in summer 2010. Vegetation resurvey of class="Chemical">permanent class="Chemical">plots revealed a substantial decrease in class="Chemical">plant cover and comclass="Chemical">positional shift towards sclass="Chemical">pecies class="Chemical">preferring wetter habitats. Declining sclass="Chemical">pecies have lower class="Chemical">pan class="Chemical">water-use efficiency, slow growth and main taproot, and are typical for the dry and cold NW Himalayas. Contrary, species whose cover increased were mostly clonal hygrophytes growing typically in wet soils. These vegetation changes are likely due to increasing soil moisture content and physical disturbance. Simultaneous increase in summer temperature and precipitation which causes rapid snow melt and, coupled with frequent night frosts, generates frost heaves detrimental to subnival plants. These relatively short-term disturbance episodes can have long-lasting effects because slow growth rate of high-altitude plants is unlikely to compensate for such severe disturbance. By combining vegetation resurvey with anatomical age determination, we conclude that climate change in the arid Himalayas produces a complex pattern of plant responses involving uphill migration due to warming, and vegetation decline due to increasing precipitation and soil disturbance. Actual effects of climate change in arid mountains are thus more complex than usually thought.

Material and Methods

The study was conducted on the southwest spur of the Tibetan class="Chemical">Plateau in Eastern Ladakh (Fig. 1A), Jammu and Kashmir State, India, on a high-altitude class="Chemical">plateau 15 km East of Tso Moriri Lake (33.1 N, 78.2 E) at elevations ranging from 5500 to 6150 m, between Chalung, Shukule and Chamser Kangri class="Chemical">pan class="Chemical">Peaks. The region is arid (Leh: 115 mm/yr, 3514 m, ca. 170 km NW of the study region, Gar: 54 mm/yr, 4232 m, ca. 160 km SE of the study region), because most of the summer monsoon precipitation is blocked by the main Himalaya Range. In the foothill desert areas around 3000–4500 m elevation, annual rainfall does not even reach 100 mm; at elevations above 4500 m precipitation tends to increase. Precipitation falling in summer above 5000 m is mostly snow. Winter precipitation is rather erratic and the snow layer is usually thin46. Soils have a coarse-grained structure, with a high percentage of large gravel, low water and organic matter contents, high pH (7–8) and relatively high concentrations of total N and P47.

In situ climate measurements

To gain detailed information about inter-annual climate variation in the region, we set up a network of automatic HOBO and TOMST® TMS microclimatic stations recording air and soil temperature (T), air relative humidity (RH) and soil moisture at hourly intervals (Fig. 2). Microclimatic stations have been in place since August 2008 to September 2014 and were placed at ~100 m vertical intervals. The climate records show that the growing season length, defined here as the number of days where the soil temperature constantly exceeded 0°C (frost-free days), shortens by 9.6 days per 100 m (from 90 to 28 frost-free days), while the mean temperature of the growing season decreases on average by 0.39°C per 100 m elevation (from 5.2 to 2.8°C) along the studied vertical gradient (5500–6150 m) encompassing the pan class="Species">alpine and subnival vegetation. Soil moisture content, averaged for the entire growing season, is generally low and varies from about 4 to 15%, with the highest values recorded at around 5700 m.

To monitor the local snow conditions, an automatic class="Chemical">PlantCamTM camera (www.wingscaclass="Chemical">pes.com) was class="Chemical">placed on a sloclass="Chemical">pe facing the class="Chemical">permanent class="Chemical">plots and climate loggers in the class="Chemical">period 2011– 2012, taking a class="Chemical">photo every day at 10 a.m. The close relationshiclass="Chemical">p between measured air relative humidity and local snow conditions allow us to calculate frequency of snow events over the whole measurement class="Chemical">period (for details see Dvorsky et al.24). Continuous snow cover usually occurs from February until Aclass="Chemical">pril but with great differences between years (Fig. 2D). There were 86 days with snow cover in a late winter class="Chemical">period (March–May) in 2009, 92 in 2010, 54 in 2011, 41 in 2012, and 89 in 2013 (Suclass="Chemical">pclass="Chemical">plementary Fig. S4). Temclass="Chemical">perature and relative humidity fluctuations were low in the Aclass="Chemical">pril–May class="Chemical">periods of 2010 and 2013 (Suclass="Chemical">pclass="Chemical">plementary Fig. S3), which coincided with unusually deeclass="Chemical">p snow covers, followed by a warm June and July. This led to high soil class="Chemical">pan class="Chemical">water content and increased the risk of frost heave (Supplementary Fig. S5), with mean daily minimum T ranging from −2.4 to −10.5 °C (absolute daily minimum T for June-July range from −4.3 to −18.1 °C), while the mean daily maxima reached +14.5 °C (Fig. 2B). The summer 2010 was snowiest in past eight years, there were 26 snowfall days in 2009 between June to August, 62 in 2010, 26 in 2011, 26 in 2012, and 27 in 2013.

Assessment of upward migration and thermophilization

To explore the possible uphill migration over the last decade, we examined the elevational distribution of subnival plant species in the study region (Supplementary Fig. S1). We compared plant distribution data collected in 2001–2003424548 with our recent distributional data collected in the same region in 2011–20142449. Depending on habitat diversity and species richness, individual localities were searched for about 30 minutes to one day. The elevation of individual floristic records was measured with an altimeter (Thommen, Switzerland) and Garmin pan class="Disease">GPS and stored in a floristic database of Ladakh4248.

To explore possible thermophilization of species assemblages, we analyzed temporal changes in species composition, cover and richness in 80 permanent plots. These plots were established within the class="Species">alpine (41 class="Chemical">plots, 5500–5800 m) and subnival zone (39 class="Chemical">plots, 5800–6030 m) in 2009, on a gentle west-facing sloclass="Chemical">pe (3–5 degrees of inclination) and similar exclass="Chemical">position. Each class="Chemical">plot was 1 m2 in size, and divided into one hundred 10 × 10 cm subclass="Chemical">plots with all sclass="Chemical">pecies of class="Chemical">pan class="Species">vascular plants and their cover being recorded. In 2012, we resurveyed all these plots following the same sampling protocol after a period of rapid climate warming (Fig. 2). The subnival zone was distinguished by a floristically-based nivality index24, the absence of steppe species and by the presence of subnival specialists. Also, the plots with cover >50% were all recorded below 5800 m (in most cases contained clonal graminoids such as Carex sagaensis and Kobresia pygmaea), marking transition between the alpine and subnival zone2445.

Plant age and growth determination

To evaluate the age and growth of plants in high elevation communities, we measured 426 individual specimens of 67 herbaceous species collected in the class="Species">alpine (43 sclass="Chemical">pecies) and subnival belt (24 sclass="Chemical">pecies) in August–Seclass="Chemical">ptember 2012–2014 (Fig. 1A). Due to dry soils in the study region class="Chemical">perennial class="Chemical">plants often have deeclass="Chemical">p taclass="Chemical">p-roots50 and distinct annual ring increments. All monocot sclass="Chemical">pecies (21%) were excluded because they are not suitable for herbochronological analysis. Three to four adult undamaged individuals of each sclass="Chemical">pecies were excavated and the washed roots class="Chemical">placed in class="Chemical">plastic bags containing class="Chemical">pan class="Chemical">ethanol until further analysis. The samples were analyzed in the Swiss Federal Research Institute WSL, in Birmensdorf, Switzerland by a method developed by Gärtner and Schweingruber51. Transverse, tangential and radial sections were cut from all studied individuals using a sliding microtome. Microscopic images of optimal ring width sequences were made for counting of the consecutive number of rings (Fig. 1D), and measuring the width of annual increments51. Plant ages were determined by ring counting in the oldest sections in the transition between the hypocotyl and the primary root (root collar). The possibility of non-transparent rings due to inappreciable growth in extremely cold or snowy summer cannot be excluded, therefore the numbers of counted rings represent minimum plant ages.

To assess the biological mechanisms indicative of thermofilization, we measured several morphological and ecophysiological traits (see Table 1 for abbreviations of the traits). These were: plant height, class="Disease">total dry biomass, number of flowers and shoots, seed weight, rooting declass="Chemical">pth, class="Chemical">pan class="Disease">leaf dry matter content, leaf carbon, nitrogen and phosphorus concentrations, leaf δ 13C, root nitrogen and phosphorus concentrations and the content and composition of nonstructural carbohydrates. Starch and fructans were analyzed using colorimetric methods. Soluble sugars, i.e. raffinose family oligosaccharides, sugar alcohols and simple sugars were quantified using high-performance anion exchange chromatography with pulsed amperometric detection (for details see Chlumská et al.52 and Supplementary Fig. S9). Leaf δC13, which measures the ratio of 13C over 12C (%), is an integrated, long-term measure of the ratio between internal and ambient CO2 concentrations (Ci/Ca) that reflects the intrinsic water use efficiency (WUE) of plants53.

The traits were measured in a minimum of 10 individuals (up to 50 for some species) from various elevations within the studied region and the average values were calculated per species. In addition to quantitative traits, each species was classified as clonal or non-clonal based on the type of belowground organs (rhizomes with adventitious roots versus primary tap-roots), and further into four space occupancy strategies, based on the rate of lateral spread (spreading–more than 10 cm per year; non-spreading–less than 10 cm per year) and persistence of connections between ramets (splitters–plants producing adventitious roots with main root decaying; integrators–plants not producing adventitious roots and/or with perennial main root50).

Habitat preference

In addition to plant traits, we included information on high-altitude species’ ecological optima in order to predict species’ response to climate change. Species’ optima on five environmental gradients were derived from vegetation composition of 369 plots (each 100 m2) sampled in a stratified design to cover major vegetation types over the study area244548. On each sampled plot, five environmental variables were estimated: (1) soil/substrate stability: 1–unstable (screes, dunes, solifluction soils), 2–partly stable (grasslands, steppes), 3–stable (rocky crevices, class="Species">Kobresia pygmaea mats); (2) light availability: 1–shaded (gorges, shaded rocky crevices, walls of stream banks), 2–class="Chemical">partially shaded (dense vegetation cover), 3–full light (sclass="Chemical">parse vegetation cover); (3) soil moisture: 1–dry (substrate usually without visible traces of class="Chemical">pan class="Chemical">water), 2–mesic, 3–wet (water level regularly but transiently above soil surface), 4–permanent surface water; (4) soil fertility: 1–low (semi-deserts, steppes), 2–medium (alpine meadows), 3–high (stables, animal resting places); (5) soil salinity: 1–no salt deposits on soil surface, 2–salt deposits scarce, 3–salts forming a continuous crust. From these data, indicator values of a species were calculated as arithmetical means of indicator values in individual plots in which that species was recorded, weighed by the logarithm of its cover.

To obtain a robust estimate of the elevational optima of the species, we calculated response curves fitted with HOF models54. Species response curves were derived from 4,150 vegetation plots (each 100 m × 100 m) sampled over the entire Ladakh between 1999 and 2014 (Supplementary Figs S1 and S6). The dataset contain more than 122,000 records of occurrence of pan class="Species">vascular class="Chemical">plant sclass="Chemical">pecies along exceclass="Chemical">ptional elevational gradient from 2800 m to 6150 m244548. class="Chemical">pan class="Chemical">Plant names follow Klimes and Dickore55.

Physico-chemical characterization of the soil

To characterize physico-chemical parameters of soils, 50 g was taken as a composite of 5 subsamples from randomly selected points within each plot in 2009. The samples were air dried for 24 h on class="Chemical">aluminum class="Chemical">plate, after than class="Chemical">placed into sterile class="Chemical">pan class="Chemical">polypropylene bags (Nasco Whirl-Pak®) and transported to the laboratory for analyses. In the laboratory, the soil samples were oven-dried at 100 °C, grinded in a mortar and sieved to 2 mm fraction after the removal of roots. Major cations (Ca+2 , Mg+2 , K+ and Na+) as well as nitrogen and phosphorus were measured in all soil samples. Cations were quantified through atomic absorption spectroscopy (AAS) using SpectrAA 640 (Varian Techtron) at the Analytical laboratory of Institute of Botany, Czech Republic. Ammonia, nitrate and total nitrogen were determined colourimetrically after Kjeldahl mineralization using automatic FIAstar 5010 Analyzer (Tecator). Phosphorus was determined colourimetrically after digestion in HClO 4 using SHIMADZU UV - 1650PC spectrophotometer. Other physico-chemical data were also measured: pH, water content, organic matter content (OM), and texture (fraction of particles >0.5 mm in diameter). The alpine soils had significantly higher total N, P, Ca concentrations and organic matter content and lower texture than subnival soils (Supplementary Fig. S12).

Statistical analysis

To assess temporal changes in the plot-level vegetation records (number of class="Species">vascular class="Chemical">plant sclass="Chemical">pecies, total class="Chemical">plant cover, individual sclass="Chemical">pecies’ frequency and cover) between the surveys (2009 vs 2012), and between the habitats (class="Chemical">pan class="Species">alpine vs subnival) and their interactions, we used generalized linear mixed-effect models (GLMM5657). The tests were based on the restricted maximum likelihood (REML) approach approximated by chi-square distribution. The statistical significance was assessed by computing Bayesian highest probability (HPD) intervals using Markov chain Monte Carlo simulations (999 permutation in each test).

To explore shifts in species composition, we performed multivariate ordination analysis RDA, separately for the pan class="Species">alpine and subnival belts. The variance class="Chemical">partitioning class="Chemical">procedure was class="Chemical">performed in RDA with exclass="Chemical">planatory variables (year) and co-variables (elevation, class="Chemical">plot identity) to remove their effects and to obtain a net effect of the year. To assess whether comclass="Chemical">positional change is driven by uniform shift among sclass="Chemical">pecies or rather by differential sclass="Chemical">pecies resclass="Chemical">ponses, RDA was class="Chemical">performed both with and without standardization by samclass="Chemical">ple norm58. Comclass="Chemical">positional differences in class="Chemical">plant assemblages were tested by 999 class="Chemical">permutations.

To explore further the ecological causes for the species temporal responses, we used plant functional traits and habitat indicator values in explaining the vegetation changes. Two approaches were employed: (1) plot-based analysis using trait averages (weighted by log-transformed species cover abundance; hereafter CWM, community weighted means), compared between two sampling periods using GLMM as described above. Because obtained results reflect mainly the features of dominant species, we also used an alternative approach which places individual taxa into focus instead of individual plots: (2) species-based analysis of the trait-environment relation was a two-step analysis, in which the temporal change of individual species was first quantified (by scores on the first ordination axis of partial RDA, with the year as the only explanatory variable and elevation and plot identity as covariables) and then those scores were related to the traits and habitat indicator values of individual species using a conditional inference trees59. Conditional inference trees has several advantages over other classification and regression trees algorithms, including the statistical testing of each split through permutation, no need for problematic pruning of over-fitted trees, and no selection bias towards variables with many possible splits or missing values.

Additional Information

How to cite this article: Dolezal, J. et al. Vegetation dynamics at the upper elevational limit of pan class="Species">vascular plants in Himalaya. Sci. Reclass="Chemical">p. 6, 24881; doi: 10.1038/sreclass="Chemical">p24881 (2016).