Accelerating upward treeline shift in the Altai Mountains under last-century climate change

Treeline shift and tree growth often respond to climatic changes and it is critical to identify and quantify their dynamics. Some regions are particularly sensitive to climate change and the Altai Mountains, located in Central and East Asia, are showing unequivocal signs. The mean annual temperature in the area has increased by 1.3–1.7 °C in the last century. As this mountain range has ancient and protected forests on alpine slopes, we focus on determining the treeline structure and dynamics. We integrated in situ fine-scale allometric data with analyses from dendrochronological samples, high-resolution 3D drone photos and new satellite images to study the dynamics and underlying causal mechanisms of any treeline movement and growth changes in a remote preserved forest at the Aktru Research Station in the Altai Mountain. We show that temperature increase has a negative effect on mountain tree growth. In contrast, only younger trees grow at higher altitudes and we document a relatively fast upward shift of the treeline. During the last 52 years, treeline moved about 150 m upward and the rate of movement accelerated until recently. Before the 1950s, it never shifted over 2150–2200 m a.s.l. We suggest that a continuous upward expansion of the treeline would be at the expense of meadow and shrub species and radically change this high-mountain ecosystem with its endemic flora. This documented treeline shift represents clear evidence of the increased velocity of climate change during the last century.


study Area
The Aktru Research Station is located in the Southeastern Altai Republic (Russia) close to the borders of Mongolia and China in the Central Eurasian Continent (50°06′03″N, 87°40′14″E). With an altitude of 2100 m a.s.l., the station is situated in the high alpine range of these mountains. There are no major anthropogenic disturbances on the surrounding ecosystems and vegetation, since the area is located in a very remote (with a 3-day drive from the closest airport) and a high-altitude part of the Asian continent. The nearest settlement (Kuray, a small village) is located downstream from the station at about 30 km southeast and the nearest town is Gorno-Altaisk, about 250 km northwest of the station. The Aktru Research Station was founded in 1956 by V.M. Tronov and has been in operation until today under the supervision of Tomsk State University. There is no evidence or documentation on any type of logging, big herbivore (wild or domesticated) grazing, or fire-related impact, at least for the last 50 years, on the vegetation of the mountainside where our study site is situated.
The study site ( Fig. 1) is located on the northwestern slope of the Severo-Chuyskiy Mountain Chain, just above the research station, at an altitude between 2150 and 2300 m a.s.l. The climate of the area is characterized by low temperatures (annual average −5.2 ± 1.9 °C, summer average +8.7 ± 2.3 °C), high average annual precipitations (500-800 mm), and high diurnal temperature variation (of 15-20 °C) 29 . The site is considered to be climatically representative of the Altai, and includes several mountain peaks and glaciers 29 . The slope of this mountainside is 43.4-54.1%.

Results
The densest timberline in the study site, composed by two species (Pinus siberica 50-60%, Larix siberica 40-50%), is located up to an altitude of ≈2150 m a.s.l. (Fig. 2). Above this level (approximately between 2150 and 2250 m a.s.l.), a treeline composed of sparse groups of trees (belonging to the same two species of the treeline) and shrubs (mostly of Juniper spp., Vaccinium spp. and Ribes spp.) were detected. From 2250 m up to 2300 m a.s.l., only isolated trees and herb-rich meadows were identified (tree species line). No trees or shrubs with a d.b.h. >1 cm (see methods) were detected above ≈2300 m a.s.l. A few landslides and groups of rocks were located on the northern side of the treeline area. The whole study site (Fig. 2) from 2150 m to 2300 m a.s.l. may be defined, according to the common scientific convention 49 , as a "high-altitude treeline ecotone".
As expected, we found a significant decline in tree size towards higher altitudes (Fig. 3). In particular, tree height had the highest negative correlation with the altitude (Spearman's rank correlation ρ = −0.76, P < 0.01), followed by the crown width (ρ = −0.66, P < 0.01) and d.b.h. (ρ = −0.54, P < 0.01). No trees higher than 6 m where found over 2250 m a.s.l., even though a few relatively tall trees (7-9 m) were able to live at an altitude of 2240-2250 m a.s.l.
The dendrochronological analyses of the 48 wood cores collected from different altitudes along the transect of the treeline ecotone, shows a slightly declining mean annual growth curve between 1902 and 2017, with some periods (1902-1904, 1914-1916, 1923, 1928-1932, 1952-1956, 1961-1962, 1968) that peaks above the annual average tendency (Fig. 4). However, no relevant peaks are detected after the year 1968 (i.e. during the last 50 years of the dendrochronology).
After normalizing the time-series data with different methods (Fig. 5), the detrended curves (Spline and Friedman's Super Smoother in comparison with Straight and Horizontal Lines) show similar peaks with respect to the raw curve (in Fig. 4) and a general decreasing growth trend, as confirmed by the Mann-Kendall statistics (τ = −0.315, P < 0.01).
Tree basal area shows an expected significant negative correlation with both the year of birth of the trees (Kendall's rank correlation τ = −0.72, P < 0.01; Fig. 6) and the altitude (τ = −0.69, P < 0.01). In contrast, mean annual growth is not significantly influenced by the altitude (τ = −0.23, P = 0.30; Fig. 6) and shows a hump-shaped tendency when plotted against the birth year of trees (quadratic function, R 2 = 0.42). In fact, the groups of the oldest and youngest trees show a similar mean annual growth, which is lower than that of mid-age trees (Fig. 6). The positive-saturating relationship between the altitude and the birth year of trees (quadratic function, R 2 = 0.85) shows detectable shifts as schematized in Fig. 7, where we find strong evidence that trees born www.nature.com/scientificreports www.nature.com/scientificreports/ before 1954 do not grow above 2,150 m and trees born before 1999 do not grow above 2,200 m a.s.l. The maximum altitude of 2,300 m a.s.l. was reached between 1994 and 2002, while after 2002 the youngest trees sampled (born around 2006) grow at a lower altitude of 2,250 m a.s.l. (Fig. 7).
To determine if the shape of the mean annual growth curve were a simple biological artifact or related to the changes in climatic conditions and establishment altitude, we compared only the first 12 years of growth of all sampled trees (e.g. we compared the first 12-years mean growth of trees born in 1900s vs. 1980s vs. 2000s). We show (Fig. 7, right plot) that the oldest trees (which are located, on average, at the lowest altitude on the treeline ecotone at ≈2150 m a.s.l.) grew more than mid-age and younger trees (located at medium [≈2200 m a.s.l.] and higher [2250-2300 m a.s.l.] altitudes, respectively).
From a comparison between the variation of the mean annual temperature 29 and our dendrochronological time-series (Fig. 8), we discovered that before 1978 the regional mean annual temperature remained below −4.3 °C. This threshold was, therefore, used to interpret the behavior of the hyperbolic curve ( Fig. 8 upper-left panel) that emerged from the relationship between the mean annual growth and mean annual temperature (R 2 = 0.34). This curve, at about −4.3 °C can be divided in two linear regressions (<−4.3 °C negative, R 2 = 0.59; 4.3°C positive, R 2 = 0.39), which allow an easier estimation of the missing annual temperature records (before 1956 and after 2009) to fill the gap of the available data in the literature for the area. The only long temperature record 29 for the study-site covers a period between 1956 and 2009 and the most recent records   31 are from close meteorological stations but situated at lower altitudes.
The estimation for the years 1902-1955 (derived from the linear regression equation at <−4 °C) shows that temperature ranged between −5 °C and −10 °C, with some lower peaks between 1952-1955, 1926-1928, 1914-1916 and 1902-1904. The estimation for the years 2010-2017, shows a slight decreasing temperature trend, with a drop of 5-10 °C in the recent years (Fig. 8).
The results of the split sample test showed that the correlation coefficient and the product means for all calibration and verification periods were significant at P < 0.01 level. The sign tests were statistically significant in calibration and in the verification periods from 1956 and 1978 and from 1979 to 2009 at P < 0.05 level. The values of RE (0.86 and 0.77, respectively) and CE (0.32 and 0.36, respectively) were always positive, indicating that the regression was stable and reliable. www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion and Conclusion
Similar to the patterns documented in other high-mountains of the world 4-6,10,11 , we found that the increase of temperature during the last century of climate change has led to an alteration of the tree growth trend and an upward shift until recently of the treeline in a remote preserved forest at the Aktru Research Station in the Altai Mountains. This can be attributed to the fact that areas with a harsh growing condition for woody species, previously covered by snow for most of the year, would become suitable for tree growth as the mean annual temperature increases, thus enabling a upward shift and expansion of the treeline ecotone 12 .
We showed that, as predictable, biometric (allometric) parameters, such as tree height, d.b.h. and crown width decrease with the altitude, with significant differences between trees in the dense treeline (at about 2,150 m a.s.l.)  www.nature.com/scientificreports www.nature.com/scientificreports/ and those growing at higher elevation (up to 2,300 m a.s.l.). These differences, however, could be due not only to an upward shift of the new saplings but also to a differential growth trend related to altitudinal effects 50,51 . In fact, trees living at the limits of the treeline may show a reduced growth when subject to harder climatic conditions  The raw time-series (highest plot) is detrended and standardized by five methods. The detrended curve shown accounts for the tree's natural biological growth trend. Except for the Ar, the other four detrended curves are plotted over the raw series and show a decreasing growth trend, particularly after 60 years. Notice that the x-axis here is the opposite to that in Fig. 4 (showing the curve increasing from the oldest to the youngest age, i.e. the abscissa represents the age of the tree while in Fig. 4 it represents the year of birth) and should be chronologically read from right to left. www.nature.com/scientificreports www.nature.com/scientificreports/ then their physiological limits (i.e. lower temperatures, frozen soil, stronger wind, etc.) 52 . To test if a variation in the tree growth trends and an upward shift of the treeline ecotone might be caused by increasing local and global temperatures, we further analyzed dendrochronological samples collected along some altitudinal transects.  From an analysis of the complete dendrochronological time series from 1902 to 2017, we found a slight decrease in the general growth trend and relevant growth peaks for many years, all dated before 1968, when the mean annual temperature was always below the identified −4.3 °C threshold. In the last 50 years, trees did not show any sign of increased growth compared to the average trend for the century. This could be due to the fact that high-mountain trees are quite sensitive to low temperatures and grow more under a certain threshold (−4.3 °C in our study-site) [53][54][55] . In fact, comparing the available mean annual temperature observations from 1956 to 2009 with our dendrochronological time series, we detected that -in this time frame -the highest growth peaks (in 1952-1956, 1961-1962 and 1968) correspond (with a few years of lag) to the lowest temperature for the period. However, after 1968 (i.e. in the last 50 years except for 1983-1988, when the temperature dropped for a while and the tree growth temporarily increased), the temperatures were often above a −4.3 °C threshold, and the tree growth was constant (without any sudden growth peak) and slightly decreasing. The fact that the precipitation trend shows no change for the whole region in the last 50 years 31 gives more weight to temperature importance and to its low extremes.
Our reconstructed temperature estimates for the period 1902-1955 derived from the dendrochronological time series, which was missing from the available observation on the local weather station 29,31 , well matched the periods of extreme cold temperatures in 1902-1906, 1913-1915, 1922-1923, 1926-1928 and 1950-1952 recorded by other biological and non-biological sources at comparable regional [56][57][58] and continental scales 59,60 .
We also checked whether last-century temperature increase would have contributed to the shift of the treeline ecotone to higher altitudes. We showed that trees have actually had an accelerating upward movement during the last century, yet a recent slight downward movement. In fact, during the last 115 years, the "residence time" Our finding that young trees (approx. 17-19 years old, sampled at the higher altitudes) grow less than old trees (93 ± 25.21 years old, sampled at lower elevations), during their first 12 years of life, confirms a reduced tree growth at higher temperatures in the high-mountains 61 . In summary, this evidence suggests that the temperature increase during the past century may pose a double effect on mountain forests: in our study-site, tree growth In general, the rapid upward treeline shift that we detected combined with the reduced tree growth, may radically change the composition and landscape of the montane ecosystems. In areas on our study site where we did not recorded any trees, we found that a high percentage of the soil was covered by highly biodiverse meadows (with approx. 14 identified herb species, data not reported here) and by both economically and culinary valuable shrubs, including Juniper spp., Vaccinium spp. and Ribes spp. A continuous upward shift tendency, as the one documented for the last century in this study, would presumably lead to an expansion of the forest at the expense of meadow and shrub species, which would have no other area to migrate to once the mountain summit is reached 62 . On the other hand, we also identified a partial tendency inversion for the most recent decade of tree expansion at high elevations, which could represent a "temporary relief " for the survival of high-altitude meadow and shrub species. Trees did not move upward after 2004 and this might reflect the recent temperature decrease observed both at the local meteorological stations 31 and from our estimates derived from the hyperbolic relationship between the mean annual growth and the 1956-2009 field-recorded temperatures. In fact, we found that as −4.3 °C represents a threshold between the negative growth tendency before the 1978 and the positive one after this year, estimated temperatures from 2010 to 2017 show a decreasing trend that, at least, could temporarily hinder the upward shift of the treeline. Further resources to reach and study these remote high-mountain areas would clarify whether this threshold of −4.3 °C applies to other treelines of the Altai Mountains.
We recognize that our study focuses on a specific area of the Russian part of the Altai Mountains and we need to be careful discussing the applicability our results to other areas of the region. However, the novel findings should not wait for publication for many years so that other sites can be analyzed but rather used to stimulate similar research at other similar ecosystems by immediate publication.
The scale of impacts of climate change in the Arctic and Alpine areas are of far more than academic importance and generalizations must be balanced by site-specific studies of relevance to local populations and local biodiversity hot-spots such as the Aktru high-mountains.
In conclusion, based on in situ treeline data and a dendrochronological analysis, we found an accelerating upward shift of treeline on a high-mountain vegetation ecotone until recently and a declining tree growth rate in this ecotone that are attributable to an increase in local air temperature in the past century. Our findings provide evidence of substantial ecological impacts of global climate change in high-mountain regions. Our study suggests that climate change, if left unabated, can cause dramatic changes in the structure, diversity, and landscape composition of high-mountain ecosystems such as the Altai Mountains.

Methods
Data sampling and wood core collection. After identifying a representative slope of the treeline ecotone in the mountain chain and checking from historical information and local people the absence of any relevant anthropogenic influence, a grid composed by 3 vertical (altitudinal) and 4 horizontal (at the same elevation) transects, at a distance of 50 m each, was drawn and placed on a georeferenced satellite map of the area in order to include the whole ecotone from the timberline, through treeline, to the isolated and highest individual trees of the tree species line (Fig. 1), whose definitions are conventions for communication and do not deserve a major scientific debate 49 . This sampling frame was subdivided in a 5 × 5 m grid in order to collect landscape data from the areas included in the frame but excluded from the sampling transects (Fig. 2). The coordinates of each transect and corners of the sampling frame were recorded with a 2-m precision.
Along each transect all trees with a diameter at the breast height (d.b.h.) >1 cm were identified at the species level, measured with a laser hypsometer-dendrometer (height, dbh and crown width) and georeferenced within the sampling grid ( Fig. 2 and Table 1). Seedlings and saplings of trees (below the minimum height of 1.35 m and d.b.h. of 1 cm) were not measured but have been mapped and included in the category "Trees" (which represents all the non-sampled trees outside the transects) in Fig. 2 and Table 1. The coverage of the other vegetation composed by shrub and herb species (identified at the genus level) and of the main landscape features (big rocks and landslide areas) included in the sampling frame among transects was visually estimated and georeferenced within each cell of the grid (Fig. 2 and Table 1).
At each intersection of the 3 vertical and 4 horizontal transects, the two closest trees were sampled with Haglöf increment borers in a cross design in order to collect two whole wood cores from each tree. This yielded to a collection of 103 individual trees from which 48 dendrochronological samples were extracted.

Aerial photographs and 3-D modelling.
In order to create a 3-D model of the study site, which could represent a reference for future studies of this remote area, 10 UHD aerial photographs taken by a professional drone from different angles where combined in order to create a three-dimensional mesh (Meshlab and Sketchup software) that allowed the reconstruction of 3-D aerial model. Then we georeferenced the 3-D model overlapping it to the 2-D Google Earth satellite image (Fig. 1). Unfortunately, no historical aerial photos of this treeline were available at both the local station and Tomsk State University.
Biometric (allometric) data analysis. We tested a possible correlation between altitude and biometric (allometric) data with the Spearman's rank correlation index at a significance level of α = 0.01 (R Development Core Team 2018). Then, we plotted a correlogram to show the correlation trends (R Development Core Team 2018).
The summary values per altitudinal ranges of the biometric data for all collected trees is shown in Table 1.
www.nature.com/scientificreports www.nature.com/scientificreports/ Dendrochronological data analysis. The 48 wood cores collected were prepared following the common dendrochronological protocols and measured with a LINTAB measuring table with an accuracy of 0.01 mm, equipped with a Leica MS5 stereoscope. The tree-ring growth series were then visually (TSAPwin software, version 0.53 63 and statistically (COFECHA 64,65 cross-dated within trees and between trees (of the same site) for avoiding dating errors in the dataset. Then, an annual mean growth trend was calculated averaging the individual tree annual growth and a variance (SD) was attached (Fig. 4). The basal area of each sampled tree was also estimated.
We then detrended the tree-ring series (Fig. 5) through the estimation and removal of the tree's natural biological growth trend, and standardized the detrended values by dividing each series by the growth trend to produce units in the dimensionless ring-width index (RWI). We used the most adopted methods available for detrending trough the package 'detrendeR' in R Development Core Team 2018. We implemented are a smoothing spline detrending via ffcsaps (method = "Spline"), a modified negative exponential curve (method = "ModNegExp"), a simple horizontal line (method = "Mean"), an "Ar" approach, and a "Friedman" approach.
The "Spline" approach uses a spline where the frequency response is 0.50 at a wavelength of 0.67 * "series length in years". This attempts to remove the low frequency variability that is due to biological or stand effects.
The "ModNegExp" approach attempts to fit a classic nonlinear model of biological growth of the form f(t) = a exp(b t) + k, where the argument of the function is time 66 .
We also checked if a suitable nonlinear model could not be fitted (function is non-decreasing or some values are not positive) and we fitted a linear model. The "Mean" approach fits a horizontal line using the mean of the series. This method is the fallback solution in cases where the "Spline" or the linear fit (also a fallback solution itself) contains zeros or negative values, which would lead to invalid ring-width indices.
The "Ar" approach is also known as prewhitening where the detrended series is the residuals of an Ar model divided by the mean of those residuals to yield a series with white noise and a mean of one. This method removes all but the high frequency variation in the series and should only be used as such 61 .
These methods are chosen because they are commonly used in dendrochronology. There is a rich literature on detrending and many researchers are particularly skeptical of the use of the classic nonlinear model of biological growth (f(t) = a exp(b t) + k) for detrending 67 .
Finally, we detrended our time series with the "Friedman" approach that uses Friedman's 'super smoother' as implemented in supsmu (R Development Core Team 2018).
We tested the significance of the time-series with a Mann-Kendall method, which is an effective test to detect the long-term change in time series 68 . In this study, this method was applied to detect the long-term trend change of the mean annual growth of the sampled trees. Moreover, we checked the possible autocorrelation on our dendrochronological time-series (with an autocorrelation function, acf implemented R Development Core Team 2018). The resulting plot (Fig. 9) shows a significant correlation at lag 1 (≈1), lag 2 (≈0.5) and lag 13 (≈0.4) that decreases after a few lags. This pattern indicates an autoregressive term. Therefore, we used the partial autocorrelation function to determine the order of the autoregressive term (Fig. 10).
After accounting for the partial autocorrelation function we derived an additional curve based on this tree ring index (Fig. 4), which, however, does not differ much from the curves detrended with different methods (particularly from the Friedman's curve). Most of the growth peaks and the general growth trend remained unaltered.
Finally, in order to evaluate the bias of our time series and the level of its accuracy (defined in terms of standard error) to sample estimates, we performed a block bootstrap analysis (tsboot package in R). This technique allowed an estimation of the bias of our distribution using random sampling methods.
We used a fixed block length of 1 to account for each lag of the time series and the results show a bias of −0.32 and a very low standard error of 0.06. We then plotted the histogram of the block bootstrapping method (Fig. 11), which shows that the mean difference of bootstrapping is well below the observed difference (dashed vertical line) in the time series of the mean annual growth (with the calculated bias of −0.32). The Q-Q plot of the 500 random bootstrapping with replacement shows that our time series has a normal distribution (Fig. 11). www.nature.com/scientificreports www.nature.com/scientificreports/ Dendrochronology vs. altitude, year of tree birth, and climate data analysis. We tested a possible correlation between our dendrochronological time series (mean annual growth and basal area) vs. altitude, year of tree birth, and climatic variations with the Kendall's correlation index at a significance level of α = 0.01 (R Development Core Team 2018).
Because the relationship between the year of tree birth vs. both the altitude and the mean annual growth was non-linear, we applied a generalized linear model (glm package in R) and we considered a quadratic function and the related determination coefficients to describe these curvilinear relationships, in the mathematical form of: To disentangle the curvilinear relationship between the mean annual growth and the year of tree birth, we analyzed the growth trends at the rarefied minimum age of all our core samples (12 years). In this way, we were able to compare the growth of the trees grouped per current age and altitude (4 groups at 2,150, 2,200, 2,250, 2,300 m a.s.l.) during their "youth", in their first 12 years of life and to evaluate their growth variability.
Although some mean precipitation and temperature observations for the whole Altai region have been recently published 31 , we used the longest climatic time series available in the literature 29 , with a time span of 53 years ). This was recorded from a meteorological station (that in Kosh-Agach) located at a comparable altitude (≈1,900 m a.s.l.) and location (linear distance ≈50 km) with respect to our study site.
Since, from the data in the literature, there is no evidence of a significant variability in the precipitation trends 29,31 , we considered only temperature records.
The mean annual temperature variability was correlated with the mean annual growth variability and a generalized linear model (quadratic function) was used to describe the relationship. Because the trend resulted in an almost hyperbolic curve (quadratic function), we inspected in any factor related to the predictor variable (temperature) could have influenced it. We discovered that a threshold of −4.3 °C was never passed before the 1978 and that this temperature approximately represents the turning point (close to the minimum point x = −b/2a of the eq. 1) of the curvilinear relationship between the mean annual growth and the mean annual temperature.
Therefore, to derive linear equations that would allow an estimation of the temperature during the missing years of observation from our dendrochronological time series, we divided the temperature-growth plot in two subplots (at −4.3 °C threshold). The parameter from the linear regression equation derived from the data points  www.nature.com/scientificreports www.nature.com/scientificreports/ at a temperature <−4.3 °C were used to estimate the temperature for the period 1902-1955 (because of we show that the temperature was always below −4.3 °C before 1978) and those derived from the data points at a temperature >−4.3 °C were used to estimate the temperature for the period 2010-2017 (because the data points >−4.3 °C exclusively represents most of the growth trends after 1978).
Following a common approach to reconstruct temperatures from tree-ring records [69][70][71] , the stability and reliability of the regression Eq. 1 were assessed using the split sample method. We performed a validation by calibrating climate data from a sub-period by dividing the data into two parts: 1956-1978 and 1979-2009) and verifying the reconstruction using the remaining data. The results were evaluated by the correlation coefficient (r), the sign test (ST), the reduction of error test (RE), the coefficient of efficiency (CE) and the product means test (t) during the verification period. The actual and estimated data in year i of the verification period and the mean of the actual data in the calibration and verification periods were compared to calculate the sign of RE and CE following refs 72,73 , respectively.
We displayed the reconstructed curves of temperature variability and compared their general trends and selected relevant periods with the available information in the literature (see the Discussion section).