Soil Enzyme Activity and Soil Nutrients Jointly Inﬂuence Post-Fire Habitat Models in Mixed-Conifer Forests of Yosemite National Park, USA

: Disentangling the relative importance of habitat ﬁltering and dispersal limitations at local scales ( < 1 km 2 ) in shaping species composition remains an important question in community ecology. Previous studies have examined the relative importance of these mechanisms using topography and selected soil properties. We examined both topography and edaphic properties from 160 locations in the recently burned 25.6 ha Yosemite Forest Dynamics Plot (YFDP) in Yosemite National Park, California, USA. In addition to eight soil chemical properties, we included phosphatases and urease enzymes in a deﬁnition of habitat niches, primarily because of their rapid changes with ﬁre (compared to soil nutrients) and also their role in ecosystem function. We applied environmental variables to the distributions of 11 species. More species–habitat associations were deﬁned by soil properties (54.5%) than topographically-deﬁned habitat (45.4%). We also examined the relative importance of spatial and environmental factors in species assemblage. Proportions explained by spatial and environmental factors di ﬀ ered among species and demographic metrics (stem abundance, basal area increment, mortality, and recruitment). Spatial factors explained more variation than environmental factors in stem abundance, mortality, and recruitment. The contributions of urease and acid phosphatase to habitat deﬁnition were signiﬁcant for species abundance and basal area increment. These results emphasize that a more complete understanding of niche parameters is needed beyond simple topographic factors to explain species habitat preference. The stronger contribution of spatial factors suggests that dispersal limitation and unmeasured environmental variables have high explanatory power for species assemblage in this coniferous forest.


Introduction
Niche-based processes and neutral processes both shape species assemblages in ecological communities. Niche theory assumes that different species have their own niche and species adaptation to specific environmental heterogeneity and biotic interactions determine the species coexistence. In contrast, neutral theory emphasizes the role of stochastic events such as dispersal-assembly in shaping community structure. Neutral theory assumes that all individuals of all species are ecologically equivalent and the environmental variables play no role in community structure [1]. Dispersal limitation controls the spread of individuals into various habitats, while habitat filtering helps multi-species remain in coexistence through interspecific competition for the same restricted resources [2].
Understanding the relationships between species demographic metrics (stem density, basal area increment, mortality, and recruitment) and habitats in forests is important not only in shaping The five most abundant species are Abies concolor (white fir), Pinus lambertiana (sugar pine), Cornus nuttallii (Pacific dogwood,) Calocedrus decurrens (incense-cedar), and Quercus kelloggii (California black oak) ( Table 1, Supplementary material, Figure S1). The YFDP is situated on two soil polygons of the Clarks Lodge-Ultic Palexeralfs complex and the Typic Dystroxerepts-Humic Dystroxerepts complex [29].   Each stem was revisited annually between 2011 and 2019, and the status (live or dead) was checked each year, with diameters remeasured in 2014 and 2019. Unburned patches ≥1 m 2 (unburned litter and duff layer) were mapped at the beginning of the growing season immediately after the fire [34]. Topographic variables (elevation, aspect, and slope) of each 20 × 20 m quadrat were calculated based on the surveyed position and elevation of the 20-m grid reference corners. Elevation was taken as the average of elevation of four corners of each quadrat, and slope was measured as the mean angle of the four panels by connecting three corners of a quadrat. Aspects between 135 • and 225 • were considered south facing, because they receive the most direct solar exposure [39]. Aspect >225 • and <135 • were considered as one group due to the lower amount of sun radiation and temperature. As aspect is a land-surface variable, we used a cosine transformation to obtain a continuous gradient, describing the north-south gradient.
Cumulative infiltration and hydraulic conductivity were calculated using mini disk infiltrometer in 56 burned and 39 unburned sites. The infiltrometer was placed on the soil surface and the water was pulled from the tube by soil suction. The volume of water was recorded at 30 s intervals and plotted (cumulative infiltration versus the square root of time) according to the methods of Zhang [40] where C 1 is the slope for the cumulative infiltration vs. the square root of time, and A is a value that relates the van Genuchten parameters for a given soil texture class to both disk radius and the suction we selected. A is computed from the below formula: where r is the disk radius, h is the suction at the disk surface, n and α are the van Genuchten parameters for the soil. The van Genuchten parameters for the 12 texture classes were obtained from Carsel and Parrish [41] (Table S1). Soil samples were collected at 160 points (98 samples from burned sites and 62 samples from unburned patches) within the YFDP in May 2017. Samples were air dried at temperature (22 • C) and sieved to remove stones (with < 2 mm sieve). The BaCl 2 method was used to determine the concentration of Ca (calcium), K (potassium), Mg (magnesium), and Mn (manganese). The Bray method was used to measure the concentration of P (phosphorus). Soil samples were extracted in 0.1 M BaCl 2 for two hours, and the concentration of Ca, K, Mg, and Mn were determined by Inductivity Coupled Plasma Analyzer [42]. Effective cation exchange capacity (ECEC) was calculated as the sum of the exchangeable cations, which are mostly Ca, Na (sodium), K, and Mg. Cation exchange capacity (CEC) was calculated as a total quantity of negative surface charges. Total exchangeable bases (TEB) was obtained from summation of exchangeable K, Ca, Mg, and Na. Base saturation (BS) was calculated by dividing TEB by CEC value and multiplying by 100. Soil samples were collected at the same locations (160 quadrats; 98 burned patches and 62 unburned patches) for measuring the alkaline phosphatase, acid phosphatase, and urease activity in 2018. We collected three soil samples per quadrat and mixed them thoroughly. The mixed samples were considered as the representative of a sample for each quadrat. Samples were sieved from quadrats and maintained at < 5 • C during transport to the lab. We allowed them to equilibrate at room temperature before starting enzymes measurements. Enzyme activity analysis was conducted using the methods developed by Dick [43]. Urease activity was assayed according to the methods of Kandeler and Gerber [44]. We used 2.5 milliliters (ml) of urea solution and 20 mL borate buffer containing disodium tetraborate for each 5 g soil sample and incubated them at 37 • C for two hours. A 30 mL potassium chloride (2 M)-hydrochloric acid (0.01 M) solution was added, and the mixtures were shaken on a shaker for 30 min. Soil suspensions were filtered and filtrates analyzed for ammonium by colorimetric procedure. Phosphatases (acid and alkaline phosphatases) were measured by the method of Tabatabai and Bremner [45,46], which includes colorimetric estimation of p-nitrophenol release (acid solution of the p-nitrophenol is colorless and the alkaline solution has yellow color) when 1 g of soil is incubated with 0.2 mL toluene and 4 mL of buffered sodium p-nitrophenyl phosphate solution (pH for buffer were considered equal to 6.5 for acid phosphatase and 11 for alkaline phosphatase) at 37 • C for 1 h. After incubation, CaCl 2 -NaOH treatment was used to extract the p-nitrophenol released by phosphatase activity.

Habitat Definition
We identified two classes of habitat predictors (topographic and soil variables) to define habitat maps. Topographic variables were comprised of elevation, aspect, and slope. Soil variables were Ca, K, Mg, Mn, total exchangeable bases (TEB), base saturation (BS), P, pH, and soil enzymes, including acid and alkaline phosphatases and urease. We calculated topographic variables (elevation, aspect, and slope) at the 1 × 1 m and 20 × 20 m scales ( Figure S2 and Figure 2) within the YFDP. The optimal number of habitats was determined by elbow and gap statistic methods using the fviz_nbclust function from factoextra package version 1.0.3 [47]. In the elbow method, a K-means clustering algorithm was run on the data set and the total within-cluster sum of square (WSS) was calculated. By plotting the WSS curve and number of clusters, the point of inflection on the curve was chosen as the optimal number of clusters. We verified the appropriate number of clusters using complementary methods (gap statistic and NbClust function). The hierarchical clustering was used to classify each quadrat within a plot into a habitat based on the environmental variables. Selective cuts across dendrogram were made to generate habitats based on the optimal number of habitats, which were determined by previous step. All analyses were performed in R version 3.4.3 [48].  We performed a species-habitat association test (torus translation) on species with ≥25 stems (stem density ≥1 stem/ha) (eleven species) ( Table 2). This threshold for local abundance was applied to differentiate rare from abundant species [39,49]. The associations of stem abundance in 2019, basal area increment from 2014 to 2019, mortality from 2014 to 2019, and recruitment from 2014 to 2019 in these eleven species were assessed within 160 quadrats (20 × 20 m). The torus translation test was conducted by following the methods of Harms et al. [50]. This test calculates the observed abundance  We performed a species-habitat association test (torus translation) on species with ≥25 stems (stem density ≥1 stem/ha) (eleven species) ( Table 2). This threshold for local abundance was applied to differentiate rare from abundant species [39,49]. The associations of stem abundance in 2019, basal area increment from 2014 to 2019, mortality from 2014 to 2019, and recruitment from 2014 to 2019 in these eleven species were assessed within 160 quadrats (20 × 20 m). The torus translation test was conducted by following the methods of Harms et al. [50]. This test calculates the observed abundance of each species in each habitat type and compares these observed values with abundance values obtained from simulated habitat maps. Simulated maps were generated by shifting the actual habitat map in four directions by 20-m increments, while the location of the stems did not change. A species was significantly positively (aggregated) or negatively (repelled) with a specific habitat type at (α = 0.05) if observed abundance was higher (lower) than at least 97.5% (or 2.5%) of the simulated abundance in simulated maps ( Figure S3).

Principal Coordinates of Neighbor Matrices
Principal coordinates of neighbor matrices (PCNM) proposed by Bocard and Legendre [51] were used to model spatial variation. Generation of spatial variables was conducted using the pcnm function from the "vegan" package version 2.5-6 [52]. The distance between spatial data was represented as a Euclidean distance matrix. This method creates a set of spatial explanatory variables and determines significant variables based on the statistical responding of the response variable [53]. Data was normalized using the Hellinger transformation before PCNM analysis. The PCNM function provides negative and positive eigenvalues as predictors, but only positive eigenvalues were selected as explanatory variables.
The number of variables was reduced by selecting variables with a statistically significant contribution on variation of species abundance (α = 0.05) using forward selection with the ordistep function (999 permutations) [54]. The variation partitioning was conducted using the varpart function from the "vegan" package [52] to partition the explained proportions of variation in species composition by environmental and spatial variables. The significance of each component was tested using anova and rda functions from the "vegan" package.

Results
Most soil properties were not significantly correlated with topography (Supplementary material, Figure S4, Table S2). Hydraulic conductivity was slightly greater in unburned sites, but the difference between burned and unburned sites was not significant five years after fire ( Figure 3).
Differences in enzyme activity (urease, phosphatase, and alkaline phosphatase) between burned and unburned sites were only significant for urease five years after fire (p-value < 0.05) ( Figure 4).
Hydraulic conductivity and alkaline phosphatase were added to our soil data as predictors, which resulted in a lower explained proportion of edaphic component in species demographic metrics, compared to those with consideration of two enzymes (acid phosphatase and urease) (Supplementary material, Figures S5 and S6 and Figure 6). The number of habitats as identified by the combination of the elbow method (Supplementary material, Figure S7), gap statistic, and the diagnostics of the NbClust package resulted in four and seven habitats based on the topographic (slope, elevation, and aspect) and eleven soil variables (eight soil chemical properties plus three soil enzyme activities) ( Figure 5, Supplementary material, Figure S8, Table S3). using anova and rda functions from the "vegan" package.

Results
Most soil properties were not significantly correlated with topography (Supplementary material, Figure S4; Table S2). Hydraulic conductivity was slightly greater in unburned sites, but the difference between burned and unburned sites was not significant five years after fire ( Figure 3). Differences in enzyme activity (urease, phosphatase, and alkaline phosphatase) between burned and unburned sites were only significant for urease five years after fire (p-value < 0.05) (Figure 4).  Most soil properties were not significantly correlated with topography (Supplementary material, Figure S4; Table S2). Hydraulic conductivity was slightly greater in unburned sites, but the difference between burned and unburned sites was not significant five years after fire ( Figure 3). Differences in enzyme activity (urease, phosphatase, and alkaline phosphatase) between burned and unburned sites were only significant for urease five years after fire (p-value < 0.05) ( Figure 4).   metrics, compared to those with consideration of two enzymes (acid phosphatase and urease) (Supplementary material, Figures S5, S6 and 6). The number of habitats as identified by the combination of the elbow method (Supplementary material, Figure S7), gap statistic, and the diagnostics of the NbClust package resulted in four and seven habitats based on the topographic (slope, elevation, and aspect) and eleven soil variables (eight soil chemical properties plus three soil enzyme activities) ( Figure 5, Supplementary material, Figure S8, Table S3). Every other quadrat was assigned to a specific habitat, and the unassigned quadrats were removed from the analysis. "HS" and "LS" indicate high and low slope in habitats. "North" and "south" show north or south facing habitats.
Among the eleven species, stem abundance of five species in 2019 (45.5% of stems) were negatively or positively associated with habitats ( Table 2). The number of significantly associated species in habitats defined by soil variables was slightly greater compared to total number of species associated with habitatsdefined by topographic factors alone (6 versus 5). The total number of demographic metrics (basal area increment, mortality, and recruitment) of species associated with habitats were smaller than number of species abundance associated with habitats (one (9.1%), two (18.2%), and two (18.2%), respectively). Among the eleven species, stem abundance of five species in 2019 (45.5% of stems) were negatively or positively associated with habitats ( Table 2). The number of significantly associated species in habitats defined by soil variables was slightly greater compared to total number of species associated with habitatsdefined by topographic factors alone (6 versus 5). The total number of demographic metrics (basal area increment, mortality, and recruitment) of species associated with habitats were smaller than number of species abundance associated with habitats (one (9.1%), two (18.2%), and two (18.2%), respectively).  Only 27 PCNMs were selected to predict the variation in community composition. The adjusted cumulative square for all 27 PCNMs was 27.9% (Supplementary material, Table S4). The proportion of variance explained by spatial and environmental variables with and without soil enzymes as a predictor for stem abundance was 45% as opposed to 41%; for species basal area the increase was 10% vs. 7%; for species mortality 53% vs. 52%; and for species recruitment, 52% vs. 51%, respectively ( Figure 6).  Table S4). The proportion of variance explained by spatial and environmental variables with and without soil enzymes as a predictor for stem abundance was 45% as opposed to 41%; for species basal area the increase was 10% vs. 7%; for species mortality 53% vs. 52%; and for species recruitment, 52% vs. 51%, respectively ( Figure 6). The variation explained by spatial variables alone was greater compared to other variables for stem abundance, mortality, and recruitment in the YFDP (Figure 6). The contributions of only the topographic component in species abundance, basal area increment, and mortality were decreased  The variation explained by spatial variables alone was greater compared to other variables for stem abundance, mortality, and recruitment in the YFDP ( Figure 6). The contributions of only the topographic component in species abundance, basal area increment, and mortality were decreased by removing soil enzymes data from edaphic predictors. Soil variables explained more variation than topographic variables in species abundance.

Associations of Different Species with Habitat Types
About half of the species were positively (six species) or negatively (seven species) associated with specific habitats. Species that are positively associated with a specific habitat may be more competitive than the species that are negatively repelled or neutrally (no association with respect to habitat) associated with the same habitat [55]. Five species were associated with habitats defined by topographic variables. Slope is an important factor, likely due to its effect on water availability, especially during the dry seasons [50]. Aspect often plays a role in species composition [56], by influencing water potential, organic matter, irradiance availability at ground level, and the creation of different microclimates [57]. Generally, low-slope north-facing sites experienced cooler temperature, a lower solar radiation and evapotranspiration rate due to the lower exposure of sunlight, greater runoff water accumulation due to the deep soil [58], and a greater amount of organic matter. Abies concolor grows in the environment with heterogenous soil conditions and shows the best growth on a moderate slopes and level ground [59]. The abundance of Abies concolor showed positive association with the low slope. Consistent with those results, mortality of Abies concolor was negatively associated with north-facing, low slopes (observed mortality number from habitat map was <2.5% of the simulated mortality value from torus-translation). The importance of water availability as a restricting factor in Abies concolor development was also found by Laacke [59].
Recruitment of Cornus sericea was positively associated with habitat 1. The levels of P concentration and K were high in these habitats. However, this positive association may be related to other factors including the high soil moisture in this habitat and the proximity to high abundances of parent plants at moist sites (considerable reproduction for this species is vegetative). Quercus kelloggii mortality was positively associated with habitat 6, where phosphorus, calcium, and urease enzyme levels were high. This association could be created as a result of higher competition in habitats with greater nutrient sources, which could result in a greater number of observed mortalities. Basal area increment of Quercus kelloggii was positively associated with habitat 7, where phosphatase enzyme activity, Ca, K, and Mg were all high. Additionally, Quercus kelloggii basal area increment was negatively associated with habitat 5 where Ca, Mg, and phosphatase levels were the lowest among all habitats and P concentration was not high. Neba et al. [60] found that the addition of Mg resulted in a better height and diameter growth due to a better root growth and greater nutrient uptake from the soil. The important effect of P in dry matter production and basal area increment was also found by another study [61]. Increase in tree growth with the availability of Ca was presented by Baribault et al. [62]. In addition, a significant effect of Mg on stem diameter growth at breast height by increasing nutrient uptake was confirmed by other studies [63].
The habitat map created by edaphic variables produced a more heterogeneous pattern than a habitat map generated by topographic variables in this study ( Figure 5). The result was a greater number of species associated with edaphically-defined habitats in comparison with the number of species associated with topographically-defined habitats. The greater number of species associated with habitats in a more complex habitat map (heterogeneous pattern) was supported by Borcard and Legendre [51].

Niche vs. Dispersal Limitation Drive Variations in Species Abundance and Basal Area Increment
The role of niche and dispersal limitation in shaping forest communities within the YFDP was investigated by partitioning the variation in species demographic metrics into different portions determined by edaphic, topographic, and spatial variables. The variance explained by purely spatial variables was attributed to dispersal-assembly and responses of species to the unmeasured environmental variation [64]. Although in general, variance partitioning analyses with observational data cannot distinguish unmeasured environmental variables and neutral processes [65], this analysis included a more comprehensive environmental dataset than that used by Legendre et al. [65], which considered topography as the principal environmental factor. We thus decreased the effect of unmeasured environmental variables in the pure spatial fraction. However, other unmeasured environmental variables (such as light availability, soil temperature, soil moisture, and competition in the local tree neighborhood) undoubtedly contribute to the community structure. Dispersal limitation has a strong potential to structure communities at fine scales, especially in species with a lower dispersal ability, whose seeds are dispersed close to their parents [66][67][68]. Spatial, topography, and soil resources (soil properties with and without enzymes) were all statistically significant in their contribution to species abundance (P = 0.01 and P = 0.03, respectively) and basal area increment (P = 0.02 and P = 0.03, respectively). Results showed that a large contribution (more than 30%) of total variation of species abundances was explained by spatial variables. The important effects of biotic processes such as dispersal, stochasticity process such as demographic stochasticity, and the weak effects of habitat filtering in structuring species composition at small scale (10 m to 20 m) were presented by Méndez-Toribio et al., Yuan et al.,and Shipley et al. [58,69,70].
The incremental contribution of spatial and environmental factors varied among species (Tables  S5 and S6). The variation in stem abundance explained by spatial variables was the highest for Pinus lambertiana, which has heavy seeds with small wings that could result in a shorter primary dispersal distances [71] (Figure 7). Abies concolor grows on a variety of soil conditions and pH levels. In addition to fire history, their abundance mostly depends on water availability and temperature [59], supporting the high contribution of topographic variables in explaining variation in Abies concolor stem abundance (Figure 7). included a more comprehensive environmental dataset than that used by Legendre et al. [65], which considered topography as the principal environmental factor. We thus decreased the effect of unmeasured environmental variables in the pure spatial fraction. However, other unmeasured environmental variables (such as light availability, soil temperature, soil moisture, and competition in the local tree neighborhood) undoubtedly contribute to the community structure. Dispersal limitation has a strong potential to structure communities at fine scales, especially in species with a lower dispersal ability, whose seeds are dispersed close to their parents [66][67][68]. Spatial, topography, and soil resources (soil properties with and without enzymes) were all statistically significant in their contribution to species abundance (P = 0.01 and P = 0.03, respectively) and basal area increment (P = 0.02 and P = 0.03, respectively). Results showed that a large contribution (more than 30%) of total variation of species abundances was explained by spatial variables. The important effects of biotic processes such as dispersal, stochasticity process such as demographic stochasticity, and the weak effects of habitat filtering in structuring species composition at small scale (10 m to 20 m) were presented by Méndez-Toribio et al., Yuan et al., and Shipley et al. [58,69,70]. The incremental contribution of spatial and environmental factors varied among species (Tables  S5 and S6). The variation in stem abundance explained by spatial variables was the highest for Pinus lambertiana, which has heavy seeds with small wings that could result in a shorter primary dispersal distances [71] (Figure 7). Abies concolor grows on a variety of soil conditions and pH levels. In addition to fire history, their abundance mostly depends on water availability and temperature [59], supporting the high contribution of topographic variables in explaining variation in Abies concolor stem abundance (Figure 7).

The Contribution of Environmental and Spatial Variables in Forming Mortality Abundances Across Species
Spatial and topographic variables were significant (P = 0.02) in their contribution to species mortality and not significant considering the effect of soil factors (soil properties with and without soil enzymes). The higher contribution of the spatial variables in explaining the variation of species mortality may be related to strong neighborhood competition in species with limited dispersal ability due to a higher density of small individuals near the parent tree [72]. As opposed to recruitment,

The Contribution of Environmental and Spatial Variables in Forming Mortality Abundances Across Species
Spatial and topographic variables were significant (P = 0.02) in their contribution to species mortality and not significant considering the effect of soil factors (soil properties with and without soil enzymes). The higher contribution of the spatial variables in explaining the variation of species mortality may be related to strong neighborhood competition in species with limited dispersal ability due to a higher density of small individuals near the parent tree [72]. As opposed to recruitment, mortality in old-growth forests is often due to insects; physical damage by wind, snow, other falling trees; disease; and intense neighborhood competition [73]. Furniss et al. [22] found that mortality following the fire was differentiated based on diameter class, and that large-diameter trees had higher survival rates than small-diameter trees. The changes in variation of species mortality explained by inclusion of soil enzymes into edaphic factors was marginal (1%). The negligible proportion of soil variables in explaining mortality indicates that soil variables are not differentiating factors for mortality in old-growth forests.
The variation in mortality explained by environmental and spatial components varied with species (Table S7). This could be related to soil nutrient availability [74,75]. The contribution of topographic variables was the highest for Cornus nuttallii, indicating the hydrological variations related to topography.

The Contribution of Environmental and Spatial Variables in Forming Ingrowth Numbers Across Species
Spatial and topographic variables were significant (P = 0.01) contributors to recruitment and not significant when considering soil factors (soil properties with and without soil enzymes) alone. The fraction of the spatial component in explaining variation of species recruitment was the highest among the other variables ( Figure 6). This showed the principal role of seed availability (or vegetative propagation) in recruitment at a local scale [76]. The low contribution of environmental heterogeneity to recruitment may be related to the importance of other factors, such as fecundity, germination rates, and initial growth rates of large-seeded species [77,78]. It is likely that other soil properties including temperature, especially during the January to March, affect the survival rate of seedlings due to the susceptibility of young seedlings to low temperature [79]. In addition, other factors include litter layer depth, which may prevent seedling emergences in small-seeded species [79].
The contribution of environmental and spatial components in explaining recruitment changed with species (Table S8). The proportion of environmental variables was the lowest for Chrysolepis sempervirens, potentially due to the hypogeal germination [80], clonal nature of this species, and low sample size.

Edaphic Effects
Compared to topography, we found that soil variables explained a greater proportion of the variance in stem abundance (14% vs. 6%) within the YFDP (Figure 6), although the total explained variance was low. Lin et al. [68] found that edaphic properties explained more variation in species distribution compared to the topographic variables by having the direct effect on the plant growth at local scales [81]. Potassium, phosphorus, calcium [82], and micronutrient deficiency [83] can limit plant growth and function. We found that the distribution of 45.5% of species was associated with edaphic properties (Table 2), consistent with results showing that 40% of species distribution was associated with soil nutrients [84]. The association of species to soil properties can be related to the direct effect of species characteristics on soil nutrients inputs and uptake, which contribute to species-soil associations as a function of species abundance [85]. We included soil enzymes in the list of soil variables due to their key role in ecosystem dynamics and biochemical functioning through the decomposition of soil organic matter and release of nutrients such as nitrogen (urease enzyme), and phosphorus (phosphatase enzyme) [12] into the soil. Soil enzymes are sensitive to small changes that occur in the environment and catalyze many essential processes necessary for soil microorganisms' life and affect the stabilization of soil structure. Their earlier response to soil disturbance compared to other soil quality indicators, made them an appropriate tool to evaluate the degree of soil alteration following fire. Soil enzyme activity showed a significant difference in urease activity between burned and unburned patches four years after fire occurrence (P = 0.01). This decrease may be related to the reduced microbial activity and biomass in the soil after fire. The decrease may also reflect the decreased soil pH in the burned microsites compared to the unburned patches (5.93 versus 7.07; P = 0.04). The long-term changes in soil acidity may affect microbial activity in burned sites and result in a higher release of urease in the unburned patches (higher pH) compared to those in the burned sites. Additionally, the reduced urease activity, which is the first hydrolytic enzyme involved in the breakdown of urea, may be related to the increase in non-hydrolysable N forms after fire [86,87].
We expected that the amount of inorganic N would have been higher (due to the activity of urease enzyme) in the unburned patches. However, there were no significant differences (P = 0.7) in NH4 + between the burned and unburned sites. This result may be related to the nutrient loss by leaching following the fire. Additionally, the availability of substrate (ammonium) to the nitrifying organisms may increase nitrification, which in turn leads to a decrease in the level of ammonium in the soil. Furthermore, the inclusion of soil enzyme activity improved (albeit by 5%) the explanatory power of soil properties in explaining variation in species stem abundance and basal area increment (Figure 6a-d). Soil enzymes (acid phosphatase and urease) alone were significant (P = 0.01) in their contribution to species abundance and basal area increment, even though the amounts of variation improvement explained by enzymes were small. The contribution of more explanatory variables (alkaline phosphatase and hydraulic conductivity shown in Figure S6) alone were not significant (P = 0.4) to species abundance and basal area increment.

Conclusions
The total number of species associated with habitats defined by soil properties was slightly greater than those associated with topographically-defined habitats. This finding suggests that niche partitioning caused by edaphic variables played a more important role compared to topographic variables in shaping species distributions. In addition, the contribution of spatial variables over topography and soil factors in explaining variation in species demographic metrics (stem abundance, mortality, and recruitment) indicates that community assembly was largely driven by spatially structured processes consistent with dispersal limitation and responses of species to the unmeasured environmental variables. Inclusion of two soil enzymes statistically improved predictions of species abundance and basal area increment, suggesting that future studies of soil enzymes may improve habitat definitions in forests. Adding soil enzymes to habitat definitions improved the explanatory power of edaphic variables to species abundance over the predictive ability of topography and soil nutrients alone. Species habitat associations and higher explanatory power of spatial factors compared to environmental variables suggest that both niche processes and dispersal limitations affect species distributions, but dispersal processes and unmeasured environmental variables were more important in the YFDP. The implication of a stronger contribution of neutral processes could reduce some concerns about the effects of increasing disturbance, decreasing habitat heterogeneity, and climate change on local species extinction in the future.