Disentangling the drivers of taxonomic and phylogenetic beta diversities in disturbed and undisturbed subtropical forests

Understanding the relative importance of dispersal limitation and environmental filtering processes in structuring the beta diversities of subtropical forests in human disturbed landscapes is still limited. Here we used taxonomic (TBD) and phylogenetic (PBD), including terminal PBD (PBDt) and basal PBD (PBDb), beta diversity indices to quantify the taxonomic and phylogenetic turnovers at different depths of evolutionary history in disturbed and undisturbed subtropical forests. Multiple linear regression model and distance-based redundancy analysis were used to disentangle the relative importance of environmental and spatial variables. Environmental variables were significantly correlated with TBD and PBDt metrics. Temperature and precipitation were major environmental drivers of beta diversity patterns, which explained 7–27% of the variance in TBD and PBDt, whereas the spatial variables independently explained less than 1% of the variation for all forests. The relative importance of environmental and spatial variables differed between disturbed and undisturbed forests (e.g., when Bray-Curtis was used as a beta diversity metric, environmental variable had a significant effect on beta diversity for disturbed forests but had no effect on undisturbed forests). We conclude that environmental filtering plays a more important role than geographical limitation and disturbance history in driving taxonomic and terminal phylogenetic beta diversity.

Phylogenetic tree. Botanical nomenclature of the species in the forest plots were standardized according to The Plant List (version 1.1; http://www.theplantlist.org/). The ten plots included 376 native tree species belonging to 61 families and 147 genera. We used genera and species present in our study plots to prune the mega-phylogeny PhytoPhylo 38 , which is an updated version of the phylogeny published in Zanne et al. 39 . This mega-phylogeny was constructed based on seven gene regions (i.e., 18S rDNA, 26S rDNA, ITS, matK, rbcL, atpB, and trnL-F), which include both slowly and quickly evolving regions 39 . The time scale for the phylogeny was based on 39 fossil calibrations 39 . Both minimum and maximum age constraints were utilized for each fossil calibration. All of the 61 families and 146 of the 147 genera in our study plots were included in PhytoPhylo. Thus, the phylogeny used in our study was completely resolved at the family level and nearly completely (99.3%) resolved at the genus level (Fig. 2). Of the 376 species present in our plots, 265 were also included in PhytoPhylo and additional 16 species belong to genera with only one species in our data and thus were represented by branches of their respective genera. Thus, the vast majority (75%) of the species present in our plots were resolved in the phylogeny extracted from PhytoPhylo. For the sole genus (i.e., Cerasus) and the remaining species of the study plots that were not in PhytoPhylo, we used the software S.PhyloMaker 38 to add them to the phylogeny, using Scenario 3, which is analogous to using Phylomatic with BLADJ to generate a phylogeny 40 .

Beta diversity indices.
Previous studies suggest that multiple beta diversity indices be used in an analysis 15,16 because different indices may emphasize on different aspects of community similarity and thus are complementary to each other to some degree. Accordingly, we used more than one beta diversity index for each of the two beta diversity categories (i.e., TBD vs. PBD) to quantify beta diversity among the ten forest plots. For TBD, we used Bray-Curtis and Jaccard indices 3,41 . PBD indices may be broadly divided into two classes: basal metric vs. terminal metric 15 . A basal metric of PDB emphasizes evolutionary divergence near the root node of the phylogenetic tree, which reflects historical events occurring in the remote past, whereas a terminal metric of PBD emphasizes evolutionary divergence near the terminal tips of the phylogenetic tree, which represent recent evolutionary events. Previous studies used PhyloSor, UniFrac and Dnn as terminal metrics and Dpw, Rao's D and Rao's H as basal metrics 15,16,42 . Because UniFrac and Rao's D are redundant with PhyloSor and Dpw, respectively 16 , we used two terminal (PhyloSor, and Dnn) and two basal (Dpw and Rao's H) metrics for PBD in the present study (Table S2). Larger values in Bray-Curtis, Jaccard, Dnn, Dpw and Rao's H represent higher beta diversity.
Geographical and environmental data. Geographical distance (GeoDist) between a pair of plots was measured as the Euclidean distance between the centers of the plots using DIVA-GIS software (http://www.divagis.org/). Environmental distance was measured from eight climatological parameters, which were extracted from the WORLDCLIM database (www.worldclim.org; for the period of 1950-2000 at a 30 arc-second resolution) using DIVA-GIS software 43 , including mean annual temperature (MAT), temperature seasonality (TS, i.e. standard deviation *100), maximum temperature of the warmest month (MTWM), minimum temperature of the coldest month (MTCM), annual precipitation (AP), precipitation of the wettest month (PWM), precipitation of the driest month (PDM) and precipitation seasonality (PS). Two additional climate-based parameters, annual actual evapotranspiration (AET) and annual potential evapotranspiration (PET), were obtained from the MODIS evapotranspiration data set (MOD16) 44 . Three topographic variables, elevation (ELE), aspect (ASP) and slope (SLOPE), were measured with a handheld global positioning system (GPS) receiver and clinometer. Aspect (ASP) is a circular variable and standardized using the formula of -cos((2πASP)/360) to make the maximum value at South and the minimum value at North so that it can be used in linear models (see Table S1).

Statistical analysis.
To disentangle the effect of environmental variables and spatial variables on beta diversity, two complementary approaches were used 45 . First, a principle component analysis (PCA) was performed on all environmental variables (i.e. standardized climatic and topographical variables) to capture differences in the environmental variables between plots (Table S1). Environmental distance (EnvDist) between plots was calculated as the Euclidean distance based on the first four principal components (PCs), which accounted for 91.07% of the total variance. Variance in disturbance history also calculated between plots (AgeDist). To estimate the relative importance of EnvDist, GeoDist and AgeDist on beta diversity metrics, we used the multiple linear regression model to assess the correlation between beta diversity and distance matrix 17,46 . We calculated the standardized regression coefficients using the "lm.beta" function in the R package QuantPsyc from the linear regression model.
Second, distance-based redundancy analysis (dbRDA) was conducted to assess the explanatory power of environmental variables, spatial variables on beta diversity of all forest plots as a whole as well as disturbed forests and undisturbed forests separately 47 . Seven spatial variables were obtained using the principal coordinates of neighboring matrices (PCNM) analysis 48 . Environmental variables included climate and topographic factors. Due to strong collinearity among some environmental factors, we removed environmental factors (e.g. ELE) that were primarily correlated with other factors (e.g., MAT) ( Table S3). As a result, the selected environmental factors were MAT, AP, PET, ASP and SLOPE. We performed forward selection method ("forward.sel" function in the R package packfor 0.0-8) to select environmental, and PCNM variables that had significantly effects on beta diversity. The disturbance history was also considered into forward selection as environmental variable for all forest plots as whole. The significant explanatory environmental and PCNM variables retained in dbRDA were used to partition beta diversity into four fractions: variation explained by space only, variation explained by environment only, variation explained jointly by space and environment, and variation explained by neither.
We tested for differences in beta diversity between disturbed forests and undisturbed forests using a nonparametric analysis of Wilcoxon signed-ranks test. Mantel test was used to detect the correlation among beta diversity metrics. All analyses were conducted using R software (version 3.2.2, R Foundation for Statistical Computing, Vienna, Austria) with the packages of vegan, ecodist and picante.
Because gymnosperms have on average much longer evolutionary histories and thus longer branch lengths in a phylogeny than do angiosperms, including both gymnosperms and angiosperms in a phylogenetic analysis may obscure overall phylogenetic patterns 32 even though gymnosperms accounted for only 3% of the species in our data set. Accordingly, we conducted two sets of analyses: one including both gymnosperms and angiosperms, the other including only angiosperms.

Data availability.
The data that support the findings of this study are available on request from the corresponding author (M.Y.). The data that we have used in this paper were collected through of joint projects of several research groups, and is the first time to be used for publishing in an international scientific journal. According to the agreement among the research groups who were involved in data collection, no research group can release the data until the projects finished. Thus, the agreement does not allow us to make the data to the public with our paper.

Results
Regardless of whether an analysis included both gymnosperms and angiosperms or only angiosperms, multiple linear regression showed that AgeDist and GeoDist had no influence on beta diversity of all forests, EnvDist significantly correlated with TBD and PBDt metrics (P < 0.05 in both cases, except for the Jaccard metrics), but no variables were significantly correlated with PBDb metrics (Table 1; Table S4).
MAT and AP were the significant environmental variables ( Table 2; Table S5), and together explained 7-27% of the variance in TBD and PBDt (Figs 3 and S1), whereas the pure spatial variable explained less than 1% of the variance in TBD and PBDt of all forests combined (Figs 3 and S1). The relative importance of environmental and spatial variables differed between disturbed and undisturbed forests ( Table 2; Table S5). Bray-Curtis, PhyloSor and Dnn of disturbed forests were significantly affected by MAT (P < 0.05), but those of the undisturbed forests were determined by AP (Dnn: R adj = 0.375, P < 0.05), topographical variable of SLOPE (Jaccard: R adj = 0.094, P < 0.05) and spatial variable (e.g. Bray-Curtis: R adj = 0.32, P < 0.05) ( Table 2; Table S5). None of the environmental and spatial variables significantly explained variance in Rao's H, and AP explained very little variance in Dpw ( Table 2; Table S5). In general, the results from the analysis including both gymnosperms and angiosperms are consistent with those including only angiosperms (compare Table 2 with Table S5, and Fig. 3 with Figure S1).
Of the six indices of TBD and PBD, only PhyloSor and Dpw marginally significantly differed between disturbed and undisturbed forests (Wilcoxon test: P = 0.06; Fig. 4c,e) when both gymnosperms and angiosperms were considered (Fig. 4). Specifically, Dpw for disturbed forests was higher than that for undisturbed forests (Fig. 4e). In contrast, the PhyloSor was lower for disturbed forests (Fig. 4c). However, when only angiosperms were considered, there were no significant differences in TBD and PBD between disturbed and undisturbed forests except for Dpw (Wilcoxon test: P < 0.05; Fig. S2). Dpw was significantly higher for undisturbed forests than for disturbed forests (Fig. S2).
We found that the beta diversity metrics that were calculated using both gymnosperms and angiosperms were significantly correlated with the beta diversity metrics that were calculated using only angiosperms (mantel test: r ≥ 0.615, P < 0.01) except for Dpw (r = − 0.164, P > 0.1). TBD metrics were significantly correlated with PBDt metrics (r ≥ 0.757, P < 0.001) ( Table 3; Table S6). The Rao's H metric was correlated with the metrics of Jaccard (r = 0.557, P < 0.01), Dnn (r = 0.552, P < 0.05) and Dpw (r = 0.610, P < 0.05); and the Dpw metric was only correlated with Rao's H metric when both gymnosperms and angiosperms were considered (Table 3).  Table 1. Results of multiple linear regression models. The significance for standardized regression coefficients and adjusted R-square (R 2 adj ) were calculated in models with environmental distance (EnvDist), geographical distance (GeoDist) and the variance of forests disturbance history (AgeDist) as explanatory variables and beta diversity matrices as the response variable. Both gymnosperms and angiosperms were included in the beta diversity matrices. Significance level: ** P ≤ 0.01, * 0.01 < P ≤ 0.05.

Discussion
The drivers of taxonomic and phylogenetic beta diversity include both dispersal limitation and environmental filtering. Understanding which drivers are more important in regulating patterns of beta diversity can provide insights into the mechanisms underlying community assembly 3,14,21,22 . Most studies on comparison of taxonomic and phylogenetic beta diversity and their relationships with environmental (i.e., niche process) and spatial (i.e., neutral process) variables for local forest communities were conducted in temperate and tropical forests. To our knowledge, the present study is the first to simultaneously compare taxonomic and phylogenetic beta diversity for forest communities at a regional extent in a subtropical region, to compare beta diversity between disturbed and undisturbed forests using a variety of beta diversity indices quantifying both shallow and deep histories of evolution, and to relate beta diversity with environmental and spatial distances. We have found some patterns, which we discuss below.
Contributions of space and environment to beta diversity. The main aim of our study was to assess which processes structure the beta diversity in subtropical forests. The result of multiple regression analysis   shows that environmental distance was significantly correlated with TBD and PBDt metrics (Table 1), whereas the geographical distance and disturbance history have no influence on TBD and PBDs ( Table 1), suggesting that environmental filtering plays a key role in structuring taxonomic and closely related species assembly in the subtropical forests. Generally, closely related species are likely to share more similar trait values than that expected under a Brownian motion model of trait evolution 16,49 . Our finding that environmental distance between forests have played a more important role in driving beta diversity than geographic distance is contrary to that of Saito, et al. 50 , who found that geographical distance, rather than environmental distance, better explained patterns of phylogenetic beta diversity in Neotropical stream meta-communities. However, our finding is consistent with the finding of Hardy et al. 11 who found that spatial distance does not have a significant effect on phylogenetic turnover of trees in tropical forests when accounting for rainfall. These studies suggest that dispersal limitation does not contribute to taxonomic and phylogenetic turnover in some study systems at a regional scale 11 . The importance of environment in driving beta diversity among the subtropical forests studied is also consistent with previous studies which found that environmental filtering has played a strong role on selecting plant traits in subtropical forests [51][52][53] . Furthermore, if distance generating geographic isolation is a strong force in causing the regional pattern of species distribution, patterns at shallow levels in the phylogeny (i.e., PBDt) might be influenced by spatial distance; alternatively, if geographic isolation causing dispersal limitation is a strong force in forming the regional species pool, patterns of PBD emphasizing deep nodes in the phylogeny (i.e., PBDb) might be influenced by spatial distance 14 . The result found in this study that spatial distance had no power in explaining PBD (Table 1; Fig. 3) suggests that dispersal limitation might not be an important factor in determining the spatial distribution of the study species in the present study. It is also possible that the relatively small sample size (i.e., 10 sites) of our study has a limited power in detecting the effect of spatial distance on beta diversity, although similar sample sizes have been used in previous studies (e.g., 12 sites were used in study of Culmsee and Leuschner 54 ).
The present study showed that disturbance history has no power in explaining variation in TBD and PBD, suggesting that local scale processes, such as succession and dispersal, might not be important in determining the taxonomic and phylogenetic turnover in community composition among regional forests, although they might be important in determining the taxonomic and phylogenetic community structures within these forests 24 .
Previous studies found that environmental variables related to temperature and/or precipitation drive community assembly of temperate forests 12,22,55 . While dispersal limitation and/or environmental filtering can both be the key determinants of beta diversity in tropical forests 10,42 , in subtropical forests, previous studies conducted at local scales found that space and environment were both related to beta diversity and their relative importance depended on spatial scale or sampling size 53,56 . According to these findings, environmental filtering seems to become more important from tropical to temperate areas, which might be due to geographic variation in temperature and precipitation 3,57 . Similar to the findings in temperate forests 22 , our study indicated that MAT and AP were the environmental variables that explained the most variation of TBD and PBDt (Table 2), which supports the notion that temperature and precipitation limit the distribution of tree species in eastern China 58 .
Moreover, there is still a large amount of the variation in TBD and PBDs unexplained in this study. We suspect this may be due to a combination of local stochastic processes (e.g., ecological drift) underlying community assembly and/or unmeasured environmental and spatial variables 48,56 . Especially, soil characteristic may be important in determining species composition in subtropical forests 53 and should be considered in future studies.
Beta diversity in disturbed and undisturbed forests. The second major goal of the present study was to compare beta diversity in disturbed and undisturbed forests and to determine the relative importance of environmental and spatial variables on beta diversity. Our results indicate that taxonomic beta diversity between disturbed and undisturbed forests did not significantly differ regardless of which beta diversity index was used and whether gymnosperms were considered. However, when phylogenetic beta diversity was quantified using PhyloSor and Dpw, phylogenetic beta diversity differed substantially between disturbed forests and undisturbed forests when both gymnosperms and angiosperms were considered. Feng et al. 24 found that past tree harvesting could affect subtropical forest structure by promoting the establishment of certain light-demanding and pioneer species. Our results support the notion that environmental filtering is the main process in driving the beta diversity patterns of the studied forests. We might find that the species composition will be more similar among disturbed forests, as the habitat environment become more homogeneous after clear-cut. Tree establishment will generally shift from shade-intolerant, short-lived pioneer trees and shrubs in disturbed forests to shade-tolerant, long-lived pioneers and wind-pollinated species in undisturbed forests 59 . The relative importance of environmental and spatial variables will change in disturbed and undisturbed forests when habitat heterogeneity in these  forests change. As we showed in this study, TBD and PBDt of disturbed forests were only explained by MAT, while in undisturbed forests they were determined by AP, slope and spatial variable ( Table 2). This may explain why PhyloSor for gymnosperms and angiosperms combined and Dpw for only angiosperms in disturbed forests were lower than those in undisturbed forests (Figs 4 and S2). However, Dpw for disturbed forests was higher than that for undisturbed forests when both gymnosperms and angiosperms were considered (Fig. 4e). This may be because gymnosperms such as Pinus massoniana are pioneer species and are often dominant trees in disturbed forests. Because gymnosperms have on average much longer branch lengths in a phylogeny than do angiosperms, phylogenetic turnover at more basal nodes will be higher than that at more terminal nodes.
Correlations among beta diversity metrics. The beta diversity metrics that reflect different evolutionary depths have different sensitivities to detect the beta diversity along environmental and spatial gradients.
In this study, we found that PBDt metrics (PhyloSor and Dnn) were significantly correlated with TBD metrics (Bray-Curtis and Jaccard) ( Table 3). This finding is consistent with those of previous studies 15,16,42 . The high correlation between PBDt and TBD may suggest that more recent radiations exist in the study region, and PBDt and TBD will be similarly tippy if taxonomic species' relationships are represented as a polytomous star-like phylogenetic tree 16 . However, because PBDb is sensitive to turnover deeper in the phylogeny 15,42 and many dominant tree species belong to the same genera (e.g., Ilex, Quercus and Prunus) (Fig. 2) 36 , the PBDb may have less power to detect the effect of environmental and spatial variables in this region. In addition, the Dpw that was calculated using both gymnosperms and angiosperms was not correlated with Dpw for which only angiosperms were considered; as a result, whether gymnosperms are included in an analysis would substantially change the result of the analysis, as shown in the present study. Overall, we investigated the drivers of taxonomic and phylogenetic beta diversity in subtropical forests in eastern China. We found that factors of environmental filtering, mainly temperature and precipitation, played a significant role in determining TBD and PBDt of the studied forest communities. In addition, we found that the climate variable of temperature played an important role in disturbed forests, whereas the climate, topographical and spatial variables all played an important role in driving beta diversity in undisturbed forests. TBD and PBDt were slightly higher in undisturbed forests than in disturbed forests. Our results would deepen our understanding of the processes that underlie the taxonomic and phylogenetic beta diversity in subtropical forests, which is vital to understanding the relative importance of niche and neutral process in driving community assembly of forests from tropical toward temperate regions.