Effects of local biotic neighbors and habitat heterogeneity on seedling survival in a spruce‐fir valley forest, northeastern China

Abstract Seedlings are vulnerable to many biotic and abiotic agents, and studying seedling dynamics helps understand mechanisms of species coexistence. In this study, the relative importance of biotic neighbors and habitat heterogeneity to seedling survival was examined by generalized linear mixed models for 33 species in a spruce‐fir valley forest in northeastern China. The results showed that the relative importance of these factors varied with species and functional groups. Conspecific negative density dependence (CNDD) was important to the survival of Abies nephrolepis and Picea koraiensis seedling, whereas phylogenetic negative density dependence (PNDD) was critical to Pinus koraiensis and Betula platyphylla, as well as functional groups of tree, deciduous, and shade‐intolerant seedlings. For shrubs and Acer ukurunduense, habitat heterogeneity was significant. Despite of the significance of CNDD, PNDD, and habitat heterogeneity on seedling survival, large proportions of the total variance were not accounted for by the studied variables, suggesting the needs to examine the influences of other factors such as pests, diseases, herbivores, forest structure, species functional traits, and microclimatic conditions on seedling survival in the future.

Species of different attributes respond differently to intrinsic and extrinsic factors (Comita & Hubbell, 2009). Growth form may be a key trait in determining the patterns of species survival. Compared with shrubs, trees have larger statures and crowns and can thus produce more seeds and gather more resources; hence trees are often thought to have stronger NDD effects than shrubs (Terborgh, Zhu, Alvarez-Loayza, & Cornejo-Valverde, 2014). Some studies have found that leaf habit is related to species survival. Evergreen species can invest a higher proportion of resources to defensive compounds (such as lignin and tannin) than deciduous species; thus deciduous species are usually assumed to suffer more damages from extrinsic factors (Coley & Barone, 1996). Shade tolerance also is an important determinant of species' reactions to their local biotic neighbors because shadetolerant species are less susceptible to enemy (such as herbivores and pathogens) attack than light-demanding species, so can survive and recover from NDD effects more easily (Coley & Barone, 1996;Kitajima & Poorter, 2010;Kobe, 1997;Myers & Kitajima, 2007). Within a species, tree size could influence the effects of biotic and abiotic variables on survival; smaller individuals may be more sensitive to neighbor density due to asymmetric competition with taller individuals (Uriarte, Condit, Canham, & Hubbell, 2004). Species responses may also vary with dispersal modes, seed production, and individual abundance (Clark, Macklin, & Wood, 1998;Johnson, Beaulieu, Bever, & Clay, 2012).
In tree communities, the transition from seedling to sapling is generally a bottleneck in tree establishment (Queenborough, Burslem, Garwood, & Valencia, 2007). Compared with saplings, seedlings may suffer more from NDD effects and abiotic factors (Wright, 2002); thus, it is essential to explore the process of maintaining diversity and species coexistence in the early stages of plant communities.
Few studies have been undertaken to examine the relative importance of biotic neighbors by consideration of both phylogenetic relatedness and habitat variables in studying seedling survival patterns. In this study, we examined the relative importance of biotic neighbors and habitat variables in seedling survival over 4 years using a dataset of 6,256 seedlings and 33 species in a spruce-fir valley forest, northeastern China. Our specific questions are as follows: (1). How do biotic neighbors and habitat variables affect seedling survival? Is the CNDD more important than PNDD and habitat variables in seedling survival?
(2). Does the importance of biotic neighbors and habitats differ with growth forms, leaf habits, and shade tolerance? More specially, is CNDD or PNDD higher for trees than for shrubs, for deciduous trees than for evergreen trees, and for shade-intolerant trees than for shade-tolerant trees?

| Study site
The study was conducted in a spruce-fir valley forest in the Liangshui National Nature Reserve (128°53′20″ E, 47°10′50″N), northeastern China ( Figure 1). The elevation is 346-352 m ( Figure 2). The mean annual temperature is −0.3°C with a mean daily maximum temperature of 7.5°C and minimum temperature of −6.6°C. The mean annual precipitation is 676 mm with a total evaporation of 805 mm. The soil is dark-brown forest soil and 60 to −80 cm in depth [by Chinese classification, see Shi and Jin (2016)].
In 2006, a 9.12-ha (380 × 240 m) plot was established in the spruce-fir valley forest, and censuses were carried out every 5 years.
All woody stems with a diameter at breast height (DBH) ≥1 cm were tagged, identified by species, measured, and mapped ( Figure 2).

| Seedling census
In 2007, a total of 912 seedling quadrats (2 × 2 m) were established in the west corner of each 10 × 10 m subplot within the 9.12-ha plot ( Figure 2). All live woody seedlings ≥30 cm tall and <1 cm dbh were tagged, measured, and identified to species within each seedling plot every 2 years, and after 2011, we added small seedlings ≥10 cm and <30 cm tall to the census.

| Biotic neighborhood variables
Four biotic variables for each focal seedling individual were calculated to quantify local neighbor effects, two for seedling neighbors F I G U R E 1 An view of a spruce-fir valley forest in Liangshui, Heilongjiang province, northeastern China. The species composition primarily includes Abies nephrolepis, Picea koraiensis, Acer ukurunduense, Pinus koraiensis, Betula costata, and Larix gmelinii (conspecific neighbor densities and heterospecific neighbor densities) of trees and shrubs using 2011 data, and two for adult tree (all woody stems ≥ 1 cm dbh) neighbors (conspecific neighbor density index and average phylogenetic dissimilarity index) using 2016 data. For seedling neighbors, conspecific seedling neighbor density (S con ) and heterospecific seedling neighbor density (S het ) were calculated. The seedlings on the 10 m edge of the 9.12-ha study plot were excluded from the calculations, which included 120 of 912 total focal seedling quadrats.
For adult tree neighbors, both tree size and phylogenetic factors were considered in the calculations of the conspecific neighbor density index (A con ) and the average phylogenetic dissimilarity index of the heterospecific neighbors (A phylo ; Chen et al., 2016). Four biotic variables for each focal seedling individual were calculated as follows: where S con is the conspecific seedling neighbor density index; S het is the heterospecific seedling neighbor density index; A con is the conspecific adult tree density index; A phylo is the heterospecific adult tree phylogenetic dissimilarity index; N con is the total numbers of conspecific neighbors within the focal individual seedling quadrat (2 × 2 m); N het is the total numbers of heterospecific neighbors within the focal individual seedling quadrat (2 × 2 m); BA con i is the conspecific neighbor basal area; BA het i is the heterospecific neighbor basal area; i is the neighboring adult tree; W i is the Gaussian weight function, which represents the influence weight of neighbor adult tree i to the focal seedling; PD i is a pairwise phylogenetic distance between focal seedling and its adult tree neighbors; Dist i is the spatial distance between focal seedling and the ith adult tree neighbor; R is the neighborhood influence radius (Table S1); and n is the total number of heterospecific neighbors within the focal individual neighborhood influence radius.
The Gaussian weight function was used instead of the former inverse distance weighted function (i.e., the reciprocal of the distance between a focal individual and its neighbors) because the former method cannot explain the variations in neighbor influence radius among species. Moreover, species of different neighbor influence radius have different Gaussian density neighbor kernel. The Gaussian weight function (W i ) decreased with the spatial distance between focal seedling and its adult tree neighbors (Dist i ), which indicates the neighbors that are more closely to focal individual plant will have stronger effects ( Figure 3).
The neighbor influence radius depended on the specific species and its topographic position. Given the relative small variation in elevation (<6 m), we chose the valley position' neighbor radius for these seven major species (Table S1, Yang, 2006), and the nearest phylogenetic-related species' neighbor influence radius for species that are not listed in Table S1. A phylogenetic tree was established using Phylomatic based on angiosperm phylogeny group (APG) III backbone phylogeny (http://phylodiversity.net/phylomatic/).

| Abiotic variables
The abiotic variables included soil properties and topography. The 9.12-ha plot was divided into square grids of 20 × 20 m. The elevation was measured at every grid point of 20 × 20 m squares. For soil properties, three soil samples were taken to a depth of 10 cm at random distance combinations of "0, 2, and 5 m"; "0, 2, and 8 m"; or "0, 5, and 8 m" in a random direction from grid points. In total, 780 soil samples were obtained for the 9.12-ha plot ( Figure 2). Ten soil properties were recorded as follows: organic carbon (C), total nitrogen (TN), hydrolyzable nitrogen (HN), total phosphorus (TP), available phosphorus (AP), rapidly available potassium (RAK), soil pH, volumetric moisture (VM), mass moisture (MM), and bulk density (BD). The values of abiotic variables measured at 20 × 20 m grid were interpolated to a 5 × 5 m grid with ordinary kriging, a popular geostatistical analysis tools. A principal components analysis (PCA) was performed to reduce the multicollinearity of soil variables. The first three principal components accounted for 71% of total variation in the ten soil variables. The first principal component (31% of the total variation) was associated with high C, TN, pH, and soil moisture, the second principal component (27% of the total variation) with low C and AP, and high TP and BD, the third principal component (13% of the total variation) with high HN and RAK (Table 1).

| Data analysis
The survival of individual seedling from 2011 to 2015 was modeled as a function of biotic neighbors and abiotic habitat variables using the lme4 package in R software (R3.1.3; http://www.r-project.org) for generalized linear mixed models (GLMMs) with binominal distribution of errors (Bolker et al., 2009). The response variable was a logistic transformation of the seedling survival status, either 1 (alive) or 0 (dead). The seedling survival included seedling height as a covariate to account for age effects and two random variables, seedling plots to account for possible spatial correlation and species to account for species differences in their responses to neighborhood effects (Chen et al., 2010;Comita & Hubbell, 2009;Lin, Comita, Zheng, & Cao, 2012). The data of continuous explanatory variables were standardized prior to the statistical analysis.
To test the relative importance of biotic and abiotic variables, four candidate models were constructed: (1)  (AIC), and models with a AIC difference <2 were judged equally valid (Burnham & Anderson, 2002). The variance explained by fixed factors was included in marginal R 2 (R 2 mar ) and that by both fixed and random was in conditional R 2 (R 2 con ) of the models (Nakagawa, Schielzeth, & O'Hara, 2013).
Seedling survival was examined at three scales. First, the seedlings of all species were included with species and plot as random effects in the model. Second, species were included in the model by different functional groups of different growth forms (tree or shrub), species leaf habits (evergreen or deciduous), shade tolerance (shade-tolerant or shade-intolerant), with species and plot as random effects. Finally, the analysis of seedling survival was conducted by individual species for five most abundant tree species (n > 99 seedlings), including Abies nephrolepis, Picea koraiensis, Acer ukurunduense, Pinus koraiensis, and Betula platyphylla, with plot as a random effect.

| RESULTS
Of the 6,256 tree and shrub seedlings tagged in 2011, 3,581, 57.24%, were still alive in 2015. The variables used in the models were summarized in Table 2. The probability of seedling survival generally increased with seedling height (Figures 4-6), and effects of other variables in the best-fit models are described below.

| All-species model
All biotic and abiotic factors were retained in the best-fit seedling survival model for all species. The fixed factors (seedling height, biotic neighbors, and habitat variables) explained 11.4% of the variance (Tables 3 and S2). Both heterospecific seedling density and soil PC1 had significant, positive effects on seedling survival, as indicated by coefficient estimates (>0). In contrast, seedling survival decreased with conspecific adult tree density and elevation (Figure 4).

| Growth form
Again, the best-fit model for trees was the full model, with fixed factors (seedling height, biotic neighbors, and habitat variables) explaining 14% of the total variance. The seedling survival was significantly increased with heterospecific seedling density and soil PC1, and marginally increased with heterospecific phylogenetic dissimilarity, but decreased with elevation. In comparison, the survival of shrub species was best described by the abiotic model, with fixed factors (seedling height and habitat variables) explaining 4.8% of the total variance. The seedling survival of shrubs decreased with soil PC3 that was primarily associated with HN and RAK (Tables 3 and S2, Figure 5a,b).

| Leaf habit
The best survival model by leaf habits was the full model, with the fixed factors (seedling height, biotic neighbors and habitat variables) explaining 14.8% of the total variance for evergreen seedlings and 11.3% for deciduous seedling (Tables 3 and S2). A positive effect of heterospecific phylogenetic dissimilarity on seedling survival was stronger for deciduous seedlings than for evergreen seedlings. The model for evergreen seedlings included positive effects of heterospecific seedling density and soil PC1, but negative effects of elevation and soil PC2. Comparatively, the survival of deciduous seedlings increased with higher soil PC1 (Figure 5c,d).

| Shade tolerance
The best-fit model for shade-tolerant seedlings was the full model, whereas that for shade-intolerant seedlings was the biotic model. The fixed factors (seedling height, biotic neighbors, and habitat variables) explained 12.7% of the total variance for shade-tolerant seedlings and 10.7% for shade-intolerant seedlings (Tables 3 and S2). For shadetolerant seedlings, their survival increased with heterospecific seedling density and soil PC1, but decreased with increasing elevation and soil PC2. For shade-intolerant seedlings, their survival significantly increased only with heterospecific phylogenetic dissimilarity ( Figure 5).

| Individual species model
The best-fit model for A. nephrolepis and P. koraiensis was the full model that explained 22.5% and 20.1% of the total variance by the fixed factors (seedling height, biotic neighbors, and habitat variables).
In comparison, the best-fit model for A. ukurunduense was the abiotic model, explaining 10.8% of the total variance by fixed factors (habitat variables) and that for P. koraiensis and B. platyphylla was the biotic model explaining 9.5% and 18.7% of the total variance, respectively, by the fixed factors (Tables 3 and S2). The seedling survival of A. nephrolepis species increased with heterospecific seedling density, conspecific adult tree density, and elevation (Figure 6a), while that, for P. koraiensis decreased with conspecific seedling density and

| Effects of density, habitat, and phylogenetic dissimilarity
At all levels of survival analyses, seedling height was positively associated with seedling survival, consistent with the findings by others (Metz et al., 2010). The possible reason is that large seedlings are likely less susceptible to biotic stresses (e.g., herbivores and pathogens) and in better positions in competition for light (Chen et al., 2010;Comita & Hubbell, 2009;Queenborough et al., 2007).
Numerous studies have indicated that NDD is an important mechanism for maintaining species coexistence at seedling stage (Bagchi, Press, & Scholes, 2010;Johnson et al., 2014;Nathan & Muller-Landau, 2000;Yan, Zhang, Wang, Zhao, & von Gadow, 2015). This is also true in our study with a strong CNDD effect on seedling survival (Figures 4   and 6a). This is likely due to the fact that conspecific adult trees have similar resource requirements, as well as share similar pests and pathogens with seedlings, which is conductive to high competition and mortality. On the other hand, seedling survival was not negatively affected by conspecific seedling density, which is similar to the findings by Paine, Harms, Schnitzer, and Carson (2008) and Svenning, Fabbro, and Wright (2008), but different from those by Chen et al. (2010) and Comita and Hubbell (2009). This inconsistency may result from the levels of competition among seedlings; the seedling of small size and low density would not be strongly influenced by CNDD (Paine et al., 2008;Svenning et al., 2008). As expected, seedling survival was positively associated with heterospecific seedling density (Figures 4, 5a,c,e, and 6a,b,d), an indication of heterospecific positive density dependence (HPDD). As suggested by the "species herd protection hypothesis" (Peters, 2003), the seedlings of more heterospecific neighbors would have less competition and therefore higher chance of survival.
Our results showed that heterospecific phylogenetic dissimilarity was positively associated with seedling survival, indicating positive effects of PNDD in the spruce-fir valley forest, northeastern China (Figures 4, 5a,c,d,f, and 6d,e). Focal seedlings would have a low survival when surrounded by conspecific tree neighbors, due to increased competition for resources and sharing of common pests and pathogens (Gilbert & Webb, 2007;Liu et al., 2012;Novotny et al., 2002).
On the other hand, conspecific tree neighbors could also be positively  total carbon, and soil total nitrogen were all positively associated with seedling survival in the spruce-fir valley forest, northeastern China.
The possible reasons would be that seedlings have shallow roots and therefore limited access to resources, and should benefit from high moisture, pH, carbon, and nitrogen in the soil (Engelbrecht, Kursar, & Tyree, 2005;Lin et al., 2012). This may also be applied to variations in topography that the increase in elevation was associated with decreased total carbon, total nitrogen, and moisture content in the soil (p < 0.001), leading to decrease in seedling survival (Figures 4, 5a,c,e, and 6c). The strong influences of abiotic factors suggest that the seedlings in the spruce-fir valley forest are restricted by soil moisture and nutrients and very sensitive to changes of resources availability, even by slight variation in elevation (<6 m).
Small proportions of the total variance (12.6%-58.5%) were accounted for by the selected fixed and random variables, indicating significant influences of other factors, such as pests, diseases, forest structure, species functional traits, and microclimatic conditions on seedling survival. Future studies should examine the influences of these factors, particularly under the context of climate change. In this study, random effects accounted for more variance than fixed effects (Table S2), indicating significant differences in seedling survival among species and spatial variations of site conditions.

| Differences among functional groups
Our results do not support the suggestions that trees produce more seeds than shrubs and are likely to affected more by CNDD (King, Wright, & Connell, 2005;Terborgh & Petren, 1991;Terborgh et al., 2014). Comparatively, tree seedlings were more influenced by biotic adult tree neighbors than shrubs, indicating greater needs for resource sharing and pressures for competition. Between evergreen and deciduous species, deciduous seedlings uptake and allocate more resources on growth, whereas evergreen seedlings tend to allocate more on defensive compounds, which would help survival (Coley & Barone, 1996;Villar, Robleto, De Jong, & Poorter, 2006). The greater growth by deciduous seedlings would require a greater level of resource sharing and therefore are more influenced by PNDD. We did not see greater influences of CNDD on survival of shade-intolerant seedlings due to their low carbohydrates and tissue density, and high vulnerability to competition and external damages (Piao, Chun, Yang, & Cheon, 2014); in facts, the survival of shade-intolerant seedlings was not affected by conspecific neighbor density. On the other hand, the survival of shade-intolerant seedlings decreased with the amount of phylogenetic-related neighbors, indicating greater importance of PNDD than CNDD in survival of shade-intolerant seedlings.

| Differences among individual species
Different species have difference physiological and functional traits and therefore different critical factors for seedling survival. In our study, the survival of A. nephrolepis seedlings was negatively correlated with conspecific adult tree density, while that of P. koraiensis seedlings was negatively with conspecific seedling density (Figure 6a,b). The CNDD influences on seedling survival, however, disappeared in the all-species model that excluded A. nephrolepis and P. koraiensis indicating that the CNDD effects were mainly on the two species ( Figure S3).
Their survival was affected by both biotic neighbors and abiotic factors according to the best-fit models. Comparatively, the survival of A. ukurunduense depended on abiotic factors such as high soil nutrient contents and moisture and therefore niche partitioning, whereas heterospecific phylogenetic dissimilarity was critical to the survival of a P. koraiensis and B. platyphylla.

| CONCLUSION
The seedling survival in a spruce-fir valley forest, northeastern China, depended on both biotic neighbors and habitat heterogeneity (elevation and soil properties). However, the relative importance of these