The Post-Fire Assembly Processes of Tree Communities Based on Spatial Analysis of a Sierra Nevada Mixed-Conifer Forest

: Understanding the mechanisms underlying tree spatial arrangements may provide signiﬁcant insights into the processes in the maintenance of species coexistence. We examined the potential role of habitat heterogeneity, dispersal limitation, negative density dependence, ﬁre history, and unilateral intraspeciﬁc and interspeciﬁc interactions of adults on juveniles in shaping the spatial patterns of four dominant tree species ( Abies concolor , Pinus lambertiana , Calocedrus decurrens , and Quercus kelloggii ) after ﬁre in the Yosemite Forest Dynamic Plot, California, USA. We used the univariate pair correlation function and implemented three point pattern processes (homogeneous Poisson process, inhomogeneous Poisson process, and homogeneous Thomas process) to evaluate the potential contributions of habitat ﬁltering and dispersal limitation. We used a bivariate null model to evaluate unilateral intraspeciﬁc and interspeciﬁc interactions of adults on juveniles. We also used the pairwise correlation function to investigate the spatial patterns of density dependence. To understand the e ﬀ ect of ﬁre, we used the univariate pair correlation function to investigate pattern changes during the six years following ﬁre. We compared spatial pattern changes in both sprouting species ( Quercus kelloggii ) and seeding species ( Abies concolor ), and also examined the changes in patterns of large-diameter individuals of Abies concolor , Pinus lambertiana , and Calocedrus decurrens in 2013 (pre-ﬁre), 2016 (two years post-ﬁre), and 2019. Comparing the contributions of the homogeneous Thomas process and the inhomogeneous Poisson process at di ﬀ erent spatial scales showed the importance of dispersal limitation and habitat heterogeneity at ﬁner scales (0 m to 5 m) and coarser scales (5 m to 60 m), respectively, which suggests that the joint e ﬀ ects of dispersal limitation and habitat heterogeneity contribute to the spatial patterns of these three dominant tree species. Furthermore, the results showed that the young individuals of Abies concolor and Pinus lambertiana were more commonly found around the conspeciﬁc adults. Juvenile regeneration to the 1 cm diameter threshold was highly aggregated following the ﬁre. Large-diameter trees of Abies concolor , Pinus lambertiana , and Calocedrus decurrens generally did not exhibit patterns di ﬀ erent from complete spatial randomness ( Calocedrus decurrens ), or displayed only slight aggregation ( Abies concolor and Pinus lambertiana ). In addition, Abies concolor and Pinus lambertiana showed positive and negative conspeciﬁc density dependence in the immediate post-ﬁre period, respectively.


Introduction
Studying processes underlying tree species spatial distributions provides insights into mechanisms that maintain species coexistence. Niche-based processes theory assumes that each species has its

Study Area
The study was carried out in the Yosemite Forest Dynamics Plot (YFDP), a 25.6 ha (320 m × 800 m) plot established according to the Smithsonian ForestGEO protocols near Crane Flat in Yosemite National Park, central Sierra Nevada, California, USA ( Figure 1) [40,41]. The YFDP is further divided into 640 quadrats of 20 m × 20 m. The climate is Mediterranean, with hot, dry summers and cool, wet winters. Annual minimum and maximum monthly temperatures are in January (−1.82 • C to 10.01 • C) and July (13.9 • C to 27.14 • C) [42]. Mean annual precipitation is 1070 mm, with most precipitation occurring between November and March [42]. Elevation ranges from 1774.1 m to 1911.3 m. The principal tree species are Abies concolor (white fir), Pinus lambertiana (sugar pine), Cornus nuttallii (Pacific dogwood,) Calocedrus decurrens (incense-cedar), and Quercus kelloggii (California black oak) [23], with distribution and abundance of these species jointly controlled by climate and fire [43]. The YFDP was burned in August 2013 by a management-ignited backfire to prevent the spread of the Rim fire (details about fire behavior and fire effects: [27,[44][45][46]).

Field Methods
Within the YFDP, each tree was revisited annually between 2011 and 2019 and the status (live or dead) was checked each year, with tree diameters remeasured in 2014 and 2019. We analyzed the four most abundant tree species; Abies concolor, Pinus lambertiana, Calocedrus decurrens, and Quercus kelloggii (2862, 857, 450, and 1282 live stems, respectively, in 2019). To evaluate the potential influence of topographic variables on species spatial patterns, elevation, aspect, and slope were calculated metrics for each 20 m × 20 m quadrat using the surveyed elevations of the quadrat corners. Each quadrat was divided into four triangular planes (each triangle is formed by connecting three corners of the quadrat) and the slope was calculated by taking the mean angular deviation of these panels. Elevation was calculated by averaging the elevations of four corners of the quadrat [7]. As aspect is a land-surface variable, we used a cosine transformation to obtain a continuous gradient.

Null Models and Spatial Point Process Models
We used the pair correlation function to analyze the spatial pattern of individuals ≥1 cm diameter at breast height (dbh) [47,48]. The pair correlation function, g(r), replaces the circle in the Ripley's k-function with a ring to analyze spatial distributions [47,49]. The g(r) is the expected density of points falling within a circular ring of radius (r) around a focal tree divided by the intensity of the pattern [47]. If g(r)~1, then the inter-tree distances are consistent with complete spatial randomness. If g(r) values are larger than one, trees are more aggregated (clustered) than expected at radius r. Values of g(r) less than one indicate that trees are more dispersed than expected at distance r. Point process models use one of several distributions as null models [50]. We used inhomogeneous Poison processes, homogeneous Thomas processes, and three null models including complete spatial randomness, antecedent condition, and random labeling to address our questions (Table 1).

Null Models and Spatial Point Process Models
We used the pair correlation function to analyze the spatial pattern of individuals ≥1 cm diameter at breast height (dbh) [47,48]. The pair correlation function, g(r), replaces the circle in the Ripley's kfunction with a ring to analyze spatial distributions [47,49]. The g(r) is the expected density of points falling within a circular ring of radius (r) around a focal tree divided by the intensity of the pattern [47]. If g(r) ~1, then the inter-tree distances are consistent with complete spatial randomness. If g(r) values are larger than one, trees are more aggregated (clustered) than expected at radius r. Values of g(r) less than one indicate that trees are more dispersed than expected at distance r. Point process models use one of several distributions as null models [50]. We used inhomogeneous Poison processes, homogeneous Thomas processes, and three null models including complete spatial randomness, antecedent condition, and random labeling to address our questions (Table 1).

Complete Spatial Randomness (CSR)
The simplest and perhaps most common null model for univariate point patterns is complete spatial randomness (CSR), which indicates that the spatial distributions of trees are not influenced by any underlying biological mechanisms at any distance [52]. For CSR patterns, each point has the same probability of occurring at any place in the study area. Gradients due to first-order environmental differences are not considered.

Complete Spatial Randomness (CSR)
The simplest and perhaps most common null model for univariate point patterns is complete spatial randomness (CSR), which indicates that the spatial distributions of trees are not influenced by any underlying biological mechanisms at any distance [52]. For CSR patterns, each point has the same probability of occurring at any place in the study area. Gradients due to first-order environmental differences are not considered.  Table 1. Overview of null models, their results, and explanations used to investigate the spatial structure of dominant tree species in the Yosemite Forest Dynamic Plot.

Analysis Null Hypothesis Results Explanation
Complete Spatial randomness Trees are randomly distributed. Aggregation of trees ( Figure 2).
Tree aggregation may be caused by seed dispersal, environmental heterogeneity, interaction, fire, and density dependence.

Inhomogeneous Poisson Process
Trees are spatially aggregated due to habitat preference.
Aggregation of all trees at scales ≤15 m ( Figure 2).
Tree aggregations are based on habitat heterogeneity.

Homogeneous Thomas process
Trees are spatially aggregated at smaller distances due to seed dispersal limitations, especially in species with known restricted seed dispersal.
Spatial aggregation at scales < 5 m especially in species with restricted seed dispersal ( Figure 2).
Tree aggregations are based on species dispersal limitations.

Antecedent condition null model
Juveniles are spatially aggregated around conspecific adults.
Facilitation due to habitat suitability, resilience promotion, mycorrhizal fungi, and reduced competition.

Random labeling
Conspecific density-dependent mortality exists among species.
Density-dependent mortality due to strong intraspecific competition, lower allocation on defense mechanisms, increased susceptibility to bark beetles and drought stress.

Inhomogeneous Poisson Process (IPP)
The inhomogeneous Poisson process (IPP) considers first order gradients in the data. In this null model, the distribution of points is assumed to be based on an intensity function (λ) that depends on the interaction between environmental covariates and tree species density [53]. Aspect, slope, and elevation strongly affect the spatial distribution of solar radiation [54]. Solar radiation influences surface temperature, evatranspiration, and soil moisture content [55], which modify the local environment and influence species composition and productivity [56]. In this study, three topographic variables including aspect, slope, and elevation were selected as environmental variables. We used the heterogeneous Poisson process to estimate the effect of topographic heterogeneity on the local tree species density.

Homogeneous Thomas Process (HTP)
This point process is a Poisson cluster process which can model the impact of dispersal limitation (independent of the habitat heterogeneity) with the aggregation of the offspring around the parent trees. This process is modeled by two steps [48,57]. First, the locations of parent trees are generated from a homogeneous Poisson point process with λ > 0 on simulation window. Second, a number of offspring with mean µ > 0 are independently generated around each parent. These offspring form the clusters and were displaced from the parent by independent Gaussian dispersal kernel with standard deviation σ. For all analyses, 999 Monte Carlo simulations were used to evaluate departures of g(r) from the null models. We estimated 95% confidence envelopes from the twenty-fifth lowest and twenty-fifth highest values obtained from 999 simulations using the "spatstat" package (version 1.64-1; [58]) in R version 3.6.3 [59]. The distance associated with the Thomas process is the distance at which g(r)~1.

Contributions of Habitat Heterogeneity and Dispersal Limitation at Different Scales
The goodness-of-fit of the homogeneous Poisson process, inhomogeneous Poisson process, and homogeneous Thomas process were calculated by Akaike's information criterion (AIC). AICs were calculated using the "nmle" package (version 3.1-149; [60]). To evaluate the contributions of models at different spatial distances, distances were divided into the six categories: 0-2 m, 2-5 m, 5-10 m, 10-20 m, 20-40 m, and 40-60 m and the AIC values were calculated over each scale range. The best model was the one that had the lowest AIC among other models at a specific scale.

Antecedent Conditions Null Model
Point pattern analysis can be employed to analyze spatial distribution and unilateral interactions of juveniles with mature individuals [61]. This bivariate null model keeps the locations of antecedent (mature trees) pattern constant and randomizes the juvenile over the study area not occupied by mature trees [52]. Randomization of individuals was performed by random shifts which combine the rotation and mirroring method, while the shape of an individual occupied in different cells needs to be maintained. We used 999 Monte Carlo simulations and used the 25 largest and 25 smallest to define a 95% confidence interval. When the observed pattern falls above the 97.5% confidence interval, there is a positive spatial correlation between juveniles and adult or clustered distribution. When the observed pattern falls below the 2.5% confidence interval, adult trees negatively affect juvenile density and yield an overdispersed spatial distribution. Juveniles were considered to be those trees with 1 cm ≤ dbh < 5 cm and adults were considered to be those trees with dbh ≥ 20 cm in 2019 (diameter cutoffs followed selection criteria of [62]). Stems ≥ 5 and <20 cm were not included in the analysis because of the uncertainty in their reproductive status based on diameter threshold alone [62]. We performed a sensitivity analysis of ±50% on the diameter threshold values (2.5 cm to 7.5 cm dbh threshold values for juveniles and threshold values ≥ 10 cm to 30 cm dbh for adults) (Figures S1-S3).

Assessing Spatial Pattern
The spatial pattern was assessed using the univariate pair correlation function g(r). Changes in spatial patterns following the fire were examined for large-diameter trees and juvenile regeneration in 2013, 2016, and 2019 (different post-fire years). Large-diameter trees were considered to be those with dbh ≥ 60 cm. We limited our analysis to tree species with more than 70 stems since we needed a sufficient number of individuals to detect patterns in point pattern statistics [53]; those species were Abies concolor, Pinus lambertiana, and Calocedrus decurrens. The changes in juveniles spatial patterns were investigated in a sprouting species (Quercus kelloggii) and a non-sprouting species (Abies concolor). We applied a 9-m radius as the local neighborhood [23]. We used the goodness-of-fit test proposed by Loosmore and Ford [57,63], and the "ecespa" package [64] was used to examine the spatial pattern against the CSR. We set α = 0.05 and applied a Bonferroni correction for multiple tests (n = 23), resulting in a threshold p-value of 0.002.

Conspecific Negative Density Dependence
Conspecific negative density dependence (CNDD) was calculated with g d, d+l (r) − g l, d+l (r), which compares the conspecific neighborhood of trees that died and survived (subscripts show tree status: dead (d), live (l)). Dead trees were individuals that survived from fire but died in 2015-2019 and live trees indicate the trees that survived through the time period (2015-2019). Positive and negative values indicate that the neighborhood of trees that died was more (less) crowded than the neighborhood of trees that survived. Random labeling null model was implemented to generate 999 simulations. The goodness-of-fit test proposed by Loosmore and Ford [63] was used to calculate p-value at inter-tree distances up to 4 m for inference.

Edge Correction
Edge effects may arise in estimating point pattern due to the falling of sample rings partially outside of the study area. In this paper, we corrected edge effects by mirroring trees within 30 m of YFDP boundary to create a simulated stem map buffer zone around the study area [20].

Results
The observed (g) functions from CSR showed that species fell above the simulation envelopes at most scales, suggesting non-random species spatial distributions under the CSR null model ( distances up to 4 m for inference.

Edge Correction
Edge effects may arise in estimating point pattern due to the falling of sample rings partially outside of the study area. In this paper, we corrected edge effects by mirroring trees within 30 m of YFDP boundary to create a simulated stem map buffer zone around the study area [20].

Results
The observed (g) functions from CSR showed that species fell above the simulation envelopes at most scales, suggesting non-random species spatial distributions under the CSR null model (Figure  2A,D,G; Figure S4A; Figure S5A; Figure S6A).

Topographic Effect (Heterogeneous Poisson Process)
The observed clustered patterns of species could be generated by environmental gradients. To obtain an estimation of the magnitude of the habitat heterogeneity effects on the abundant species density, tree species distributions were generated based on topographic heterogeneity with the IPP null model. Abies concolor, Pinus lambertiana, and Calocedrus decurrens were aggregated at scales ≤ 15 m ( Figure 2B,D,H; Figure S4B; Figure S5B; Figure S6B). The percentage of aggregated species decreased with the increase in the spatial distance. In addition, the species densities were predicted based on the covariates data (topographic) in the plot ( Figure 3). The results show that the predicted patterns missed some clusters of trees. This indicates that point patterns analysis with topographic variables will be not enough to examine the species distributions.

Topographic Effect (Heterogeneous Poisson Process)
The observed clustered patterns of species could be generated by environmental gradients. To obtain an estimation of the magnitude of the habitat heterogeneity effects on the abundant species density, tree species distributions were generated based on topographic heterogeneity with the IPP null model. Abies concolor, Pinus lambertiana, and Calocedrus decurrens were aggregated at scales ≤15 m ( Figure 2B,D,H; Figures S4B-S6B). The percentage of aggregated species decreased with the increase in the spatial distance. In addition, the species densities were predicted based on the covariates data (topographic) in the plot ( Figure 3). The results show that the predicted patterns missed some clusters of trees. This indicates that point patterns analysis with topographic variables will be not enough to examine the species distributions.

Dispersal Limitation Effect (Homogeneous Thomas Process)
In this study, the effect of dispersal limitation on species distributions was tested separately to examine the importance of neutral-based (dispersal limitation) process in generating tree species spatial patterns. Results display that Abies concolor, Pinus lambertiana, and Calocedrus decurrens showed aggregated patterns at scales lower than 5 m, as with the inference of topographic effect ( Figure 2C,F,I; Figures S4C-S6C).

Contributions of Habitat Heterogeneity and Dispersal Limitation
Results indicated that spatial scale was likely to affect the suitability of dispersal limitation and habitat heterogeneity in explaining tree species spatial patterns ( Table 2). The homogeneous Thomas process was generally a better model in explaining species spatial patterns at finer scales. The inhomogeneous Poisson process performed better at coarser scales and the relative importance of the homogeneous Thomas process decreased with distance. At finer scales, individuals experienced relatively homogeneous environmental conditions and the species spatial distributions were severely limited by dispersal limitations. As the scales increased, the increased environmental heterogeneity appeared to be overbalanced by the lower effect of dispersal limitation. The results of the goodness-of-fit obtained from Monte Carlo simulation were relatively consistent with the results from AIC.

Dispersal Limitation Effect (Homogeneous Thomas Process)
In this study, the effect of dispersal limitation on species distributions was tested separately to examine the importance of neutral-based (dispersal limitation) process in generating tree species spatial patterns. Results display that Abies concolor, Pinus lambertiana, and Calocedrus decurrens showed aggregated patterns at scales lower than 5 m, as with the inference of topographic effect ( Figure  2C,F,I; Figure S4C; Figure S5C; Figure S6C).

Contributions of Habitat Heterogeneity and Dispersal Limitation
Results indicated that spatial scale was likely to affect the suitability of dispersal limitation and habitat heterogeneity in explaining tree species spatial patterns ( Table 2). The homogeneous Thomas process was generally a better model in explaining species spatial patterns at finer scales. The inhomogeneous Poisson process performed better at coarser scales and the relative importance of the

Spatial Patterns of Juvenile and Adult Trees
Juveniles of Abies concolor, Pinus lambertiana, and Quercus kelloggii were clustered around conspecific adults at 0-10 m, 0-7.5 m, and 0-4.5 m scales, respectively ( Figure 4A,C,E,G; see the location of the observed values in relation to simulation envelope in Figure S7A,C,E,G; see species totals in Table S1 and species spatial distribution in Figure S8A,C,E,G). Juveniles of Abies concolor and Quercus kelloggii exhibited a significant attraction with adults of other species at 0-6.5 m and 0-2.5 m scales respectively ( Figure 4B,D,F,H; see the location of the observed values in relation to simulation envelope in Figure S7B,D,F,H; see species spatial distribution in Figure S8B,D,F,H). The juveniles of Calocedrus decurrens were randomly distributed around conspecific and heterospecific adults at most distances.

Spatial Patterns of Juvenile and Adult Trees
Juveniles of Abies concolor, Pinus lambertiana, and Quercus kelloggii were clustered around conspecific adults at 0-10 m, 0-7.5 m, and 0-4.5 m scales, respectively ( Figure 4A,C,E,G; see the location of the observed values in relation to simulation envelope in Figure S7A,C,E,G; see species totals in Table S1 and species spatial distribution in Figure S8A,C,E,G). Juveniles of Abies concolor and Quercus kelloggii exhibited a significant attraction with adults of other species at 0-6.5 m and 0-2.5 m scales respectively ( Figure 4B,D,F,H; see the location of the observed values in relation to simulation envelope in Figure S7B,D,F,H; see species spatial distribution in Figure S8B,D,F,H). The juveniles of Calocedrus decurrens were randomly distributed around conspecific and heterospecific adults at most distances.

Overall Changes in Tree Spatial Patterns
The changes in the spatial pattern of regeneration (1 cm ≤ dbh < 5 cm) were investigated for A. concolor and Q. kelloggii in 2013 (pre-fire), 2016 (little post-fire), and 2019 ( Figure 5; Figure S9; Figure  S10; Table S2). Most of the post-fire recruitments (2014-2015) of these species were from small seedlings that survived the fire. The regeneration in A. concolor displayed a significant aggregation

Overall Changes in Tree Spatial Patterns
The changes in the spatial pattern of regeneration (1 cm ≤ dbh < 5 cm) were investigated for  Carlo goodness-of-fit test; 2016 and 2019: P = 0.001). The g(r) value for A. concolor in 2016 rose significantly at small scale, reached to the aggregation peak at about 8 m and decreased gradually at coarse scale. The A. concolor spatial pattern in 2019 was approximately the same as 2016, but displayed aggregation at greater scales. The spatial arrangement was different from complete spatial randomness from 0-9 m for Q. kelloggii in 2016 (Monte Carlo goodness-of-fit test; 2016: P = 0.001). The g(r) value decreased from 18 to 45 m and showed spatial inhibition at around 45 m.

Conspecific Negative Density Dependence
Diameter distributions ( Figure S14) showed that the plot was characterized by abundant trees in the smaller diameter classes in 2019. The number of living Quercus kelloggii was relatively greater than of dead Quercus kelloggii. Dead stems were in the smaller diameter classes ( Figure S10). The results (Figure 7) showed that live A. concolor had more conspecific neighbors (living and dead) at 1-

Conspecific Negative Density Dependence
Diameter distributions ( Figure S14) showed that the plot was characterized by abundant trees in the smaller diameter classes in 2019. The number of living Quercus kelloggii was relatively greater than of dead Quercus kelloggii. Dead stems were in the smaller diameter classes ( Figure S10). The results (Figure 7) showed that live A. concolor had more conspecific neighbors (living and dead) at 1-

Conspecific Negative Density Dependence
Diameter distributions ( Figure S14) showed that the plot was characterized by abundant trees in the smaller diameter classes in 2019. The number of living Quercus kelloggii was relatively greater than of dead Quercus kelloggii. Dead stems were in the smaller diameter classes ( Figure S10). The results (Figure 7) showed that live A. concolor had more conspecific neighbors (living and dead) at 1-4 m than did dead A. concolor which indicates the conspecific positive density dependence. In addition, Pinus lambertiana trees that died during the study had significantly more crowded neighborhoods than trees that survived. 4 m than did dead A. concolor which indicates the conspecific positive density dependence. In addition, Pinus lambertiana trees that died during the study had significantly more crowded neighborhoods than trees that survived.

Effect of Random Process, Habitat Heterogeneity, and Dispersal Limitation on the Formation of the Spatial Patterns of Abundant Species
The homogeneous Thomas process was the best model in describing the spatial pattern of three species compared to the inhomogeneous Poisson process at finer scales, consistent with a greater contribution of neutral processes and unmeasured environmental variables at those scales. Dispersal limitation introduced a mechanism in which a parent defines the location of an individual. Poor seed dispersal is related to a greater clumping pattern [65]. In species with animal dispersal mode, it is expected to have a higher dispersal ability and less aggregated pattern compared to species with wind or gravity dispersal mode. Seeds in P. lambertiana are not often dispersed great distances by wind, because they are large and heavy [66]. Furthermore, secondary mechanisms of dispersal (foodstoring rodents and birds) contribute substantially in P. lambertiana seed dispersion by dispersing seeds away from parent trees [17]. In addition, seeds of A. concolor disperse shorter distances due to their short and broad wing relative to their weight [18]. Calocedrus decurrens is a wind-dispersed species that releases seeds in late autumn when winds can be high, increasing dispersal distances. The better contribution of dispersal limitation in comparison with the environmental heterogeneity in shaping species spatial aggregation relative to environmental heterogeneity was consistent with several studies [67][68][69].

Effect of Random Process, Habitat Heterogeneity, and Dispersal Limitation on the Formation of the Spatial Patterns of Abundant Species
The homogeneous Thomas process was the best model in describing the spatial pattern of three species compared to the inhomogeneous Poisson process at finer scales, consistent with a greater contribution of neutral processes and unmeasured environmental variables at those scales. Dispersal limitation introduced a mechanism in which a parent defines the location of an individual. Poor seed dispersal is related to a greater clumping pattern [65]. In species with animal dispersal mode, it is expected to have a higher dispersal ability and less aggregated pattern compared to species with wind or gravity dispersal mode. Seeds in P. lambertiana are not often dispersed great distances by wind, because they are large and heavy [66]. Furthermore, secondary mechanisms of dispersal (food-storing rodents and birds) contribute substantially in P. lambertiana seed dispersion by dispersing seeds away from parent trees [17]. In addition, seeds of A. concolor disperse shorter distances due to their short and broad wing relative to their weight [18]. Calocedrus decurrens is a wind-dispersed species that releases seeds in late autumn when winds can be high, increasing dispersal distances. The better contribution of dispersal limitation in comparison with the environmental heterogeneity in shaping species spatial aggregation relative to environmental heterogeneity was consistent with several studies [67][68][69].
However, the results displayed the contribution of habitat heterogeneity in forming tree species spatial patterns at coarser scales (>5 m), indicating the joint effects of the dispersal limitation and habitat heterogeneity in forming the species spatial pattern [70]. Species clumping patterns at fine scales were determined by dispersal limitations, whereas the impacts of habitat heterogeneity occurred at larger scales. These findings indicate that the effects of dispersal limitation and environmental heterogeneity did not overlap at the same scales. Dispersal limitation first generates the template for tree distribution (seed aggregation close to the parent trees), but species habitat requirements may result in uneven seedling survival (if species survival depends on the habitat characteristics) and modify the sapling and adult distribution. Thus, the spatial pattern of seedlings and young trees may display a combination of both habitat heterogeneity and dispersal limitation. Therefore, it is better to evaluate the potential contribution of two processes simultaneously while examining spatial distribution.

Biotic Interactions of the Four Abundant Species
The analysis of juveniles around adults can be considered as an approach to examine the facilitation or microhabitat selection. According to this analysis, in three out of four species including A. concolor, P. lambertiana, Q. kelloggii, we found juveniles aggregated around conspecific adults. This also can be related to dispersal limitation in these species as the regeneration occurred near the seed source. California experienced a severe drought from 2012 to 2016 [36]. The aggregation may be attributed to microhabitats less affected by the drought [71] and higher juveniles resilience after the intense drought [37] around the conspecific adults. Positive interactions with mycorrhizal fungi may facilitate the establishment of conspecifics. The vast majority of tree species are colonized by either arbuscular (AM) or ectomycorrhizal (ECM) fungi in the most temperate forests. AM trees do not have fungal external structures that protect root cells from entering pathogens. In contrast, ECM trees have a fungal structure that prevents pathogens entry. A. concolor and P. lambertiana juveniles were clustered around the conspecific adults at <10 m and <7.5 m scales, respectively. These aggregated patterns can be related to the mycorrhizal association (both are ECM trees and have ectomycorrhizal fungi), which caused positive interactions among their juveniles with conspecific trees [72,73]. Furthermore, the aggregation of A. concolor juveniles around the conspecific adults at fine scales may have arisen from a lower mortality risk with proximity to conspecific neighbors in old-growth conifer forests [74]. The juveniles of Q. kelloggii exhibited attractions with conspecific adults at the short distances due to the quick regeneration by sprouting in a response to fire [75,76]. Additionally, the juveniles of A. concolor and Q. kelloggii were aggregated around adults of other species at small scales, respectively. These interspecific aggregations suggest that the intraspecific self-thinning process may reduce competition, promote interactions, and create appropriate conditions for heterospecific recruitment. Results showed that young individuals were more commonly found around conspecific adults in comparison to heterospecific adults (three species vs. two species). The observed random distributions of Calocedrus decurrens juveniles around adults indicate that juveniles did not display association with local adults' neighborhood. This suggests that the density-dependence effect was not striking to blur juveniles, and asymmetric competition between adults and juveniles represses the juveniles growth rather than eliminate them. The fact that more than 50% of these abundant species displayed intraspecific aggregation between juveniles and adults suggests the more significant effect of seed dispersal limitation rather than facilitation. The results of the antecedent condition null model would be better if we applied it to specific pairs of species, but due to the lack of enough sample size in particular species, we conducted it in all species in the plot.

Effect of Disturbance on the Spatial Pattern of Juveniles Regeneration and Large-Diameter Trees
Fire can alter plant community composition by generating post-fire conditions that favor in the establishment of some species and change the species composition in the remaining forest [43]. Fire was suppressed for a century in the YFDP which resulted in a nearly closed canopy forest, high stem density, and significant increase in surface and ladder fuels accumulations. The YFDP was burned in 2013 and the fire severity was classified as low-and-moderate severity fire regime with small unburned surface fuels and high severity patches within the YFDP [27,77]. The fire removed >90% of the fire surface fuels and >61% of the coarse woody debris [46,78]. Most post-fire mortality occurred in the first year, especially for small trees [79]. The proportional area of unburned fire refugia was 4.9% (12,597 m 2 ) with the lowest number in high-severity pixels (derived from Landsat) within the YFDP [27]. When lowand mixed-severity fires burn through sugar pine-white forest, they can create a mosaic with patches including open patches, patches with spaced single trees, and patches with aggregated trees in the forest [80][81][82][83]. Open patches of mineral soil are favorable sites for seedling establishment and increased soil moisture and reduced competition in these patches [84] may promote seedling and sapling survival. Furthermore, burned patches may show a moderate rise in the availability of resources including soil nutrients, light, and water increasing the chances of establishment and survival [70]. The fact that light was not a dominant driver of regeneration success (unsurprising in this post-fire forest) was confirmed by the decreased attraction between juveniles and adults at greater distances. Oak's traits and growth strategies allow it to take advantage of post-fire conditions. The better adaptation of oak can be traced to the ability to make root reserves by allocating carbohydrates to root growth [85], producing deep post-fire root systems from top-killed trees [86]. These characteristics allow oak to have a post-fire competitive advantage over seeders due to their shorter time to contribute to measured recruitment, i.e., 1 cm dbh [87]. Increasing the proportion of fire refugia may be considered as an important factor for post-fire species regeneration. As a nearby seed source is a principal requirement for A. concolor regeneration, there is likely to be higher juvenile regeneration near surviving mature trees in unburned patches. In addition, the residuals overstory in unburned patches produces shade which can limit shrub competition and promote seedling establishment. Bloomdahl et al. [27] reported that the proportion of fire refugia was higher in low-and moderate-severity burns and tree survival rates were higher in unburned patches for diameter classes ≤60 cm four years after the fire. In our study, most of the regeneration (three and six years following the fire) in A. concolor was from small seedlings that survived the fire because of their location in small fire refugia or completely unchanged patches (as inferred from Landsat). Results displayed the aggregated pattern in post-fire regeneration and increases in observed clustered patterns at finer scales following the fire. Additionally, other factors including favorable local microhabitat, vicinity to standing logs and stumps, which reduce summer drought by decreasing the solar radiation and increase soil moisture, could affect the regeneration's aggregated pattern [88,89]. All of the recruitment for Q. kelloggii was from sprouting and resulted in a strongly aggregated post-fire pattern of at scale about 1 m. Our results suggest that post-fire regeneration had a tendency to aggregate at distances up to 20 m, consistent with Fulé and Covington [90]. Given continued fire exclusion in the plot, it is likely that the shade-tolerant and fire-intolerant A. concolor may have the competitive advantage and will emerge as dominant conifer regeneration within the plot [91].
The degree of disturbance affects the spatial pattern of large-diameter trees in different ways [92,93]. Ground fuel accumulations and increased tree density can result in mortality for most sub-canopy trees while most of the larger trees survive [94]. Although low-to-moderate intensity fire often does not have a direct effect on large-diameter tree mortality it can result in injuries that reduce large-diameter tree physiological functionality [95]. Small patches of high-severity fire can kill most or all large-diameter trees, but lower severity patches leave even medium-sized trees alive [91]. Spatial patterns of A. concolor were clumped at scales of 0-9 m in 2013 and 2016. The survival of the large-diameter of A. concolor in refugia [4,90] would result in a clustering pattern following the fire. Random spatial distributions of large-diameter trees were observed in most scales following the fire. The spatial pattern of large-diameter P. lambertiana did not change much from 2013 to 2019, most likely due to its fire resistance. Clustering patterns at 0-9 m scales were observed in P. lambertiana and A. concolor following the fire occurrence. This can be related to the survived stems in scattered fire refugia [23,96] within the plot. The post-fire spatial pattern in C. decurrens was not different from CSR which can be related to the fire-resistant characteristic in large trees (fire-intolerant when they are small).

Effect of Conspecific Negative Density Dependence in Regulating Dominant Tree Species Spatial Pattern
The conspecific negative density dependence (CNDD) was examined to assess the possible effect of density-dependent mortality on species spatial patterns. The observed positive density dependence in A. concolor at <4 m scales may be associated with environmental heterogeneity within the YFDP. Some residual trees may be found in the fire-created patches which may also be good habitat for soil fungi and persistent soil nutrients. Furthermore, this positive density effect is possibly related to the fungal structure in the roots of some species where saplings clustered around conspecific trees experience higher performance due to shared ectomycorrhizal protection against pathogens [73]. Additionally, CNDD regulation was only detected at scales <4 m in Pinus lambertiana. This may be related to the lack of separation based on neighboring tree diameter size as the intensity and tendency of density dependence differ according to the tree diameter [20]. The increase in post-fire conspecific density mortality in Pinus lambertiana may be related to the reduction in trees allocation to resin defense (resin duct area) due to the higher intraspecific competition for water, which increases their vulnerability against bark beetles attack and drought stress [97].

Conclusions
Understanding processes underlying spatial arrangement reveals important clues for mechanisms of species coexistence. Investigating the effect of environmental heterogeneity, dispersal limitation, biotic interactions, and disturbance in spatial patterns development would help us to predict plant community responses to environmental changes as a result of global change. This study showed that both seed dispersal limitation and habitat heterogeneity likely shape the spatial distribution of abundant species. We found that fire can play a key role in forest composition by creating opportunities for tree species to be reorganized following fire, although a longer time period is needed to study the regeneration spatial pattern as most of the recruits analyzed here originated from small seedlings which survived the fire rather than recruits developing from seed post-fire. Additionally, a more comprehensive study of the effect of fire on the tree mortality spatial pattern is needed to examine fire's influence on species spatial pattern. The results of conspecific density dependence indicate the role of CNDD and CPDD in shaping abundant spatial distribution. Each of these processes may have acted individually or jointly with other processes to influence spatial patterns. Explaining the coexistence of species remains a challenge and future work could investigate the contribution of other ecological factors including bark beetles, pathogens, climate and drought, post-fire mortality pattern, physical mortality agent (including wind and crushing), and other environmental heterogeneity (including light, soil moisture, temperature) in species spatial pattern.
Supplementary Materials: The following are available online at http://www.mdpi.com/2571-6255/3/4/72/s1, Figure S1. Distribution and density of the four most abundant species in the Yosemite Forest Dynamic Plot in 2019; Figure S2. Results of the sensitivity analysis for diameter cutoff values chosen for the grouping of adult and juvenile stems in 2019 in the Yosemite Forest Dynamic Plot; Figure S3. Results of the sensitivity analysis for diameter cutoff values were chosen for grouping the adult and juvenile stems in 2019 in the Yosemite Forest Dynamic Plot; Figure S4. Panels display the relationship between the univariate pair-correlation function g(r) and distances for Abies concolor in 2019 in the Yosemite Forest Dynamic Plot; Figure S5. Panels display the relationship between the univariate pair-correlation function g(r) and distances for Calocedrus decurrens in 2019 in the Yosemite Forest Dynamic Plot; Figure S6. Panels display the relationship between the univariate pair-correlation function g(r) and distances for Pinus lambertiana in 2019 in the Yosemite Forest Dynamic Plot; Figure S7. Left panels show the relationship between the pair-correlation function g(r) and scales for the Abies concolor (A), Calocedrus decurrens (C), Pinus lambertiana (E), Quercus kelloggii (G), to assess the interaction between juveniles (1 cm ≤ dbh < 5 cm dbh) and conspecific adults (individuals ≥ 20 cm dbh) in 2019 in the Yosemite Forest Dynamic Plot; Figure S8. Left panels display the distributions of juveniles (1 cm ≤ dbh < 5 cm dbh) and conspecific adults (dbh ≥ 20 cm) in Abies concolor (A), Calocedrus decurrens (C), Pinus lambertiana (E), and Quercus kelloggii (F), in 2019 in the Yosemite Forest Dynamic Plot; Figure S9. Comparison of juveniles regeneration spatial patterns (1 cm ≤ dbh < 5 cm) in seeder species (Abies concolor) and sprouting species (Quercus kelloggii) in 2013 (pre-fire), 2016 (post-fire), and 2019 (postfire) in the Yosemite Forest Dynamic Plot. Black lines display the observed g(r), values above (below) simulation envelopes indicate aggregated (dispersed) pattern; Figure S10