Ecological Dynamics and Regeneration Expansion of Treeline Ecotones in Response to Climate Change in Northern Bhutan Himalayas

: The alpine treeline ecotones are an early indicator of vegetation’s response to changes in climate, and the advancement of diffuse treeline ecotones has been associated with mean annual warming temperatures. However, the knowledge of how tree demographic size, age and population distribution, and regeneration decrease with increasing elevation and mean annual temperature remain fragmentary in Bhutan. There was no explanation of how treelines migrate in response to the climate. Therefore, the objectives of this study were to investigate tree demographic size and age and population distribution, as well as the regeneration expansion of treeline ecotones of Abies densa trees in response to climate change. Demographic data from thirty transect bands from treeline ecotones and reconstructed mean annual temperatures from tree-rings were used. Regression analysis was used to establish a relationship between elevation/temperature and demographic tree size and age, as well as to determine recruitment frequency distributions and whether these could be driven by climate change. The tree demography indicated that the treeline ecotone in our sampling site is temperature limited. Hence, cooler temperatures at higher elevations should drive decreases in basal diameter, age and recruitment frequencies. From the dendroecological analysis, the diffuse treeline ecotones appear to be climbing on average 1.00 m per year in Northern Bhutan. We also found that the recruitment frequency has increased over recent years (1850–2017), as temperatures continue to rise. The thermal treeline ecotones will be likely to serve as a line of bioclimatic reference against which other zones of bioclimate can be deﬁned. With documented responses of treeline ecotones toward mean annual temperatures, the expectation is that additional warming will continue to inﬂuence regeneration expansion in the future. This dynamic response of treeline ecotones towards the climate acts as an indicator of climate change. Information about climbing treelines and altered ecotones should be a vital part of the material for decision makers to consider, to assess impacts and threats to Himalayan alpine biota.


Introduction
The alpine treeline ecotone in high mountainous regions is one of the important altitudinal boundaries of vegetation; it marks the transition of subalpine forest to the alpine zone [1][2][3][4]. This is an ecotone built up of various microsites in which elements of both ecosystems, the alpine zone and subalpine forest, are nested [1,5]. Treelines are early indicators of vegetation's response to changes in climate [2][3][4][5][6][7][8][9][10][11]. Across the world, studies of arctic or alpine treelines have been conducted with the aim to understand the responses to variations in temperatures [8]. The studies were undertaken to detect changes and assess the threats to arctic and alpine biota in response to the upward movement of treelines in high-altitude and -latitude communities [8]. However, few studies focus on the impacts of environmental changes and climate changes on treeline ecotones, forests and biodiversity aspects [12] as reported earlier [12] in the Himalayan region.
Globally, evidence reveals that tree species generally migrate towards higher latitudes and altitudes [13]. This response and the sensitivity of treelines to climate change is further modulated by regional and local topographic conditions and their treeline forms [5,12]. Though the shift of treelines is not uniform and universal, the shift made in response to changes of climate is widely observed throughout the world [10,12]. Observations of migrating treelines have previously been made in various regions of the Himalayas [10,12,14]. A majority (80%) of diffuse treeline studies have also indicated accelerating advancement under warming conditions of climate, i.e., mean annual, winter and summer temperatures [6,8,15]. The advance of diffuse treelines has been associated with mean annual warming temperatures, which also envelope the growing season as an integral component [8,15,16]. However, the existence of a temperature association with regard to the advancement of the diffuse treeline has not yet been explored in Himalaya and Bhutan. Thus, intensified efforts and dendroclimatic analysis are needed to establish if this condition is also valid in this region.
The structure, position and physiognomy of the treeline ecotones have responded to change in climate, particularly to temperature, in recent decades [5,6,8,12,15,16]. One of the consistent patterns with all treelines of alpine and arctic regions is the gradual decrease of height in trees when approaching the upper limits of treelines [8]. Thus, a significant negative correlation of elevation (masl) against basal diameter (cm), density and age (year) are expected characteristics and patterns of diffused treeline ecotones [6,8,15,17]. the elevation of a treeline is affected by low temperatures [3,4,7] that limit growth and survival and hence evidence should show decrements in diameter with increments in elevation [8,15,18]. In line with the general global responses [6,11], the reports from studied treelines in eastern Nepal Himalaya reveal that temperature is the dominant climatic factors that control tree regeneration and growth dynamics [12,19]. The diffuse treelines are formed and maintained mainly by limitation of growth under low temperature of growing season within 5-8 • C as the primary stressor [8,11,15]. However, the literature on how decreasing air temperatures with altitude influence the growth of trees and their regeneration in mountainous system of the Himalaya still remains fragmentary [20].
The upper natural treeline ecotones of Bhutan are mainly formed by Abies densa Griff. Juniperus species, Betula utilis and Rhododendron species growing at 4000 m on north-facing slopes [20,21]. Scientifically, a majority of these species at the treeline ecotones of Himalaya have proven potential for dendrochronological cross-dating and thus provide opportunities to find out changes in their structure and position of treelines in their response to climate change [10,12,20,[22][23][24]. However, dendrochronological studies of the temperature related to treeline ecotones are still outstanding in Bhutan. Perhaps, there was no straightforward explanation of how and if these ecosystems and ranges shift act as indicator of climate change [12,25,26] from Bhutan. Moreover, the dendrochronological research in seedling survival, mortality, establishment and dynamics of forests could elucidate the structure of age and size and perhaps provide information of the persistence and longevity of forests and tree species [27]. Over the last three decades, dendrochronology has made progress in identifying the responses of treeline ecotones [7,28,29]. This explosion in studies clearly reveals the lack of knowledge and experience in dendrochronological methods [12,26] in Bhutan Himalaya. We therefore established a network of sampling sites at treeline ecotones of Abies densa trees at remote high elevation (3800-4500 masl) regions of northern Bhutan (Figures 1 and 2). The objectives of this study were to investigate the tree population size and age distribution structures, ecological dynamics, and regeneration expansion of treeline ecotones of Abies densa tree species in response to climate change in Northern Bhutan Himalaya. We mainly used methods of dendrochronology to investigate the ecological dynamics of treeline ecotones under the influence of a changing climate [6,7,24,29].   Treeline ecotones are transition zones that bear the evidence of the responses of rapidly shifting vegetation [6][7][8][9]. However, a definition of a treeline as single convention does not exist. [6,8,15]. Within the context of the tempo-spatial structure and growth form of vegetation in mountainous regions [6] of Bhutan Himalaya, we clarify and modify the definition in the following. Here, we define the treeline ecotone as a wide belt between the upper boundary of the closed forest (i.e., timberline) and the upper boundary of the solitary trees (i.e., treeline). Thus, this study focused on the diffuse type of treeline that exhibits a gradual decline in tree density and height along the transition zone [6,8,10]. Trees at the upper range frequently display stunted growth forms (Krummholz) of the relevant tree species (tree-species line) [20] regardless of the typical size, form and age of the species [8,30,31]. In this zone, we sampled Abies densa at (3800-4600 masl) from the timberline to treeless alpine vegetation on north-exposed slopes in the North of Bhutan (Figures 1 and 2).

Methods and Materials
Specifically, the mountain treeline ecotones of Lunana, Laya and Singye were targeted in the central, western and eastern parts, respectively. Here, Abies densa Griff. forms almost pure stands [21,32]. However, the stands are interlaced with understory species including Rhododendron hodgsonii Hook, R. campylocarpum, R. lanatum, and R. campanulatum are found in front-range areas of precipitation higher than 3000 mm annually [20,21,[32][33][34]. This extends from 3600 to 4100 mals in regions of somewhat drier interiors (precipitation around 1200 mm) [21,32,33]. These conditions are conducive to old growth forests and ancient trees which make the forested mountainous region ideal for dendrochronological research [35]. In all three orographically regional sites, the geological formations (soil types, slope gradient and aspects) and climatic conditions, vegetation types and management regimes sharing similarity in their history were identified and used for sampling [7].

Sampling Transects-Enumeration, Measurement and Coring of Trees
The transect bands of long rectangular shape along elevation gradients were used for sampling across the gradients in terms of measures such as aspect, degree of slope, altitude or moisture [7,10,12,29]. The transect bands were oriented parallel to the ground slope, covering the full length of the ecotone, defined above [7,10,12,29,36]. This approach was used to accurately capture the variation in stand age structure and succession along ecotones. Ten transect bands were used at each site, where each band was 24 m width and 300-700 metres in length depending on the extension of the ecotone ( Figure 2). The altitude for both the lower and upper parts of the transect bands were recorded using GPS [7,10,12]. All the ground measurements were made using clinometers and measuring tapes. The distance between bands was set to 300-500 metres.
The enumeration and measurements of trees, including standing, stumps and dead trees, and the tree distance to seed source zone/timberline were carried out. Along the transect band all trees (>2 m), saplings (0.5-2 m) and seedlings (<0.5 m) were enumerated. The basal diameter (BD) of all trees were measured [12,19]. The age of seedlings and saplings were ascertained using their branch whorls and the bud scars on the stem [12,19,30]. Measurements of diameters were made including the bark. Trees unsuitable for coring were measured for their BD and distance from timberline/seed source zone line [36]. The distance of each individual tree from timberline was measured along the length of transect band. The data on physiographic such as altitude, slope and aspect were incorporated across all band transects [7].
Dendrochronological field techniques were used in coring tree-ring samples from Abies densa tree species. Across all transects, all individual trees suitable for dendrochronology were cored for tree-ring samples using increment borers at a direction parallel to counter lines, intended for measurements of calendar dated ring-widths [7,28,[37][38][39][40]. In this study, we cored two cores per tree to facilitate cross-dating within a tree and to represent overall tree growth [7,29]. The trees were cored near to the ground surface [12] to establish the relationships among germination year, height, BD [10,12], annual tree ring width and treeline migration. For investigation of successional processes in stand-age, the exact germination date is important, which was why the trees were cored at the ground base [7]. Depending on availability, the stem discs were collected and prepared using a chainsaw from the base of the dead trees and wind-fallen trees at several location within the transects [28,37,39,40]. The tree ages for uncored tree were estimated using regressions between BD and age.
The same layout of transect band sampling and measurements of treeline ecotones were repeated and used in all three treeline ecotone regions ( Figure 2). The figures were generated using ArcGIS10.8. The ArcGIS10.8, DEM and GIS climate data were sourced from the National Land Commission (NLC), Ministry of Agriculture and Forest (MOAF); Bhutan and GPS data and field photos were sourced from research field grounds.
The analysis of cores and stem discs was performed using techniques of cross-dating, i.e., skeleton plotting [7,37,41]. The principle and technique of cross-dating were used to ensure precisely assigning the correct calendar year to each tree ring and to take into account problematic rings such missing rings, false rings, etc. [37,45]. A complete master chronology was developed with correctly assigned calendar years to represent the overarching signal of the stand level. The statistical tool-the COFECHA computer program [46,47]-was used as second check for the dating of the tree-ring width measurements [43].
Tree-age estimation for trees missing pith and un-cored trees for samples was completed using a geometrical technique. The non-destructive method of obtaining the ages of trees, the counting of tree-rings from increment core samples, was used [48]. If the increment core failed to reach the pith due to heart rot or too short increment corer, the pith year was estimated geometrically [48]. The geometrical method is commonly used to determine the position of the pith using the curvature of the innermost rings [48] in studies of the tree population and dynamics [48,49]. The Pith Indicator method [7] was also used for some cores that missed the pith (James, 2010) to compare and evaluate the above method.
Transfer functions from regressions between dendrochronologically established tree ages and BD were used to establish the ages of the trees not cored for tree-ring samples in the field. In total, 186 tree-ring samples from the Lunana, 131 from the Laya and 302 from the Singye Dzong treelines were cored, amounting to 619 core samples. The total number of enumerated and measured trees was 1070 from Lunana, 1419 from Laya and 891 from the Singye Dzong treeline region, amounting to 3380 trees. For regression analysis, the XLSTAT-Ecology 2020 program [50] was used to carry out the linear regression of the basal diameter against tree age. If the regression revealed a significant relationship, p < 0.001, an estimation of the tree ages was carried out [12].  [51,52] was used in this study [35,53]. As per an earlier report [53], this CRU 4.01 [51] has shown significant relationship with local climate data (DHMS, 2016, DGM, 2016) and tree-ring data of those higher elevation treeline regions of Northern Bhutan [53]. Thus, the average of mean annual temperatures of six CRU grids-19, 27, 20, 28, 21 and 29 located over all three treeline regions of Northern Bhutan was used as climate variability [53]. Firstly, the six grids of CRU TS 4.01 data (hereafter NBGridTEM-6) were used for the comparison and evaluation of earlier mean annual temperatures reconstructed (NB-TEMR) from the same sampled trees and sites used in Northern Bhutan [53]. Secondly, the CRU TS 4.01 (32 grid stations) was used in GIS analysis of spatial interpolation and modelling of climate and elevation, delimiting the distribution and of vegetation (Fir) and biological boundary (. Thirdly, the NBGridTEM-6 was used along with reconstructed mean temperature (hereafter NB-TEMR) to assess the response of treeline regeneration and ecological dynamics to climate change.
The Point-by-Point-Regression (PPR) correlation ®and regression (RSQ) statistics of calibration (CR = 0.67, CRSQ = 0.46; p < 0.001) and validation (VR = 0.50, VRSQ = 0.28, p < 0.001) between NBGridTEM-6 (predictand) and mean composite tree-ring index chronology (predictor) was shown to be highly significant in modelling NB-TEMR in an earlier study [53] from those same sampled tree species of the same sites. The full model f-ratio = 49.57, df = 1 and 55 (p < 0.0001) explained skilful reconstruction model. Moreover, the strong positive correlation (R = 0.61 and RSQ = 0.37) between NB-TEMR and NBGridTEM-6 data significant at p < 0.001; the NB-TEMR has a longer time series that fits for use as climate variability ( Figure 3). As such, the NB-TEMR [53] was found to be appropriate and used as climate variability in this study ( Figure 3). Using the Pearson coefficient correlation (hereafter R) and the regression coefficient of determination (hereafter RSQ) for NB-TEMR and NBGridTEM-6, respectively, against tree demographic population size, age and recruitment frequency significant at p < 0.001 level (2-tailed), a regression model and relationship were established. The XLSTAT-Ecology 2020 program [50] was used in carrying out regression analysis. Figures were generated using Kaleida Graph [53,54].

Tree Demographic Size, Age and Recruitment Frequency Tempo-Spatial Distribution in Relation to Elevation and Climate
Through the distribution of tree size, age, recruitment and frequency [6,12], the regeneration statuses (trees, saplings and seedlings) were established in relation to their distance from the timberline to the tree species limits. The classes of the elevations, tree size and age, and frequency were formed in finding how the tree size, age, recruitment and frequency changed as their distances (elevation) increased from the timberline [36]. Through the divisions of classes based on size, age, recruitment and frequency (for instance, 5-year bin), the stake columns, scatter plots and polygon lines were produced to portray their regeneration distribution with elevation and climate. Thus, the dynamics of the treeline ecotones were analysed using the tempo-spatial distribution of tree recruitments, basal diameter, age and tree frequency against elevation gradients [12]. The XLSTAT-Ecology 20,202 program [50] was used in carrying out regression of elevation (masl) against basal diameter (cm), age (year) recruitment (year) and their respective frequencies. We used R and RSQ at significant [p < 0.001 level, (2-tailed)] to validate the regressions establish for all treeline ecotones in Lunana, Laya and Singye.

Dendroecological Analysis-Tree Regeneration, Delimiting Distribution and Upward Expansion of Treeline Ecotones Driven under Changes in Climate and Elevations
The dynamics of treeline ecotones were analysed through the assessment of the tempospatial distribution patterns of tree size, age, recruitment year, frequency against elevation and climate change trends [12]. The NB-TEMR [53] and NBGridTEM-6 were used for analysis and confirming the significant relationship of climate with tree regeneration expansion, and recruitment frequency at treeline ecotones ( Figure 3). Using the relationship revealed by R and RSQ of NB-TEMR and NBGridTEM-6 against the trees in the Singye treeline recruitment distributions, which were significant at p < 0.001 level (2-tailed), the relationship between distributions and temperature over the past years was confirmed.
All tree size measurements, age and their frequency data from all transect bands, were pooled together to represent respective treeline ecotones. These pooled data were used in analysis of variance (One-Way-ANOVA) to determine site-specific dynamic response of treeline ecotones. Statistical analysis included variance, homogeneity and equality of means of tree size, and age distributions. IBM SPSS, Version 23 [55] was used for the statistical analysis of the demographic population distribution of the treeline ecotones.
In delimiting the upward expansion of the treeline ecotones driven under climate change, the dendroecological analysis was carried out. The expansion of the limit of upper tree species was analysed through observations of the ages of individual tree in whole transect bands [12,14,30]. The rate (m/year) of shifting in the treelines or species limit was calculated by means of dividing the changes in elevational positions of the limit of trees by elapses in time [10,12]. Using the dates of tree establishment, the range of maximum of shift in treeline was calculated [10] as detailed.
The best transect bands for each treeline ecotone in Lunana, Laya and Singye were selected that contained complete tree age tree-ring samples (FYOG to LYOG) and tree stands meeting the characters of treeline ecotone in calculating the regeneration expansion rate. The calendar years of the cross-dated chronologies of the tree-ring samples were used to fix all recruitment years (RY)-first year of growth (FYOG), last year of growth (LYOG)-and years of sample collection and tree ages (TAge) at elevations of various transect bands within respective treeline ecotone. Elevation lapses/differences and tree age lapses between trees at the uppermost limit elevation (UE) and trees at the lowermost elevation (LE) of the transect bands were used [10,12]. The expansion rate in every treeline ecotone was the average of the expansion rates of both the complete tree age based on cross-dated chronologies from tree-ring samples and the estimated tree age based on the regression method to maximize the accuracy of expansion rate.

GIS Analysis-Spatial Interpolation and Modelling of Climate and Elevation Delimiting the Distribution of Fir Forests
In Geographical Information System (GIS) analysis, the ArcGIS10.8 [56] and Global Position System (GPS) were used as tools for the mapping and interpolation of gridded CRU TS 4.01 temperature (32 grid stations) and DEM to evaluate influence of climate and elevation on Abies densa forests distribution over Bhutan [38,[56][57][58]. Detail is provided in Figure 10. The digital elevation model (DEM) and CRU TS 4.01 (32 grid stations) mean annual temperature were used as data in running Kriging for interpolation and modelling analysis. The results of these two analyses were compared and used to evaluate the delimiting distribution of Fir forests influenced by lowering temperature and increasing elevation (masl). GPSMap 60CSx was used to collect geocoordinates of transect bands and treeline ecotones.

Tree Demographic Population Size and Age Class Frequency in Relation to Tempo-Spatial Distribution
Temporal and spatial continuity in regeneration and mortality were portrayed through these age and size structures of population in three treeline ecotones (Figures 4 and 5). The treeline ecotones experienced continuity and increment of growth over the periodic span of 1515 to 2017 in all treeline ecotones ( Figure 4). The basal diameter-frequency distribution and tree age-frequency distribution were decreasing in trend with increasing elevation towards upper bin limit [8] showing right-skewed unimodal (reverse J-shape). These are shown in their basal diameter classes (0-160 cm), tree age classes (1-500 year) and elevation upper-bin limits (3930-4327 masl).  Detail is shown in stake column graphs (Figures 4 and 5). This frequency in diameter distribution decreases as tree size increases with every increase of elevation upper-bin limit ( Figure 4). The unimodal (reverse J-shape) of both basal diameter and age frequency were firstly right skewed with increasing frequency of juveniles in every diameter size and age (year) classes. And secondly, the frequency of both diameter (cm) classes and age (year) classes decreased with successive increments in size and age with increasing elevation upper-bin limits in all treeline ecotones (Figures 4 and 5). As the only exception, the Lunana treeline ecotone followed a bimodal (bell shape) basal diameter-frequency distribution and tree age-frequency distribution in the upper elevation bin limits of 3930 and 3967 masl at lower elevation transect bands (Figures 4 and 5). Normal tree structure distribution at the timber line or the seed source zone was perhaps attributed to this bimodal bell shape.
The basal diameter and age class frequency counts were maximal at lower elevation in all size classes, and, with the increasing elevation, the older-aged and bigger-sized tree frequencies decreased successively in every class leading to disappearance. Subsequently, the juvenile smaller-sized and younger-aged individuals successively dominated the uppermost elevation of the treeline ecotones (Figures 4 and 5).
In a diameter and age bin of 0-20 cm and 1-50 year at an elevation bin limit of 3997 masl and 4027 masl, the frequency counts of juvenile size and age trees dominated the Singye treeline ecotone as compared to others. This is attributed perhaps to profuse juvenile sizes and young-aged tree regeneration under wider opening crowns of scattered old-age trees at lower elevation (Figures 4 and 5). This was perhaps caused by earlier past regime opening gaps under force of unknown events.
This study presented good regeneration of tree species, for example, the establishment of seedlings of Abies densa extended beyond the upper limit of adult trees and at elevations of the uppermost tree species' limit, i.e., 4327 masl (Figures 4 and 5). Both age and size distribution detail the history, pattern and process of dynamics and recruitment of the Abies densa tree species in diffused treeline ecotones. The elevation of the treeline is affected by low temperatures that limit the growth and survival, and hence evidence should be with decrement in diameter and age with increment in elevation [8]. These present high significant spatial and temporal matches with the interpolated model, showing decreased trends and patterns influenced under lowering temperatures and increasing elevations and latitudes across Bhutan (Figure 10).  Based on the regression models of recruitment frequency (polygon) against elevation for Lunana (y = 449.1 − 0.1025x, R = −0.64, p < 0.001) and Laya (y = 1838 − 0.4238x, R = −0.94, p < 0.001), Singye (y = 2883 − 0.6978x, R = −0.83, p < 0.001), the strong significant negative correlations were yielded ( Figure 6). These regression statistics of elevation (masl) against basal diameter (cm) and recruitment frequency are significant at p < 0.001 (2-tailed) and indicated strong trends and characteristics of diffused treeline ecotones in their gradual decline from the timberline to limited tree species [6,8]. This regression model showing a statistically significant relationship of recruitment frequency (polygon) with elevation of treeline ecotones (Figure 6 supports tree-regeneration distributions, trends and patterns of expansion along elevation gradients of treeline ecotones. Moreover, the analysis results portrayed sporadic and cohort distributions at base elevations of 3980 to 4020 masl in the Singye ecotone (Figures 6 and 8). As per filed observations, the causes for these cohorts and sporadic distributions are perhaps attributed to the heterogeneity in conditions of microtopography and microclimate such as slope gradients, rocks, scree, moisture and nutrients [8,59,60]. However, in the Lunana and Laya treeline ecotones, the distribution was uniform and in continuity of gradual frequency decline with increasing elevation.

Tree Size and Age Distribution as Site Specific Dynamic Response of Treeline Ecotones
The tree size (BD) and age (year) declined gradually with increase of elevation; however, the evidence of site-specific dynamics of tree size and age distribution patterns were revealed in the statistical analysis of variance.
All measurements of tree basal diameter (cm) and age (year) revealed a significant difference (p < 0.001) in their variance, homogeneity of variances, equality of mean and normality. The statistics for basal diameter are as follows: Thus, each variable is not identical to the others, and, as such, all samples come from different populations. All these statistical analysis test results proved differences in means of variances and homogeneity between the three sites significant at p < 0.001. The differences were revealed between the treeline ecotones of three different regions in their means of basal diameters (cm) and age (year) distributions of the Abies densa trees species. The causes for these differences in the means of basal diameter and age distribution between regions are perhaps variation in their spatial heterogenetic conditions of microclimate and microtopography such as ground slope gradients, rock cliffs, moraines, moisture, nutrients [8,59,60] and undergrowth alpine vegetation.

Tree Recruitment Frequency Distribution in Relation to Climate Change
In regression analysis, the R and RSQ statistics of NBGridTEM-6 with Singye recruitment frequency (R = 0.46, RSQ = 0.21, p < 0.005) was significant. The correlation of NB-TEMR with Singye tree recruitment frequency (R = 0.47, RSQ = 0.22, p < 0.005) was also significant, revealing their strong relationship. Both NB-TEMR and NBGridTEM-6 have indicated strong relation with recruitment frequency from the Singye treeline ecotone and the so-proven influence of mean annual temperature on tree regeneration and growth in all treeline regions (Figure 7). Moreover, the graphical positive slope (stake columns) reveals its strong relationship with the temperature curve line. These regression statistics for recruitment frequencies were employed as a reference to prove and validate the strong significant relationship between mean annual temperature and tree recruitment frequency in all three regions (Figure 7). This graphical analysis showcase increases of recruitment frequency with increasing temperature over recent years from past years (Figure 7). The graph depicts degradation in past years beyond 1800 perhaps attributed to the closure of crowns that limit undergrowth regeneration and the domination by larger-sized and older-aged trees in lower elevations ( Figure 8). This mean annual temperature of Northern Bhutan treeline regions has shown temperature as a dominant climatic factor that influences tree regeneration, growth and dynamics. This strong relationship is significant at p < 0.001 and revealed favour of climatic conditions during the time of recruitment that serve as crucial for regulating dynamics of treelines. The strong positive correlation significant at p < 0.001 between recruitment frequency and mean annual temperature proves that conditions of climate that influence regeneration are similar to those favourable to tree radial growth in treeline regions [12]. Therefore, it was found that ascending annual temperatures could increase the rates of tree growth and survival of seedlings, leading to more rapid recruitment beyond the diffuse treelines.  With the increasing elevation, the growth rates of trees are decreasing as per basal diameter (Figure 4). Further, the production and viability of seeds are perhaps declining with increasing elevation as per the recruitment frequency that is positively related to temperature [8] (Figures 7, 8 and 10). These present significant spatial and temporal matches, with interpolated models showing decreased trends and patterns influenced by lowering temperatures and increasing elevations and latitudes across Bhutan (Figure 10).

Tree Recruitment Year Distribution in Relation to Climate Change and Elevation
This graphical analysis also shows the decrease of the recruitment frequency distribution with the increase of elevation and increase of recruitment downward in the recent years from past years in all treeline ecotones (Figure 8).
Temporal and spatial continuity in regeneration and mortality were portrayed through these recruitment population structures (Figures 7 and 8). Numerous saplings, seedlings and trees of Abies densa were established adjacent to limits of tree species in upper elevations of the established position of treeline ecotones (Figures 7 and 8). This indicates that shortage of seed dispersal is not the limiting factors for shifts of treelines in recent times [12,26,58,61]. The recruitment trend decreases with the increase of elevation (m) and recruitment increased over recent years with the increase of mean annual temperature. These depict characteristics of treeline ecotones and the influence of temperature in all three treeline regions (Figure 8). All three treeline ecotones portray shifting regimes as younger tree regeneration distributes and moves upward toward higher elevation mountain tops from older regeneration at lower elevations ( Figure 8). In these spatial and temporal scales, the recruitment trends and processes indicate that the advance of the treeline is finally the function of recruitment explained in the context of the micro-climate and macro-climate in addition to the conditions of microsites [5,8]. The increase in the abundance of trees and seedlings in the higher elevation and beyond the treelines is the evidence that indicates the increased of likelihood/probability to occur and to establish beyond the position of the current treeline [6] under a warming climate (Figure 8).

Tree Regeneration and Age Distribution in the Upward Expansion of Treeline Ecotones Driven under Climate Change
This model (Figure 9) presents the tree age distribution in relation to changes in elevation and temperature. This model showcases the decrease of tree age (year) with the increase of elevation and the increase of tree frequency upward younger age from side of older age with increasing temperature in recent years.
Thus, all three treeline ecotones portray shifting regimes as younger-aged trees distribute and move upward toward higher elevation mountain tops from older-aged trees at lower elevation.
Plus, the distribution of a greater number of younger regeneration-aged trees in recent time with increasing mean annual temperatures is an indication of a treeline ecotone shifting regime under the warming of the climate. This reveals the treeline expansion and upward shifting as tree age decreases with the increase of elevation and tree frequency increases with the increasing mean annual temperature (Figures 6 and 9). Based on the above model (tree measurement data of thirty transect bands), the upward shifting rates of the treeline ecotones were computed as detailed below.   (2000-2021). Regenerations and growth of Abies densa tree species are enhancing adequately at and below this upper limited elevation of 4159.30 masl that continuously building adequate stand stock capable for initiating further upward expansion. The Singye Dzong treeline expansion rate was comparatively lower than Lunana and Laya treeline ecotone. The causes for such variations are perhaps attributed to denser Krummholz Rhododendrons; undergrowth alpine vegetation; microtopographic conditions such as rock cliffs, topographical slope gradients; moisture, and nutrients [59,60] As per the field evidence, the slowing down of the recent advance of the treeline has perhaps been interrupted by the competition and dominance of Krummholz Rhododendron trees at several locations in subalpine regions. Dense Krummholz Rhododendron trees appear to operate as an effective barrier for the upslope migration of tree species in some other Himalayan regions [12,26]. If these barriers are not present, the treeline ecotones are perhaps likely to expand at higher rates of 20 to 80 m every decade [25] covering more upward alpine areas per the rising rate of temperature. All in all, the recruitments have been established before decades in all treelines and occurring below the position of treeline ecotones, indicating stable dynamics in the long-term context for influencing the advance of the treeline position and in response to the recent warming of the climate (Figures 6, 8 and 9). All these population demographic structures, growths, and distribution analyses presented tree size, age, recruitment and their frequencies indicated formation of diffuse treelines and maintenance of treelines are under limitation/control of mean annual temperature (Figures 6 and 8-10).

Geographical Position of Vegetation Margins in Response to Changes in Climate and Elevation
Distribution of Abies densa tree species (fir) in response to elevation (DEM) and temperature (CRU 4.01) is detailed below. Thirty-two gridded mean annual temperatures (CRU TS 4.01) represent local instrument climate data of Bhutan. As per this temperature interpolation and elevation model, the fir forests are distributed commonly in colder regions/areas of temperature operating at approximately −2.267-11.23 • C only in high elevation beyond 2700 masl in mountains. In addition, consequently, the fir forests are not found in warmer regions depicting higher than 11.23 • C and lower than 2700 masl ( Figure 10). This spatial analysis indicated fir forest ecosystem coverage in areas of colder regions exhibited their linkages to the limitation by the climate and limitation caused by increment in elevation [18] for their growth. The Figure 10 model supports trends and patterns of how tree size, age, recruitment and their frequencies (Figures 4, 5 and 7-10) decrease under the influence of lowering temperatures and increasing elevations (altitude) and latitudes across Bhutan. The results revealed treelines are in equilibrium with climatic conditions ultimately limited by low mean annual temperatures and increased elevations and showcases their positive response to the warming of the climate ( Figure 10). Figure 10 showcases increase of elevation (masl) decreases mean annual temperature ( • C) ultimately limiting growth and distribution of forests over Bhutan Himalaya. The increase of elevation forces decreases of temperature that predominantly influence vegetation [6,62,63] and, most strikingly, the formation of treeline ecotones under such global thermocline caused by decline of temperature [6,63]. All these treelines appear to serve as a line of global bioclimatic reference against which other zones of bioclimate can be defined. This study provides areal statistics for various climatic belts in the mountains. Although the interpolation and distribution of temperature and elevation (DEM) indicates their dominant influence on distribution of fir forests, it reveals the influence of these two factors on distribution of agro-ecological zones and agricultural crop plants under changing climate in Bhutan.

Tree Demographic Size, Age and Population Distribution Triggered by Changes in Elevation
One of the consistent patterns with all treelines of alpine and arctic regions is the gradual decrease of height in trees as one approaches the upper limits of treelines [8,15]. These agreed with two potential causes for the declining of tree growth with proximity to upper elevation limit of treeline. These are limitations in essential nutrients because of the low supply or uptake especially the carbon and low temperature impairing the formation of tissue [6,8,15]. Thus, the decreasing distribution trends of basal diameter and tree age with increasing elevation strongly agreed with other reports [8] in diffused treeline ecotones [8,15]. The negative Pearson coefficient correlation (R) and regression coefficient of determination (RSQ) of elevation (masl) against basal diameter (cm) and age (year) significant at p < 0.001 (2-tailed) proved strong characteristics and patterns of diffuse treeline ecotones of alpine and artic regions in their gradual decrease of size and age from timberline to treeline/limited tree species [6,8,15]. The diffuse treelines are formed and maintained mainly by limitation of growth with low temperature of growing season within 5-8 • C as primary stressor identified in earlier reports [8,11,15]. However, this study's results proved that the primary stressor in forming and regulating these diffuse treelines through limitation of growth is perhaps the lowering of mean annual temperature. The cause of formation of treeline is the growth limitation due to low temperature [16]. As such, the elevation of treeline is affected by low temperatures that limit growth and survival, and hence evidence should be gathered of decrement in basal diameter and age with increment in elevation.

Tree Regeneration and Recruitment Dynamics Driven by Warming Mean Annual Temperatures
The significant correlation and regression statistics (p < 0.001) of recruitment frequencies against temperature revealed strong response to warming climate and expanding upward mountain top (Figure 9). In line with worldwide common responses [6,11,16] and reports from studied treelines in Nepal Himalaya [12], this mean annual temperature of Northern Bhutan treeline regions have shown temperature as a dominant climatic factor that control tree regeneration, growth and dynamics.
Recruitment of Abies spectabilis D. Don is sensitive to temperature, moisture and/or conditions of drought as reported in other studies [12,14]. High temperatures of winter and summer favoured the regeneration of Abies spectabilis with a correlation of positively significant [12] with provided condition of no limitation to moisture in Nepal Himalaya [12]. The results are consistent with expectations [8] in diffuse treelines that responded to warming of annual temperatures, whereas Krummholz, abrupt and island treelines responded to warming of winter [12,64]. The diffuse treelines are limited by conditions of the growing season and are most responsive to the warming of the climate [8,15]. In this study, the annual mean temperature, comprised of high temperatures of winter and summer favoured the regeneration of Abies densa, provided adequate moisture in treeline regions. The moisture is adequate as reported earlier that stands of Abies densa trees are found in front range areas of precipitation higher than 3000 mm annually in Bhutan [20,21,[32][33][34].
With respect to climate-influenced spatial displacement, the several other researchers reported changes of structural features of the vegetation at position of treeline. These structural changes include increased rates of the growth of trees and increase in tree density of the populations [9,[65][66][67][68]. Numerous saplings and seedlings of Abies densa were established adjacent to limits of tree species in upper elevation of the stabilized position of treeline ecotones (Figures 6 and 8). This indicates that the shortage of seed dispersal is not the limiting factor for shifts of treelines in recent times of rising temperature [12,26,58].
The recruitment trend decreases with the increase of elevation (m) and recruitment increases towards recent years with increase of mean annual temperature. These depict characteristics of treeline ecotones and influence of temperature in all three treeline regions (Figures 9 and 10. The diffuse treelines are formed and maintained mainly by limitation of growth with low temperature of growing season within 5-8 • C as primary stressor identified in earlier reports [8,11,15]. However, this study indicates that these diffuse treelines were formed and maintained primarily by limitation of growth by low temperature of mean annual temperature as primary stressor. This means the temperature of growing season within 5-8 • C was also responsible as the associates of mean annual temperature for dynamics of treeline ecotones and act as stressor [8,11,15]. Thus, this study's results proved that primary stressor in forming and regulating these diffuse treelines through limitation of growth is perhaps mean annual temperature.

Regeneration Dynamics in Upward Expansion of Treeline Ecotones Driven by Climate Change
Regeneration and growth of Abies densa tree species are enhancing in adequate at and below upper limited elevations of 4159.30 masl (treeline and tree species limits) that continuously building adequate stand stock capable (ecological and growth conditions) for initiating further upward expansion Abies densa tree species (Figures 6, 8 and 10). The diffuse treelines regulated by limitation of growth due to influence of the low temperature of growing season are making advance [8,15]. The distributional margin of tree species would shift to higher latitudes or altitudes if their growth were limited by low temperatures of growing season [6,8]. Thus, with no doubt that temperature-related phenomena exert controls of dominant over the position of treeline at global scale [6,16]. Though the shift of treeline is not uniform and universal, the shift made in response to change of climate is widely observed throughout the world. The earlier observations made regarding the shifting of treelines were in consistent and similar among various reports in Himalaya [10,12,14]. The average rate of upward shifting was 0.93 m per year was reported in Abies spectabilis tree species with stabilization in recent time at majority of the sites in Nepal Himalaya [12]. There were site-and species-specific patterns of shifting in most of the study sites with little shift of treeline in recent years [10,12,14]. In association with the changing climate and environmental conditions in Himalaya, the upward shifting of the treelines is not a unidirectional response. The treeline shifts vary depending on the history of the climate of the region involved, the inclusion of species and site, and also the analysis approach [10][11][12]14,26]. However, in this study, the advance of diffuse treeline ecotone (Abies densa tree species) was found in all three different regions on average 1.289 m per year in Lunana, 1.022 m per year in Laya and 0.503 m per year in Singye.
The knowledge of how susceptible the trees at elevational timberlines and their seedlings to frost drought is very fragmentary in Himalaya [20], and studies have not been carried out till now in Bhutan Himalaya. The dieback happened to suffer during dormant season particularly during winter [8,15]. The dieback is also highly sensitive to extremes of low temperature and other stressors [15]. However, this dieback is not common in Bhutan and not observed in research sites of this study. The top dying Abies densa trees were observed in several trees at limits of tree species perhaps caused by damages from snow, wind, frosts etc. These are coppicing well in maintaining their growth in high elevation mountains of Bhutan. Thus, the field evidence is adequate to explain the impacts of spring frosts cannot prevent migration of treeline ecotones. The growth of Abies spectabilis close to treelines in central [14] and western Nepal [62] has shown negative correlation with temperatures of spring season (MAM) and positive correlation with precipitations belonging to same months [63]. All these show that sudden late frost events are not likely to prevent the upward migration of treelines in Bhutan Himalaya. As per this research data, the strong positive correlation significant at p < 0.001 between recruitment frequency and mean annual temperature prove climatic conditions that influence regeneration are similar to those favourable to tree radial growth in treeline regions [12,64]. Therefore, these prove only that the warming annual temperature could increase the rates of the tree growth and survival of seedlings, leading to more rapid recruitment beyond the diffuse treelines.
The results (expansion rates) of this study indicated ongoing warming mean annual temperatures comprised of both summer and winter season prove to dominate as the primary factor in expanding treeline ecotones upward higher elevations. As per the results of this study, the presence of adequate regeneration of saplings and trees beyond the age of more than 17 years displayed secured establishment and stability at their upper most elevation limits for upward expansion under ongoing warming climate in Northern Bhutan (Figures 6, 8 and 10). The recruitments have been established over previous decades and as occurring below the position of treelines indicating stability dynamics in long-term context [8,15] for influencing the advance of the treeline position in response to recent warming of climate.

Conclusions
The dendrochronological statistics of tree demography indicated the decrement of tree size, age, recruitment and growth with increasing elevation gradients portraying temperature as a limiting factor. As such, the elevation of the treeline ecotone is affected by low temperatures that limit growth and survival, and hence evidence should be gathered with regard to the decrement in basal diameter, age and recruitment frequency with the increment in elevation. The recruitment frequency has increased over recent years (1850-2017) with the increase of mean annual temperature depicting influence of temperature in all three treeline regions. The highly significant correlation of tree age and recruitments and their frequencies with temperature revealed strong response of treeline ecotones to warming climate as a dominant factor that controls tree regeneration, growth and dynamics.
As per dendrochronological analysis results, the advance of diffuse treeline ecotone (Abies densa tree species) was found in all three different regions on average 1.289 m per year in Lunana, 1.022 m per year in Laya and 0.503 m per year in Singye. In this study, the mean annual temperature, comprising both the summer and winter seasons, proved to dominate as the primary factor in expanding treeline ecotones. Furthermore, the presence of adequate regeneration of saplings and trees beyond the age of more than 17 years displayed secured establishment and stability at their uppermost elevation limits (4313 masl) for upward expansion under ongoing warming climate in Northern Bhutan. With the responses of treelines and trees to mean annual temperatures over more recent and earlier periods in the past, the expectation is that the warming of the climate will exert such influences on regeneration expansion in future as well. The changes in climate and geographical elevation have driven treeline ecotones and Fir Forest distributions in Bhutan. The thermal treeline ecotone distribution is likely to serve as a line of bioclimatic reference against which other zones of bioclimate can be defined. Though elevation and temperature exert dominant control on fir forest distribution, it reveals these two factors perhaps influence distribution of agro-ecological zones and agricultural crop plants under changing climate. The existence of these relationships of tree species with geographical elevation and climate is likely to be relevant and useful for planners in assessment of development of plans.
The future research on dendroclimatology and dendroecology can be useful for additional discovery of climatic sensitive diffuse treeline ecotones to warming climate in other biogeographic regions. Thus, similar trends of research will aid in detecting changes and assessing the threats to alpine biota in response to the movement of treeline ecotones upward higher elevation limits. The dynamic response of treeline ecotones to to changes in climate act as an indicator of climate change, and, as such, the setting up of these treeline ecotones areas as research core zones will further help to predict changes in climate and its consequences.