Timberline structure and woody taxa regeneration towards treeline along latitudinal gradients in Khangchendzonga National Park, Eastern Himalaya

With main purpose of developing a coherent baseline information of timberline forests of Sikkim, Eastern Himalaya, we investigated, (i) forest structure and composition, (ii) regeneration status of dominant timberline tree species between timberline and treeline, (iii) influence of environmental variables in species composition, and (iv) relationship between environmental variables and ecological attributes. The study was conducted along the timberline stretch of Dzongri landscape nested within the Khangchendzonga National Park (UNESCO’s World Heritage Site, 2016), a core zone of Khangchendzonga Biosphere Reserve (in UNESCO WNBRs, 2018), Sikkim, India. We employed quadrat method to sample 9 contiguous sites to capture all possible variations in timberline composition. Transect method was used to study the regeneration of woody taxa between timberline to treeline. In total, 20 woody species belonging to 10 genera and 6 families were recorded. Among these, Abies densa, Rhododendron lanatum and Sorbus microphylla exhibited higher dominance, comprising of 50% of the total importance value index (IVI) weightage. Betula utilis the common treeline species in much of the western and central parts of Himalaya was absent here. Tree density in studied timberline was significantly higher than its western Himalayan counterparts of Indian Himalayan region. Environmental variables viz., elevation, slope, and humus were observed determinants of species composition across the study area. The species dominance correlated negatively (p<0.01; n = 9) with species diversity and richness. We observed an irregular spatial pattern of timberline across the 9 study sites, and the extent of timberline elements (seedling; sapling; live tree or dead tree) ranged between 5.3m to 187.7m higher than the current timberline at different sites. The present trends suggest that upslope advancement is unlikely to occur in near future, while treeline densification is anticipated. Further investigations are suggested to develop a holistic understating of these timberline patterns across the Indian Himalayan region.


Introduction
Ecotones are the transition zones between the adjacent ecological systems or vegetation types [1,2]. Across the globe in the high mountain areas, alpine timberline which marks the upper limit of closed canopy forests represents a major ecotone between forested and non-forested vegetation [3]. In recent times rising temperature has led to accelerated changes in mountain ecosystems [4]. It is assumed that these climate-induced-changes will often occur first at the margin of a plant's occurrence, rather than in the centre of distribution, and can be detected rapidly through ecological monitoring [5]. Understanding of the environmental factors responsible for determining the geographical and ecological extent of species is important [6]. Studies have indicated that the timberline vegetation responds to changing environment through changes in their growth, growth forms, regeneration and shift in their habitat (spatial) [7,8]. The heterogeneity of timberline vegetation increases from global to local scale [9], so the factors influencing the plant diversity at local scale are not the same as those acting at global scale [10,11]. Thus, the understanding on the dynamics of these forest ecosystems, coupled with other influencing abiotic components of the ecosystem, is important [12]. Temperature is a principal influencing factor for prompting snowmelt and the length of growing season which are critical in determining the nature of plant communities in the alpine timberline ecotones [13,14]. Further, it directly affects the productivity, composition, and diversity of an ecosystem [15]. However, seedling recruitment and colonization, stand densification and regeneration can vary greatly between the slopes and aspects, due to the differential presence of permafrost [16], and factors such as stand history, dispersal ability, habitat suitability, disturbance and ecological interactions [17]. Timberline species are observed at the limits of their ranges because at least one environmental factor is at a critical minimum [18]. However, the cumulative effects of several factors restrict the upper limit of timberlines and the factor that primarily determines timberline can be used to define the timberline type. The four types of upper timberline boundaries were recognized by Plesník [19], i.e., climatic, orographic, avalanche and edaphic. These timberline boundaries were more precisely described by Jodłowski [20] as, Orographic, Morphological, Edaphic, Mechanical and anthropogenic.
Himalaya is a biologically, culturally, and ecologically important [21][22][23][24]. Owing to its unique physiogeography with an average elevation above 4000m asl., Himalaya is highly sensitive and vulnerable to climate change [25][26][27][28]. Compared to other parts of Himalaya, the eastern Himalayan region is generally more moist and biologically and culturally more diverse [29]. Himalayan timberlines have been less investigated for vegetation changes under the influence of changing climate than its European counterparts [30], and because of their remoteness and inaccessibility most of the studies have largely been confined to monitor their spatial distribution [31,32]. Across the Indian Himalayan Region (IHR), the vegetation surveys mainly focused on the typical forests or the vegetation zone [21,[33][34][35][36]. Few studies on ecological assessment within the timberline of the Himalayas are available and majority of them concern with Western Himalaya in Uttarakhand [37][38][39]. A need for a long term research to ensure conservation and management of timberline areas of Eastern Himalaya has been stressed [40]. To fill this gap and to strengthen the understanding of timberline, the present study was undertaken along the timberline ecotone of the Dzongri landscape of the Khangchendzonga National Park of Sikkim, Eastern Himalaya. To capture representative picture of the region's timberline ecotone we sampled 9 contiguous transects in Khangchendzonga National Park.
The study focuses on (1) the structure and composition of timberline forests, based on 9 contiguous transects, and their regeneration status of important species, (2) structural and positional extent of upward/or downward shift of important timberline tree species, to examine whether there is upslope advancement in treeline or densification of margin between timberline to treeline, (3) to understand the influence of environmental variables in species composition, and their relationship with ecological attributes.

Study area
The study was conducted in Sikkim, the eastern Himalayan state of India, which extends between 27˚03'41'' and 28˚7'34"N latitude and 88˚03'40" and 88˚57'19" E longitude, and 300m asl (foothills) and 8586m asl (Mount Khangchendzonga). Much of the geographical area of Sikkim (7096 km 2 ) forms the catchment of river Teesta. The state shares borders with Bhutan in East, Nepal in West, Tibetan autonomous Region of China in North and West Bengal state of India in South. Due to complex topography, wide altitudinal and climatic variations, different ecological zones have developed, which uphold a rich diversity of species (flora and fauna) and diversified landscapes [24,41]. It is a monsoon (June to September) drenched region, with an annual precipitation exceeding 2500 mm in many areas.
For the timberline investigations, we carried out intensive field study in Khangchendzonga National Park (KNP). The KNP, a core zone of Khangchendzonga Biosphere Reserve, is a part of the Himalaya global biodiversity hot spot. Inscribed as the World heritage site by UNESCO in July 2016 (http://whc.unesco.org/en/newproperties); the KNP lies along 1,829 to 8,586 m altitude and 27.60˚latitude and 88.19˚longitude and covers an area of 1784 km 2 , approximately 25.14% land area of Sikkim state. Khangchendzonga Biosphere Reserve is now included in the UNESCO's World Network of Biosphere Reserves in July 2018 (https://en.unesco.org/ news/twenty-four-new-sites-join-unesco-s-world-network-biosphere-reserves; accessed on 4.8.2018), as a result of strong conservation policies of Sikkim Government. These conservation leads through organic policies of Sikkim government get further affirmed by the UN's Food and Agriculture Organization's (with World Future Council and IFOAM-Organics International) 'Future Policy Gold Award' on 15 th October 2018, conferred on to Dr Pawan Chamling, Hon'ble Chief Minister of Sikkim. Despite remote and difficult terrain with frequent barriers to movement (walking being the only transport mode) from one side to another, we managed to sample nine timberline sites across the latitudes, in subalpine zone of Dzongri landscape in West Sikkim (Fig 1 and Table 1). These nine sites were located between the conflux of two snow fed rivers Prekchu and Rathongchu originating from the Mount Khangchendzonga range and constitute the Rangit basin. The Rangit is a tributary of Teesta, the biggest river of state Sikkim.

Study area characteristics
The Khangchendzonga National Park is known for (1) habitat contiguity from sub-tropical forests to cold desert, (2) diverse Rhododendron krummholz elements between sub-alpine and alpine zone, (3) high mountain ranges of the total geographical area 90% above 3000m, 70% above 4000m asl and 34% area as glaciers, ice-sheets or perpetual snow [41], and (4) Over 1580 species of vascular plants occurring [42] under three major vegetation zones i.e. the sub-tropical forests (ca 1800m), temperate forests (1800-3500m), and alpine forests (ca3500-5000m), which are further divisible in to10 distinct forest types and subtypes [43] (S1 Table). In the study area the mean annual temperature is 3.6˚C with maximum temperature of warmest month 14.16˚C and, minimum temperature of coldest month -12.16˚C. The area receives 765.29 mm annual precipitation, and it remains snow covered for almost 6 months. The 1 km 2 resolution climate data of the study area was procured from WorldClim online data repository hosted at http://www.worldclim.org and processed using Arc GIS 10.2 special analyst tool [44]  Timberline structure along latitudes in Khangchendzonga National Park, Eastern Himalaya

Field sampling
In this study, we treated timberline as the upper limit of timber yielding trees with over 30% crown cover and their continuity with the subalpine forest below. We carried out field investigations (largely woody species) during September, October, November of 2016, and April of 2017. The sites occurred within 20 Km horizontal extent of Dzongri landscape. The nine contiguous sites roughly covered the entire range of timberline communities of Dzongri landscape (for latitudes of sites, please see Table 1).
At each study site, we randomly laid three adjacent plots (20m×50m) covering the total cover area of 3000m 2 per site. Within each of three plots, we sampled 5 quadrats (10m×10m) for tree/ large sapling/ seedling enumeration and nested within each 10×10 m quadrat, we laid one 5m×5m random quadrat for the sampling of shrubs/saplings. Thus, 15 quadrats (n = 125) were laid for the sampling of trees, and 15 quadrats (n = 125) for the sampling of shrubs/saplings in each study site. Woody taxa having the �30 cm CBH (circumference at breast height, 1.3 cm above ground level) were considered as tree, taxa having CBH range between � 10cm to 30 cm were counted as saplings and those with circumference <10 cm/ or collar diameter <3.14 cm at 1 cm above the ground level were counted as seedlings. We measured the circumference using measuring tape, plant height and slope of each 10m×10m quadrat using clinometer, seedling diameter with the help of digital Vernier calliper (Model: CD-8˝CS, Mituyoto, Japan), and canopy cover using spherical-crown-densiometer (Convex model-A, Jafri Survey Instruments, India) [45]. In each plot we measured the humus depth by inserting steel scale on the ground randomly at each 10m×10m quadrat. The permission to work in KNP/KBR was given by the Chief Wildlife Warden of FEWMD department, Sikkim. And, our field studies did not involve endangered or protected species.

Transect above timberline to treeline
At most sites timberline gradually gave away to treelines. To observe the extent of upward presence of dominant timberline tree species, and regeneration status in this ecotone, we laid out three vertical transects (2m wide) from each site up to the last individuals (seedling; sapling; tree or dead tree) present in vertical transects. We recorded the number of individuals (seedling; sapling; tree or dead tree) and measured their distance from timberline plots ( Table 1 and Fig 1). The advancement of timberline tree species was estimated by subtracting the elevation of adult trees in timberline study plots with the elevation of the last individual (seedling, sapling, live tree/dead tree) at tree line in vertical transect [46].

Data analysis
Tree and shrub layers were analysed for density, frequency, abundance, relative density, relative frequency and regeneration of tree species using the method described by Phillips [47] and Misra [48]. We determined the species richness as the number of species per unit area [49]. We calculated the total basal area (TBA) of tree species in m 2 ha -1 and used it to express the relative dominance of a species. We calculated the importance value index (IVI) by following Curtis [50], i.e. sum of relative frequency, relative density and relative dominance. We computed the diversity Index [51], Simpson's index of Dominance, Margalef's Index of Species richness and Shannon Index of Species evenness following Magurran [52]. Briefly, we pooled the frequencies of all woody taxa of study area and classified them in five frequency classes i.e.  [53].
To analyse the distribution patterns within study sites we used an abundance-to-frequency ratio (A/F) [54]. The A/F value, 0.05 indicated contagious distribution, 0.025 to 0.05 indicated random distribution and 0.025 indicated regular distribution [55]. We checked the sampling adequacy to cover species diversity of the study area by using two non-parametric estimators (Chao 2 and Jack 1) using the statistical software Estimate S Version 9.1.0 [56,57] and estimated values were found in close connection with the tested estimators (S1 Fig). All the tested estimators Chao 2 (20.18 species; 99.11%), Jack 1 (21.78 species; 91.83%) and Cole rarefaction (20 species; 100%) authenticated our sampling efforts (20 species). We used Duncan's multiple range test (DMRT) to compare the significant differences (p<0.05) in mean values of important value index (IVI) and density values of woody taxa between different sites. We used Pearson's correlation to assess the relationship between ecological and environmental variables. To study microclimatic similarities among studied sites, we used an unweighted pair group method with arithmetic mean (UPGMA) agglomerative hierarchical clustering, based on presence and absence data [58] of 20 species of woody taxa from all studied nine sites, and the dendrogram was prepared based on the Jaccard's similarity matrix. Furthermore, we analysed the species abundance and environmental variables of all nine sites to study the interactions between environmental variables and vegetation composition, using canonical correspondence analysis (CCA) [34]. The CCA was performed using XLSTAT trial version. Prior to perform CCA, collinearity between the independent variables was checked through variance inflation factors (VIF) and independent factors having the VIF values less than 3 (slope, aspect, humus depth, distance of treeline and elevation) were selected to avoid the multicollinearity issue. All the statistical analysis (unless otherwise stated) were done using SPSS statistical package for Windows (Version 20; SPSS Inc., Chicago, USA) statistical software package.

Species composition
We recorded 20 woody species (11 tree and 9 shrub), belonging to 10 genera and 6 families from the timberline forest of KNP, of which Ericaceae with 13 species, was the most prominent family, followed by Rosaceae (3 species) and remaining families each was represented one species each (S2 Table). The canopy tree layer was dominated by Abies densa and undercanopy layer consisted of Sorbus microphylla, Rhododendron lanatum and Rhododendron wightii, while Rosa sericea, Ribes glaciale and Juniperus recurva were dominant in the shrub layer (S2 Table). In several sites understory tree species were more prominent than A. densa or were sole tree species (site 9). The undercanopy tree species had about 180-190 IVI at sites 2,3,4 and 6, about 231 at site 5, 255 at site 7 and 250 at site 8. While A. densa trees were up to 40 m tall, the undercanopy trees could grow only up to 8-12 m tall. The tree species number across the study sites ranged from 9 to 14. We identified five forest communities using Duncan's multiple range test based on important value index (IVI) of tree species as following: Abies densa, (site 1,6); Abies densa-Rhododendron lanatum (site 2,4); Abies densa-Sorbus microphylla (Site 3); mixed broadleaved-Abies densa (site 5,7,8); and Rhododendron wightii (site 9). Rosa sericia was observed as a dominant shrub species in most of the communities (S3 Table, Fig 2). At other sites species of Rhododendron and Gaultheria shared importance values.
Across the sties density ranged from 120 individual ha -1 to 374 ha -1 (p>0.05). Undercanopy species like, Sorbus microphylla (104 individual ha -1 ) and Rhododendron lanatum (93 individual ha -1 ) had higher densities than the canopy species Abies densa (64 individual ha -1 ), while Pieris villosa (<1 individual ha -1 ) indicated the lowest values. The average sampling area exhibited overall sapling density of 339 ha -1 , which varied from 138 to 576 individual ha -1 ) across the sites. Sorbus microphylla (94 individual ha -1 ) and Rhododendron lanatum (91 individual ha -1 ) were most common in sapling layer, whereas Rhododendron arboreum (4 individual ha -1 ) indicated the least value. There were fewer seedlings than saplings, the average being varied from 63 to 491 individual ha -1 , across the sites. A. densa had the highest seedling density (121 individual ha -1 ), indicating its ability to grow in its own shade. The entire shrub layer recorded the total average density of 162 individual ha -1 , and across the sites the density varied from 57 to 304 individual ha -1 . Rosa sericea and Rhododendron lepidotum were common shrub species ( Table 2 and  S4 Table). Relatively lower standard error values indicate uniform distribution of species while the higher values compared to mean value indicated the speckled distribution i.e., species was found only in one or two plots of a site (S4-S7 Tables).

Frequency and dominance distribution
The Raunkiaer's frequency analysis revealed that the overall frequency distribution of woody taxa across the forest communities of timberline is slightly heterogeneous and individuals of each frequency class were present in the study area. Sorbus microphylla and Rhododendron lanatum were among the most frequent angiospermic tree species in the landscape followed by Abies densa, the only gymnosperm tree species. Pieris villosa, Rhododendron thomsonii, Viburnum cordifolium and Rhododendron arboreum were among the least frequent species. This pattern of frequency class distribution suggests the heterogeneous composition of forest communities in the study area (S3 Table). The comparative histogram of standard values and observed values is depicted in Fig 3. Alternatively, the trees and saplings demonstrated a homogeneous distribution of frequency classes on separate analysis, while seedlings displayed decreasing trend of frequencies towards higher frequency classes that indicates heterogeneous occurrence of seedlings in study area.

Population structure
For tree species, the girth class wise distribution pattern gradually declined in density with increasing girth class. The lower girth class (CBH: < 30cm) including seedling and sapling together accounted for 56% of the total standing individuals followed by young tree individuals (CBH: 30-90cm) which contributed 37%. Similarly, lower girth class (CBH: 30-60 cm) contributed to the maximum total basal area (30%), then the CBH concurrently declined up to 180 cm CBH class. Despite, having only 2% contribution to density the trees of highest girth class (CBH: > 180 cm) contributed 28% to the total basal area (Fig 4). We encountered high girth classes of Abies densa trees in timberlines, even up to CBH of 530 cm.
The Simpson's index of dominance varied from 0.25 (site 6) to 0.49 (site 9) for trees and 0.32 (site 9) to 0.68 (site 4) for shrubs across the sites. Shannon-Weiner diversity index was between 0.77 (site 9) and 1.47 (site 6) for trees and from 0.51 (site 4) to 1.29 (site 9) for shrubs. In trees, the site 8 recorded maximum value for the Margalef's index of species richness (5.8)  Table).

Cluster analysis and species composition
The UPGMA agglomerative clustering based on the Jaccard's similarity matrix data is depicted in S2 Fig. Four clusters were defined based on the similarity observed in the presence-absence  matrix of 20 woody taxa among 9 sites. Clustering separated the sites of higher species richness (site 7 and site 8) from the remaining sites at about 43% of similarity level.

Interaction of environmental variables and vegetation composition
Canonical correspondence analysis diagram (CCA-tri-plot) results showed a total inertia of 1.302. The first axis alone explained 44.52% of total unexplained variance. Together, the first and second axes of the CCA data set explained 68.9% of the total inertia, which indicates a good correlation between woody taxa composition and environmental variables. The CCA ordination indicated that the axis 1 was primarily correlated with elevation and aspect; the axis 2 was partially correlated with slope; while the axis 3 was mainly correlated with humus and treeline distance (Fig 5). Of these, the most influential factors were elevation, and humus. The two relatively higher elevation sites, site 7 and site 9, and their respective vegetation are clustered to the right of the CCA plot, clearly separated from the low elevation site1 and site 6 ( Fig  5).  Timberline structure along latitudes in Khangchendzonga National Park, Eastern Himalaya subalpine forests which reaches up to the timberline and sometime up to the treeline. The species like G. pyroloides, P. villosa, R. hodgsonii and R. arboreum, showed a strong relation with the humus (Tables 1 and 2 and Fig 5).

Transect above timberline to tree line
It was observed that most of the 9 sites had distinct timberline and diffuse types of treelines. The plot-wise elevation trend of timberline and treeline of study area is depicted in S3 Fig. We observed the maximum distance (187.7 m) between timberline and treeline for site 8, at site 9 the treeline and timberline were located at the same elevation. Similarly, the location distribution of individuals varied across the sites. We recorded maximum number of individuals from the site 6 (seedling: 24; sapling: 8; tree: 4; dead tree: 0), followed by site 1 (seedling: 12; sapling: 3; tree: 14; dead tree: 6) (Fig 6A). A. densa regenerated well between timberline and treeline, though the number of seedlings declined with the increasing distance from the timberline. Interestingly, in each size-class we recorded the presence of adult trees and seedlings, but conversion rate of seedlings into saplings was poor along the transect, and above 150m distance from the timberline no sapling was observed (Fig 6B). Further, the presence of dead standing trees at higher distance (>180m) from the current timberline indicates elevation decline of current timberline in certain sites (Fig 6A and 6B).

Discussion
In Dzongri timberline zone, we observed Abies densa is the only canopy tree species, but it did not dominate at every site, and was absent from the highest elevation site (3989 m), where Rhododendron wightii was most prominent tree species. This Rhododendron occurs in eastern Nepal, Sikkim, Bhutan, S.E. Tibet to N. Burma, primarily in timberline ecotone between 3500 and 4500m asl. In Sikkim, it grows between the elevations of 4000-4500m asl [59]. A. densa grows in the eastern Himalayas between 2450-4000m asl on rocky steep slopes [60]. So, it is an subalpine forest species which goes up to timberline ecotone. Betula utilis, which predominates the Western Himalayan timberlines [33,39], was absent in Dzongri timberline, although we observed its few solitary trees on the western aspect below the timberline plots. Whereas, the Rhododendron spp. are prevalent throughout the eastern Himalayas [61]. In the present study area, at some sites no tree individuals were observed beyond timberline (forest with a least 30% tree crown density), so treeline and timberline were the same. In the others, a few individual of tree species adult trees, sapling and seedlings, mostly of A. densa were able to establish above treeline. The site 9 had abrupt timberline having no trees beyond timberline, perhaps due to growth limiting factors such as soil texture, topography [20]. In our study A. densa form the diffused treelines owing to fewer trees (Fig 6 and S3 Fig). On comparison, our study found low species richness than the other Himalayan forest types [34,36,62], which can be attributed to site characteristics and position of the study sites, as the present study was focused on ecotonal region viz. timberline, where harsh climatic conditions prevails, supporting to the survival of only few species. Yet, woody taxa can be compared with timberline communities of Kedarnath Wildlife Sanctuary, Uttarakhand, which ranged from 9 to15 species [39]. Across the longitudinal extent of the IHR, the proportion of trees increases by half and the shrubs increases by one-third, from west to east [63]. Similar trend in woody taxa is observed in other studies in North America [64], South Africa [65] and Italy [66].
The species distribution is considered homogenous when the numbers of individuals are the same in all parts of the community [34]. The woody taxa in Dzongri timberline had contagious distribution, which indicates harsh conditions in study area. The contagious distribution is more common in nature than random distribution (occurs only in a very uniform environment), and regular distribution (occurs when severe competition prevails between individuals) [67]. Here treeline occurred up to 4011 m, which is higher than treeline elevation of Kashmir and Uttarakhand in the western IHR (< 3500m). We recorded comparatively higher species diversity for woody taxa in KNP, which sustains varied moisture regime, supporting a diversity of trees, shrubs, and herbs [68]. However, the fewer number of trees decreases the inter-specific competition, resulting in to high degree of tree evenness [69], which we observed in Dzongri timberline. In our study, tree dominance increased with the elevation but tree diversity declined with increasing elevation; (at site 9, treeline and timberline are at the same Timberline structure along latitudes in Khangchendzonga National Park, Eastern Himalaya altitude, thus can be referred to as abrupt timberline. Thus, the nature of treeline can vary, depending upon micro-site factor within the same area. Nevertheless, within similar altitudinal range various co-factors such as aspect, topography and exposure can affect the forest composition [70]. The presence of a large number of seedlings and saplings in the timberline ecotone (Fig 4A  and 4B), indicates that in near future the treeline ecotone will become more dense [71]. This is also supported by occurrence of the higher weightage of juveniles and younger trees, especially Abies densa (> 59%), Sorbus microphylla (53%) and of Rhododendron lanatum (55%)s [72]. The seedlings of Abies densa contributed to highest regeneration potential (46%) under different size classes but their recruitment rate in to saplings was only 13% (S4 Fig). This might be attributed to the lower climatic threshold of seedlings than adults [73] and their greater sensitivity to other microclimatic factors. The poor regeneration along timberline ecotone and poor conversion rate of seedling to sapling has also been observed in Western Himalaya [39].
Density and frequency distribution of trees contributes to the structure of forests [36]. Tree density in our study are comparable to different forest types across the Himalaya [32,[34][35][36]39] and elsewhere in the world [74,75]. Such variations across the sites can be attributed to their microclimatic conditions, as monsoon and westerlies strongly influence the climate [76]. The relative abundance of species is one of the fundamental aspects of community structure [77] which are expressed in our the canonical correspondent analysis. Despite lower elevation range difference in our study (3787-3989 m asl) the overall woody taxa have shown a mid-elevation peak in species richness, which is commonly observed pattern in variety of ecosystems [78], and governed by a series of interactions between biological, climatic, historical and spatial factors [79,80].
Worldwide, timberlines do not show consistent directional shifts in response to warming climate [81]. In our study, elevation along with slope and humus was found to be the most determining factors for upward movement of tree species. The maximum average distance of treeline from timberline (>100m) was recorded from the sites having lower elevation range viz., site 8, site3, and site 1. Although, the site 6 was located at the lowest elevation range the 56˚average slope restricted the movement of tree species within 30.56 m. Notably, in high elevation sites, the presence of tree species beyond the timberline was limited to 5.3m in site 5 and no tree species was recorded above timberline in site 9. The available humus was also reduced significantly along the elevational gradients. Therefore, elevation may signify its importance in the determination of tree species distribution by influencing the microhabitats [33,34,82]. Apart from elevation, the position of the alpine timberline results from the interaction of several factors including climatic variables and disturbances along with forest succession [7,8,83]. Further, across the alpine timberlines, the competition between plant species is not usually as great as amongst lowland communities but the severity of other environmental factors might be much greater, viz. deep snow, short flowering season, heavy winds, lack of soil and unstable topography, etc. [84]. The formation of dense patches of small shrubs like Rhododendron setosum, R. anthopogon along with crooked, bent growth of R. campanulatum and Juniperus recurva as common across the krumholtz zone above timberline.
Out of five timberline boundaries' types of Jodlowski [20], we observed orographic boundaries in most of our study sites where trees were growing in the ridges besides the morphological timberline boundaries in site 6, where upper limit of the timberline got restricted due to steep slope or footsteps of rock faces.

Conclusion
The main conclusions of this study are as following: 1. Overall, Abies densa dominated the timberline vegetation, however, at several small undercanopy tree species (<10 m) namely, Sorbus microphylla, Rhododendron lanatum and Rhododendron wightii, far exceeded A. densa in cover.
2. The elevation of timberline varies between 3787-3989m asl and environmental factors such as humus, elevation and slope played an important role in shaping the vegetation composition as well as timberline boundaries of the landscape.
3. Presence of dead standing trees of A. densa above the current timberline indicates, that timberline was higher. It could be remnant treeline. Their higher position in the past can also be related with the prehistoric anthropogenic disturbance, as the Yak herders were using this area before governmental ban on grazing.
4. Fair regeneration of A. densa between the timberline and treeline, but its absence beyond the adult trees indicates that, the treeline and timberline ecotone can get densified in near future, but their upslope advancement is unlikely to occur.
5. Timberline, being one of the most sensitive ecotones and an indicator of warming climate, requires further comprehensive studies in order to understand its dynamics. Our study demonstrates a methodology for monitoring the status of timberlines in KNP; we recommend, replication of this/similar studies throughout Himalaya for comparative analysis of timberline structures.
Supporting information S1  Plots (P1 to P25) were sampled (each of 50x20 m) across 9 sites (mostly 3 plots a site) over about 20 km. In timberline, tree crown cover is at least 30%, treeline connects the highest trees of each transect, distances between timberline and tree ranges from 0 to 155 m. At many sites treeline and timberline are the same, thus distance being 0. None of the seedlings and saplings extent beyond treeline, so treeline is unlikely to advance, however, the treeline and timberline ecotone can get densified.