Geographic patterns and environmental correlates of taxonomic and phylogenetic beta diversity for large-scale angiosperm assemblages in China

A full understanding of the origin and maintenance of β -diversity patterns in a region requires exploring the relationships of both taxonomic and phylogenetic β -diversity (TBD and PBD, respectively), and their respective turnover and nestedness components, with geographic and environmental distances. Here, we simultaneously investigated all these aspects of β -diversity for angiosperms in China. Specifically, we evaluated the relative importance of environmental filtering vs dispersal limitation processes in shaping β -diversity patterns. We found that TBD and PBD as quantified using a moving window approach decreased towards higher latitudes across the whole of China, and their turnover components were correlated with latitude more strongly than their nestedness components. When quantifying β -diversity as pairwise distances, geographic and climatic distances across China together explained 60 and 53% of the variation in TBD and PBD, respectively. After the variation in β -diversity explained by climatic distance was accounted for, geographic distance independently explained about 23 and 12% of the variation in TBD and PBD, respectively, across China. Overall, our results suggest that environmental filtering based on climatic toler-ance conserved across lineages is the main force shaping β -diversity patterns for angiosperms in China.


Introduction
Understanding the mechanisms that shape geographic variation in species diversity remains a challenging theme in ecology and biogeography (Mittelbach et al. 2007, Svenning et al. 2011. Gradients of species diversity are not only associated with variation in species number at local (i.e. α-diversity) and regional (i.e. γ-diversity) scales but also with variation in community composition among sites (i.e. β-diversity) (Gaston et al. 2007, Baselga 2010. Large-scale β-diversity patterns can be explained by geographic distance and environmental difference between sites (Harrison et al. 1992, Qian and Ricklefs 2007, Soininen et al. 2007, which reflect dispersal limitation and contemporary environmental filtering, respectively (Nekola and White 1999, Condit et al. 2002, Qian and Ricklefs 2007, Svenning et al. 2011. Because β-diversity provides a key link for understanding the relationships between species and their environment (König et al. 2017), understanding geographic patterns of β-diversity and the relationship between β-diversity and environment across a region is essential for understanding mechanisms generating variation in species diversity across the region.
β-diversity is commonly measured with indices representing pairwise compositional dissimilarity between sites, such as the Jaccard and Sørensen dissimilarity indices (Baselga 2010). These indices measure total β-diversity between sites, which results from two components: turnover and nestedness. The turnover of species between sites due to species replacement is primarily caused by dispersal limitations and niche constraints of the species in the sites; the nestedness of species between sites represents the degree to which the species in the site with fewer species is a subset of the species in the site with more species (Baselga 2010).
Most areas with a warm climate today showed stable climatic conditions during the Quaternary, which preserved species from extinction and provided more opportunities for speciation, and consequently generated more small-ranged species than areas with cold climate (Jansson and Dynesius 2002). As a result, the contribution of species turnover to β-diversity is expected to be higher in areas with a warm climate (Dobrovolski et al. 2012). In contrast, the nestedness component represents species loss or gain. Moving along a temperature gradient from warm to cold, such as from the equator to the Arctic, increasing numbers of species would have gone extinct during glaciation cycles and more species in modern communities are recent recolonizers from neighbouring areas after glaciations. This would result in stronger β-diversity patterns caused by nestedness in cold and harsh environments (Baselga 2010, Dobrovolski et al. 2012. As a result, the contribution of species nestedness to overall beta diversity is expected to increase with decreasing temperature (Dobrovolski et al. 2012). A better understanding of the origin and maintenance of geographic patterns of β-diversity in a region, therefore, requires investigating the relative importance of each of the two components to total β-diversity , Leprieur et al. 2012, Peixoto et al. 2017.
While a taxonomic measure of β-diversity (taxonomic β-diversity, hereafter TBD) and its turnover and nestedness components can effectively provide information about the degrees of overlap and distinctness of species between sites, it does not allow associating species distribution to their evolutionary history (Graham and Fine 2008). Incorporating information regarding the phylogenetic relatedness between species in the study of β-diversity may substantially advance our understanding of the ecological and evolutionary mechanisms structuring communities (Ives and Helmus 2010). Analogous to taxonomic β-diversity, which measures change in species composition across space, phylogenetic β-diversity measures the extent to which assemblages differ in terms of the evolutionary relationships of its members (Graham and Fine 2008). Phylogenetic β-diversity (PBD hereafter) can also be partitioned into two components accounting for 'true' phylogenetic turnover and phylogenetic diversity gradients (Leprieur et al. 2012). Examining TBD and PBD simultaneously can provide insights into both contemporary ecological and historical evolutionary mechanisms shaping variation in species diversity and composition among local and regional assemblages (Graham and Fine 2008, Peixoto et al. 2017, Cássia-Silva et al. 2019.
A full understanding of the origin and maintenance of β-diversity patterns in a region requires the understanding of the relationships of both TBD and PBD and their respective turnover and nestedness components with geographic and environmental distances, and the relationships between the turnover and nestedness components of β-diversity. However, there are, to our knowledge, few regions in the world for which all these aspects of β-diversity have been examined for a particular group of organisms.
In the present study, we simultaneously explore these aspects of β-diversity using two complementary approaches, namely the neighborhood approach (Peixoto et al. 2017, Pinto-Ledezma et al. 2018) and pairwise approach Ricklefs 2007, Qian et al. 2013), and we test relevant hypotheses for large-scale angiosperm (flowering plant) assemblages (in 10 000 km 2 ) across China, which has a broad geographic extent and a great variation in climate. Specifically, we test the 'dispersal limitation' and 'environmental filtering' hypotheses, which are not mutually exclusive.
H1: if dispersal limitation processes (e.g. due to the emergence of historical barriers or differences in species dispersal ability) have left a strong imprint on current patterns of angiosperm β-diversity in China, we would expect a minor role of contemporary environmental conditions in explaining spatial patterns of both TBD and PBD. According to this hypothesis (Nekola andWhite 1999, Hubbell 2001), assemblages that experience similar climatic conditions in geographically distant areas would exhibit distinct species compositions, as a result of dispersal limitation.
H2: if environmental filtering based on climatic tolerances conserved across lineages (i.e. phylogenetic niche conservatism, see Qian and Ricklefs 2016) is the main force structuring regional assemblages of angiosperms in China, we would expect a strong spatial congruence between TBD and PBD along climatic gradients, with greater TBD than PBD (Jin et al. 2015). Indeed, if species have tracked their ancestral climatic niches, we would expect to observe predominantly the spatial turnover of closely related species rather than distantly related species (Graham et al. 2009, Qian et al. 2013. Then, we can expect decreasing phylogenetic turnover towards cold temperatures at higher latitudes . Assemblages in regions with more stressful climates (colder and drier) should be composed of more closely related species (i.e. phylogenetically clustered, see Qian et al. 2019, Mienna et al. 2020 for example) because most tropical clades failed to disperse into extratropical regions, lacking adaptations to survive winter temperatures below freezing (Ricklefs 2006). Consequently, the turnover component of PBD should contribute less to overall PBD in cold environments (i.e. assemblages are more phylogenetically similar) than the phylogenetic diversity component (see also Peixoto et al. 2017). Similarly, we can expect that pairwise dissimilarity values for both TBD and PBD would be mainly explained by climatic distances, as shown for angiosperms in North America (Qian et al. 2013).

Species distribution data
The data set with species distributions in 100 × 100 km grid cells reported by Lu et al. (2018) was used in our study as a primary data source. We supplemented the data with two additional sources: 1) county-level species distributions derived from specimen records with the National Specimen Information Infrastructure (NSII; <www.nsii.org.cn>) and the Global Biodiversity Information Facility (GBIF; <www. gbif.org>), which were described in Qian et al. (2018), and 2) plant checklists for local floras (e.g. counties, nature reserves and national forest parks) in the literature (Zhao et al. 2006, Qian and Chen 2016, Qian et al. 2017. Botanical nomenclature of species was standardized according to The Plant List (ver. 1.1, <www.theplantlist.org>). Infraspecific taxa were combined with their respective species. Non-native distributions were excluded based on Wu et al. (1994Wu et al. ( -2013 and provincial floras (Liu et al. 2007). Species distributions obtained from NSII and GBIF and from plant checklists for local floras were converted to distributions in 100 × 100 km grid cells as in Lu et al. (2018). The species assemblage of each grid cell was defined as all species with distributions falling within the grid cell. We combined species distributions from the above-mentioned three data sources. Only those grid cells which have 75% or more of their area on land were included in this study.
Because there are few inventory-based species checklists of vascular plants at the spatial scale of 10 000 km 2 in China, we are not able to assess the degree of completeness of the species lists used in this study for 100 × 100 km grid cells. Considering that the species lists of 100 × 100 km grid cells reported by Lu et al. (2018) and by others for Chinese plants have been commonly used in macroecological studies (Feng et al. 2016, Lu et al. 2018, Xu et al. 2018, Zhao et al. 2018, Ye et al. 2019, the use of occurrence data compiled for Chinese plants in 100 × 100 km grid cells in the present study is consistent with the current literature. Furthermore, because the plant data used in this study is an updated version of the data reported by Lu et al. (2018), with 10.5 millions of specimen records from NSII and GBIF (Qian et al. 2018) and over 200 local and regional floras (Qian et al. 2017(Qian et al. , 2018 having been used to update Lu et al.'s data, we believe the species lists used in this study are more complete than those used in most, if not all, of previous macroecological studies on Chinese plants at the spatial scale of 10 000 km 2 .

Phylogeny reconstruction
There were 27 905 species and 2822 genera of angiosperms in our data set. We used the R package V.PhyloMaker (build. nodes.1; Jin and Qian 2019) to generate a phylogeny for the Chinese angiosperm species. Specifically, V.PhyloMaker uses an updated and expanded version of the dated megaphylogeny GBOTB reported by Smith and Brown (2018) as a backbone to generate phylogenies. All families of Chinese angiosperms were resolved in this megaphylogeny. Of the 2822 genera, 2468 were included in the megaphylogeny. Of those genera that occurred in our data set but were not included in the megaphylogeny, 214 have closely related (often sister) genera in the megaphylogeny based on Zanne et al. (2014) and Lu et al. (2018). We treated them as sisters to their respective closely related genera in the megaphylogeny. As a result, phylogenetic relationships of 2682 (95%) of the genera in our data set were resolved. At the species level, 11 132 species of the angiosperm species in our data set were included in the megaphylogeny. For the genera and species in our data set that are absent from the megaphylogeny, we added them to their respective genera (in the case of species) and families (in the case of genera) using the Phylomatic and BLADJ approaches (Webb et al. 2008) implemented in V.PhyloMaker (scenario 3; Jin and Qian 2019). Specifically, V.PhyloMaker sets branch lengths of added taxa of a family by placing the nodes evenly between dated nodes and tips within the family, and adds a missing genus at the crown node of its family and a missing species at the crown node of its genus. This approach of adding taxa to a megaphylogeny is broadly used in community phylogenetics (Webb andDonoghue 2005, Shiono et al. 2018). Recent evaluations on phylogenies generated in this or a similar way have shown that such phylogenies are appropriate for macroecological and biogeographic studies (Qian and Jin 2016, Qian and Zhang 2016, Li et al. 2019). Finally, we pruned the megaphylogeny to retain only the 27 905 species occurring in our data set.

Metrics of taxonomic and phylogenetic β-diversity
We used the Sørensen dissimilarity index (β sor ) to measure both taxonomic and phylogenetic β-diversity. Following Baselga (2010), we partitioned β sor into two components, β sim and β nes . β sim , which is the Simpson dissimilarity index, quantifies β-diversity due to turnover (or replacement) of species between sites (i.e. β-diversity resulting from true substitution of species between sites); thus β sim quantifies true turnover of species without the influence of difference in species richness between sites. β nes quantifies β-diversity resulting from nestedness of species between sites. Nestedness of species between sites occurs when the species of the site with the lower number of species are a subset of the species of the richer site (Baselga 2010). β sim and β nes result from two antithetic processes, i.e. species replacement and species loss (or gain), respectively (Baselga 2010). β sor and its two components are defined as follows (Baselga 2010): where, when taxonomic β-diversity is concerned, a is the number of species shared by the two sites, b is the number of species unique to one site and c is the number of species unique to the other site (Baselga 2010). When applied to phylogenetic β-diversity, shared and unique species are replaced with shared and unique branch lengths, respectively (Leprieur et al. 2012). We indicate β-diversity and its two components, respectively, as β sor.tax , β sim.tax and β nes.tax for taxonomic β-diversity, and as β sor.phy , β sim.phy and β nes.phy for phylogenetic β-diversity. Both β sor.tax and β sor.phy vary from 0 (all species or branch lengths shared by the two sites) to 1 (no species or branch lengths shared by the two sites). For both taxonomic and phylogenetic β-diversity, we calculated the ratio of β nes to β sor (i.e. β ratio = β nes /β sor ; Dobrovolski et al. 2012; see also Peixoto et al. 2017, Pinto-Ledezma et al. 2018. β ratio indicates the relative influence of turnover and nestedness on total β-diversity (Dobrovolski et al. 2012, Peixoto et al. 2017. A value of β ratio smaller than 0.5 indicates that β-diversity is determined mainly by turnover whereas a value of β ratio greater than 0.5 indicates that nestedness is the more important component than turnover in driving taxonomic β-diversity (Dobrovolski et al. 2012). For phylogenetic relatedness, a value of β ratio smaller than 0.5 indicates that the turnover of lineages is more important than differences in phylogenetic diversity in shaping phylogenetic β-diversity. We differentiated β ratio for taxonomic and phylogenetic β-diversity as β ratio.tax and β ratio.phy , respectively.
As noted by Jin et al. (2015) and Cadotte and Davies (2016, p. 161), taxonomic β-diversity will commonly be higher than phylogenetic β-diversity. This pattern has been observed in previous studies (Meynard et al. 2011, Qian et al. 2013, Guevara et al. 2016, McFadden et al. 2019. Assessing patterns of deviations of phylogenetic β-diversity from taxonomic β-diversity across different geographic areas can provide insights into the mechanisms that generate phylogenetic β-diversity across a broad region (Peixoto et al. 2017, Pinto-Ledezma et al. 2018. We therefore calculated phylogenetic β-diversity deviation (β dev ) from taxonomic β-diversity using the following formula: β dev = (β sor.tax − β sor.phy )/β sor.tax . A value of β dev would reflect how much more novel information is captured by phylogenetic β-diversity over taxonomic β-diversity. When values of β dev are compared across a geographic region, β dev identifies places where turnover of phylogenetic lineages is high or low with respect to turnover of species. When β sor. phy is the same as β sor.tax (i.e. in the case of a star phylogeny), β dev would be zero, indicating that phylogenetic β-diversity is maximized. In contrast, if β sor.phy approaches zero, β dev would approach 1, indicating that there is little or no lineage exchange between sites. Accordingly, a smaller value of β dev indicates a greater turnover of phylogenetic lineages and thus a greater effect of phylogenetic β-diversity between sites. β dev reflects the degree of lineage exchanges between areas over time (Peixoto et al. 2017, Pinto-Ledezma et al. 2018.

Climate data
Mean annual temperature and annual precipitation are two key climate variables that determine plant distributions at a broad spatial extent (Whittaker 1975). Stressful climates (including extreme cold and drought) and seasonal variation in climate (temperature and precipitation seasonality) are also considered as constraints of the distributions of species (Kamilar et al. 2015, Weigelt et al. 2015, Patrick and Stevens 2016. Accordingly, we used six variables to characterize climate in each grid cell: mean annual temperature, annual precipitation, minimum temperature of the coldest month, precipitation during the driest month, temperature seasonality and precipitation seasonality. We obtained values for these variables from the WorldClim database (Hijmans et al. 2005; <www.worldclim.org>: variables bio1, bio12, bio6, bio14, bio4 and bio15, respectively). The mean value of each of the six climate variables was calculated for each site using 30-arcsecond resolution data.
In addition to using the original climate variables in the analyses of β-diversity quantified based on the neighbourhood approach, we also used principal components analysis (PCA) to generate a synthetic climatic variable (i.e. the first principal component, PC1) that represents the six original climatic variables for each grid cell and used the synthetic variable in some analyses. PC1 accounted for 69% of the variance in the six original climatic variables.

Beta diversity analyses
We applied two complementary approaches to investigate patterns of β-diversity across China: a neighborhood approach (also called moving window approach), which is for investigating β-diversity between sites separated by short geographic distances (< 200 km), and a pairwise comparison approach, which is for investigating β-diversity between sites separated by varying geographic distances (up to 3700 km).
With the neighborhood approach, we first calculated the average β-diversity between a focal grid cell and each of its first-order (i.e. immediate) neighboring grid cells (n = 8 at maximum; Supplementary material Appendix 1); second, we calculated average β-diversity between the focal grid cell and each of its second-order neighboring grid cells (n = 16 at maximum; Supplementary material Appendix 1). We excluded those focal grid cells for which the number of actual cell pairs was equal to or smaller than 50% of the number of overall cell pairs (i.e. the number of first-order neighboring grid cells is 4 or fewer, or the number of second-order neighboring grid cells is 8 or fewer; Supplementary material Appendix 1). We visually compared maps of β-diversity between the two frameworks for each of β-diversity metrics and correlated the two sets of β-diversity values. We found that geographic patterns of β-diversity were highly congruent between the two frameworks (Spearman's rank correlation coefficient r s = 0.737 and 0.760, respectively, for β sor.tax and β sor.phy ). Accordingly, we only map and analyze β-diversity between a focal cell and its first-order neighboring cells. This spatial framework has been commonly used in recent studies (Peixoto et al. 2017, Pinto-Ledezma et al. 2018, McFadden et al. 2019. With the pairwise approach, we explored how pairwise compositional dissimilarity (or pairwise β-diversity) was related to environmental and geographical distances between sites. For each metric of β-diversity, the 903 grid cells in our data set would result in 407 253 pairs of grid cells. Due to limitations on workspace with some of the software packages that we used, we were not able to compute β-diversity metrics for such a great number of pairwise grid cells based on such a large data matrix of species by grid cell (i.e. 27 905 species by 903 grid cells). As a result, we randomly selected one third of the 903 grid cells, and computed each metric of β-diversity for the selected data, which resulted in 45 150 pairs of sites for each metric of β-diversity. The climatic distance between each pair of sites was measured as the Euclidean distances of the six climate variables, as stated above. The geographic distance between each pair of sites was measured as the Euclidean distance between mid-points of the sites. Following previous studies, we log 10 -transformed geographic distance, because distance-decay patterns tend to be exponential (Nekola and White 1999, Tuomisto et al. 2003, Rodrigues and Diniz-Filho 2017. Difference in elevational range between each pair of sites was used as a measure of difference in topographical heterogeneity between the sites (Rodrigues and Diniz-Filho 2017).

Statistical analyses
When considering the neighbourhood approach, we performed simple and multiple regression analyses to explore the relationships of β-diversity with latitude and the considered environmental variables. To account for spatial autocorrelation, we used spatial autoregressive (SAR) models (Kissling and Carl 2008).
When considering the pairwise approach, we applied simple and multiple regressions on distance matrices, with statistical significance being determined based on permutations (Legendre et al. 1994, Legendre andLegendre 2012), to assess the relationships between each β-diversity metric and 1) geographical distance, 2) difference in contemporary climatic conditions (climatic distance) and 3) difference in elevational ranges. We also conducted regression on distance matrices (Legendre et al. 1994) from which we extracted the standardized regression coefficient (i.e. equivalent to the Pearson's correlation coefficient if the two matrices are normalized) and the coefficient of determination (R 2 ).
Because climatic distance between cells is expected to decline with increasing geographical distance, disentangling the relative roles of environmental filtering and dispersal limitation in shaping variation in β-diversity is not straightforward (Gilbert and Lechowicz, 2004). We therefore conducted multiple regression on distance matrices (MRM, Legendre et al. 1994) including geographic distance, climatic distance and elevation range difference, and calculated the standardized partial regression coefficients that enabled us to compare the per-unit effect of a given explanatory variable on pairwise β-diversity, while controlling for the effect of the other variables (Baselga and Leprieur 2015). We also reported the coefficients of multiple determination (R 2 ) for these complete models. To overcome the problem of lack of independence between paired sites, the significance of the regression coefficients and the coefficients of multiple determination (R 2 ) was assessed using a permutation test (n = 999).
To further distinguish between the effect of variables describing dispersal limitation (geographic distance) versus environmental conditions (climatic distance), we applied a variance partitioning approach (Legendre and Legendre 2012) based on MRM to separate the total explained variance (R 2 ) into three parts: explained uniquely by geographic distance, explained uniquely by climate distance, and explained by geographically structured climate variation (i.e. variance jointly explained by geographic and climate distances). Following previous studies (Qian et al. 2013, Mazel et al. 2017), we considered the variance explained by climate distance uniquely and geographically structured climate variation as the variance explained by climate distance, because climatic variables are strongly structured geographically at large scales and thus effects of geographically structured climate variation should be considered as indirect effects of climate (Mazel et al. 2017).

Geographic patterns of β-diversity across China
Geographic patterns of total TBD and PBD (β sor.tax and β sor.phy , respectively) and their components of turnover (β sim.tax and β sim.phy ) and nestedness (β nes.tax and β nes.phy ) for angiosperm assemblages among neighboring grid cells were highly congruent across China (Fig. 1). The correlation was strong in all the three cases (Spearman's rank correlation: r s = 0.963 for β sor.tax versus β sor.phy , 0.975 for β sim.tax versus β sim.phy and 0.959 for β nes.tax versus β nes.phy (n = 853 in all cases). The turnover component was negatively and weakly correlated with the nestedness component for both TBD and PBD (Spearman's rank correlation: r s = −0.227 and −0.219, respectively).
Across China, both total TBD and PBD were found to be negatively and weakly associated with latitude (SAR: R 2 = 0.067 and 0.101 for β sor.tax and β sor.phy , respectively, Supplementary material Appendix 2). When the turnover components of TBD and PBD were considered, the negative association of β-diversity with latitude became stronger (SAR: R 2 = 0.188 and 0.250 for β sor.tax and β sor.phy , respectively; Supplementary material Appendix 2).
Geographical patterns of β-diversity across China appear to be more apparent by distinguishing between the turnover and nestedness components of TBD and PBD (β sor.tax and β sor.phy ). In particular, the southeastern half of China appears to be in contrast with the northwestern half of China with regard to both turnover and nestedness components. In general, turnover was high and nestedness was low in the southeastern half of China; in contrast, turnover was low and nestedness was high in the northwestern half of China. Across overall China, the turnover component was as high as the nestedness for TBD (mean ± SD, 0.124 ± 0.065 for β sim.tax versus 0.121 ± 0.064 for β nes.tax , n = 853, paired t-test, p = 0.397), but the former was lower than the latter for PBD (0.078 ± 0.042 for β sim.phy versus 0.094 ± 0.053 for β nes.phy , n = 853, paired t-test, p < 0.001).

Geographic patterns of β ratio and β dev across China
Geographic patterns of β ratio , which quantifies how much the nestedness components contribute to total β-diversity, were highly congruent between taxonomic and phylogenetic β-diversity (Spearman's rank correlation: r s = 0.972). In general, the ratio was higher in the northwestern half of China, which is with a cold and dry climate, than in the southeast, which is under warm and moist climate (Fig. 2). The vast majority of the grid cells in the southeastern half of China had values of β ratio.tax and β ratio.phy smaller than 0.5 whereas the vast majority of the grid cells in the southeastern half of China had values of β ratio.tax and β ratio.phy greater than 0.5. Figure 1. Patterns of taxonomic and phylogenetic β-diversity (β sor.tax and β sor.phy ) and their components of turnover (β sim.tax and β sim.phy ) and nestedness (β nes.tax and β nes.phy ) for angiosperms in China. A value of β-diversity for a cell represents the average of β-diversity values resulting from all pairwise comparisons between the cell and each of its first-order neighboring cells (i.e. using the neighborhood approach; see Supplementary material Appendix 1 for details). β dev was positive in all angiosperm assemblages across China (Fig. 2), confirming that taxonomic β-diversity was higher than phylogenetic β-diversity in every assemblage. β dev tended to increase with latitude, reflecting that phylogenetic β-diversity decreased faster with increasing latitude than expected with respect to taxonomic β-diversity. This indicates a greater turnover of phylogenetic lineages and thus a greater effect of phylogenetic β-diversity between sites in lower latitudes. Regions with the least differences between taxonomic and phylogenetic β-diversity were primarily located in the Himalayas and the Hengduan Mountains in southwest China. These areas are mountainous regions with great environmental differences produced by elevational differences. We observed that these regions have high turnover of both taxonomic and phylogenetic β-diversity (Fig. 1).

Relationship between β-diversity and climate
Total TBD and PBD were found to be weakly but significantly associated with the climatic variables (SAR: R 2 < 0.05 for both β sor.tax and β sor.phy when considering either the mean annual temperature or the synthetic climate PC1 variable, Supplementary material Appendix 2). Stronger relationships with the mean annual temperature were found when analyzing the turnover component of both TBD and PBD (R 2 = 0.132 and 0.123 for β sim.tax and β sim.phy , respectively, see Supplementary material Appendix 2).

Relating pairwise β-diversity to geographic and climatic distances and elevational range difference
TBD and PBD were strongly to moderately associated with geographic distance (R 2 = 0.504 and 0.369, respectively, for β sor.tax and β sor.phy ), and moderately associated with climatic distance in both cases (R 2 = 0.364 and 0.407, respectively, for β sor.tax and β sor.phy , p < 0.001 in all cases based on 999 permutations, Table 1). Similar results were found when only considering the turnover components of TBD and PBD (Table 1). In contrast, the nestedness components of TBD and PBD were found to be weakly associated with geographic and climatic distances (Table 1). Multiple regression models testing the combined effects of geographic and environmental distance-based variables on TBD and PBD and their turnover components accounted for more than 48% of the variation in dependent variables in all cases (Supplementary material Appendix 3). Standardized partial regression coefficients indicated that climatic distance, after having controlled for the influence of both geographical distance and elevation range difference, was a strong predictor of TBD and PBD and their turnover components (Supplementary material Appendix 3). Figure 2. Patterns of the relative importance of the nestedness component of β-diversity (β ratio.tax for taxonomic β-diversity; β ratio.phy for phylogenetic β-diversity) and the deviation of phylogenetic β-diversity with respect to taxonomic β-diversity (β dev ) for angiosperms in China. A value for a cell represents the average of β-diversity values resulting from all pairwise comparisons between the cell and each of its firstorder adjoining neighboring cells (i.e. using the neighborhood approach; see Supplementary material Appendix 1 for details). Table 1. Results of simple linear regression based on distance matrices relating each β-diversity metric to geographic distance (Geog.dist), climatic distance (Clim.dist) and elevational range difference (Elev.range.diff). We report the standardized regression coefficients (Coeff.) and the coefficient of determination (R 2 ). Significance was tested using a permutation test (999 permutations): ** p < 0.001, * p < 0.05, ns p > 0.05.

Partitioning of variance in β-diversity explained by geographic and climatic distances
Geographic and climatic distances together explained 59.8 and 52.9% of the variance in taxonomic and phylogenetic β-diversity, respectively (Table 2, Supplementary material Appendix 4). They explained much more variance in the turnover component than in the nestedness component of β-diversity (64.4% versus 12.5% for taxonomic β-diversity; 48.4% versus 4.0% for phylogenetic β-diversity). Geographic distance uniquely explained about 23.4 and 12.2% of the variance in taxonomic and phylogenetic β-diversity (Table 2), and explained much less variance in the nestedness component than in the turnover component of β-diversity (5.7% versus 26.2% for taxonomic β-diversity; 0.3% versus 14.0% for phylogenetic β-diversity; Table 2). When differences in elevational ranges were added to the multiple regression on distance matrices, including both geographic and climatic distances (Supplementary material Appendix 5), they only explained additional < 1% of the variance in β-diversity. This suggests that differences in elevational ranges between sampling sites did not play an additional role in driving taxonomic and phylogenetic β-diversity after geographic and climatic distances were accounted for.

Discussion
Understanding why species composition varies along spatial and environmental gradients is a long-standing issue in the fields of ecology and biogeography (Gaston et al. 2007). Numerous hypotheses have been proposed to explain species richness gradients at large spatial scales, but much less attention has been given to the determinants of spatial patterns of β-diversity, especially with consideration of both the taxonomic and phylogenetic dimensions of β-diversity. This study analyzed a comprehensive data set to address various issues on β-diversity for angiosperms in China. We discuss some of the important findings of this study below.

Taxonomic versus phylogenetic β-diversity
Geographic variation in taxonomic versus phylogenetic β-diversity (i.e. β sor.tax versus β sor.phy ) across China is highly consistent (Fig. 1). This is true also for their turnover and nestedness components (i.e. β sim.tax versus β sim.phy ; β nes.tax versus β nes.phy ; Fig. 1). As found for regional assemblages of angiosperm in North American (Qian et al. 2013), taxonomic β-diversity is higher than phylogenetic β-diversity for each angiosperm assemblage across China (β dev > 0). This suggests that β-diversity patterns predominantly consist in the exchange of phylogenetically close species across space (i.e. species belonging to the same lineage) rather than distant ones (i.e. species belonging to different lineages). We also found that PBD and TBD, and more particularly their turnover components, were significantly related to the gradient of mean annual temperature across China (Supplementary material Appendix 2). Altogether, these results suggest that species sorting related to differential temperature tolerances conserved across lineages (i.e. phylogenetic niche conservatism) might have played a major role in structuring regional assemblages of angiosperm in China, as also showed in previous regional-scale studies (Qian et al. 2013, McFadden et al. 2019). However, a large proportion of variation in PBD and TBD remained unexplained (Supplementary material Appendix 2), as also found by McFadden et al. (2019) for the New World's plant assemblages, and consequently this suggests that other factors may be important in shaping geographic patterns of β-diversity such as the strength of biotic interactions (Ohlmann et al. 2018) and/or dispersal limitation processes (Qian 2009).
We observed that β dev is much lower in the southern part of Tibet (including the northern part of the Himalayas) and the broad region of the Hengduan Mountains (Fig. 2c). This indicates that phylogenetic β-diversity is higher in these regions than elsewhere in China when taxonomic β-diversity is accounted for. Mountain ranges in the Himalayas and Hengduan Mountains are high and steep, with many peaks being over 7000 m above sea level. Mountain building in these regions has generated steep environmental gradients along elevational gradients and high environmental heterogeneity, which would result in both high taxonomic β-diversity and high phylogenetic β-diversity, as we observed in this study (Fig. 1). Strong environmental gradients across elevation, along which members from distinct lineages are segregated among different areas, would produce higher than expected phylogenetic β-diversity (Graham et al. 2009).

Opposing patterns for the turnover and nestedness of β-diversity
Patterns of variation in the turnover and nestedness across China tend to be opposite for both taxonomic and phylogenetic β-diversity. Specifically, turnover tends to be higher in the southeastern half than in the northwestern half of China whereas nestedness tends to be higher in the northwestern half than in the southeastern half of China (Fig. 1). This finding is generally consistent with those of previous studies (Baselga 2010, Dobrovolski et al. 2012, Soininen et al. 2018, which found that turnover tends to decreases with latitude whereas nestedness tends to increase with latitude. These authors attributed the greater role of the nestedness component of β-diversity in determining β-diversity at higher latitudes to the coverage of ice-sheets at high latitudes during the Pleistocene. Specifically, they reasoned that the extinction and slow recolonization of species in areas covered by ices during the Last Glacial Maximum would result in the higher prevalence of the nestedness component of β-diversity whereas the extinction events in areas without ice coverage during the Last Glacial Maximum should have been less severe, which allows the maintenance of a greater number of species and the emergence of endemism, which in turn would result in turnover of species. However, because eastern China was not covered by ice-sheet during the Last Glacial Maximum, the glaciation-hypothesis proposed by previous authors may not explain the pattern of increasing nestedness (and β ratio.tax and β ratio.phy ) with increasing latitude in China. In our view, the pattern of increasing nestedness with latitude might not be directly caused by coverage of ice-sheets during the last glaciations and recolonization after the glaciations, although stronger decrease in temperature at higher latitudes during Pleistocene glaciations might have played a role. We believe that the latitudinal nestedness gradient is caused by the well-known pattern of decreasing temperature with latitude (Morueta-Holme et al. 2013). While still being debated (Igea and Tanentzap 2020), speciation rate is expected to be high at low latitudes and species ranges tend to be small, leading to different small-ranged endemic species present in different areas. This would in turn lead to high turnover of species in regions at low latitudes. In contrast, at high latitudes, speciation rate is thought to be low and extinction rate is high, the latter of which would lead to extinction of small-ranged species during Pleistocene glaciations, which would decrease turnover but increase nestedness of species among sites. Furthermore, the gradient of increasing temperature seasonality with latitude would result in the gradient of increasing species range with latitude (i.e. Rapoport's rule; Stevens 1989). Areas with more large-ranged species might show more importance of nestedness, compared to turnover (Tomašových et al. 2016). Thus, the emergence of the opposite latitudinal pattern between the turnover and nestedness components of β-diversity might have been caused by 1) species sorting processes along the well-known latitudinal gradients of mean annual temperature and temperature seasonality and 2) by the tendency of lineages to maintain similar environmental limits to their range extents over time (phylogenetic niche conservatism, see H2).

Relations of β-diversity with geographic and climatic distances
Our study found that pairwise β-diversity is associated with climatic distance more strongly than with geographic distance, and more particularly when considering PBD (Table 1, 2; Supplementary material Appendix 3). For example, after accounting for the effect of climatic distance on β-diversity, geographic distance alone explained little variation in both PBD and its turnover component across China (less than 15%, Table 2). In addition, the per unit effect of climatic distance on both TBD and PBD and their turnover components remained relatively high after having accounted for the effects of geographic distance and elevational range difference (Supplementary material Appendix 3). Overall, our results suggest that patterns of taxonomic and phylogenetic β-diversity across China were primarily driven by species sorting along climatic gradients, as proposed by the environmental filtering hypothesis (H2), rather than by dispersal limitation (H1), a result that is consistent with previous regional studies on angiosperms (Qian et al. 2013, Jin et al. 2015, Pinto-Ledezma et al. 2018.

Conclusions
The fact that PBD was found to be lower than TBD along climatic gradients suggests that species sorting related to differential climatic tolerances conserved across lineages (i.e. phylogenetic niche conservatism) might be a major force structuring regional assemblages of angiosperms in China (see also Qian et al. 2013 for North America). When analyzing TBD and PBD as pairwise distances, we found that regional assemblages of angiosperms in China were primarily explained by geographically structured climatic conditions, which provides a further support for the 'environmental filtering' hypothesis (H2). However, this does not imply that dispersal limitation processes (H1) were unimportant in shaping modern assemblages of angiosperms in China. Because angiosperm lineages may differ in their dispersal abilities (Qian et al. 2013), future research may compare the decay in compositional similarity with geographic distance and, independently, with climatic distance, for diverse lineages with contrasting dispersal abilities and climatic niches (Gómez-Rodríguez and Baselga 2018). In addition, geographic distance alone does not reflect the presence or absence of geographic barriers (i.e. spatially close locations may be separated by strong dispersal barriers such as mountain ranges) and the development of metrics of 'barrier distance' is a promising avenue of research to unravel the role of dispersal limitation in shaping large-scale β-diversity patterns (Eiserhardt et al. 2013).

Data availability statement
All data used in this study have been published and are accessible to readers from the cited sources.