Spatial Pattern and Ecological Process Difference Analyses of the Boundary Habitats of a Treeline Patch: A Case Study from the Li Mountain, North China

: Treeline patches are among Earth’s most sensitive and are important model ecosystems for assessing climate change trends. To explore ecological factors that limit the species’ survival in treelines, the treeline patch of Li Mountain National Nature Reserve was selected as the research site. Pinus armandii ( P. armandii ), Betula albo-sinensis ( B. albo-sinensis ), and Betula utilis ( B. utilis ) were selected as research species based on their dominance. Two 50 m × 50 m plots were established separately from the upper and lower limits of the highest treeline patch for point pattern analysis. Five 10 m × 10 m quadrats per plot were sampled to investigate the ﬂora and environmental factors. The results showed that: (1) Slope and community composition at tree layer in quadrates had signiﬁcant differences between upper and lower limits. Pinus armandii had a greater population size at the upper limit. Seedling recruitment restricted population development for B. albo-sinensis at the lower limit and B. utilis at the upper limit and less regeneration of B. albo-sinensis at the upper limit. (2) More aggregation scales occurred at the upper limit, and P. armandii had more aggregation scales than the other two species at 0–25 m. The heterogeneity caused by density distribution affected P. armandii pattern at the upper limit, and heterogeneity of seed dispersal could explain species patterns in both limits. Distinctness of size difference may have an inﬂuence on inter-speciﬁc species correlations.


Introduction
Since the last century, global warming has been a significant problem, especially regarding vegetation composition and spatial distribution of an ecosystem in a high-altitude region [1,2]. Alpine treeline dynamics have proved important for predicting climate change patterns and for studying species response strategies [3]. The term "treeline" is typically defined as a transition zone extending from the uppermost part of a closed canopy forest to the scattered individuals [4]. The treeline patch extends from the lowermost of <50% coverage forest to an open meadow [5]. Climate research on high-elevation treelines has investigated environmental factors, such as high sunlight [6,7], low temperature [8,9], topography [10,11], water stress [12], fire disturbance [13], snow cover [14]. Species response research has focused specifically on vegetation growth form [15], vegetation growth rate [16,17], soil microbiomes [18,19] and plant physiology [12], as well as seedling recruitment and regeneration [20,21].
The responses of treelines on temporal and spatial scales are different. Climate factors determine elevation limits for treeline species on a global scale [6,8,9,22]. At a local scale, terrain and edaphic factors, such as slope, slope aspect, soil fertility, and depth, significantly influence treeline shift and vegetation structure in treeline patches [23]. Research on the relationship between the treeline pattern and topography in Changbai Mountain shows that factors such as high solar radiation and low snow cover created a high coverage rate of trees; the maximum value of tree coverage rates occur on low slopes at a lower elevation and occurs on medium slopes at higher elevations [11]. The edaphic environment's heterogeneity can determine treeline structure in the Krummholz zone due to decreased nutrient availability at higher elevations [24][25][26].
Additionally, seed dispersal and seedling survival are key processes in shaping species distribution patterns; inadequate accumulated temperature in summer decreases seeds' productivity and viability [27]. Some research shows that density-dependent mortality often occurs at the seedling stage [22,28], and the seedling establishment rates in alpine treeline ecotones may be low and may vary among species [20]. However, other research indicates that the increase in seedling facilitation is due to the environmental severity of treelines [29], which is accorded to the stress hypothesis [30], and the intraspecific aggregated intensity will be higher at the harsh upper treeline limit [31]. Hence, multiple ecological processes shape species patterns at various scales.
The treeline patch closest to the top of the mountain has a compassionate response to climate change [9]. Vegetation changes are indicative of environmental changes [1]; hence, community differences and ecological processes differences between the upper limit and lower limits in the highest treeline patch might be significant. In ecological studies, spatial point pattern analysis is an informative method that has mainly been applied to explaining tree distributions [32] and structures [33], as well as interactions between trees and shrubs at the treeline [32]. A study site in the Shunwangping district of Li Mountain was set; P. armandii, B. albo-sinensis, and B. utilis as common species in the treeline area were selected to explore it as target species.
Our hypotheses are as follows: firstly, treeline patches can be formed at an altitude of 2300 m, which is not a particularly high altitude), we speculate that even if the area of the treeline patch is small, there are significant differences in the environmental factors between the upper and lower limits of the treeline patches and significant differences in community structure and species spatial patterns. Secondly, there are differences in the population structure of the three species at the upper and lower limits of the treeline patches, which contain large trees, medium trees, and small trees; hence, we surmise that the spatial patterns of the three species are affected by different ecological processes, the population composed of smaller individuals is significantly affected by environmental heterogeneity. Thirdly, the three species have the characteristics of seed dispersal and propagation; we hypothesize that their patterns are all related to seed dispersal; in addition, the population structure of the three species differences will influence the interspecific correlation between species by the distinctness of tree size differences.

Study Area and Field Sampling
The study area is located in Li Mountain National Nature Reserve (111 • 51 -112 • 5 35 E, 35 • 16 30 -35 • 27 20 N), which belongs to a part of Zhongtiao Mountains and is located in Shanxi Province, northern China. It is one of the most biodiverse regions in North China, with complete natural vegetation ( Figure 1). Sampling plots were established in the Shunwangping Scenic spot, the only primordial forest in North China. Peaks and valleys of Shunwangping were formed by the Himalayan mountains folding orogeny. The highest peak of Shunwangping has an elevation of 2358 m. The annual precipitation is 600-800 mm. The study area has a temperate continental monsoon climate. According to the data from the Yuanqu County Meteorological Bureau, the temperature in the hottest month is 27 • C~28 • C, the temperature in the coldest month is −1 • C~2 • C, and the extremely hottest temperature is 32 • C [34]. The accumulated temperature greater than 10 • C is 4160 • C. Most soils are subalpine meadow soil; a few are brown forest soil at >1800 m. From the top to the base of Li Mountain, soils transition from mountain meadow Land 2022, 11, 2064 3 of 18 soils to brown forest soils, then to montane-leached cinnamon soil, and finally to montane cinnamon soils [35]. Li Mountain is a warm temperate deciduous broadleaf forest belt, with Betula platyphylla, P. armandii, Pinus tabulaeformis, Platycladus orientalis, Quercus variabilis, and Quercus wutaishanica as the predominant tree vegetation [36]. Among them, P. armandii, B. albo-sinensis, and B. utilis are the most abundant trees in the treeline patch. month is 27 °C~28 °C, the temperature in the coldest month is −1 °C~2 °C, and the extremely hottest temperature is 32 °C [34]. The accumulated temperature greater than 10 °C is 4160 °C. Most soils are subalpine meadow soil; a few are brown forest soil at >1800 m. From the top to the base of Li Mountain, soils transition from mountain meadow soils to brown forest soils, then to montane-leached cinnamon soil, and finally to montane cinnamon soils [35]. Li Mountain is a warm temperate deciduous broadleaf forest belt, with Betula platyphylla, P. armandii, Pinus tabulaeformis, Platycladus orientalis, Quercus variabilis, and Quercus wutaishanica as the predominant tree vegetation [36]. Among them, P. armandii, B. albo-sinensis, and B. utilis are the most abundant trees in the treeline patch. In early July 2017, two sampling plots of 0.25 ha (50 m × 50 m) from the highest treeline patch of Shunwangping, were set at the highest peak in Li Mountain. One sampling plot (N 35°24.999, E 111°57.857) at an elevation of 2300 m was at the upper limit of the treeline patch, and the other plot (N 35°25.227, E 111°57.729), at an elevation of 2250 m, was located at the lower limit of the treeline patch. Differences in environmental factors, community dissimilarity, population development of treeline patches and spatial point patterns were separately analyzed in each sampling plot.
Five 10 m × 10 m quadrats were established in each sampling plot, with four quadrats located at midpoints from the midpoint of the sampling plot to four corners, and one quadrat was the midpoint of the sampling plot ( Figure 2). Elevation was measured by Handy GPS in two plots; slope and slope aspect were measured by compass in each quadrate; temperature and pH of surface soil layer (0-10 cm) were measured by a handheld pH meter; an earth borer measured soil layer thickness, a steel ruler measured litter layer thickness, and tree coverage was estimated by canopy projection percentage at noon. In the case of a community survey, the interest objects of the community were tree and shrub layers. Vegetation below 5 m with no obvious main stem was defined as a shrub (Technical guidelines for biodiversity monitoring-terrestrial vascular plants from National Environmental Protection Standards of the People's Republic of China). Therefore, tree In early July 2017, two sampling plots of 0.25 ha (50 m × 50 m) from the highest treeline patch of Shunwangping, were set at the highest peak in Li Mountain. One sampling plot (N 35 • 24.999, E 111 • 57.857) at an elevation of 2300 m was at the upper limit of the treeline patch, and the other plot (N 35 • 25.227, E 111 • 57.729), at an elevation of 2250 m, was located at the lower limit of the treeline patch. Differences in environmental factors, community dissimilarity, population development of treeline patches and spatial point patterns were separately analyzed in each sampling plot.
Five 10 m × 10 m quadrats were established in each sampling plot, with four quadrats located at midpoints from the midpoint of the sampling plot to four corners, and one quadrat was the midpoint of the sampling plot ( Figure 2). Elevation was measured by Handy GPS in two plots; slope and slope aspect were measured by compass in each quadrate; temperature and pH of surface soil layer (0-10 cm) were measured by a handheld pH meter; an earth borer measured soil layer thickness, a steel ruler measured litter layer thickness, and tree coverage was estimated by canopy projection percentage at noon. In the case of a community survey, the interest objects of the community were tree and shrub layers. Vegetation below 5 m with no obvious main stem was defined as a shrub (Technical guidelines for biodiversity monitoring-terrestrial vascular plants from National Environmental Protection Standards of the People's Republic of China). Therefore, tree coverage, quantity, and DBH (diameter at breast height) were recorded in the tree layer, and the height, coverage, and quantity of shrubs were recorded in the shrub layer. Land 2022, 11, x FOR PEER REVIEW coverage, quantity, and DBH (diameter at breast height) were recorded in the tr and the height, coverage, and quantity of shrubs were recorded in the shrub laye More than 90% of the tree layer was constituted of P. armandii, B. albo-sinens B. utilis; therefore, the three species were chosen as target species for population d ment and spatial point pattern analysis. It identified and recorded DBH and coo of all independent trees with a diameter of at least 1 cm at breast height (DBH) in t species by measuring tape [37]. DBH was for population development analysis, a dinates of the three species were obtained for the spatial point pattern analysis vided the sampling plot (50 m × 50 m) into grids (10 m × 10 m) and measured each in the grids with measuring tape. Finally, individual locations in 10 m × 10 m w verted to coordinates of individuals in 50 m × 50 m ( Figure 2).

Environmental Factors Analysis
Slope, soil layer thickness, litter layer thickness, pH, the temperature of the soil layer, and tree coverage were surveyed in the five quadrats of each sampling unpaired t-test was used to test differences for the above factors between the two the p < 0.05 level. The slope aspect and elevation were measured in the centre of × 50 m sampling plot. The number of dead trees in the 50 m × 50 m sampling plot recorded. More than 90% of the tree layer was constituted of P. armandii, B. albo-sinensisis, and B. utilis; therefore, the three species were chosen as target species for population development and spatial point pattern analysis. It identified and recorded DBH and coordinates of all independent trees with a diameter of at least 1 cm at breast height (DBH) in the three species by measuring tape [37]. DBH was for population development analysis, and coordinates of the three species were obtained for the spatial point pattern analysis. We divided the sampling plot (50 m × 50 m) into grids (10 m × 10 m) and measured each location in the grids with measuring tape. Finally, individual locations in 10 m × 10 m were converted to coordinates of individuals in 50 m × 50 m ( Figure 2).

Environmental Factors Analysis
Slope, soil layer thickness, litter layer thickness, pH, the temperature of the surface soil layer, and tree coverage were surveyed in the five quadrats of each sampling plot. An unpaired t-test was used to test differences for the above factors between the two limits at the p < 0.05 level. The slope aspect and elevation were measured in the centre of the 50 m × 50 m sampling plot. The number of dead trees in the 50 m × 50 m sampling plot was also recorded.

Importance Value Calculation
Species in different layers (tree and shrub) have different ecological functions, and the importance value was used to evaluate species status, with importance value equations as follows [38]: Importance value (tree) = (relative coverage + relative abundance + relative dominance)/3 Importance value (shrub) = (relative height + relative coverage + relative abundance)/3 Mean importance value = Sum of importance values from all five quadrats/5 Relative abundance = Quantity of individuals of a given species/Total quantity of individuals of all species Relative dominance = Dominance of a species/Total dominance of all species Relative coverage = Coverage of a species/Coverage of total species Relative height = Average height of a species/Average height of all species.

Community Similarity Analysis
PERMANOVA was used to carry out a community similarity test separately in the tree and shrub layers between two sampling plots. First, the Bray-Curtis dissimilarity was calculated by using the importance value of every species in the quadrats, which is a good tool for community composition analysis. Then, the "Pseudo-F statistic" method was used to test the similarity difference in two sampling plots based on the Bray-Curtis dissimilarity values.
The Bray-Curtis dissimilarity [39] (Anderson 2001): where i and j are two quadrat codes, y is the species importance value, and k is the total species richness of two quadrats, and different layers have different k values.
where SSB is the sum of squares between groups, SSW is the sum of squares within groups, t is the number of groups, and N is the number of total samples.

DBH Distribution Analysis
To assess the status of population development in the survey area, DBH data from the sampling plots were used. The DBH of every individual that met diameter criteria was recorded, with categories designated at diameters of 5 cm, 10 cm, 20 cm, and 30 cm to divide the DBH values into five bins. The quantity of individuals of each species for each DBH bin was also recorded.

Point Pattern Analysis
The theory of the spatial point process has been widely used to fit the specific real distribution of corresponding species by adopting many ecological hypothesis models. Considering the community traits in this study, we decided to choose four null models (complete spatial randomness, heterogeneous Poisson process, and homogeneous Thomas process, inhomogeneous Thomas process) for analysis. The pair-correlation function g(r) was calculated for each analysis. Univariate pair correlation is defined as the expected point density at distance r from a typical point of the pattern, relative to the intensity λ of the pattern [40]. It is related to the K function, whereby: where K(r) is defined as the ratio of the expected number of the other points to the point density of the sample within a circle with an arbitrary point as the centre and radius r [41]. The univariate g(r) was used in the analysis of a single species. When the value of g(r) is Land 2022, 11, 2064 6 of 18 larger than 1, the individuals in the sample have an aggregated distribution at distance r. When the value of g(r) is smaller than 1, individuals are regularly distributed at distance r. When the value of g(r) is equal to 1, individuals are randomly distributed at distance r. g(r) could extend to bivariate functions g 12 (r), which describes the association (interaction) between two different types of species; values of the g 12 (r) function that significantly deviates from 1 indicate a significant attraction (or repulsion) between the two types of points [42] (Dixon 2002). When g 12 (r) > 1, it shows a positive correlation; when g 12 (r) = 1, it shows no correlation, and when g 12 (r) < 1, it shows a negative correlation.
A is the area of the sampling plot, n indicates the total point in the sampling plot, r is the distance between two points, u ij is the distance between i and j when u ij < r, I r u ij = 1, or u ij > r, I r u ij = 0; w ij indicates the ratio of a circle with i as the centre point and u ij as the radius to A.
For a single species or pair species, the pair correlation function in the null model for ecological process analysis was performed. The Monte-Carlo approach was applied to test significant g(r) values at the r scale. Each of the n = 199 simulations of a point process underlying the null model generated a summary statistic (e.g., a pair correlation function g 12 (r); simulation envelopes with α = 0.01 were calculated in the 199 simulations) [43]. A significant departure from the null model occurred at scale r if the test statistic was outside the simulation envelopes. This approach allowed us to approximately assess scale effects for illustrative purposes and to determine the type of significant effect. The complete spatial randomness process, Heterogeneous Poisson process, Homogeneous Thomas process, and Inhomogeneous Thomas process were used as null models [37,[44][45][46].

Model Test
To avoid problems caused by simultaneous inference [47], we evaluated the overall ability of a given null model to describe the data by means of a goodness-of-fit test (GoF) [48]. A GoF was used to test whether a given null model fit a given summary statistic of the observed data over a given distance interval (r min , r max ) [48]. The u i statistics were calculated for the observed data (i = 0) and simulated data (i = 1, . . . , 199), and the rank of u 0 among all u i was determined. If the rank of u 0 was larger than199, then the data showed a significant departure from the null model (across the distances of interest) with the error rate p < 0.05.

Software Analysis and Figure Processing
The "adonis2" function in the "vegan" package was applied for the PERMANOVA test for determining community composition differences. The "t.test" function in the "stats" package was applied for the analysis of the environmental factors between sites. The "Spatstat" package was applied for the analysis of spatial point pattern analysis. All analyses and the figures were plotted in R Project 3.4.3 (R Core Development Team 2017).

Environmental Conditions within Sampling Plots
The lope difference was significant (p < 0.05), and other environmental factors differences were not significant; soil layer thickness, litter layer thickness, pH of the surface soil layer, and tree coverage were greater at the upper limit than at the lower limit. Still, the surface soil layer temperature was greater at the lower limit (Table 1).

Composition and Similarity Differences in Tree Layer and Shrub Layer
We detected significant differences in community composition in the tree layer, but no significant differences in the shrub layer between the upper and lower limits were found. There were five tree species and six shrub species at both limits of the treeline patch, and the average importance values of P. armandii and B. utilis were separately at the largest at the upper limit and lower limit, followed by B. albo-sinensis at both limits. In addition, Cerasus clarofolia, B. albo-sinensis, and Salix wallichiana were observed in both tree and shrub phenotypes. The B. albo-sinensis shrubs had a greater importance value than the shrubs of most other species, especially at the upper limit (Table 2).

DBH Numerical Distribution
The individual distributions in Figure 3 show that the population of P. armandii was greatest at the upper limit, and that B. utilis had the largest quantity at the lower limit. The interval of 10 cm ≤ D ≤ 20 cm was the common DBH range for the three species. Pinus armandii was concentrated in a smaller DBH interval; most individual DBH values of P. armandii were distributed within the range of 10 cm ≤ D ≤ 20 cm at the upper and lower limits. The DBH values of B. albo-sinensis appeared in every DBH bin, the largest quantity of individual DBHs having D ≤ 5 cm at the upper limit, DBH values of B. albo-sinensis distributed within the range of 5 cm ≤ D ≤ 30 cm at the lower limit. Betula utilis was concentrated in larger DBH intervals (10 cm ≤ D ≤ 30 cm) at both limits, with a larger number at the lower limit than at the upper limit and individual DBHs having 5 cm ≤ D ≤ 10 cm only at the lower limit.

DBH Numerical Distribution
The individual distributions in Figure 3 show that the population of P. armandii w greatest at the upper limit, and that B. utilis had the largest quantity at the lower limit. T interval of 10 cm ≤ D ≤ 20 cm was the common DBH range for the three species. Pin armandii was concentrated in a smaller DBH interval; most individual DBH values of armandii were distributed within the range of 10 cm ≤ D ≤ 20 cm at the upper and low limits. The DBH values of B. albo-sinensis appeared in every DBH bin, the largest quant of individual DBHs having D ≤ 5 cm at the upper limit, DBH values of B. albo-sinen distributed within the range of 5 cm ≤ D ≤ 30 cm at the lower limit. Betula utilis was co centrated in larger DBH intervals (10 cm ≤ D ≤ 30 cm) at both limits, with a larger numb at the lower limit than at the upper limit and individual DBHs having 5 cm ≤ D ≤ 10 only at the lower limit.

Completely Spatially Random Null Model Analysis
The spatial distributions of three dominant species were as shown in Figure 4. Pin armandii distributed more on the upper limit of the treeline patch, and its aggregation ar was larger. The distribution of B. albo-sinensis and B. utilis on the upper and lower lim of the treeline patch were scattered, and both their aggregation area were small. The d tribution of the three species in the sampling plot did not indicate much overlap.

Completely Spatially Random Null Model Analysis
The spatial distributions of three dominant species were as shown in Figure 4. Pinus armandii distributed more on the upper limit of the treeline patch, and its aggregation area was larger. The distribution of B. albo-sinensis and B. utilis on the upper and lower limits of the treeline patch were scattered, and both their aggregation area were small. The distribution of the three species in the sampling plot did not indicate much overlap.  In comparing the real g(r) values with the envelope g(r) values and GoF test (Table 3), the results were as follows: There are more real g(r) values beyond the envelope at the upper limit than those at the lower limit for P. armandii and B. albo-s while the g(r) values of B. utilis beyond the upper envelope were less. There ar points at r < 10 m on the upper limit and at r < 15 m on the lower limit. The aggrega P. armandii occurred on the broadest scales, Betula albo-sinensis aggregated at the s broadest scale, and aggregation of B. utilis occurred at the smallest scale ( Figure  CSR could not fit the real patterns of the three species well (p = 0.005). Note: "CSR", "IP", "HT" and "IT" are abbreviations of Complete Spatial Randomness, Inh neous Poisson, Homogeneous Thomas and Inhomogeneous Thomas, respectively. * means cant differences between the expected line simulated by the null model and the observed lin 0.05, ** means extremely significant differences between the expected line and the observe In comparing the real g(r) values with the envelope g(r) values and GoF test results (Table 3), the results were as follows: There are more real g(r) values beyond the upper envelope at the upper limit than those at the lower limit for P. armandii and B. albo-sinensis, while the g(r) values of B. utilis beyond the upper envelope were less. There are a few points at r < 10 m on the upper limit and at r < 15 m on the lower limit. The aggregation of P. armandii occurred on the broadest scales, Betula albo-sinensis aggregated at the secondbroadest scale, and aggregation of B. utilis occurred at the smallest scale ( Figure 5). The CSR could not fit the real patterns of the three species well (p = 0.005).

Heterogeneous Poisson (HP) Null Model Analysis
Environmental factors produced no significant differences between the upper and lower limits (Table 1). In reference to density in Figure 6, the three species' densities were selected as background information to simulate heterogeneity.
At the upper limit, the GoF test results in Table 3 and comparisons between the real values with expected values in Figure 7 show that B. albo-sinensis and B. utilis have many g-values beyond the upper envelope, and the result of the GoF test (p < 0.05) showed their expected lines were quite different from the actual lines at 0-25 cm. While the total actual line of P. armandii was surrounded by envelopes and had a similar tendency at 0-25 m at the upper limit, the GoF test value (p = 0.115) also showed that the HP simulated pattern was not significantly different from the real pattern of P. armandii at 0-25 m on the upper limit. However, at the lower limit, the actual lines of the three species were significantly different from simulated lines at 0-25 m, respectively (p < 0.05).

Heterogeneous Poisson (HP) Null Model Analysis
Environmental factors produced no significant differences between the upper and lower limits (Table 1). In reference to density in Figure 6, the three species' densities were selected as background information to simulate heterogeneity.  At the upper limit, the GoF test results in Table 3 and comparisons between the real values with expected values in Figure 7 show that B. albo-sinensis and B. utilis have many g-values beyond the upper envelope, and the result of the GoF test (p < 0.05) showed their expected lines were quite different from the actual lines at 0-25 cm. While the total actual line of P. armandii was surrounded by envelopes and had a similar tendency at 0-25 m at the upper limit, the GoF test value (p = 0.115) also showed that the HP simulated pattern was not significantly different from the real pattern of P. armandii at 0-25 m on the upper limit. However, at the lower limit, the actual lines of the three species were significantly different from simulated lines at 0-25 m, respectively (p < 0.05).

Homogeneous Thomas (HT) Null Model Analysis
Observations from g(r) values and GoF test values (Table 3) showed the difference between the observed line and fitted expected line was significant (p < 0.05) only for P. armandii at the upper limit, where many of the g(r) values of P. armandii were greater than g-values of the upper envelopes at r < 15 m scales. The results of the GoF test at the lower Figure 7. HP null model analysis of spatial patterns of P. armandii, B. albo-sinensis, and B. utilis at the upper limit and the lower limit of the treeline patch. Note: the red points are the g-values beyond the upper envelope at a corresponding scale, the blue points are the g-values under the lower envelope at a corresponding scale.

Homogeneous Thomas (HT) Null Model Analysis
Observations from g(r) values and GoF test values (Table 3) showed the difference between the observed line and fitted expected line was significant (p < 0.05) only for P. armandii at the upper limit, where many of the g(r) values of P. armandii were greater than g-values of the upper envelopes at r < 15 m scales. The results of the GoF test at the lower limit were not significant for P. armandii (p = 0.675). The observed lines of B. albo-sinensis and B. utilis were all located inside the envelopes in both sampling plots, and the differences between their observed lines with corresponding expected lines were not significant from 0-25 m with the GoF test values (p > 0.05) ( Figure 8).
limit were not significant for P. armandii (p = 0.675). The observed lines of B. albo-sinensis and B. utilis were all located inside the envelopes in both sampling plots, and the differences between their observed lines with corresponding expected lines were not significant from 0-25 m with the GoF test values (p > 0.05) (Figure 8). .

Inhomogeneous Thomas (IT) Null Model Analysis
Most observed g(r) values were included in the 99% confidence interval for the three species, the real g(r) values of B. utilis were smaller than g(r) values of lower envelope at r = 2 m and real g(r) values of P. armandii were greater than g(r) values of the upper envelope at r = 16 m (Figure 9). GoF test values in Table 3 showed that the differences between their observed lines with the corresponding expected lines for three species were not significant at the 0-25 m scale at both the upper and lower limits (p > 0.05).

Inhomogeneous Thomas (IT) Null Model Analysis
Most observed g(r) values were included in the 99% confidence interval for the three species, the real g(r) values of B. utilis were smaller than g(r) values of lower envelope at r = 2 m and real g(r) values of P. armandii were greater than g(r) values of the upper envelope at r = 16 m (Figure 9). GoF test values in Table 3 showed that the differences between their observed lines with the corresponding expected lines for three species were not significant at the 0-25 m scale at both the upper and lower limits (p > 0.05).

Homogeneous Poisson Null Model Analysis Inter-Specific Correlation
As for the upper limit, the real g 12 (r) values for P armandii and B. albo-sinensis were greater than the g 12 (r) values on the upper envelope at 10-17 m. The real g 12 (r) values for P. armandii and B. utilis were smaller than the g 12 (r) values on the lower envelope at a total 0-25 m. The real g 12 (r) values for B. albo-sinensis and B. utilis were smaller than the g 12 (r) values on the lower envelope only at a few scales of 10 m < r < 15 m ( Figure 10).
As for the lower limit, the real g 12 (r) values for P. armandii and B. albo-sinensis were smaller than the g 12 (r) values on the lower envelope at r < 10 m and r = 25 m. The real g 12 (r) values for P. armandii and B. utilis were smaller than the g 12 (r) values on the lower envelope at 0-20 cm, and the significant g 12 (r) values were fewer than those at the upper limit of treeline patch. The real g 12 (r) values for B. albo-sinensis and B. utilis were greater than the g 12 (r) values on the upper envelope at r = 15 m ( Figure 10).
The results showed that P. armandii and B. albo-sinensis had significantly positive correlations at the upper limit but a significantly negative correlation at the lower limit. Pinus armandii and Betula utilis had significantly negative correlations at both limits, and B. albo-sinensis and B. utilis had significant correlations on few scales.

Homogeneous Poisson Null Model Analysis Inter-Specific Correlation
As for the upper limit, the real g12(r) values for P armandii and B. albo-sinensis were greater than the g12(r) values on the upper envelope at 10-17 m. The real g12(r) values for P. armandii and B. utilis were smaller than the g12(r) values on the lower envelope at a total 0-25 m. The real g12(r) values for B. albo-sinensis and B. utilis were smaller than the g12(r) values on the lower envelope only at a few scales of 10 m < r < 15 m ( Figure 10).
As for the lower limit, the real g12(r) values for P. armandii and B. albo-sinensis were smaller than the g12(r) values on the lower envelope at r < 10 m and r = 25 m. The real g12(r) values for P. armandii and B. utilis were smaller than the g12(r) values on the lower envelope at 0-20 cm, and the significant g12(r) values were fewer than those at the upper limit of treeline patch. The real g12(r) values for B. albo-sinensis and B. utilis were greater than the g12(r) values on the upper envelope at r = 15 m ( Figure 10).
The results showed that P. armandii and B. albo-sinensis had significantly positive correlations at the upper limit but a significantly negative correlation at the lower limit. Pinus armandii and Betula utilis had significantly negative correlations at both limits, and B. albosinensis and B. utilis had significant correlations on few scales.

The Differences of Community Structure and Species Spatial Pattern at the Upper and Lower Limits of Treeline Patch
Slope differences were significant between the upper and lower limits of the treeline patch (Table 1, p < 0.05); there was significant differences in the community structure in the tree layer, but not in the shrub layer; it was not consistent with the evidence that slope had significant effects on the shrub and herbage [49,50]. In the tree layer, community structure was mainly attributed to three dominant species (P. armandii, B. albo-sinensis, and B. utilis)-they were at different statuses and stages of population development. Pinus armandii was a developing population-most of its DBH was distributed on a scale of smaller than 20 cm; the abundance of P. armandii was greater at the upper limit than that at the lower limit, which showed that population size was greater at the upper limit than that at the lower limit. The reason might be that the highly viable seeds of P. armandii help it survive in occasionally benign conditions [51]; in addition, by comparing the environmental factors of upper limit and lower limit. We also noticed that even the differences in soil thickness, litter thickness, and tree coverage were not significant, but they were superior at the upper limit to that at the lower limit ( Table 1). The DBH of B. albo-sinensis distributed on a scale of smaller than 5 cm at the upper limit than that at the lower limit, population structure tended to be younger individuals. More individuals with shrub lifeform were at the upper limit, the shrub phenotype of B. albo-sinensis was a harm avoidance strategy in a stressful habitat [52]. All DBH of B. utilis was greater than 5 cm, and few shrub individuals were at the lower limit; hence, the reason might be that the recruitment and regeneration in recent years were scarce at the upper and lower limits [53]. Hence, the lifeforms and population development of the three species affected community structure in tree layers.
The occurrence scales of the three species aggregation patterns were all greater at the upper limit than those at the lower limit. Pinus armandii had the broadest aggregation scale over that of B. albo-sinensis and B. utilis. Most of the P. armandii and B. albo-sinensis were both young individuals at the upper limit, where crowded young individuals initially promoted each other's survival and development [54]. The mutual benefits in their own populations were important for surviving in harsh environments [55]. Most individuals of the B. utilis were old trees, and it had a tolerance for low temperatures and poor habitat; hence, its pattern was a random distribution at most scales [56].
The above results accorded to our first hypothesis, and slope as an environmental factor differed significantly. Community structure and spatial pattern also had significant differences between the upper and lower limits.

Effect of Community Density Heterogeneity on Seedling Pattern
Environmental heterogeneity simulated by the density of the three species had a significant effect on the pattern of P. armandii at the upper limit; the population of P. armandii consisted mostly of medium individuals at the two limits, they were more sensitive to their surroundings [44], such as soil temperature, summer drought, and competition between the shrubs affected their survival in the treeline [57,58]. More small individuals of P. armandii were at the lower limit, and more small individuals of B. albo-sinensis were at the upper limit, but the community heterogeneity had insignificant effects on their patterns; hence, we speculated that medium tree patterns were more sensitive to density heterogeneity, which was inconsistent with the second hypothesis.

Seed Dispersal Effect and Size Effect on Species Patterns
Simulated seed dispersal of P. armandii in homogeneous conditions differed significantly from its real pattern at the upper limit; high vegetation density was an important restriction factor for seed regeneration [57]. Simulated inhomogeneity seed dispersal could well explain the patterns of three species. The pattern of the P. armandii simulated by heterogeneity and seed dispersal simultaneously could especially explain the real pattern. Pinus armandii only had seed propagation in its life history; seed productivity and dispersal were important ecological processes for population development [59], which were impeded by species distributions and the habitat conditions around it, as death density dependency [60,61] and environmental stress theory [62] stated. Hence, integrating heterogeneity characteristics and seed dispersal processes could have fitted the three species' patterns well in our study.
The larger individuals of B. utilis distributed in the limits and competitive relationships among the canopy trees played a key role in generating stand structure [63]. The significant negative correlations between P. armandii and B. utilis occurred at both limits; this was also the case with B. albo-sinensis and B. utilis. Previous studies in the mixed forests of the Qinling Mountains suggested that less intense bamboo could impede Betula regeneration [64]; hence, both P. armandii and B. albo-sinensis had an impact on the regeneration of B. utilis.
While the correlation between B. albo-sinensis and P. armandii was positive at the upper limit and negative at the lower limit. At the upper limit, most individuals of B. albo-sinensis were smaller than P. armandii; the inverse situation happened at the lower limit. The size-adjusted inter-specific correlation, when the neighbouring species were larger than the target species, had a positive effect on it at first, then a negative effect when the neighbour species' size reached a certain scale [62]. In addition, the DBH differences between P. armandii and B. utilis were obvious in both limits-the size effect leads to negative correlations between them [62].
As the third hypothesis stated that seed dispersal had significant effects on species patterns, the differences in size affected the inter-specific correlation, with a negative correlation often occurring between two species with a greater size difference.

Conclusions
Community composition in the tree layer varied significantly. The environmental factors we measured demonstrated that the slope significantly differed between the upper limit and the lower limit of the highest treeline patch. Pinus armandii had a greater population size at the upper limit, and B. utilis had a greater population size at the lower limit. The aggregation scales were greater at the upper limit for all three species. The density heterogeneity impacted the P. armandii spatial pattern at the upper limit with most of the medium individuals in its population. Inhomogeneity seed dispersal processes significantly affected the patterns of B. albo-sinensis, B. utilis and P. armandii at both limits. Pinus armandii and Betula albo-sinensis had negative interspecific correlations, which often occurred between the species with obviously different sizes. Density, seed dispersal, and unbalanced size structure played important roles in treeline patch development.
For the non-altitudinal tree line area, the slope is an important environmental factor that restricts the forest line community, and there are differences between the upper and lower limits for tree communities in the tree line. At the upper limit of the tree line, P. armandii regenerated quickly, and its population is in its development stage. Due to the population density and interspecific relationship affecting its survival, it can be considered to properly decrease the number of individuals. The young trees of B. utilis regenerate slowly-according to the existing literature, B. utilis are common constructive species in a tree line, and they could create a better living environment for other species. However, larger individuals restrict the survival of other species, so they can be planted under harsh conditions with less vegetation. The scale and size had significant effects on the relationships between the species by considering the differences in individual size and planting distance among the species. In addition, interspecific competition at some scales and some growth stages can be avoided. Seed dispersal is an important ecological process for the natural breeding of offspring, and even in harsh environments, seed dispersal characteristics should be considered for population development.