Phylogenetic and Functional Traits Verify the Combined Effect of Deterministic and Stochastic Processes in the Community Assembly of Temperate Forests along an Elevational Gradient

: Explaining community assembly mechanisms along elevational gradients dominated by deterministic processes or stochastic processes is a pressing challenge. Many studies suggest that phylogenetic and functional diversity are signiﬁcant indicators of the process. In this study, we analyzed the structure and beta diversity of phylogenetic and functional traits along an elevational gradient and discussed the effects of environmental and spatial factors. We found that the phylogenetic and functional traits showed inconsistent changes, and their variations were closely related to the abiotic environment. The results suggested that the community assembly of woody plants was obviously affected by the combined effect of deterministic processes and the stochastic hypothesis (primarily by the latter). Phylogenetic and functional traits had a certain relationship but changed according to different rules. These results enhance our understanding of the assembly mechanism of forest communities by considering both phylogenetic and functional traits.


Introduction
The assembly mechanism of plant communities has always been an important topic in ecological research [1]. Understanding what drives the variation in a community can provide insights into the mechanisms of community assembly [2]. There are two main hypotheses used to explain the mechanism of community assembly: the deterministic hypothesis based on niche theory and the stochastic hypothesis based on neutral theory. The deterministic hypothesis states that the restriction of environmental conditions and the interaction between organisms will screen the species entering the community and result in a stable community. The main drivers of this hypothesis are habitat filtering and competitive exclusion [3]. In contrast, the stochastic hypothesis states that the probability of entering a community is equal among all species; niche differences and competition do not affect species coexistence, community assembly is a random ecological drift process of equivalent individuals, and the main drivers are dispersion limitation and random effect [4,5]. Both hypotheses could solve some problems; however, some existing phenomena remain unexplained [6]. Therefore, it is generally believed that the two hypotheses should be incorporated to explain the mechanism of community assembly [7,8]. Some studies found that both hypotheses reflected different aspects of community construction, and the more important hypothesis may depend on the current environmental conditions [9]. For example, the deterministic hypothesis may play a greater role in adverse environments [10]. bly mechanisms along elevational gradients. To date, many studies have shown that the deterministic process dominates community assembly along elevational gradients, and the change patterns of phylogenetic and functional traits are also different [23,45,46]. In this study, based on elevational gradient transect data and other scholars' conclusions, we propose theoretical expectations: the changes in phylogenetic and functional traits along the elevational gradient are consistent in mountain forest communities, the changes in abiotic environmental gradients have certain effects on the beta diversity of phylogenetic and functional traits, and the combined effect of deterministic processes and stochastic hypotheses guide community assembly.
To test these hypotheses, we analyzed the structure and beta diversity of phylogenetic and functional traits, tested the phylogenetic signals of functional traits, and discussed the interpretive ability of environmental factors to explain the phylogenetic and functional traits.

Study Site
Donglingshan Mountain is an extension of the Xiaowutaishan Mountains and belongs to the broader Taihangshan Mountains, and is located 100 km northwest of Beijing city, China. The study area, the Beijing Forest Ecosystem Research Station of the Chinese Academy of Sciences, is located at 40 • 00 -40 • 03 N and 115 • 26 -115 • 30 E. The soil type of the area is brown soil classified as Eutric cambisol. The area has a typical warm temperate continental monsoon climate with an average annual precipitation of 500-650 mm. The mean annual temperature is 5-10 • C. The elevation of most of the area is more than 1000 m above sea level, and the elevation of the highest peak is 2303 m.

Data Collection
We chose the well-developed Quercus wutaishanica Blume. forest on the west slope, and ten transects were established according to the elevational gradient. Several 10 m × 10 m quadrats were separated along each transect, with a total of 119 quadrats (Appendix A). The coordinates, elevation, slope exposition, and slope were recorded, and standard plant community surveys were carried out. We investigated the species identity, diameter at breast height (DBH), and height of all tree species, divided the quadrat into four 5 m × 5 m quadrats, and investigated the species identity, quantity, height, and coverage of all shrub species. Additionally, soil samples were collected by a 5 cm diameter earth drill, and five sites were selected randomly from each 10 m × 10 m quadrat to collect the surface soil (0-10 cm); these were mixed into one composite sample and used for soil physiochemical property analysis after passing through a 2 mm diameter soil sieve and being placed indoors to naturally air dry [47]. Soil temperature and moisture were measured by a soil moisture tester (TDR 350, Spectrum); soil pH was measured by a pH meter (PB-10, Sartorius); soil total carbon, organic carbon, and total nitrogen content were determined by an element analyzer (Vario EL III, Elementar); and soil total phosphorus content was determined by a continuous flow analyzer (SKALAR SAN++, SKALAR) [48].
Leaves of all mature woody plants were collected in August 2019. We collected leaves from each species. At least five individuals in each 10 m × 10 m quadrat and 10-20 intact, mature, and healthy leaves were collected on the dayside near the canopy for each individual, and the functional properties were determined after returning to the laboratory. We selected 13 traits that reflect the characteristics of plants: leaf area (LA), specific leaf area (SLA), leaf length (LL), leaf width (LW), leaf thickness (LT), and leaf dry matter content (LDMC) were selected for the chemical properties of leaves; leaf total carbon content (LTC), leaf total nitrogen content (LTN), leaf total phosphorus content (LTP), leaf total potassium content (LTK), the ratio of carbon to nitrogen of a leaf (L C/N), and the ratio of nitrogen to phosphorus of a leaf (L N/P) were selected for the physical properties of leaves; finally, the relative content of chlorophyll (RCC) was selected for the photosynthetic capacity of a leaf (Appendix B). LL, LW, and LT were measured by a Vernier caliper (DL 3941, Deli, Ningbo, China); LA was collected by a CanoScan LiDE 220 scanner (Canon, Tokyo, Japan) and measured by ImageJ 1.49v software (Wayne Rasband, National Institutes of Health, Stapleton, NY, USA); the relative chlorophyll content was measured by a SPAD-502plus chlorophyll meter (KONICA MINOLTA, Tokyo, Japan); leaf fresh weight was measured by an electronic balance; leaf dry weight was also measured by an electronic balance after drying in an oven at 80 • C for 48 h; and LTC, LTN, LTP, and LTK were measured by a J200 LIBS laser spectral element analyzer (ASI, Sacramento, CA, USA).

Phylogenetic Tree
According to the investigation, 30 species of woody plants were recorded in 10 transects, and there were 15 tree species, such as Quercus wutaishanica Blume, Fraxinus chinensis Roxb., and Acer mono Maxim. The shrub species included Spiraea trilobata L., Abelia biflora Turcz., Rhododendron micranthum Turcz., and 12 other species (Appendix C). We counted the family and genus information of all woody plants, and a phylogenetic tree of woody plants was reconstructed by using phylomatic (V3) software based on the Zanne 2014 phylogenetic hypothesis (Appendix E) [49]. Quercus wutaishanica Blume and Spiraea trilobata L. were the absolute dominant species in the tree layer and shrub layer, respectively, and the two populations had high abundances and were present in each transect.

Data Analysis
To test the correlation between functional traits and species evolutionary history, phylogenetic signals of functional traits were calculated by Blomberg's K method [50]. If the P value was not significant, there was no phylogenetic signal. Significant values lower than 0.8 indicate a low/weak phylogenetic signal, values within 0.8-1.1 indicate the trait follows a Brownian motion, and values > 1.1 indicate a strong phylogenetic signal. Phylogenetic and functional trait structures were calculated by the mean phylogenetic distance between species (MPD) and the mean nearest taxon distance (MNTD). Phylogenetic and functional trait beta diversities were calculated by the mean pairwise distance (Dpw) and the nearest neighbor distance (Dnn), respectively. We transformed the empirical value to a standardized effect size (SES value) by subtracting the mean of the null distribution and dividing by the standard deviation of the null distribution, and species richness was considered in the calculation process [51,52]. The null model "richness" was used, and the analysis was completed in R 4.0.3 [53]. Structure overdispersion occurred if the SES value > 0, clustering occurred if the SES value < 0, the structure was completely random if the SES value = 0, and the SES value was significant when it was greater than 1.96 or less than −1.96.
Environmental factors and spatial factors were investigated as variables for redundancy analysis (RDA) of the structure of phylogenetic and functional traits. Elevation, slope exposition, slope, soil temperature, soil moisture content, soil pH, soil total carbon content, soil total nitrogen content, soil total phosphorus content, soil organic carbon content, ratio of soil carbon to nitrogen, and the ratio of soil nitrogen to phosphorus were chosen as environmental factors. Six principal coordinates of neighboring matrices (PCNM) spatial variables were generated according to the PCNM based on the longitude and latitude of transects and used as spatial factors. Monte Carlo permutation (n = 999) and forward selection were used to screen the significant environmental factors and spatial factors (p < 0.05) as explanatory variables. Then, we carried out RDA and adopted the adjusted R 2 value. Variation partitioning was used to understand the explanatory ability of each factor to the structure of phylogeny and functional traits. Multiple regressions on distance matrices (MRM) analysis was used to test the effects of environmental factor distance and spatial factor distance on phylogenetic and functional trait beta diversity. MRM entails a multiple regression of a response distance matrix on two or more environmental, spatial, or other explanatory distance matrices, each unfolded into a distance vector. The significance of an MRM model and its regression coefficients are tested by permuting matrices while holding the explanatory matrices constant [54]. Phylogenetic Dpw and Dnn indices were analyzed as response matrices, while environmental factor and spatial factor matrices were regarded as explanatory distance matrices. Furthermore, phylogenetic Dpw and Dnn indices were added as explanatory distance matrices when functional traits' Dpw and Dnn indices were analyzed. MRM analysis could test the correlation between response matrices and explanatory distance matrices, and the Bray-Curtis distance was used as the matrix distance index.

Phylogenetic Signal
We analyzed the phylogenetic signal intensity of functional traits by using Blomberg's K index ( Table 1). The results showed that LTC, LTN, LTP, and L C/N and L N/P had significant p values, but Blomberg's K values were lower than 0.8, which meant that they showed weak phylogenetic signals. Meanwhile, no phylogenetic signals were detected from other traits. Therefore, we found that the accumulation of nutrient elements in leaves had a more obvious phylogenetic signal than other functional traits.

Phylogenetic and Functional Trait Structures along Elevational Gradients
We found that the phylogenetic SES.MPD showed a trend of lower elevations greater than zero and higher elevations less than zero (Figure 1a), and this was also found for the phylogenetic SES.MNTD (Figure 1b), but the SES values were significant at high elevations only (<−1.96). This result suggested that the phylogenetic structure of woody plants changed from a stochastic process to clustering with increasing elevation. However, the trait SES.MPD showed a unimodal curve pattern that changed from less than zero to greater than zero and then back to less than zero (Figure 1c), and the trait SES.MNTD was significant (<−1.96) at most elevations ( Figure 1d); the trend was different from that of the phylogenetic structure. This result suggested that the functional trait structures of woody plants were stochastic, and the trait distributions of relative species showed obvious clustering. tions only (<−1.96). This result suggested that the phylogenetic structure of woody plants changed from a stochastic process to clustering with increasing elevation. However, the trait SES.MPD showed a unimodal curve pattern that changed from less than zero to greater than zero and then back to less than zero (Figure 1c), and the trait SES.MNTD was significant (<−1.96) at most elevations ( Figure 1d); the trend was different from that of the phylogenetic structure. This result suggested that the functional trait structures of woody plants were stochastic, and the trait distributions of relative species showed obvious clustering.

Phylogenetic and Functional Trait Beta Diversity along an Elevational Gradient
The phylogenetic SES.Dpw of woody plants at adjacent elevations showed a trend of lower elevations greater than zero and higher elevations lower than zero (Figure 2a), while the beta diversity of the functional traits showed the opposite trend (Figure 2c), and only few SES values were significant (>1.96 or <−1.96). The phylogenetic SES.Dnn showed that the extent of evolution was lower than the expected value of a random process and showed random changes along the elevational gradient, and similar results were found for the functional trait SES.Dnn (Figure 2b,d).

Phylogenetic and Functional Trait Beta Diversity along an Elevational Gradient
The phylogenetic SES.Dpw of woody plants at adjacent elevations showed a trend of lower elevations greater than zero and higher elevations lower than zero (Figure 2a), while the beta diversity of the functional traits showed the opposite trend (Figure 2c), and only few SES values were significant (>1.96 or <−1.96). The phylogenetic SES.Dnn showed that the extent of evolution was lower than the expected value of a random process and showed random changes along the elevational gradient, and similar results were found for the functional trait SES.Dnn (Figure 2b,d).

Environmental Interpretation of Phylogenetic and Functional Traits
Significant environmental factors were selected as variables through forward selection (Appendix D). The results showed that environmental factors and spatial factors explained 5.4% and 4.5% of the variation in the phylogenetic SES.MPD, respectively. The

Environmental Interpretation of Phylogenetic and Functional Traits
Significant environmental factors were selected as variables through forward selection (Appendix D). The results showed that environmental factors and spatial factors explained 5.4% and 4.5% of the variation in the phylogenetic SES.MPD, respectively. The interaction between the environmental and spatial factors explained 35.6% of the variation, and the unexplained part remained at 54.5% (Figure 3a). Additionally, both the single interpretation rate and the common interpretation rate were low for the phylogenetic SES.MNTD (residual = 0.723) (Figure 3b).    To explore the impact of phylogeny on functional traits, we added the phylogenetic structure into the functional traits' RDA as an explanatory variable. The screening method for significant variables was the same as that described above (Appendix D). The results showed that environmental factors explained 4.3% of the variation in the functional traits' SES.MPD, while the spatial factors explained 11.6% of the variation alone. However, the individual explanatory ability of phylogenetic factors was relatively low. The interaction of the three factors explained 5.8% of the variation, and the unexplained part remained at 77.6% (Figure 3c). The functional traits' SES.MNTD and environmental factors explained 14.1% alone, and the interpretation rates of the spatial and phylogenetic factors were similar (9.3% and 8.6%, respectively). The interaction of the three factors explained 15.8% of the variation, and the unexplained part remained at 50.5% (Figure 3d).

Discussion
We used MRM to test the correlation between the phylogenetic and functional trait distance matrices. The results showed that the spatial and environmental factor distances were both extremely significantly correlated with the phylogenetic beta diversity distance, but the spatial factor distance had a greater influence. The interaction of the two factors also significantly affected the phylogenetic beta diversity between transects. Similarly, to test the impact of phylogeny on functional traits, we added phylogenetic beta diversity to the MRM analysis of functional traits as a factor. The results suggested that the phylogenetic beta diversity distance was extremely significantly correlated with the functional trait beta diversity distance. The spatial factor distance had a significant correlation with the beta diversity distance of functional traits, while the environmental factor distance had no significant correlation. Their impact on the beta diversity of functional traits also decreased progressively. Furthermore, the interaction of the three factors significantly affected the beta diversity of the functional traits between transects ( Table 2).

Discussion
Phylogenetics could reflect the complex adaptive strategies or some unmeasured conservative functional traits of species, which are closely related to community assembly mechanisms. Therefore, the method combined with the analysis of phylogenetic and functional traits has become an effective research tool that can be used to evaluate the ecological process of community assembly [36,46,56]. Furthermore, this method could complement studies on the relationships between conventional forest stand attributes and structural indices in order to help resolve this issue, which is of high practical importance for forest managers [57,58]. In this study, we analyzed the phylogenetic and functional traits of 30 woody plants along an elevational gradient on Mt. Dongling. The lowest elevation we selected represents the typical low-elevation forest vegetation in this area, while the highest elevation is near the tree line; thus, this elevational gradient is ideal for investigating the phylogenetic patterns and functional traits of woody plant assemblages in the study area [5]. We tested the mechanism of community assembly by the changes in phylogenetic and functional traits along the elevational gradient. The main results showed that the phylogenetic and functional traits exhibited general random component and less nonrandom component transformation along the elevational gradient, indicating that the stochastic hypothesis played the more important role in controlling community transformation.
The effects of elevation on phylogenetic and functional traits have been proven in many studies. In general, the structure of phylogenetic and functional traits changed regularly along the elevational gradient, including consistent and inconsistent changes [46,59]. Many studies found that competitive exclusion leads to the overdispersion of phylogenetic (functional) structure, while environmental filtering leads to clustering [60][61][62]. Some studies have suggested that the relative importance of competitive exclusion in high elevational areas is greater than that of environmental filtering [63,64], but our research showed the opposite result. We found that the significant phylogenetic SES values (less than −1.96) appeared in the high-elevation area; they showed that the woody plants were clustered and indicated that the environmental filtering significantly affected the distribution of woody plants. The clustering of relative species at high elevations also confirmed this conclusion. However, the traits' SES values showed considerable randomness, which suggested that the change in elevation in this area did not cause significant differences in the functional traits of woody species, perhaps because the change in elevation was not correlated with great changes to the environment or because the degree of environmental change has not high enough to cause a different response among woody species. In addition, functional traits' SES.MNTD values were less than zero and greater than −1.96 along most of the elevational gradient, indicating that the response of the traits of relative species to environmental changes were similar, and the change in elevation had no obvious effect on the traits of relative species.
The beta diversity also differed between phylogenetic and functional traits. Most SES values of structure and beta diversity were not significant (between −1.96 and 1.96), indicating that only a few elevations had deterministic processes, and there were many random effects in the community assembly of this region. Although most of the SES values were not significant (indicating the dominant role of randomness), the SES.Dpw showed a completely opposite change trend, and the SES.Dnn showed a similar trend; furthermore, the results were not consistent with our hypotheses. Several studies suggested that the changes between phylogenetic and functional traits were relatively independent, but others showed a high correlation [65,66]. A study also suggested that the beta diversity of functional traits increases with elevation [52]. It is worth noting that the beta diversity was significant in both some low-elevation and high-elevation regions. Combined with the overall trend, we considered that elevation gradient had a certain impact on phylogenetic and functional traits, which led to the opposite changes, but most of the changes did not reach a significant degree of difference. Furthermore, the SES.Dnn showed a similar negative curve, which indicated that the variation along the elevation gradient had few impacts on the phylogenetic and functional traits of relative species. This is also consistent with the results of structure analysis. One study indicated that phylogenetic diversity may be used as a crude surrogate measure of functional diversity [67], but we concluded that phylogenetic diversity should not be used as a surrogate of functional diversity since the change trends in functional traits differed with phylogeny.
Stronger environmental filtering leads to phylogenetic clustering that is easier than competition exclusion when functional traits with significant phylogenetic signals participate in the community assembly process [63]. Another study demonstrated that the competition pressure of species was a crucial factor driving and regulating the distribution of biomass [68]. We detected the phylogenetic signals of 13 kinds of leaf functional traits, and the elemental contents showed some weak phylogenetic signals. The signals of elemental contents were relatively significant, which was similar to the results of other studies [46,67,69]. Combined with the changes in phylogenetic and functional traits along the elevational gradient, the results also verified our study. The accumulation of nutrient elements in plant leaves more easily showed obvious signals of phylogeny than other functional traits, which indicated that the influences of phylogeny on the functional traits of woody plant leaves might originate from changing the accumulation of nutrient substances. In addition, the insensitivity of functional traits to phylogenetic signals shown in our results might be related to the selection of traits directly.
The change in species along the environmental gradient was not only due to diffusion limitations but also due to the underlying environmental gradient that plays a pivotal role in species replacement and community assembly [52,70,71]. Several studies have confirmed that an important part of the abiotic environment cannot always be clearly quantified due to strong spatial autocorrelation [15,34,72]. To verify the effect of the deterministic process, we screened and analyzed the interpretive abilities of environmental factors to phylogenetic and functional traits. The mountainous environment provided a natural laboratory for exploring the relationship among environmental variables, spatial variables, and elevational gradient diversity [73]. Our studies indicated that both environmental and spatial variables could explain phylogenetic and functional traits to a certain extent, which also verified our hypothesis, but phylogeny had a low explanation rate for the functional trait structure. One study showed that spatial and environmental variables determined the structure of community functional traits, while phylogenetic information had little influence [74]; the results were consistent with ours. Furthermore, we found that the trait variation of relative species was related to environmental conditions, spatial changes, and phylogenetic structure of species through the high explanation rate for the SES.MNTD. In general, the variation in phylogenetic and functional trait structure was clearly related to the abiotic environment, which also strengthened the support of our hypothesis. Additionally, abiotic environmental filtering could not explain the community assembly process entirely. We suggested that the phylogeny showed a lesser impact on functional traits, and the high residuals in variation partitioning might also confirm the existence of environmental factors and random effects that have not been considered.
It has been found that the interaction between abiotic environmental gradients dominates the variation in beta diversity in phylogenetic and functional traits [24]. Different species have their own distribution range along the elevational gradient, and spatial location, seed diffusion, environmental conditions, and other factors limit the expansion of species in the new environment [75,76]. Plants at different elevations will be affected by their own genes and spatial conditions, and with changes in their surroundings, they will experience corresponding changes in functional traits to adapt to the current habitats. The results of MRM analysis also showed that there were close correlations among environmental variables, spatial variables, and phylogenetic and functional traits, which confirmed our hypothesis. However, how phylogeny affects functional traits remains to be further analyzed from the perspective of molecular biology.

Conclusions
In this study, we used the structure and beta diversity of phylogenetic and functional traits to test the hypotheses of their relationship. The results suggested that the change patterns of phylogenetic and functional traits varied due to the different main drivers of community evolution, but strong random effects were found at most elevational gradients, and the combined effect of deterministic processes and stochastic hypotheses (with the primary role of stochastic processes) guided community assembly along the elevational gradient. Additionally, the weak phylogenetic signals indicated the slight effect of phylogeny on the selected functional traits along the elevational gradient. The variation in phylogenetic and functional trait structure had a certain correlation to the abiotic environment, but abiotic environment filtering could not explain the community assembly process entirely. Elevation gradient had a certain impact on phylogenetic and functional traits, but most of the changes did not reach a significant degree of difference. Our study emphasized the importance of joint research on phylogenetic and functional traits when exploring the complex mechanisms driving community assembly.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy restriction.

Acknowledgments:
We are grateful to all the students from Tianjin Chengjian University and Beijing Forestry University for their help with field data collection and laboratory experimental work. We also thank Ran Wu from Zhejiang University for her help with data analysis.

Conflicts of Interest:
The authors declare no conflict of interest.  Note: " " indicates that the species is distributed in the transect; " " indicates that the species is not distributed in the transect.