Tree-to-tree interactions slow down Himalayan treeline shifts as inferred from tree spatial patterns

17 rectangular plots at treeline ecotone with a size 30 m × 150 m (except M1, 30 m × 100 m) were established on gentle slopes ranging from 12º to 22º, and placed away from previously affected areas by landslides or avalanches. All plots included the upper treeline and the forest limit which was defined as the elevation of closed forest with at least 30 % canopy cover (Camarero & Gutiérrez, 2004; Liang et al., 2011; Liang et al., 2016). The plots had its longer side (y axis) located parallel to the elevation gradient. For all plots, coordinates (x, y) = (0, 0) were situated in the left, bottom corner (Camarero & Gutiérrez, 2004). The geographic information of all four corners was noted by using a Global Positioning System (GPS) with ± 4 m accuracy. The location of each individual of targeted species was based on the x, y reference points. Coordinates were measured to the nearest 0.1 m (Liang et al., 2016). To quantify the spatio-temporal patterns of Himalayan treelines, the changes in tree density and treeline elevation were calculated for 50-year age-classes (1-50 yrs., 51-100 yrs., and 101-150 yrs.). We used the “spatstat” package (Baddeley, Rubak, & Turner, 2015) in R program (R Development Core Team 2018) to calculate nearest-neighbour (NN) distances for all trees grouped in the three age classes. To summarize the spatial distribution patterns of each age-class using a simple measure, we used the R statistic (clustering intensity) calculated in PAST software (Hammer, Harper, & Ryan, 2001). If a Poisson random process and an exponential NN distribution are assumed, μ may be defined as: (1) where A is the area and n the number of trees. Then, the R statistic is calculated as: (2) A spatial pattern may be described as clustered, random (Poisson) or overdispersed for R < 1, R ~ 1, and R > 1, respectively (Clark & Evans, 1955).