Spatial Distribution Patterns in the Very Rare and Species-Rich Picea chihuahuana Tree Community (Mexico)

The very rare Mexican Picea chihuahuana tree community covers an area of no more than 300 ha in the Sierra Madre Occidental. This special tree community has been the subject of several studies aimed at learning more about the genetic structure and ecology of the species and the potential effects of climate change. The spatial distribution of trees is a result of many ecological processes and can affect the degree of competition between neighbouring trees, tree density, variability in size and distribution, regeneration, survival, growth, mortality, crown formation and the biological diversity within forest communities. Numerous scale-dependent measures have been established in order to describe spatial forest structure. The overall aim of most of these studies has been to obtain data to help design preservation and conservation strategies. In this study, we examined the spatial distribution pattern of trees in the P. chihuahuana tree community in 12 localities, in relation to i) tree stand density, ii) diameter distribution (vertical structure), iii) tree species diversity, iv) geographical latitude and v) tree dominance at a fine scale (in 0.25 ha plots), with the aim of obtaining a better understanding of the complex ecosystem processes and biological diversity. Because of the strongly mixed nature of this tree community, which often produces low population densities of each tree species and random tree fall gaps caused by tree death, we expect aggregated patterns in individual Picea chihuahuana trees and in the P. chihuahuana tree community, repulsive Picea patterns to other tree species and repulsive patterns of young to adult trees. Each location was represented by one plot of 50 x 50 m (0.25 ha) established in the centre of the tree community. The findings demonstrate that the hypothesis of aggregated tree pattern is not applicable to the mean pattern measured by Clark-Evans index, Uniform Angle index and Mean Directional index of the uneven-aged P. chihuahuana trees and P. chihuahuana tree community and but to specific spatial scales measured by the univariate L-function. The spatial distribution pattern of P. chihuahuana trees was found to be independent of patches of other tree species measured by the bivariate L-function. The spatial distribution was not significantly related to tree density, diameter distribution or tree species diversity. The index of Clark and Evans decreased significantly from the southern to northern plots containing all tree species. Self-thinning due to intra and inter-specific competition-induced mortality is probably the main cause of the decrease in aggregation intensity during the course of population development in this tree community. We recommend the use of larger sampling plots (> 0.25 ha) in uneven-aged and species-rich forest ecosystems to detect less obvious, but important, relationships between spatial tree pattern and functioning and diversity in these forests.

In this study, we examined the spatial distribution pattern of trees in the P. chihuahuana tree community in 12 localities, in relation to i) tree stand density, ii) diameter distribution (vertical structure), iii) tree species diversity, iv) geographical latitude and v) tree dominance at a fine scale (in 0.25 ha plots), with the aim of obtaining a better understanding of the complex ecosystem processes and biological diversity [21]. Because of the strongly mixed nature of this tree community, which often produces low population densities of each tree species and random tree fall gaps caused by tree death, we expect aggregated patterns in individual P. chihuahuana trees and in the P. chihuahuana tree community [26], repulsive Picea patterns to other tree species and repulsive patterns of young to adult trees [41]. We also assumed no differences between i) the spatial distribution of the northern and southern populations of the P. chihuahuana tree community, and ii) the spatial distribution of suppressed and dominant trees, because of similar degrees of competition-induced mortality [42].

Material and Methods
We confirm that the field studies provide the specific location of study (Fig 1, S1 Dataset). No vertebrate studies were carried out. Field permit was granted by SEMARNAT, Mexico (http:// www.semarnat.gob.mx/).

Study area
Chihuahua spruce grows in areas characterized by an average temperature of between 9 and 12°C [1], precipitation ranging from 600 mm to 1,300 mm [43] and soil pH of 5.3-6.3 [4]. In order to determine the spatial tree structure of the P. chihuahuana tree community, 12 locations where the community occurs in the State of Durango and Chihuahua (north-western Mexico) were considered (Fig 1). Each location was represented by one plot of 50 x 50 m (0.25 ha) established in the centre of the tree community. Trees of all species of diameter at breast height (DBH) ! 7 cm were fully scored. The DBH, height and x, y coordinates were also recorded. The stem number per hectare (N), stand basal area (G), quadratic mean diameter (dg), mean breast height diameter (d), mean total height (h), maximum diameter (d max ), and maximum height (d max ) of all tree species (total) and Picea chihuahuana M. (Pch) were computed ( Table 1). The total numbers of tree species within each of the populations in the P. chihuahuana tree community have been reported by Quiñones-Pérez et al. [44].The DBH structures (as parameters of vertical structure) in the 12 plots considering all tree species showed a reverse J-shaped form (Fig 2) typical of uneven-aged forests. Fig 2 also demonstrates that the minimum balanced structure area of this tree community is very small (< 3ha) [45]. In total, 15 tree species were foundin the 12 plots: Abies durangensis Martínez, Cupressus lindleyi Klotzsch ex Endl., Juniperus deppeana Steud., P. chihuahuana, Pinus arizonica Engelm., Pinus strobiformis Engelm., Pinus cooperi Martínez, Pinus durangensis Martínez, Pinus leiophylla Schl. & Cham., Pinus teocote Schiede ex Schltdl. & Cham., Populus tremuloides Michx., Prunus serotina Ehrh., Pseudotsuga menziesii (Mirb.) Franco, Quercus sideroxyla Humb. andQuercus crassifolia Humb. In each plot, P. chihuahuana grew along with three to eight other tree species [44]. To represent the diversity profile (v sp,a ) of the tree species, we selected the described diversity for each location. Thus, each location of the P. chihuahuana tree community was characterized by the total number of tree species (species richness, (υ sp,0 )), effective number of tree species (Simpson index, (υ sp,2 )) and the number of prevalent tree species (υ sp,1 ), as Hill numbers [46] in each plot. The diversity values were taken from [13].

Spatial Structural Analysis
Clark-Evans index (CE), Uniform Angle index (W ) and Mean Directional index (R). The Clark-Evans index (CE) [31], the Uniform Angle index (W ) [36] and the Mean Directional index (R) [47] were used to describe the spatial distribution of the trees in each study plot, on the basis of the spatial distribution of the n trees nearest to a reference tree i. The CE was estimated using one neighbour (n = 1), while W and R were calculated using four neighbours (n = 4) [48], [49], [40]. A Poisson distribution pattern was characterised by a CE value of 1, cluster tendency by CE < 1 and a tendency of regular distribution of trees by CE > 1, with a maximum of 2.1491 for a hexagonal arrangement of trees.
For calculation of W , W i must first be calculated. The angle α 0 was set at 72 degrees, which yielded a mean value of W = 0.5 [37]. For each tree, the value of W i was determined and the average W for all trees was calculated. W i and W values close to 0 were associated with a regular neighbourhood of tree i, while values of W i and W close to 1 corresponded to irregularity of the spatial distribution in the neighbourhood of tree i.
Finally, calculation of R also requires calculation of R i . The exact value of R for a Poisson distribution in each plot in a 4-tree sample was 1.799, as obtained by a simulation based on 10 6 trees. Values of R i and plot mean R close to 0 were associated with a regular tendency of the neighbourhood of tree i, while values of R i and R larger than 1.799 were associated with a tendency of the spatial distribution in the neighbourhood of tree i to be irregular (see more in [40]).
To exclude the edge effect, and therefore to enhance the accuracy of the estimates, the nearest-neighbour edge-correction concept (NN1) was applied, as proposed by [29], for calculating CE, W and R. Table 1. Summary of important stand parameters calculated from the tree data: stem number per hectare (N), stand basal area (G), quadratic mean diameter (dg), mean breast height diameter (d), mean total height (h), maximum diameter (dmax), and maximum height (dmax) of the all tree species (total) and Picea chihuahuana M. (Pch) in the 50 x 50 m plotsin the 12 study locations and minimum (min), mean and maximum (max) parameter values for the stands.  The hypothesis of complete spatial randomness (CSR) for the mean values of CE, W and R for each plot was tested by atwo-sided permutation test (here 10,000 permutations) If 1-P(Z ! CE), P(Z ! W ), and P(Z ! R) are non-significantly small or non-significantly high (i.e. 0.01 < P < 0.99, at the 1% acceptance level), we can expect random effects and otherwise, directed forces. If the observed 1-P(Z ! CE), P(Z ! W ) or P(Z ! R) values are smaller than 0.01, we can assume non-randomly acting diversifying forces (e.g. seed dispersal pattern, association with nutrient-rich patches) that will produce a clustered distribution. If the observed 1-P(Z ! CE), P(Z ! W ) or P(Z !R) values are larger than 0.99, we assume that non-randomly acting homogenizing forces (e.g. competition for light, water and nutrients) will yield a regular distribution [50], [51], [40]. After Bonferroni correction [52], the new (modified) critical P value (significance level Ã = 0.0002) was calculated by dividing the critical P value (here the significance level = 0.05) by the number of comparisons (hypotheses) (m = 216). Spatial Structural Analysis by Ripley's K(t)-function. Ripley's K(t)-function is used to characterize completely mapped spatial point process data. The mapped data are usually recorded in two or three dimensions and include the locations of all events in a predefined study area. Unlike other functions (e.g. mean nearest-neighbour distance or the cumulative distribution function of distance from random points to their nearest neighbours), this function preserves information about distances between all points in the pattern, thus enabling visualization of how point pattern distributions vary with scale. Ripley's K(t)function is useful for summarizing point patterns, testing hypotheses about the patterns, estimating parameters and fitting models. Bivariate or multivariate K(t) functions can be used to describe relationships between two or more point patterns [53].
Ripley's K-function [33], [35] was used to determine the scales at which the tree pattern in each plot tends to be regular, clumped or random. The function was used to describe the relationship between the spatial pattern of Picea chihuahuana and the spatial structure of the other tree species inside the 12 plots.
The univariate Ripley's K-function [53] can be estimated as where A is the area of the study region, n is the number of observed points, w ij (r) is an edge effect correction factor, δ(r) is an indicator function and d ij is the distance between the i-th and j-th points. Because of its hyperbolic behavior, the interpretation of the K-function is not straightforward and a modification, called the L-function, was proposed by Besag (1977) in order to normalize the function (Besag, 1977): Now, the expected value of the univariate L-function under CSR is 0 for all r, positive when the pattern tends to be clustered and negative when the pattern tends to be regular.
In order to test the deviation from randomness of the point pattern using the univariate Lfunction, a 99% simulation envelope of L(r) was computed, using the Monte Carlo Method [54], from 1,000 simulated CSR patterns with the same number of points contained inside a region with the same geometry.
The bivariate Ripley's K-function [53] is estimated as where n i and n j are the numbers of the points of type i and j respectively, w ik,jl (r) is an edge effect correction factor and δ(r) is an indicator function and d ik,jl is the distance between the kth point of type i and the l-th point of type j. Due to the edge effect, K ij and K ji are correlated, but not identical, and therefore the following means of estimating K B is recommended: Its associated bivariate L-function is defined as The expected value of the bivariate L-function under spatial independence is 0 for all r, positive when the two point processes tend to be aggregated and negative when the two point processes tend to be repulsive.
In order to generate the simulation envelope that corresponds to the hypothesis of spatial independence, the method holds the point pattern of type 1 and type 2 constant and then randomizes their relative position in each simulation. For more details, see Lotwick and Silverman [55].
The trees were grouped in two diameter classes to check i) whether clustering at small scales was caused by a high degree of aggregation of smaller trees (using the univariate L-function) and ii) whether young trees were clustered around the adults (using the bivariate L-function). To find the d cut (cut-off) that divides the population into smaller and larger individuals (determined by their DBH), we obtained the bivariate L-function for various values within the range 7 < d cut,DBH < 40 cm and we therefore chose the d cut where the aggregation patterns between the groups were visually more significant. d cut,DBH,all corresponded to 23.2 cm, d cut,DBH,Pch to 29.3 cm. All analyses were performed using the "Spatstat" package implemented in the free statistical application R[56] [57].
The statistical tests for spatial tree pattern, null hypothesis, interpretation and related ecological questions are summarised in Table 2.

Covariation analysis
The relationships between stand densities (N and G), relative frequency of 10 cm DBH class (fcd) and tree species diversities (the Hill numbers υsp,0, υsp,2, andυsp,1 [46]), degree of latitude (lat) and spatial pattern indices (CE, W and R) were measured by the covariation (C) described by Gregorius et al. [50]. This method can detect types of covariation that are monotonous but not necessarily linear. C ranges between -1 and 1, where C = 1 indicates an entirely positive covariation and C = -1 a strictly negative covariation. If the denominator is zero, C is undefined [50]. Formally, In order to test the possibility that the observed degrees of covariation were only produced by random events rather than directed forces, a one-sided permutation test was performed (here 10,000 permutations) [58]. After Bonferroni correction, the new critical P value was 0.002.
In order to test whether the observed differences in the average spatial pattern indices (Diff) (CE, W and R) between i) P. chihuahuana trees and all other tree species, ii) suppressed and dominant P. chihuahuana trees and iii) suppressed and dominant trees of all species in the plots were produced solely by random events rather than directed forces, a permutation test based on 10,000 randomly chosen reassignments was performed [58]. Loosely based on the BAL competition Index [59] [60], the dominant tree class used in the present study included all larger trees that together included 50% of the stand basal area. The suppressed trees included all smaller trees that together included the other 50% of the stand basal area. This permutation test constitutes a non-parametric approach, which among other uses enables comparison of two groups in terms of the mean values of some variable; however, unlike with the t test, the assumptions of normality and equality of variances do not need to be satisfied by the data [51]. After Bonferroni correction of the data, the new critical P value was 0.007.

Results
Most of the P. chihuahuana trees were randomly distributed (92% of cases), as confirmed by all measures used (and considering a significance level of 1%): the Clark-Evans index (CE), uniform angle index (W ), mean directional index (R) and univariate L-function. Based on CE, W and R, 8% of the plots display clusteringat the 1% significance level. The univariate L-function indicated CSR in all plots, mostly due to the low effective tree number (repetitions) for calculating the values ( Table 3). The bivariate L-function showed that the number of Picea trees smaller than 29.3 cm DBH in the neighborhood of larger Picea trees (or equivalently the number of larger Picea trees in the neighbourhood of smaller Picea trees) was only larger than expected in TN and CV (Table 4). However, after Bonferroni correction of the data, all indices indicated random distribution of the P. chihuahuana trees in all plots.
Considering all tree species in each plot, the CE indicated CSR in 67% of the plots. The W and R indicate CSR in 92% of the plots at the 1% significance level (  Fig 3 (above and centre)). The bivariate L-function demonstrated that the number of smaller trees in the neighbourhood of larger trees (or equivalently the number of larger trees in the neighbourhood of smaller trees) was larger than expected (Table 4). However, after Bonferroni correction, the data indicated that the trees in 92% of the plots were randomly distributed (P ! 0.0002).
The Picea and other tree species were not spatially segregated, i.e. Picea tended to be found in patches of other tree species excepting the location San Jose de las Causas (SJ) ( Table 4, Fig 3  at the bottom left and right).
Analysis of the suppressed or immature trees versus the dominant or mature trees revealed that at the 1% significance level and according to CE, suppressed P. chihuahuana M. trees Table 2. Summary of statistical tests of spatial tree structure, null hypothesis, interpretation and related ecological questions; CE = neighbourhood-based Clark-Evans index, W = Uniform Angle index, R = Mean Directional index, LU(r) = univariate L-function, and LB(r) = bivariate Lfunction.  (Tables 6  and 7). The W and R values indicated CSR for all the plots where they were obtained. However, the mean W and R values for all plots of the suppressed trees (of all species) were statistically significantly higher than the dominant trees (of all species) (P = 0.0009 and P = 0.0006, respectively), also after Bonferroni correction (P < 0.007). Therefore, the tendency for clustering was significantly higher in the suppressed trees than in the dominant trees. W and R failed in some plots because of an insufficient number of trees(repetitions) for the calculations. were not statistically significant. The strongest covariation (C) between spatial pattern indices and tree density was C[R x N] with + 0.61 (P = 0.07) (i.e. clustering tended to be associated, but not significantly, with high stand density). The strongest covariation (C) between spatial pattern indices and tree species diversity was C[W x υ sp,2 ] with + 0.51 (P = 0.12) (i.e. clustering tended to be associated, but not significantly, with high tree species diversity). None of the 12 uneven-aged forest plots showed a statistically significant regular spatial tree pattern.
No statistically significant differences (Diff) (P < 0.01) between the average spatial pattern indices (CE, W and R) were observed for i) P. chihuahuana and all tree species or ii) suppressed and dominant P. chihuahuana trees.
When latitude was analyzed, we found that after Bonferroni correction, CE, but not W and R, decreased significantly from the southern to northern plots containing all tree species (C = -0.97, P = 0.0008). Table 3. Spatial structure of Picea chihuahuana M. in the 50 x 50 m plots in the 12 study locations, based on the neighbourhood-based Clark-Evans index (CE), Uniform Angle index (W), and Mean Directional index (R) (P values estimated with 10,000 permutations) and univariate L-function . The 99% simulation envelope (dashed red lines) for the CSR hypothesis was calculated by the Monte Carlo Method (Besag 1977), with 1,000 simulations (distance interval equals 0-12 m). N equals the number of Picea chihuahuana M. trees in the plot. 1

Discussion and Conclusions
In this study, we analysed the fine-scale spatial tree patterns in a special forest tree community of P. chihuahuana M. in Mexico. We examined the spatial tree pattern and its relationships to tree stand density, vertical structure, tree species diversity, geographical latitude and tree dominance. The findings demonstrate that the hypothesis of aggregated tree pattern is not applicable to the mean pattern measured by CE, W and R of the uneven-aged P. chihuahuana trees and P. chihuahuana tree community (Tables 3 and 5) and but to specific spatial scales measured by the univariate L-function, because 58% of the plots showed clustering at small (42%), intermediate and larger scales (17%) ( Table 4). Frequent clustering at small scales was mainly caused by a high degree of aggregation of trees smaller than 23.2 cm DBH (Table 4) [42]. We also found that the mean W and R values for immature (suppressed) were significantly larger than the corresponding values for mature (dominant) trees (Table 7). We also observed that immature P. chihuahuana trees showed clustering in 25% of the plots, according to the CE index. As in many primeval forests, small (young) individuals are almost always located in groups (Table 4)-often as a result of tree fall gaps in the canopy as the dominant type of disturbance [61] [62]. The common clustering at small scales in the P. chihuahuana tree community indicates that the forest patches were often created only by one canopy tree falling, as typically observed in species-rich tropical rain forests [63] and occasionally in old-growth (subalpine) Table 4. Analysis of spatial tree structure in 50 x 50 m plots in the twelve locations including all tree species using the univariate (for all trees, smaller trees of all species [44](< 23.2 cm diameter at the breast height (DBH)), larger trees (! 23.2 m DBH), and bivariate versions of the L-function (spatial pattern of Picea chihuahuana (Pch) vs. other tree species, pattern of smaller vs. larger trees and of smaller Pch (< 29.3 cm DBH) vs. larger Pchtrees (! 29.3 cm DBH)). The 99% simulation envelope (dashed red lines) for the CSR hypothesis (for the univariate L-function) and for spatial independence hypothesis (for the bivariate L-function) was calculated by the Monte Carlo Method (Besag 1977), with 1,000 simulations (distance interval equals 0-12 m). 1

Location
Univariate spruce-fir forests. These forests tend to be horizontally structured, mainly because an initiating disturbance is followed by long periods when small-scale, low-intensity disturbances control tree regeneration [8]. Therefore, fire disturbance, which almost always homogenizes stands, is a rare event in the P. chihuahuana tree community and should not be necessary for or beneficial to the community dynamics. Moreover, fires may bring an end to this fragmented tree community in very small and isolated locations [1] [17]. In contrast, the clustering by smallscale disturbances may be mainly caused by insect attack, disease or windthrow, which may create patchiness and spatial heterogeneity within locations [8].
The rare clustering at larger scales was mainly affected by the low tendency of aggregation of canopy trees (Table 4), as also reported by Malik et al. [64], Christensen [65] and Whipple [66] for uneven-aged populations. Therefore, the overall random patterns were a result of shift from initial aggregation to a random distribution [67]). Hence, Lepš and Kindlmann [42] postulated that i) an observed random pattern does not represent evidence of independence of individuals and ii) the intensity of spatial pattern should not be considered a measure of community organization.
Self-thinning due to intra and inter-specific competition-induced mortality was probably the main cause of the decrease in aggregation intensity [42], [67] during the course of population development in the P. chihuahuana tree community. However, environmental heterogeneity, uneven-age distributions, insufficient competition, limited seed dispersal and random germination may have prevented the presence of a significantly regular pattern of the mature trees in the tree community under study [68]. The aforementioned factors, particularly insufficient competition in the plot may also have favoured clustering in the northern locations (plots).
The number of P. chihuahuana trees in the neighbourhood of other tree species (or the number of trees from other species in the neighbourhood of P. chihuahuana) was not expected. The spatial distribution of P. chihuahuana trees was independent of the patches of other tree species, except in the San Jose de las Causas (SJ) location (Table 4, Fig 3). In SJ, spruces had a Table 5. Analysis of spatial tree structure in 50 x 50 m plots in the 12 locations including all tree species (species shown [44]) and based on the neighbourhood-based Clark  repulsive pattern to other species, similarly to a study in an old growth spruce-fir forest in Changbaishan Natural Reserve, China [41]. We therefore assume that there was a similar but weaker inter-and intraspecific competition between the trees at the species level [68] and that P. chihuahuana can tolerate partial shade conditions. The bivariate L-function showed that smaller trees (of all species) often grew in the neighbourhood of larger trees (of all species) ( Table 4), typically in uneven-aged forests (Fig 2) and probably as a result of some shade-  tolerant frequent tree species under mature canopy (such as Abies durangensis, Cupressus lindleyi and Juniperus deppeana) and the slender shaped crowns of the mature canopy trees in this community [69]. The smaller P. chihuahuana trees, did not, however, generally grow more frequently in the surroundings of larger Picea trees (Table 4), perhaps only for statistical reasons (such as few tree/repetitions) or due to no special preferences of the young (small) Picea trees for mature canopy Picea trees. In this study, P. chihuahuana regeneration was occasionally found on horizontal dead trees. The plots in which a clustered structure was observed tended to be associated (not significantly) with a large number of trees because of the presence of a greater number of understory trees. Understory trees often displayed a tendency to grouping (Table 4).   Although the aggregation indices were not associated with the diameter distribution, the results showed that none of the 12 uneven-aged forest plots under study displayed a statistically significant regular spatial tree pattern; however, 58% showed clustering in specific spatial scales at the 0.1% significant level (Table 4). No covariation (C) between aggregation indices and diameter distributions was observed, because the 12 diameter distributions scarcely varied in their reverse J-shaped form (Table 2) The cluster structure was weakly positively related to higher tree species diversity, probably due to a combination of the accumulation effect [13] [70] and increasing competition in denser plots [71]. While the accumulation effect resulted in higher diversity, the self-thinning processes led to saturation in tree species diversity [72].The high tree species diversity in the P. chihuahuana community [13] may also provoke clustering at smaller scales (in small gaps) Table 7. Analysis of spatial structure of the suppressed and dominant trees in 50 x 50 m plots containing all tree species(species shown [44]),in the 12 study locations, based on the neighbourhood-based Clark-Evans index (CE), Uniform Angle index (W), and Mean Directional index (R). P values estimated with 10,000 permutations. N equals the tree number in the plot. 1) La Tinaja (TN), 2) El Ranchito (RC), 3) El Cuervo (CV), 4) Talayote (TY), 5) Las Trojas (TR), 6) El Venado (VN), 7) La Quebrada (LQ), 8) Paraje Piedra Rayada (PPR), 9) Quebrada de los Durán (Arroyo del Indio Ignacio) (QD), 10) Cebollitas (CB), 11) San José de las Causas (SJ), and 12) La Pista (Arroyo de La Pista) (LP). because the lifespan and dimension of each tree species are often different. The probability that two or more trees of different species would fall at the same time and create a gap is lower than the probability of the same happening with trees of the same species. We conclude that satisfactory understanding of spatial forest structure is essential for the sustainable conservation of this unique mixed uneven-aged Picea forest [20]. Our measures of spatial tree structure, particularly W and R failed in several plots because of an insufficient number of trees (repetitions) for the calculations. Therefore, we recommend use of larger sample plot sizes (> 0.25 ha) in uneven-aged and species-rich forest ecosystems to detect less obvious, but important, relationships between spatial tree pattern and functioning and diversity in these forests.
Supporting Information S1 Dataset. Data set used in this study. (XLS)