Upward Treeline Shifts in Two Regions of Subarctic Russia Are Governed by Summer Thermal and Winter Snow Conditions

: Climate warming impacts on alpine treeline dynamics. However, we still lack robust assessments of the long-term impacts of climate on tree recruitment at the treeline, particularly in remote areas such as the subarctic regions of Russia subjected to different climate inﬂuences. We expected that the treelines in two regions may have different features and dynamics patterns. We analyzed climate variables and assessed treeline dynamics by quantifying recruitment using the tree rings of ca. 7000 trees of four species ( Betula pubescens Ehrh. ssp. tortuosa , Pinus sylvestris L., Picea abies Ledeb. ssp. obovata , Larix gmelinii Rupr.) along 14 altitudinal transects (series of study plots). We compared the Khibiny Massif (Kola Peninsula) and the western Putorana Plateau, subjected to oceanic and continental inﬂuences, respectively. In both regions, summers became warmer, and winters became snowier during the past century. At the low part of the treeline ecotone, tree recruitment has slowly increased since the mid-18th century at the Putorana Plateau and the mid-19th century at the Khibiny but accelerated in the early 20th century at both regions and reached a maximum peak in the second half of the past century. Treeline encroachment intensiﬁed in the 1930s at the Khibiny and the 1950s at the Putorana Plateau. Trees encroached in the tundra leading to upward treeline shifts in the late 20th century. The slope exposure affected the rates of treeline shift with higher upward advances on southern-oriented slopes. Tree recruitment and early-winter precipitation were positively correlated. The differences in species composition, treeline altitude and inﬂuences of slope orientation on treeline dynamics can be explained primarily by differences in the degree of continentality. The abundance of saplings in both regions allows the future encroachment of trees into tundra and further treeline upward shifts to be forecast.


Introduction
The warming of the climate system, especially since the 1950s, is unprecedented on a scale of decades to millennia and has particularly impacted cold regions across the Arctic [1]. It is known that with the observed climate change, many plant species would shift polewards or upwards [2,3]. One of the most notable manifestations of the terrestrial biota response to climate warming is the advancement of woody vegetation into the alpine [4] and Arctic [5,6] ecosystems, which constitutes an example of ecotone dynamics. Boundaries between different plant communities such as treeline ecotones are subjected to constant changes in space and time [3,7]. This is due to the fact that the altitudinal and latitudinal position of the treeline depends primarily on temperature and, therefore, its position is highly sensitive to such climate variables [4]. It is known that repeated fluctuations of the alpine and polar treelines and forest limits occurred during the Holocene in response to temperature fluctuations [8][9][10][11].
Nowadays, treeline shifts have been observed across almost all continents [12,13], but they did not linearly respond to climate warming and some were uncoupled from temperature changes. Here we address this issue by investigating how treelines in two climatically contrasting Russian subarctic regions (maritime vs. sharply continental) respond to climate conditions during the cold season, particularly to winter snowfall, which could represent a major, but underestimated, driver of treeline dynamics. Such high-altitude, remote territories located in subarctic regions have been little disturbed by local anthropogenic use (grazing, logging) and there, tree growth and regeneration are very sensitive to temperature changes [6,14,15]. For instance, in subarctic regions of Alaska, northern Canada and the Scandinavian Alps, treelines have positively responded to warmer growing-season conditions by showing enhanced tree growth and recruitment [16][17][18][19]. In subarctic Russia, similar studies were mostly carried out in the Polar Urals [20][21][22][23][24], with a few studies done on the Putorana [25,26] and the Anabar Plateau [27]. Recently, in subarctic Russian regions such as the Khibiny mountain range [28] and the Putorana Plateau [29], an intense upward shift of the treeline was observed during the 20th century which was modulated by slope exposure. We also investigated this aspect in more detail in these two subarctic regions subjected to lower (Khibiny mountain range) and higher continentality (Putorana Plateau). In this work, we perform a comparative analysis of treeline ecotone dynamics along altitudinal gradients located on slopes of various exposures. We reconstruct these dynamics over the past three centuries to infer the underlying climatic drivers of the recent and intense treeline shifts. We compared those responses in two remote regions subjected to contrasting climatic influences, on the western side-the Khibiny Mountains located in the Kola Peninsula, and the eastern side-the Putorana Plateau situated in the northwest of Central Siberia. The two study regions are separated by ca. 2200 km.
Our aims were: (1) to reconstruct the tree stands formation within the forest-tundra transition by determining the size and age of about 7000 trees of four species (Betula pubescens Ehrh. ssp. tortuosa, Pinus sylvestris L., Picea abies Ledeb. ssp. bovate, Larix gmelinii Rupr.) on 14 altitudinal transects on the slopes, differing in exposure and dominant tree species; (2) to link these changes with climatic data from weather stations situated near the study areas; (3) to assess the links between changes in winter snow cover depth and air temperature and treeline dynamics (recruitment) along the altitudinal gradient; and (4) to compare the responses among regions, tree species, altitudes and slope exposures. Our hypotheses are: (1) the species composition, altitude position and rate of treeline shift in the two study regions would differ and depend on the degree of continentality; (2) within each region, the treeline shift would depend on altitude and slope exposure which determine temperature and solar radiation, i.e., heat supply; and (3) upward tree expansion would be closely related to climate change but also to site microclimatic conditions.

Study Sites
The Khibiny Mountains are located in the central part of the Kola Peninsula, its length is about 45 km from north to south and about 50 km from west to east; they are situated about 150 km north of the Arctic Circle ( Figure 1, Table 1). The territory of the mountain range is included in the Atlantic-Arctic zone, which is characterized by frequent influxes of warm air masses from the Atlantic Ocean and more intense Atlantic activity in winter and during transitional seasons than in summer. The climate in the Khibiny is mitigated by the proximity of the Barents Sea, warmed by the warm current of the Gulf Stream. The average annual temperature is 0.3 • C (Figure 2a). The average annual precipitation is 538 mm, with snow precipitation being 34% of the total. The region is impacted by western winds with average speeds of 8.5 m s −1 and 6.5 m s −1 in winter and summer, respectively. The Khibiny Massif has a complex geological structure with prevailing plutonic rock and dissected terrain with numerous glacial cirques and steep slopes (9-34 • ) [30]. Permafrost is not formed in this region. Soils are acidic, stony, loam Podzolized Rankers, which stay mainly unfrozen (0.1-0.6 • C) during the cold season on the treeline [28]. In the Khibiny Massif, birch (Betula pubescens Ehrh. ssp. tortuosa) and Scots pine (Pinus sylvestris L.) are the most widespread species at the forest limit, and Norway spruce (Picea abies Ledeb. ssp. obovata) is less common.   The Putorana Plateau is located in northwestern Central Siberia (Table 1). This is the largest monolithic mountain range situated in the Russian Arctic, almost entirely located to the north of the Arctic Circle. In geological and geomorphological terms, it is a basalt crystalline massif (plateau) with flat tops, raised in average of 900-1200 m above sea level [31]. With repeated uplifts, deep radial tectonic splits created gorges and canyons. The Putorana Plateau is located in the subarctic climatic belt at the junction of the Atlantic and Siberian regions, and it is characterized by a sharply continental climate [32]. The average annual temperature is subzero (−9.8 • C) (Figure 2b). The average annual precipitation is 635 mm, snow being 56% of the total. The region is characterized by a complex wind regime with a predominance of eastern winds with average speeds of 1.2 m s −1 . Although the Putorana Plateau is surrounded by terrains with continuous permafrost, on its steep (24-35 • ) and well-drained slopes permafrost is not found. Stony, acidic, sandy loam Skeletic Leptosols (Turbic) are found at the treeline sites, which may be frozen during the cold season [29]. In this region, larch (Larix gmelinii Rupr.) is the dominant tree species at the forest limit.

Studied Treelines
The present studies were carried out in the ecotone of the elevational treeline. It is a complex element of vegetation cover unique to each slope [20]. We studied it taking into account the following classification of vegetation along the altitudinal gradient considering four levels (following [7]): (1) the treeline or sparse tree stands line, which is a theoretical line connecting the uppermost individual trees or tree groups separated by 20-60 m on average; (2) the open forest line, which is the theoretical line connecting the uppermost forest patches separated by 7-30 m on average (Figure 3a); (3) the closed forest line or upper boundary of closed forests, where forest patches merge, forming more or less continuous stands separated by 2-7 m; and (4) the closed forest (50-80%) at 15-70 m altitude below the closed forest line (which, in many cases, coincides with the altitudinal position of the closed forest line in the 1950-1960s, depicted on historical topographic maps). Forest patches were defined as vegetation areas where the height of the trees exceeds the height of the shrub layer by 2-3 times, the crown cover is at least 10%, and the average diameter of the forest patches is at least 5 times the average height of the tree layer. The treeline usually corresponds as is often defined as the highest elevation where trees at least 2 m tall can be found [3]. The undergrowth includes trees no more than 1.5 m high and no more than 40 years old. If the specimen under study was less than 1.5 m, but over 40 years old, then it was considered a tree.

Field Sampling
From 2017 to 2019, sampling was carried out in three areas of the Khibiny Mountains: (1) Imandra (northwestern macroslope), (2) Maly Vudayvr (central part), and (3) Kitchepakh (southeastern macroslope). In total, 9 altitudinal gradients or transects were installed ( Table 2). On the Putorana plateau, studies were carried out in only one areaon the Sukhiye Gory massif. This is due to the extreme inaccessibility and remoteness of the Putorana Plateau. Four altitudinal transects were set up, two of which were situated on slopes with southern exposure (see also [29]). Before setting the transects, we tried to

Field Sampling
From 2017 to 2019, sampling was carried out in three areas of the Khibiny Mountains: (1) Imandra (northwestern macroslope), (2) Maly Vudayvr (central part), and (3) Kitchepakh (southeastern macroslope). In total, 9 altitudinal gradients or transects were installed ( Table 2). On the Putorana plateau, studies were carried out in only one area-on the Sukhiye Gory massif. This is due to the extreme inaccessibility and remoteness of the Putorana Plateau. Four altitudinal transects were set up, two of which were situated on slopes with southern exposure (see also [29]). Before setting the transects, we tried to select slopes that were approximately equal in steepness and stone cover. At each altitudinal level, 2-5 square plots with 20 m sides were laid down, on which we determined the exact location of each tree trunk, the diameter at the base and breast height, and the horizontal projections of the crown diameter in two perpendicular directions. To determine the age of trees, cores were taken near ground level and up to 30 cm from the soil level using Pressler increment borers (Haglof, Langsele, Sweden). The age of small or shrubby trees (<3 cm) was determined by taking a section from the base of the trunk from every third specimen using a hand saw.
In the laboratory, all cores taken were fixed to wooden strips. The cores and wood sections were cleaned with a paper knife and a razor blade, and they were carefully sanded to inspect the rings. All samples were visually cross-dated, and their ring widths were measured on a semiautomatic device Lintab 5 (RinnTech, Heidelberg, Germany) with 0.01 mm accuracy. Then, the visual cross-dating was checked using the Cofecha software [33].
To identify false or missing rings, a tree-ring width chronology or mean series was built using at least 40 cores per site, specially taken from old trees in the study area. If cores did not reach the pith, then the missing rings were calculated using a theoretical pith finder based on a transparent film with lines of concentric circles of different sizes which was matched to the innermost curved rings of cores. Using the saplings' height and age, a regression equation was calculated between these two variables, with the help of which corrections were calculated to determine the age of trees with diameters >3 cm. In total, we determined the age of about 7000 trees (see details in Table 1). It should be noted that trees growing on all surveyed altitudinal gradients have not been exposed to fires or other disturbances, since we did not find traces of fires (fire scars, soil charcoal). In addition, when setting study plots, we did not observe many dead trees or indications of massive damage by insects or other pests.

Climate Data and Statistical Analysis
For the Khibiny Massif, the data of the Kandalaksha station (67 • 10 09 N, 32 • 21 15 E, 25 m a.s.l., located 80 km away from the study area) were analyzed considering average monthly air temperature (data from 1912 to 2016), total monthly precipitation (data from 1936 to 2015), the snow depth (data from 1938 to 2013) and the total monthly sunshine (data from 1961 to 2015). For the area of the Putorana Plateau, the data of the Dudinka meteorological station (69 • 24 00 N, 86 • 10 00 E, 19 m a.s.l., located 170 km away from the study area) were analyzed considering average monthly air temperature (data from 1906 to 2012), total monthly precipitation (data from 1942 to 2010), the snow depth (data from 1937 to 2011) and the total monthly sunshine (data from 1957 to 2012). All station data were taken from the specialized Roshydromet database (http://meteo.ru/data, accessed on 2 November 2021). The data on the total amount of precipitation were adjusted with corrections for wetting and changing instruments [34,35]. The depth of the snow cover was analyzed from the beginning (31 December) to the end (31 March) of the cold period.
For the analysis of climatic data, the periods of warm and cold seasons were defined as July and October-April for Dudinka, and July and November-April for Kandalaksha. The definition of the warm period corresponded to the phase of the most active growth of vegetation in the study areas. The cold period included months with an average air temperature below 0 • C and stable snow cover. Anomalies of climatic parameters in the warm and cold periods of each year were determined through the calculation of deviations of the current value with respect to the average of the reference period 1961-1990. Linear regressions were fitted to assess trends in data on changes in climatic anomalies, and their coefficients of determination (R 2 ) were calculated.

In Situ Air Temperature and Snow Measurements
The measurement of the near-surface air temperature along the altitudinal gradients was carried out using dataloggers (Maxim iButttons, model DS1921G; San Jose, CA, USA). Two or three loggers were placed in the crown of trees at a height of at least 2.0 m. Temperature registration was carried out every 4 h. In total, 81 loggers were installed in the Khibiny sites and 48 loggers on the Putorana Plateau sites.
The snow cover depth was measured in the months of maximum snow accumulation for the study area: at the end of March in the Khibiny mountain range and at the end of April on the Putorana Plateau. The snow cover height was measured using a special metal avalanche probe tool with marks every 1 cm. On each test plot, the snow height was measured in at least 40 points.

Statistical Analyses
We analyzed the relationships between the number of recruited trees, considering five-year age classes, and the corresponding five-year average values of climatic parameters for the current and previous periods at the warm and cold periods, and also at the beginning of the cold period for each study area. We calculated average air temperature, average snow cover depth, average sunshine duration and total precipitation. Spearman correlation coefficients (Rs) were calculated since the data distribution differed from a normal distribution according to Shapiro-Wilk tests.
The usage of five-year classes allows us to consider uncertainties in age estimation and also the fact that masting is frequent in some of the study species, such as spruce, which could lead to episodic recruitment peaks [36]. In addition, when a tree forms a multi-stemmed form of growth, such as in the case of mountain birch, difficulties arise in determining the age due to the fact that most of the trunks or sprouts begin to form much later than the time of the appearance of the oldest one forming the basal stump [37].

Microclimatic Conditions
The data indicate that microclimatic conditions and snow cover depth differ both within the same slope and between regions (Tables S1 and S2). A slight decrease in air temperature was observed as altitude increased. In the Khibiny mountain range, the lowest temperatures were recorded in the central part on Maly Vudayvr, where the forest boundary was located at the highest elevation (Table 1). Comparisons of the two regions showed that the temperature on the Putorana Plateau in the forest-tundra ecotone was 1.5-2 times lower than in the Khibiny Massif, indicating a higher continentality in the eastern region. Snow cover depth in both study areas generally decreased upwards. The smallest snowpack depth within the forest-tundra ecotone was recorded on the eastern slope of the Sukhie Gory massif (from 132 ± 3 cm at the lowest-altitude level to 34 ± 2 cm at the highestaltitude level) and on the southern slope of the Imandra massif (from 62 ± 20 cm at the lowest-altitude level to 24 ± 17 cm at the highest-altitude level). In general, the snowpack depth was approximately 1.5 times higher in the Khibiny Mountains than in the Putorana Plateau sites, which confirms a higher winter precipitation in the western region. In the Khibiny Mountains, snowdrifts up to 3 m long were observed in some places.

Morphological Features and Age Structures
We found a mean 2-5-fold decrease in the average and maximum characteristics of measured variables (stem diameter and height, crown diameter and cover, stand basal area, tree age) as we moved from closed forests to the uppermost individual trees located in the tundra (Tables S3-S6). The tree density and crown area further up decreased 5-50-fold. Thus, the higher the location, the smaller, younger and sparser the trees were.
The study of age structures on Imandra mountain (Khibiny Massif) shows that the first trees were recruited at the beginning of the 19th century ( Figure 4). On southern-exposed slopes, most Scots pine recruitment occurred during the second half of the 20th century, and this happened at all altitudinal levels. At the highest levels of 1, 2, and 3 (from the treeline to the closed forest limit), this process took place mainly after the 1970s. On northern-exposed slopes, recruitment was dominated by birch which began to appear earlier than on the southern-exposed slopes, in the 1850s. Birch recruitment was progressive, without obvious peaks in the number of trees. At the highest altitudinal level (species line), birch trees appeared after the 1930s. the snowpack depth was approximately 1.5 times higher in the Khibiny Mountains than in the Putorana Plateau sites, which confirms a higher winter precipitation in the western region. In the Khibiny Mountains, snowdrifts up to 3 m long were observed in some places.

Morphological Features and Age Structures
We found a mean 2-5-fold decrease in the average and maximum characteristics of measured variables (stem diameter and height, crown diameter and cover, stand basal area, tree age) as we moved from closed forests to the uppermost individual trees located in the tundra (Tables S3-S6). The tree density and crown area further up decreased 5-50fold. Thus, the higher the location, the smaller, younger and sparser the trees were.
The study of age structures on Imandra mountain (Khibiny Massif) shows that the first trees were recruited at the beginning of the 19th century ( Figure 4). On southernexposed slopes, most Scots pine recruitment occurred during the second half of the 20th century, and this happened at all altitudinal levels. At the highest levels of 1, 2, and 3 (from the treeline to the closed forest limit), this process took place mainly after the 1970s. On northern-exposed slopes, recruitment was dominated by birch which began to appear earlier than on the southern-exposed slopes, in the 1850s. Birch recruitment was progressive, without obvious peaks in the number of trees. At the highest altitudinal level (species line), birch trees appeared after the 1930s. In the Maly Vudayvr (central part of the Khibiny Massif), the process of tree recruitment followed a different scenario ( Figure 5). The oldest birch trees found were recruited in the middle of the 19th century, on the lower-altitude slopes with eastern, southeastern and southwestern exposures. On the northeastern slope, the settlement of trees began at all altitudes in the 20th century. On the slopes with southeastern and southwestern exposures, birch trees actively recruited in the second half of the 20th century. In general, the data show that on all the studied slopes of the Maly Vudayvr mountains the most active birch expansion across the treeline ecotone occurred in the second half of the 20th century. In the Maly Vudayvr (central part of the Khibiny Massif), the process of tree recruitment followed a different scenario ( Figure 5). The oldest birch trees found were recruited in the middle of the 19th century, on the lower-altitude slopes with eastern, southeastern and southwestern exposures. On the northeastern slope, the settlement of trees began at all altitudes in the 20th century. On the slopes with southeastern and southwestern exposures, birch trees actively recruited in the second half of the 20th century. In general, the data show that on all the studied slopes of the Maly Vudayvr mountains the most active birch expansion across the treeline ecotone occurred in the second half of the 20th century. The age structures of trees growing on the slopes of the Kitchepakh mountain, in the eastern Khibiny Massif, showed that birch also actively recruited in the 20th century, although the oldest individuals were dated to the late 19th century ( Figure 6). The birch mostly populated the eastern-exposed slopes in the second half of the 20th century. Norway spruce was mainly found at low-altitude levels and appeared at the beginning of the 19th century as single trees. This recruitment process of Norway spruce expansion took place mostly throughout the 20th century, mainly on the eastern and southeastern slopes, with some individuals reaching the highest altitude in the treeline (level 1). The age structures of trees growing on the slopes of the Kitchepakh mountain, in the eastern Khibiny Massif, showed that birch also actively recruited in the 20th century, although the oldest individuals were dated to the late 19th century ( Figure 6). The birch mostly populated the eastern-exposed slopes in the second half of the 20th century. Norway spruce was mainly found at low-altitude levels and appeared at the beginning of the 19th century as single trees. This recruitment process of Norway spruce expansion took place mostly throughout the 20th century, mainly on the eastern and southeastern slopes, with some individuals reaching the highest altitude in the treeline (level 1). On the Putorana Plateau, the first larch trees were recruited at lower altitudes throughout the 17th (data not shown) and 18th centuries (Figure 7). Their most massive expansion took place on the slopes of the western, southern and eastern exposures in the 19th and 20th centuries. On the northern-exposed, low-altitude slopes, the main recruitment episode of larch occurred in the second half of the 19th century. At higher altitudes, the most massive larch regeneration phase took place after the 1800s. At the higher altitudes on the slopes of southern, western and eastern exposures, the main regeneration peak was observed in the second half of the 20th century. On the Putorana Plateau, the first larch trees were recruited at lower altitudes throughout the 17th (data not shown) and 18th centuries (Figure 7). Their most massive expansion took place on the slopes of the western, southern and eastern exposures in the 19th and 20th centuries. On the northern-exposed, low-altitude slopes, the main recruitment episode of larch occurred in the second half of the 19th century. At higher altitudes, the most massive larch regeneration phase took place after the 1800s. At the higher altitudes on the slopes of southern, western and eastern exposures, the main regeneration peak was observed in the second half of the 20th century.

Analysis of Changes in Climatic Conditions
The analysis of climate data from meteorological stations shows that in both areas there was a long-term change in climatic conditions ( Figure 8). The largest temperature drops occurred during the warm season, and the accumulation of precipitation was higher during the cold season. In the warm period, significant increases in precipitation were noted for the Kandalaksha weather station (Khibiny Massif) by +4.9 mm decade −1 (R 2 = 0.54, p < 0.001) and also in the average air temperature for the Dudinka meteorological station (Putorana Plateau region) by +0.16 °C decade −1 (R 2 = 0.21, p = 0.04). The cold-period precipitation anomaly significantly increased in the Kandalaksha weather station by +14.6 mm decade −1 (R 2 = 0.56, p < 0.001) and also in the Dudinka meteorological station by +19.2 mm decade −1 (R 2 = 0.25, p = 0.07). The snow cover depth at the beginning of the cold period increased significantly at both meteorological stations, in particular for the Kandalaksha meteorological station by +3.3 cm decade −1 (R 2 = 0.56, p < 0.001) and also for the Dudinka meteorological station by +5.2 cm decade −1 (R 2 = 0.32, p = 0.03).

Analysis of Changes in Climatic Conditions
The analysis of climate data from meteorological stations shows that in both areas there was a long-term change in climatic conditions ( Figure 8). The largest temperature drops occurred during the warm season, and the accumulation of precipitation was higher during the cold season. In the warm period, significant increases in precipitation were noted for the Kandalaksha weather station (Khibiny Massif) by +4.9 mm decade −1 (R 2 = 0.54, p < 0.001) and also in the average air temperature for the Dudinka meteorological station (Putorana Plateau region) by +0.16 • C decade −1 (R 2 = 0.21, p = 0.04). The coldperiod precipitation anomaly significantly increased in the Kandalaksha weather station by +14.6 mm decade −1 (R 2 = 0.56, p < 0.001) and also in the Dudinka meteorological station by +19.2 mm decade −1 (R 2 = 0.25, p = 0.07). The snow cover depth at the beginning of the cold period increased significantly at both meteorological stations, in particular for the Kandalaksha meteorological station by +3.3 cm decade −1 (R 2 = 0.56, p < 0.001) and also for the Dudinka meteorological station by +5.2 cm decade −1 (R 2 = 0.32, p = 0.03). Overall, the most significant expansion of trees took place in the second half of the 20th century in both regions, mainly at high-altitude sites (Figure 9). Overall, the most significant expansion of trees took place in the second half of the 20th century in both regions, mainly at high-altitude sites (Figure 9).

Relationship of Climate Variability with Tree Expansion: Roles of Altitude and Slope Exposure
For the Khibiny Massif, high positive correlation coefficients were found between tree recruitment (birch and Scots pine) and precipitation in July and during the beginning of the cold period (November-December), and also with the snow depth at the beginning of the cold period (see Table S7). The number of recruited Scots pine individuals from 1910 to 2010 on the Imandra mountain was positively related to the depth of snow at the end of December for the five years preceding tree recruitment (Rs = 0.56, p = 0.04). The relationship between the number of recruited birch individuals in the Maly Vudayvr site and the depth of snow at the beginning of the cold period of the corresponding five years also showed a high correlation coefficient (Rs = 0.72, p < 0.01). On these two sites of the Khibiny Massif (Imandra mountain and Maly Vudayvr sites) there was a general tendency of stronger responses of tree recruitment to the aforementioned climate variables at the two high-altitude levels, i.e., in the treeline and the open forest line. Indeed, the highest correlation coefficients of tree recruitment with snow depth at the beginning of the cold period were found at the highest level of the Imandra mountains (Rs = 0.77, p < 0.01) for the previous five-year period, and also at the highest level of the Maly Vudayvr (Rs = 0.83, p < 0.001) for the corresponding five-year period. We also found negative relationships of recruitment with the depth of snow cover at the beginning of the cold period at the lowest

Relationship of Climate Variability with Tree Expansion: Roles of Altitude and Slope Exposure
For the Khibiny Massif, high positive correlation coefficients were found between tree recruitment (birch and Scots pine) and precipitation in July and during the beginning of the cold period (November-December), and also with the snow depth at the beginning of the cold period (see Table S7). The number of recruited Scots pine individuals from 1910 to 2010 on the Imandra mountain was positively related to the depth of snow at the end of December for the five years preceding tree recruitment (Rs = 0.56, p = 0.04). The relationship between the number of recruited birch individuals in the Maly Vudayvr site and the depth of snow at the beginning of the cold period of the corresponding five years also showed a high correlation coefficient (Rs = 0.72, p < 0.01). On these two sites of the Khibiny Massif (Imandra mountain and Maly Vudayvr sites) there was a general tendency of stronger responses of tree recruitment to the aforementioned climate variables at the two high-altitude levels, i.e., in the treeline and the open forest line. Indeed, the highest correlation coefficients of tree recruitment with snow depth at the beginning of the cold period were found at the highest level of the Imandra mountains (Rs = 0.77, p < 0.01) for the previous five-year period, and also at the highest level of the Maly Vudayvr (Rs = 0.83, p < 0.001) for the corresponding five-year period. We also found negative relationships of recruitment with the depth of snow cover at the beginning of the cold period at the lowest altitude level, both for the corresponding five-year period (Rs = −0.72, p < 0.01) and for the previous five-year period (Rs = −0.84, p < 0.001).
An analysis of slopes with different exposures in Khibiny Massif sites such as Mount Imandra showed a significant relationship between the number of recruited pine trees and the depth of snow at the beginning of the cold period for the previous five years on the southern slopes (Rs = 0.56, p = 0.04). For the same indicator, the highest correlation coefficients between the number of birch trees recruited on the Maly Vudayvr site were Rs = 0.82 (p < 0.001) and Rs = 0.68 (p < 0.01) for the southwestern and southeastern slopes, respectively, and considering the corresponding five-year period. On Mount Kitchepakh, predominantly negative relationships (particularly for precipitation in July) with birch recruitment were noted, but they were not significant.
Correlation analysis showed that on the Putorana Plateau (Table S8)  In the Putorana Plateau, larch recruitment was related to precipitation at the beginning of the cold period (October-December) of the previous 5-year period (Rs = 0.78, p < 0.01). At the upper and lower altitudinal levels, larch recruitment correlated with precipitation at the beginning of the cold period of the previous five-year period (Rs = 0.84 and R = 0.86, respectively; p < 0.001 in both cases).
Regarding different exposures in this region, the highest correlations between larch recruitment and July temperature in the previous five-year period were found for the western (Rs = 0.60, p = 0.01) and eastern (Rs = 0.43, p = 0.07) slopes. The strongest relationships were observed between larch recruitment and precipitation in the first months of the cold period of the previous five-year period for the eastern (Rs = 0.64, p = 0.02), southern (Rs = 0.67, p = 0.01) and western (Rs = 0.47, p = 0.11) slopes. Similar associations were observed for the snow cover depth at the beginning of the cold period (determined at the end of December) and considering the corresponding five-year periods, with correlation coefficients of Rs = 0.54 (p = 0.04), Rs = 0.72 (p < 0.01) and Rs = 0.59 (p = 0.02) for the eastern, southern and western slopes, respectively.

Treeline Characteristics and Slope Exposure
Two subarctic regions of Russia, subjected to contrasting degrees of continentality, presented upward treeline shifts during the 20th century. Analyses of the age structures of tree stands on the upper limit of their growth indicated a shift at the upper boundary of the sparse stands and open and closed forest lines on the more oceanic mountains of the Kola Peninsula (Khibiny) and also on the more continental massifs of the Putorana Plateau.
At each site and along the altitudinal gradient, tree size and age decreased depending on the exposure of the slope, which also affected the maximum altitude limit of stands and the rates of treeline expansion. The highest treeline positions and the largest increase in tree recruitment occurred on southern-oriented slopes. For example, in the Khibiny Massif, the modern open forest treeline reaches its highest altitude on slopes of southwestern exposure, where it is 123 m higher than on northeastern-oriented slopes [28]. In this area, on the slopes of southern exposures, its maximum shift (99-107 m) has been recognized over the past 60 years. On the Putorana Plateau (Sukhie Gory massif), the most significant altitudinal shifts of the open forest boundary over the past 60 years were found on the slopes of southern (111 ± 74 m) and western exposures (86 ± 62 m) [29].
In some areas, depending on the exposure of the slope, there are differences in tree species composition; for example, Scots pine dominates on the drier and warmer southern slopes of Mt. Imandra of the Khibiny Massif, and birch dominates on the wetter and colder northern slopes. This agrees with the importance of the heat supply of the slopes in the cold climate of these cold-stressed subarctic sites. On the slopes oriented in a more southerly direction or closer to this direction, the amount of incoming total solar energy is greater compared to northern-oriented slopes. For example, from the Yukspor meteorological station (Khibiny Massif, 910 m a.s.l.), 34.8 kcal cm −2 are measured on the surface of southern-oriented slopes with an inclination angle of 10 • from June to August, whereas 32.9 and 30.6 kcal cm −2 are measured on western/eastern-and northern-oriented slopes, respectively. With an increase in the angle of inclination of the surface, the difference between the northern and southern slopes also increases up to 6.5, 8.8 and 17.7 kcal cm −2 at 20 • , 30 • and 50 • degrees of inclination, respectively [38]. Direct solar radiation leads to significant heating of the crowns and trunks of trees, and their temperature becomes higher than the air temperature by +7-10 • in the middle of the day and on average by 0.8-2.9 • C per month [29].
Slope exposure can affect: the wind regime [39], the thermal regime during the vegetation period [40], the start of the vegetation period due to greater warming up of tree trunks in spring [41], the start of snow melting in spring [3], as well as the degree of nitrogen mineralization in soil [42]. More solar radiation and warmer conditions can enhance the survival of tree seedlings and saplings and tree development in general [3]. Differences in the average altitude of the upper forest boundary on the eastern and western slopes, which are illuminated by the sun in summer months and warmed up at a relatively similar rate, may be explained because, on the eastern slopes, due to the prevalence of southwestern and western winds during the winter, a deeper and bigger snowpack accumulates. This can lead to a delayed snow melting and thereby to a reduction in the growing season in comparison with western-oriented slopes [29]. In the southern Rocky Mountains (USA), the expansion of trees was significantly greater on the wetter slopes with northern exposure during dry periods [43], where soil moisture was higher [44]. In the eastern Alp treelines, the cambial activity of treeline trees lasted longer on southern than on northern slopes [41]. This was also observed in dry-warm mountains such as the Snake Range (Nevada, USA), where Pinus monophylla grew less and formed smaller tracheids in the driest treeline sites [45]. In mid-latitude forests from northwest British Columbia, shaded northern-oriented slopes also showed lower productivity than those located on the sunnier and warmer southern-oriented slopes [46]. Further, in the southwestern Yukon (Canada), Picea glauca treelines advanced faster on southern-than on northern-oriented slopes [47]. Therefore, it is important to take into account the slope exposure and its effects on tree growth and regeneration when considering models of treeline dynamics [48].

Drivers of Treeline Dynamics
Given that forest boundary shifts have occurred in two distant regions with similar tree expansion periods in the last century, we assumed that the main driver was changes in climatic conditions. Analysis of data from nearby meteorological stations in the study areas showed that in both areas the climate changed in a similar way, namely summers became warmer and winters snowier. Moreover, on the Putorana Plateau, the change in climatic parameters was more pronounced. In the Khibiny Massif, there was an increase in the duration of the vegetation period, especially due to an earlier onset. Between 1961Between -1990Between and 1991Between -2010 there was an increase in the periods with temperatures above 5 • C and 10 • C by 4 and 9 days, respectively [49]. On the Putorana Plateau, the duration of the vegetation period increased by 4-7 days [29].
A number of researchers have linked treeline shifts with climate warming during the warm season because it is the time when trees grow [50,51]. These conclusions are confirmed by the high correlations we identified between the number of recruited larches on the Putorana Plateau site and the average monthly temperatures in July of the preceding 5-year period. It is known that the favorable or unfavorable conditions in the period preceding the appearance of trees determine the number of viable seeds, and the conditions during the period of their appearance and the first years of life affect the number of seedlings and the survival rate of recruits [52]. Globally, upward treeline shifts often occurred in regions where there was a change in the snow regime during winter [53]. Previously, we showed [20] that the main climatic driver of treeline shifts was a change in winter precipitation across the Ural Mountains. Similar conclusions were also obtained for the expansion of Juniperus sibirica Burgsd. [54]. It was also established that ongoing climatic changes lead to the upward expansion of tree stands on the Avam Mountains, located in the northwestern Putorana Plateau [25]. According to data presented by [25], the most massive larch recruitment episodes occurred in the middle and upper parts of the forest-tundra ecotone during the second half of the 20th century, findings which agree with our results. Earlier, based on the comparison of satellite images of the central Putorana Plateau, it was found that the average annual temperature during the last three decades of the 20th century promoted an upward shift of vegetation of about 15 m [26]. Based on a comparative analysis of satellite images, these authors [26] [26].
Analysis of the influence of precipitation on larch recruitment on the Putorana Plateau over a five-year period showed that the greatest influence was exerted by the total precipitation of the cold period of the year (from October to April), especially during the first months (from October to December). Moreover, the strongest connection was established on the eastern slope, where the lowest values of snow cover depth were recorded in the winter season (Table S2). Similar results were obtained from the Khibiny Massif. The strongest associations were found between tree recruitment on all slopes and on all transects with early winter precipitation. The highest correlations were observed on the southern slope of Imandra mountain (Khibiny Massif), where the depth of snow cover is lowest (Table S1). In this site, a high early-winter snow precipitation enhanced pine recruitment (Table S7).
Snow cover and its formation at the beginning of winter are one of the most important factors influencing the survival of tree recruits in the forest-tundra ecotone [3,15]. Tree seedlings and saplings, but also shrubby individuals, are exposed to snow abrasion, needle abrasion, shoot desiccation and frost-induced xylem embolism; therefore, snow protection allows them to mitigate these damages [7,15]. With the formation of a uniform and sufficiently deep snow cover at the beginning of winter, more favorable microclimatic conditions for the understory and tree recruits are created [55]. The deeper the snow cover, the lower the soil freezing temperature and the less damage to the root system [56,57].

Regional Climate and Treeline Features
Differences in the species composition of the study areas, the maximum elevations of the forest boundary and the nature of the slope settlement, can be primarily due to the most pronounced differences in climatic conditions. In the Khibiny Massif, the climate is maritime with oceanic influence since the study sites are located 120-170 km away from the coast of the Barents Sea, which is warmed by the warm Gulf Stream current. On the Putorana Plateau, the climate is sharply continental and study sites are situated 450 km from the coast of the cold Kara Sea. In the Khibiny Massif, the average rate of advancement of the open forest was 6.9 m year −1 , whereas in the Putorana Plateau it was 1.0 m year −1 . Despite the more severe climatic conditions on the Putorana Plateau, the average altitude of the upper forest boundary is slightly higher in contrast to the Khibiny mountain range. In the Khibiny Massif, the forest boundary occupies higher positions in the central areas (Maly Vudayvr), where, according to our assumptions, a warming effect occurs due to radiation from the opposite slopes of the valleys. On the Putorana Plateau (Sukhie Gory massif), the highest altitude of the forest boundary is 700 m a.s.l., and in the northwestern Putorana Plateau (valley of the Bolishoy Avam river) it is located up to 390 m a.s.l. [25]. The differences in the altitude of the upper forest boundary in these two regions of the Putorana Plateau are due to the large influence of Atlantic air masses in the western part of this region and, accordingly, a large amount of precipitation falling here, which also happens in other subarctic territories such as the Polar Urals [20]. On the northern Putorana Plateau, the climate is determined by the significant influence of cold Arctic air masses, whilst in southern and more continental sites of this region the forest limit and treeline are located at higher altitudes in agreement with the so-called mass elevation effect [3,15].
In fact, the forest limit and treeline on the Putorana Plateau is dominated by larch, which supports intense continentality [58], whereas in the Khibiny Massif, the forest limits and treeline are dominated by birch, with Scots pine and Norway spruce as secondary species. The dominance of birch in the Khibiny massif is a clear indicator of climate mildness and oceanic influence [4]. In our opinion, such significant differences in treeline species composition are due to: (1) differences in the temperature regime of the studied territories (the Putorana Plateau has a more severe and continental climate), (2) the duration of the vegetation period which is probably shorter in the Putorana Plateau sites, (3) the snowpack depth in the forest-tundra ecotone, which was higher in the Khibiny Massif sites (global mean of 1.4 m) than in the Putorana Plateau sites (global mean of 1.0 m; Tables S1 and S2), and (4) the number of sunshine hours, which, apparently, is not high enough for larch to grow in the Khibiny Massif. The presence of Scots pine at the Khibiny forest limit also indicates a milder climate than in the Putorana Plateau [59].

Conclusions
The results obtained confirm our hypotheses: (1) in two regions of subarctic Russia, differing in terms of continentality and tree composition, upward treeline shifts were observed; (2) treeline dynamics were driven by changes in heat supply modulated by topography and microrelief; thus, the most favorable conditions for tree settlement were found on southern-exposed slopes and in the upper and middle parts of the treeline ecotone; (3) the upward tree expansion was closely related to summer warming and increased winter precipitation, but also depended on snow accumulation.
In general, the results of this study are important for understanding the patterns and dynamics of subarctic forests and treeline ecotones under cold-limiting conditions. These forest limits and treelines are positively responding to climate warming in terms of recruitment leading to tree encroachment into tundra. However, those responses are driven by snow accumulation and melting and by slope exposure since tree regeneration depended on sufficiently deep snowpacks and enough radiation and temperature in southern-oriented slopes. We demonstrate that considering these local processes is mandatory for projecting treeline dynamics in a warmer future. Treeline dynamics models need to consider the growing season but also winter climate conditions and topographical factors (slope exposure, altitude) to make reliable and realistic forecasts.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/f13020174/s1, Table S1: The average air temperature, the height of the snow cover and the time of its descent at dif-ferent altitude levels of profiles 1 and 2 on the Imandra mountain and profiles 1, 2, 4 in the valley of Lake Small Vudayvr and profiles 1,2,3 Kitchepakh in 2018-2019. All sites are located in the Khibiny massif. Table S2: Average air temperature, the height of the snow cover and the time of its melting at different altitude levels of the profiles on the slopes of the Sukhie Gory massif (Putorana Plateau region) in 2018-2019. Table S3: Average morphometric parameters of trees (mean value ± standard deviation, SD) and areal characteristics of forest stands at different alti-tude levels of altitude profiles within the forest-tundra ecotone (1-3 levels) and the upper part of the mountain-forest belt (level 4) on the slopes of the mountains on the slopes of Mount Imandra (Khibiny massif). Species' abbreviations: BP, birch; PS, Scots pine; PO, Norway spruce. Table S4: Average morphometric parameters of birch trees (mean value ± SD) and areal characteristics of forest stands at different altitude levels of altitude profiles within the forest-tundra ecotone (1-3 levels) and the upper part of the mountain-forest belt (level 4) on the slopes of the mountains in the valley of the lake. Small Vudayvr (Khibiny massif). Table S5: Average morphometric parame-ters of trees (mean value ± SD) and areal characteristics of forest stands at various altitude levels of altitude profiles within the forest-tundra ecotone (1-3 levels) and the upper part of the moun-tain-forest belt (4-5 levels) on the slopes of Mount Kitchepakh (Khibiny massif). Species' abbrevia-tions: BP, birch; PO, Norway spruce. Table S6: Average morphometric parameters of larch trees (mean value ± SD) and areal characteristics of forest stands at different altitude levels of altitude profiles within the forest-tundra ecotone (levels 1-3) and the upper part of the mountain-forest belt (level 4) on the slopes of the Sukhie Gory massif (Putorana Plateau). Table S7: Spearman cor-relation coefficients calculated between the number of trees recruited in the Khibiny massif treeline ecotones over five years and climatic parameters corresponding to previous (first value in each cell) and current (second value in each cell) five-year periods (values with a significance level of p < 0.05 are highlighted in red). Table S8: Spearman correlation coefficients calculated be-tween the number of recruited larch trees on the Putorana Plateau treeline ecotones considering five-year age classes and climatic parameters corresponding to previous (first value in each cell) and current (second value in each cell) five-year periods (values with a significance level of p < 0.05 are highlighted in red).