Moisture‐mediated responsiveness of treeline shifts to global warming in the Himalayas

Among forest ecosystems, the alpine treeline ecotone can be considered to be a simplified model to study global ecology and climate change. Alpine treelines are expected to shift upwards in response to global warming given that tree recruitment and growth are assumed to be mainly limited by low temperatures. However, little is known whether precipitation and temperature interact to drive long‐term Himalayan treeline dynamics. Tree growth is affected by spring rainfall in the central Himalayan treelines, being good locations for testing if, in addition to temperature, precipitation mediates treeline dynamics. To test this hypothesis, we reconstructed spatiotemporal variations in treeline dynamics in 20 plots located at six alpine treeline sites, dominated by two tree species (birch, fir), and situated along an east–west precipitation gradient in the central Himalayas. Our reconstructions evidenced that treelines shifted upward in response to recent climate warming, but their shift rates were primarily mediated by spring precipitation. The rate of upward shift was higher in the wettest eastern Himalayas, suggesting that its ascent rate was facilitated by spring precipitation. The drying tendency in association with the recent warming trends observed in the central Himalayas, however, will likely hinder an upslope advancement of alpine treelines and promote downward treeline shifts if moisture availability crosses a critical minimum threshold. Our study highlights the complexity of plant responses to climate and the need to consider multiple climate factors when analyzing treeline dynamics.


INTRODUCTION
Climate fundamentally controls the ranges of species and the rates of ecological processes in natural ecosystems (IPCC, 2014). Multiple evidences indicate that recent climatic changes have large effects on the structure and dynamics of plant terrestrial ecosystems (Lenoir, Gégout, Marquet, De Ruffray, & Brisse, 2008;Peñuelas et al., 2013, Piao et al., 2014. Alpine plant communities and subalpine forests are potentially among the most affected terrestrial ecosystems, so treeline ecotones, which separate them, are a potential monitor of the impacts of climatic warming (Beckage et al., 2008;Harsch, Hulme, Mcglone, & Duncan, 2009).
The Himalayas and Tibetan Plateau hold the world's highest mountain systems, and the highest treelines in the northern Hemisphere (Miehe, Miehe, Vogel, Co, & La, 2007), and have a broad gradient of bioclimatic zones along their altitudinal gradient. Himalayan birch (Betula utilis D. Don) and Himalayan fir (Abies spectabilis (D. Don) Spach) form the natural treelines near 4000 m a.s.l. in the central Himalayas. Himalayan alpine treelines thus constitute suitable locations for determining the impacts of climate change on forest ecosystems. A recent study reported that the growth of Himalayan birch at the treeline was primarily limited by pre-monsoon precipitation (Liang, Dawadi, Pederson, & Eckstein, 2014) and its spatial distribution depends on soil moisture availability (Shrestha, Ghimire, Lekhak, & Jha, 2007). Precipitation is expected to control the rates of shifts of such moisture-driven alpine treelines. Local-scale studies (Chhetri & Cairns, 2018;Gaire, Koirala, Bhuju, & Borgaonkar, 2014;Shrestha, Hofgaard, & Vandvik, 2014;Tiwari, Fan, Jump, Li, & Zhou, 2017) have shown treeline-shift rates differ among specific areas; however, these studies may not be able to disentangle the effects of precipitation on treeline dynamics from other local-scale interactions due to a lack of approaches considering a macroscale precipitation gradient, i.e. the longitudinal range 80-88 E and approximately 800 km across the central Himalayas. Macroscale studies encompassing east-west precipitation gradients are therefore necessary to determine the responses of alpine treelines to temperature and precipitation changes across the central Himalayas.
The specific objectives of this study were to: (i) reconstruct the spatiotemporal patterns of treeline shifts during the last 150 years across the central Himalayas, and (ii) determine the role played by precipitation on the dynamics of Himalayan treelines. Based on the evidence that precipitation decreases westwards across the central Himalayas, we hypothesized that warminginduced treeline shifts would be mediated by regional precipitation gradients, and would be more substantial in the east. We also hypothesized that the components of these shifts (tree growth, tree recruitment, and treeline position) would be uncoupled and shaped by different climatic factors along this precipitation gradient.

Study area
We established a network of 20 treeline plots at six sites located across an east-west gradient in Nepal to study the macroscale dynamics of treelines of two major treeline species, Himalayan birch (Betula utilis D. Don) and Himalayan fir (Abies spectabilis (D. Don) Spach), across the central Himalayas (Figure 1). The climate of this region is influenced by two weather systems. The eastern region is highly influenced by the southern Indian monsoon during summer, and the western region is influenced by westerly winds during winter (Yao et al., 2012). The total annual rainfall ranges from about 3500 mm in the east to 800 mm in the west. About 80% of the precipitation falls during the monsoon season. The mean annual temperatures along the treeline ecotone range from 3.7º to 4.9 °C, and the mean soil temperature at the treeline during the growing season is about 7.5 ± 0.6 °C (Müller et al., 2016), slightly higher than that recorded on the southeastern Tibetan plateau (6.0 ± 0.3 °C)  and the global average (6.7 ± 0.8 °C) (Körner & Paulsen, 2004).

Research plots and sampling design
The position of the upper treeline was defined by the presence of upright trees with a minimum height of 2 m at the maximum altitude and a continuous distribution above the timberline (forest coverage >30%) in the plots (Camarero & Gutiérrez, 2004;Liang, Wang, Eckstein, & Luo, 2011. Treeline altitude varies from 3800 to 4250 m a.s.l. All plots included the upper limits of the tree species in closed forests following previous field protocols (Camarero & Gutiérrez, 2004;Liang, Wang, Eckstein, & Luo, 2011;.  (Camarero & Gutiérrez, 2004). The x, y origin (0, 0) for all plots was in the lower left corner facing upslope. The geographic information of all four corners was noted using GPS. The location of each tree of the target species was noted relative to the origin (x, y = 0, 0). The diameter at breast height (DBH) of each tree in the plot was measured at 1.3 m using diameter tapes. Tree height was measured using clinometers for trees with height >2 m and a measuring stick for trees with height ≤ 2 m (Wang, Zhu, Liang, & Camarero, 2016). The average height and coverage of shrubs and grasses were measured along altitudinal transects from the current treeline to 10 m above the treeline and the vegetation thickness index (TI) was calculated as average plant height × coverage .
A core was collected using an increment borer as close as possible to the ground from the main stem of each living tree with a DBH >5 cm for estimating the germination dates of the trees (in total, n=1594 birches and n=414 firs were sampled). A core was similarly collected at a height of 2 m for each site and species (n=110 birches and n=37 firs) for calculating the age at 2 m (Table   S1), which was then used to identify the treeline position (altitude of trees 2 m tall). The ages of seedlings and saplings (individuals with DBHs <5 cm) were estimated by counting terminal branch whorls and scars along the main stem (Camarero & Gutiérrez, 2004;Liang, Wang, Eckstein, & Luo, 2011). Then, wood sections were collected from the root collar (n=105 trees) to validate the estimates. The estimated ages of the firs were very similar (± 1.5 yrs.) to the estimated germination dates, but the estimated ages of birches were four years older than the germination dates. Age-DBH relationships were then determined for each site and species using linear regression models ( Figure S1). Age and DBH were highly correlated (P < 0.001) for all treeline sites and species.
These relationships were used to estimate the age of trees with a rotten center or missing pith.

Data analyses
The dynamics of recruitment on a decadal scale were reconstructed from the tree ages for each site and compared with past climatic records. An integrated ice-core δ 18 O series from the Tibetan Plateau was used as a proxy of annual and summer temperature (Thompson et al., 2006), and chronologies derived from the width of tree rings of the moisture-sensitive Himalayan birch from the central Himalayas (Liang, Dawadi, Pederson, & Eckstein, 2014) were used as a proxy of pre-monsoon precipitation. Similarly, the spatiotemporal variability of tree density and treeline position were characterized for 50-year long intervals (Camarero & Gutiérrez, 2004;Liang, Wang, Eckstein, & Luo, 2011;) ( Figure S2).
We developed a comprehensive data set of several potential explanatory variables of treeline shifts: aspect, slope, altitude, tree species, shrub and herb coverage above the treeline, TI, distance between a treeline and a species line (TL-SL), mean/total annual, spring, summer and winter temperatures (AT, SPT, SUT and WT) and precipitation (AP, SPP, SUP and WP) and changes in annual, spring, summer and winter temperature (CAT, CSPT, CSUT and CWT) and precipitation (CAP, CSPP, CSUP and CWP). We also considered annual mean maximum (ATmax) and minimum temperatures (ATmin).
To obtain climate data, we used high-resolution interpolated monthly climatic data (temperature and precipitation) due to the paucity of instrumental climatic data for high Himalayan altitudes. Temperature data were retrieved from the Climate Research Unit (CRU) TS3.23 database, which has a spatial resolution of 0.5° (Harris, Jones, Osborn, & Lister, 2014). We obtained these data from the 0.5° grids located nearest to each sampling site and applied a published standard rate of temperature lapse to calculate site-specific temperature after adjusting altitudinal effect (Kattel et al., 2013). Temperature data after 1970 were used for the analyses because the climatic data were most reliable for this period. In our early studies in the Langtang valley and Mt. Everest, we tested the representative of variations of CRU gridded temperature (Dawadi, Liang, Tian, Devkota, & Yao, 2013;Liang, Dawadi, Pederson, & Eckstein, 2014), and showed that the variations of the gridded monthly temperature appear to represent well those thermal conditions in the treeline sampling areas.
High-resolution (0.25° grid) precipitation data were retrieved from the Tropical Rainfall Measuring Mission (TRMM) satellite, available online (http://giovanni.gsfc.nasa.gov/giovanni/) since 1998. These data were not corrected for the effects of altitude, because few local data were available for local precipitation. TRMM radar data have been widely used to show the relationship between precipitation and topography in the Himalayas (Bookhagen & Burbank, 2006;Shrestha, Singh, & Nakamura, 2012). For example, TRMM data capture a clear zonal distribution pattern of rainfall along the southern periphery of the central Himalayas (Liang, Dawadi, Pederson, & Eckstein, 2014). We also retrieved evapotranspiration data 0.25° spatial resolution) from the Global Land Data Assimilation System (GLDAS) available online (http://giovanni.gsfc.nasa.gov/giovanni/). A generalized linear model (GLM) was used to identify the variables responsible for the rates of treeline shift during the last 150 years using R (R Development Core Team, 2017) and we used the R "relaimpo" package to evaluate the relative importance of the variables in linear regression models based on the successive-sweep method (Grömping, 2006). Partial correlation analyses were also used to measure the linear relationship between two variables, such as treeline shift and summer monsoon precipitation, after having excluded the effect of a third variable, such as premonsoon precipitation.

RESULTS
All plots showed a common demographic structure, characterized by increasing tree recruitment in recent decades, particularly after the 1940s (Figure 2). The decadal rate of tree recruitment was significantly correlated (P<0.05; Table S2) with several temperature proxies (δ 18 O series from ice cores from the Tibetan Plateau). Tree recruitment was mainly linked to the rapid temperature increase during the second half of the 20 th century. Tree recruitment above the current treeline was higher for birch than fir.
We reconstructed the spatiotemporal variation in tree density and treeline position based on the age and height of individual trees for 20 treeline plots situated at six Himalayan sites ( Figure   S2). The average upward rates of treeline shift during the last 150 years ranged from 0 to 0.37 m y -1 with an average of 0.11 ± 0.10 m y -1 , i.e. the ascent varied from 0 to 56 m. Treeline positions in eight of the 20 plots have been static or have changed little during the last 150 years. No dead trees or stumps were found above the current treelines, indicating that the sampled trees were the first generation recruited above the treeline during at least the last two centuries and suggesting that the subalpine forest was established during the last 2-3 centuries ( Figure S2). The rate of treeline advance in our study, however, was not high but it increased in eastern Nepal (0.20-0.37 m yr -1 ; see Table 1).
The rates of treeline shifts varied significantly along the east-west precipitation gradient, suggesting that precipitation and temperature interacted to control the regional treeline dynamics.
Generalized linear models indicated that spring precipitation and mean annual maximum temperature together explained 60.5% of the total variance in treeline shift (Table 2 and Figure 3).
Nevertheless, spring precipitation contributed more than mean annual maximum temperature (53% vs. 47%) to explain treeline-shift rates and it was the most significant driver (P<0.001), followed by mean annual maximum temperature (P<0.01). As showed by the partial correlation analysis, summer monsoon precipitation is not a critical factor in controlling treeline shifts (Table   S3). Spring precipitation is significantly more important than summer precipitation in affecting treeline dynamics.
We assessed the role of biotic interactions by calculating a vegetation thickness index (TI) which is the product of cover and height of alpine plant vegetation (shrubs, grasses) situated just above the treeline . The TI did not have a significant effect on treeline-shift rates, which may be explained because dwarf shrubs were sparsely distributed above the current treelines across all sites, resulting in a low TI (0.12-0.82; see Table 1).

DISCUSSION
Increasing tree recruitment at treeline ecotones in recent decades is likely linked to ongoing climatic warming. Warmer and stable spring and summer conditions can improve seed germination and increase the survival of seedlings at treelines, thus promoting tree recruitment (Camarero & Gutiérrez, 2004;Körner, 2012). Higher recruitment of birch above current treeline may be due to the fact that light seeds can be easily dispersed by wind long distances beyond the current treeline (Nathan & Muller-Landau, 2000). Himalayan fir, however, produces seeds within a cone that probably are not dispersed such long distances by wind (Coates, 2002;Greene & Johnson, 1995).
The treeline positions and their shift rate in the central Himalayas were spatially heterogeneous. Indeed, such results essentially match those from the most recent meta-analysis (Harsch, Hulme, Mcglone, & Duncan, 2009), with 52% of sites showing treeline advance upslope worldwide. The average rate of treeline shift was lower than those observed on the Tibetan Plateau

Jumla (J) and Humla (HM) sites than at the eastern study sites. Precipitation in the central
Himalayas is regulated by two weather systems, a summer monsoon during June to September and westerly winds during November to May (Yao et al., 2012). Likewise, western Nepal is situated between the influences of a cyclonic cell to the west and an anti-cyclonic cell to the east, which causes substantial dry conditions and decreasing trends of precipitation westward (Wang, Yoon, Gillies, & Cho, 2013). We observed a higher treeline-shift rate in the eastern study area which supports the idea that moisture drives treeline-shift rates across the studied E-W precipitation gradient. These results are supported by the facts that tree growth is also constrained by low spring Temperature may not drive treeline shift rates per se but through interactions with precipitation and water available to trees. Moisture stress can have significant and negative influences on seed production and viability of seeds across the circumarctic forest tundra ecotone, being a bottleneck to treeline advance (Brown et al., 2018). As shown in western North America, dry conditions can limit seedling establishment and hence tree species shifts (Elliott, 2012;Elliott & Petruccelli, 2018;Kueppers et al., 2017;Lloyd & Graumlich, 1997;Moyes, Germino, & Kueppers, 2015). Warmer temperatures could lead to a deficit in soil moisture driven by increased evapotranspiration at alpine treeline ecotones (Trujillo, Molotch, Goulden, Kelly, & Bales, 2012).
Satellite-derived estimates suggest that the winter rate of evapotranspiration has increased significantly during the last three decades ( Figure S3) at four of the six study regions, which might negatively affect moisture availability during spring. Our analysis of spring precipitation and temperature of nearby treeline sites found that precipitation has decreased significantly at some sites, whereas air temperature has increased at all sites ( Figure S4), in agreement with observations of shrinking Himalayan glaciers linked to rising temperatures (Yao et al., 2012). Decreasing trends in the snow accumulation (Thompson et al., 2000) and increasing values of drought severity reflected by changes in the Palmer Drought Severity Index (Sano, Ramesh, Sheshshayee, & Sukumar, 2011) during the last 150 years suggested weakening South Asian Monsoon ( Figure S5).
Furthermore, warming will enhance transpiration, and intensify photoinhibition and hydraulic failure in drought-prone sites (Allen, Breshears, & Mcdowell, 2015;Choat et al., 2018). Hence, interactions between spring precipitation and the annual maximum temperature are likely driving treeline shifts across the dry Himalayas.
Tree recruitment, growth and treeline shifts were regulated by distinct climatic factors, supporting our second hypothesis. Tree growth at the upper timberline of Himalayan birch and Himalayan fir was mainly limited by low pre-monsoon precipitation in the central Himalayas (Liang, Dawadi, Pederson, & Eckstein, 2014;Tiwari, Fan, Jump, Li, & Zhou, 2017), being different with thermal treelines and cold biomes of the Northern Hemisphere (Li et al., 2017;Rossi et al., 2016). Drought stress delays the onset of cambial cell division and reduces the growth rate at drought-prone forest sites (Ren et al., 2018). A high percentage of locally missing rings coincided with dry and warm pre-monsoon seasons at the upper timberline of Himalayan birch in the central Himalayas, providing additional evidence that moisture controls its growth and survival (Liang, Dawadi, Pederson, & Eckstein, 2014). In contrast, long-term tree recruitment was constrained by cold temperatures, so the treeline shifts were jointly controlled by both precipitation and temperature. These results differed from those observed at temperature-limited treelines, where tree growth, tree recruitment, and treeline position usually responded consistently to temperature (Harte & Shaw, 1995), but were in agreement with observations on Nothofagus pumilio treelines from Patagonia, on Larix sibirica forest-steppe ecotone in northern Mongolia and several conifer species in western North America whose dynamics depend on precipitation and temperature (Daniels & Veblen, 2004;Dulamsuren, Hauck, & Leuschner, 2010;Elliott & Cowell, 2015;Elliott & Petruccelli, 2018;Fajardo & McIntire, 2012;Kueppers et al., 2017;Morgan, Losey, & Trout, 2014;Moyes, Germino, & Kueppers, 2015).
To conclude, treeline shifts in the Himalayas are better explained by interactions between precipitation and temperature than to either alone, and their interaction may lead to changes in the availability of soil moisture to trees. The upslope treeline shift rates over the last 150 years were low and its pace was primarily controlled by spring precipitation in association with changes in mean annual maximum temperature. Warm and moist spring climatic conditions led to increased tree recruitment, thereby facilitating their ascent. The drying tendency in association with the recent warming trends in the central Himalayas, however, will likely hinder an upslope    Triangles represent major mountain peaks near the sites. Blue squares represent plots with Himalaya birch, and yellow squares represent plots with Himalaya fir. (six sites, 20 plots; see Figure 1). Decadal recruitment was significantly and positively correlated with mean summer temperature (r = 0.32 -0.86) reconstructed from the ice cores (see Table S1 in Supporting Information). The red line in 'b' represents the 10-year moving average.

FIGURE 3
A three-dimensional diagram depicting the relationships between changes in treeline position during the last 150 years, total spring precipitation and mean annual temperature across the central Himalayas. The surface is based on predictions of fitted generalized linear models (see Table S2 in Supporting Information). Points show actual treeline-shift rate.     Table 1 and Fig. 1). Closed symbols in the respective 50-year intervals represent trees that established during that period (indicated at the top of each window).
Open symbols represent trees that established in previous 50-year periods. The shapes of the symbols represent the establishment periods (e.g. upward triangles represent trees established during 1711-1760).