1. Introduction

Biogeographical regionalization aims to classify global or regional biota into meaningful geographical units for analysis (Mackey et al., 2007). Thus, it can provide an indispensable background for research on biodiversity and conservation (Crisp et al., 2009). Takhtajan (1978) divided the global flora into six kingdoms, 34 regions and 147 provinces, which are widely used in historical and ecological biogeography (Linder et al., 2012; Meyer et al., 2012a, 2012b).

The increasing availability of species distribution data, combined with the tree of life reconstructions, have set new dimensions for studying biogeographical regionalization (Kreft and Jetz, 2010; Holt et al., 2013; Li et al., 2015, 2018; Li and Sun, 2017). Recently, important progress has been made in biogeographical regionalization using both qualitative approaches with endemic taxa (Wu et al., 2010) and quantitative approach with distribution and evolutionary data (Zhang et al., 2016a, 2016b; Ye et al., 2019). Compared with qualitative approaches, regionalization based on quantitative approaches are much more objective and reproducible. To arrive at a reasonable quantitative biogeographical regionalization scheme, it is critical to choose an appropriate dissimilarity index that measures pairwise distances between grid cells assemblages. Several metrics have been used in biogeographical regionalization, including Euclidean distance (Zhang et al., 2016a), Jaccard (Linder et al., 2005), the mean nearest taxon distance (MNTD) (Webb, 2000), Sørensen (Xie et al., 2004) and the turnover component of the Sørensen dissimilarity index (β_sim; Baselga, 2012) (see Table 1). Some of these metrics are independent of differences in species richness among grid cells (Baselga, 2012), while some are strongly affected by differences in species richness (Lennon et al., 2001). For instance, Slik et al. (2018) suggested that the MNTD is easily affected by differences in taxonomic richness, and cannot give the correct measure of taxon turnover between two assemblages. Therefore, in studies of biogeographical regionalization, researchers must use those metrics (e.g., β_sim and β_jtu) that are least affected by variation in richness (Leprieur and Oikonomou, 2014), especially for those biotas with a large richness gradient (Mouillot et al., 2013), such as the Chinese flora (Lu et al., 2018).

Because some previous studies used unreasonable metrics, in this study, we first explored the effect of metric choice on biogeographical regionalization, then compared the regionalizations of the Chinese flora resulting from different metrics by using species distribution records and a genus level phylogenetic tree of the Chinese flora. We then analyzed the boundaries, the characteristics, and origins of the chosen regions of the Chinese flora. We anticipate this work will help identify reasonable metrics that can be used to propose a phylogenetic floristic regionalization of China.

2. Materials and methods 2.1. Distribution data and phylogenetic tree

We used the same data set as Ye et al. (2019).

2.2. Distance metrics and clustering analysis

To compare the effect of richness-dependent and richnessindependent metrics on biogeographical regionalization, we used MNTD and pβ_sim to measure phylogenetic turnover (pβ_sim), respectively. MNTD is affected by differences in species richness whereas pβ_sim is independent of differences in species richness among grid cells (Slik et al., 2018). In this study, MNTD was obtained by the 'comdistnt' function of the 'picante' package (Kembel et al., 2010) in R (R Core Team, 2019) and pβ_sim was calculated according to Ye et al. (2019).

To identify spatial clusters, we compared the output from eight hierarchical clustering algorithms based on the MNTD and pβ_sim pairwise distance matrices:(1) single linkage (SL), (2) complete linkage (CL), (3) unweighted pair-group method using arithmetic averages (UPGMA), (4) unweighted pair-group method using centroids (UPGMC), (5) weighted pair-group method using arithmetic averages (WPGMA), (6) weighted pair-group method using centroids (WPGMC), and (7) and (8), two variants of Ward's minimum variance (ward.D and ward.D2). We used the Kelley-Gardner-Sutcliffe penalty function (Kelley et al., 1996) to get a reasonable number of clusters for the phytogeographical regions. The clustering algorithms and reasonable number of clusters selected for this study helped to identify appropriate phytogeographical regions, and show the differences of the floristic regions obtained by using MNTD and pβ_sim respectively.

2.3. Characterized clade and indicator genus analyses

MNTD is sensitive to differences in taxon richness; therefore, we adopted floristic regionalization using pβ_sim. To identify clades that characterize a phytogeographical region (Ye et al., 2019) or contribute significantly to phylogenetic structure, we used the procedure nodesig in Phylocom 4.2 (Webb et al., 2008). The overall phylogenetic tree nodesig tests were used to further compare the genus pool of each phytogeographical region, and examine whether a node in the phylogeny has significantly more (sig. more) or significantly less (sig. less) descendent taxa in a sample than a null model predicts. The null model is a random draw of n taxa from the phylogeny, where n is the number of taxa in the sample (Webb et al., 2008).

We used the Dufrêne- Legendre indicator species analysis (Dufrêne and Legendre, 1997) to characterize the genus composition of each identified phytogeographical region. The indicator value (IndVal) was defined as the product of two quantities, specificity (A) and fidelity (B). A represents the probability that a surveyed site belongs to the target site group on the basis of a target genus that has been found (calculated as the relative frequency of the occurrence of each genus in the target site group, divided by the sum of relative frequencies over all groups). B is the probability of finding the target genus in a site that belongs to the site group (i.e., the relative frequency of a genus inside the target site group; De Câceres and Legendre, 2009). The significance of the IndVal was assessed by a permutation test (9, 999 permutations). A genus that had significant IndVals (P < 0.001) was ranked as a suitable indicator genus of its relevant phytogeographical region. The analysis was conducted using the 'indicspecies' package (De Câceres and Legendre, 2009).

3. Results 3.1. Floristic regionalization using MNTD and pβ_sim

The five clustering algorithms varied in performance for the MNTD pairwise distance analysis (Table 2). Ward's method was the best performing clustering algorithm (cophenetic coefficient = 0.99), followed by CL (complete linkage) (cophenetic coefficient = 0.97), while UPGMA (unweighted pair-group method using arithmetic averages) performed worst (cophenetic coefficient = 0.94). Therefore, we chose Ward's method. According to the Kelley-Gardner-Sutcliffe penalty function, eight floristic regions belonging to two groups were identified in China, with southern regions 6 and 8 forming one group (South Cluster), and the other six regions forming the other group (North Cluster) (Fig. 1a). The geographical dividing line between North and South Cluster in eastern China was roughly consistent with the Qinling Mountain-Huaihe River Line. The North Cluster further splits into two subclusters, with regions 1, 3 and 7 forming one subcluster and regions 2, 4 and 5 the other subcluster.

For the floristic regionalization using pβ_sim, five floristic regions were identified (Fig. 2; Fig. 4 in Ye et al., 2019).

3.2. Characterized clades and indicator genera

The nodesig algorithm identified differently characterized clades for different phytogeographical regions. The Tethyan region had the largest number (301 clades) of sig. more clades, which included 99.12% of genera (904 genera) in this region, but the mean number of genera contributed by each clade was the least (3.00 genera). In contrast, the East Asiatic region had relatively few sig. more clades (57 clades), which only included 42.88% of genera (988 genera) in this region, but each clade contributed the largest number (25.61 genera) of genera on average. Details of sig. more/sig. less nodes and the descendant taxa are shown in Appendix S1:Fig. S1.

Among the 2, 591 genera analyzed, 34.3% (889 genera) had a significant indicator value (Table 3). This was a substantial number of genera that qualified as indicator taxa and indicated that different phytogeographical regions had significantly different floras. The Paleotropic region had the largest proportion (40.6%, 361 genera) of total indicator genera (889 genera), whereas the Holarctic region had the least proportion of indicator genera (5.8%, 52 genera). The indicator genera identified for each phytogeographical region are provided in Appendix S1:Tables S1-S5.

4. Discussion 4.1. Choice of the right index for bioregionalization studies

Beta diversity is composed of two different components:spatial species turnover and nestedness of assemblages, which result from two antithetic processes, namely species replacement and species loss, respectively (Lennon et al., 2001; Koleff et al., 2003; Baselga et al., 2007; Tuomisto, 2010). Using European longhorn beetles as a case, Baselga (2010) exemplified the relevance of disentangling spatial turnover and nestedness patterns. Baselga and Orme (2012) published an R package "betapart" to compute these dissimilarity measures. In the studies of biogeographical regionalization (bioregionalization), we should use the turnover component to measure composition dissimilarity between grid cells to avoid the impact of diversity differences (Leprieur and Oikonomou, 2014). Kreft and Jetz (2010) recommended the beta-sim index (β_sim) and it has been applied in most of the subsequent bioregionalization studies (Table 1; Linder et al., 2012; lez-Orozco et al., 2013; Holt et al., 2013). However, some studies still used inappropriate or incorrect indices (e.g., Bray-Curtis, mean nearest phylogenetic distance), which have proved to be influenced by diversity differences (Slik et al., 2018). In this study, the floristic regions obtained by MNTD distance matrices were also influenced by diversity differences between grid cells, because the Chinese flora has a conspicuous east-west gradient in plant diversity (Lu et al., 2018).

4.2. Comparison between two regionalizations

The cluster analysis based on MNTD distance matrices presented a picture of eight regions for the Chinese flora, whereas there were five regions based on pβ_sim identified by or found by Ye et al. (2019). The floristic regions obtained by MNTD distance matrices were similar with the patterns of generic richness in China (Extended Data Figure 8a in Lu et al., 2018). Southern regions 6 and 8 have the highest richness, regions 1, 3 and 7 have the moderate richness, and regions 2, 4 and 5 have the lowest richness. As a result, regionalization based on MNTD distance matrices was strongly influenced by the nestedness component of the metric (Baselga, 2012), which was also found in a study of Indo-Pacific coral reef fishes (Mouillot et al., 2013). Furthermore, the pβ_sim was demonstrated to be able to reflect the true phylogenetic turnover for the Chinese angiosperm flora (Fig. 2b in Ye et al., 2019) and the five floristic regions-the Paleotropic, Holarctic, East Asiatic, Tethyan, and QTP regions-which clearly reflected the evolutionary history of the flora (Ye et al., 2019).

In our study, the characterized clades identified by the nodesig algorithm for each phytogeographical region differed greatly. For example, the clades in magnoliids (Myristiaceae, Annonaceae, Magnoliaceae, Calycanthaceae, Lauraceae, Hernandiaceae, Piperaceae, and Aristolochiaceae) were significantly associated with the Paleotropic region, while clades in Asteraceae occurred less often than expected in this region but were overabundant in the Holarctic region (Appendix S1:Fig. S2). We found that in the East Asiatic region, the sig. more nodes were deeper nodes with the most descendant genera on average in the five regions. This suggests the East Asiatic region has older lineages in China, while the flora of the Tethyan region with shallower sig. more nodes is younger (Table 3).

We concluded the discussion with a description of the five regions (Paleotropic, Holarctic, East Asiatic, Tethyan, and QTP regions) identified by our analysis (Fig. 2).

4.3. The five identified regions based on the pβ_sim metric

The Paleotropic region consists of tropical evergreen broadleaved forest vegetation in the southeastern coast of China and Hainan Island. Its northern boundary is roughly the Tropic of Cancer (Fig. 2). Within the five regions, this region has the least amount of areas and lowest elevation, the highest mean annual temperature (MAT) and mean annual precipitation (MAP), and the richest indicator clades and genera (Table 1 in Ye et al., 2019, Table S1 in Appendix S1). The Paleotropic region has the most evolutionarily distinct flora in China (Table 1 in Ye et al., 2019) and is characterized by the presence of thermophilic Malesian taxa (e.g., Bruguiera, Chrysophyllum, Endospermum, Stylidium, Tetracera), which were immigrated from southern and southeastern Asia (Axelrod et al., 1998). These tropical and subtropical plant clades and genera have mostly not migrated to the subtropical and temperate areas in central and northern China and alpine areas in southwestern China, implying the prevalence of Niche Conservatism Hypothesis (Wiens and Donoghue, 2004), thus resulting in the main split of the Chinese flora between the southern Paleotropic and northern Chinese floras, as shown in our clustering dendrogram (Fig. 1b). The flora of the Paleotropic region is the oldest and represents a floristic museum (Lu et al., 2018).

The Holarctic region consists of Inner Mongolia steppe in the west and deciduous broad-leaved forest vegetation in Northeast and North China. Within North China, it is roughly bounded by the Helan Mountains in the west and the Taihang Mountains in the east (Fig. 2). The mean elevation is approximately 700 m, MAT is 4.1℃, and MAP is nearly 500 mm (Table 1 in Ye et al., 2019). Flora of this region includes 6, 878 species that belong to 1, 254 genera and 170 families. It represents the Boreal-Tertiary flora in China that preserved Northern Hemisphere temperate taxa with origins in the late Cretaceous and early Paleogene, such as Acer, Alnus, Corylus, Eleutherococcus, Rhododendron, and Tilia (Sun, 2002b). This region has the fewest indicator genera among the five regions that are mainly composed of north temperate taxa (e.g., Filifolium, Hypochaeris, Saposhnikovia, Scabiosa, and Trientalis) (Table 1 in Ye et al., 2019, Table S2 in Appendix S1).

The East Asiatic region includes deciduous broad-leaved forest vegetation in the east of the Taihang Mountains (Fig. 2), subtropical evergreen broad-leaved forest areas in most of southern China, and seasonal rain forest vegetation in southern Yunnan. Its western border extends along the western margin of the Sichuan Basin to the eastern and southern margin of the Yunnan Plateau. Flora of the East Asiatic region is closely related to that of the Holarctic region. The mean elevation is approximately 600 m a.s.l.; MAT is 16℃; MAP is nearly 1, 200 mm (Table 1 in Ye et al., 2019). The East Asiatic region is characterized by the presence of genera that are endemic to East Asia (e.g., Akebia, Emmenopterys, Pinellia) and genera with the distribution pattern of the East Asiae-North American disjunction patterns (e.g., Liriodendron, Shortia) (Table S3 in Appendix S1). The East Asiatic region hosts the highest plant species richness (20, 317 species) among the five regions (Table 3). First, as a stable and ancient area, the East Asiatic region preserved the flora from the early Cretaceous to early Paleogene (Wu et al., 2010). Second, with the uplift of the QTP, and intensification of both the Pacific and Indian Ocean monsoons in this area, new habitats have appeared during the Neogene, which might have promoted the rapid radiation of some plant taxa, such as Begoniaceae, Gesneriaceae, Orchidaceae, and Rubiaceae (Favre et al., 2015). Third, the East Asiatic region preserved Boreotropical plant taxa that were once common in the Northern Hemisphere but extinct from other areas during Quaternary glaciation, such as Cyclocarya, Davidia, Emmenopterys, and Eucommia (Wolfe, 1975; Wen, 1999; López-Pujol and Ren, 2010). Therefore, the East Asiatic region harbors both early and recently diverged taxa in China and East Asia.

The Tethyan region mainly includes the desert areas in North-west China (Fig. 2). The mean elevation is approximately 1, 800 m; MAT is 5.6℃; and MAP is 111.2 mm. The Tethyan region has the lowest MAP of the five regions. As a temperate desert flora, plant diversity in the Tethyan region is low, with only 5, 412 species that belong to 912 genera and 128 families (Table 1 in Ye et al., 2019). This region has 147 indicator genera, which are mainly temperate desert xerophytic shrubs and herbs (e.g., Ammopiptanthus, Glaucium, Lagochilus, Potaninia, Tugarinovia, and Zygophyllum) (Table S4 in Appendix S1). However, the Tethyan region had a high evolutionary distinctiveness in comparison with other regions (Table 1 in Ye et al., 2019). After the Oligocene, the inland Tethys Sea retreated and was accompanied by an increasingly dry climate, which resulted in diversification of the xerophytic Tethyan flora (Palamarev, 1989). Consequently, this region is the diversity center for Amaranthaceae (Chenopodioideae), Apiaceae, Brassicaceae, and Lamiaceae (Sun, 2013). Moreover, the Tethyan region was also the most important glacial refugium for Tertiary relics throughout Quaternary glacial cycles (Qiu et al., 2011); for example, Tetraena mongolica is considered a relic of the Tethys sea tropical phytogeographical region (Wu et al., 2010).

The QTP region is bordered to the south by the Himalayan Mountains, to the east by the Hengduan Mountains, to the north by the Kunlun Mountains and Qilian Mountains, and to the west by the Karakorum Mountains (Fig. 2). The QTP region is the largest region and has the highest average elevation (4, 248 m a.s.l.) and lowest MAT of the five regions. The QTP region is characterized by the presence of alpine genera (e.g., Cyananthus, Diapensia, Lancea, Meconopsis, Morina, and Solms-laubachia) (Table S5 in Appendix S1). This region possesses one of the world's richest flora, with 16, 220 species that belong to 1, 897 genera and 205 families (Table 1 in Ye et al., 2019), and numerous Chinese and local endemic species (Huang et al., 2016; Zhang et al., 2016a). Because of the uplift of the QTP, this region is also an important origin and diversification center of north temperate and alpine floras, such as the genera Aconitum, Corydalis, Gentiana, Rhododendron, Pedicularis, and Primula (Sun, 2002a; Zhang et al., 2016a).

5. Conclusions

We demonstrated that the choice of metrics could have a dramatic impact on the biogeographical regionalization. Those regionalizations based on metrics that are sensitive to taxa richness can reflect the pattern of turnover and nestedness together, but cannot reflect the true pattern of taxonomic and phylogenetic turnover (Baselga, 2012). The regionalization based on the MNTD metric is greatly influenced by the richness differences in the Chinese flora. Our phylogeny-based regionalization using beta-sim metric provides a foundation for future biogeographical and biodiversity studies and can also be used to develop strategies to conserve taxa and phytogeographical regions with different origins and evolutionary histories more effectively in China.

Declaration of Competing Interest

The authors declare that there is no conflict of interest regarding this manuscript. All the authors agreed to submit this manuscript.

Acknowledgements

This work was supported by grants from the National Natural Science Foundation of China (NNSF 31800178, 31590822, 31500179 and 31461123001), the Strategic Priority Research Program of Chinese Academy of Sciences (XDA19050103 and XDB31000000), the National Key Basic Research Program of China (2014CB954101).

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2019.12.003.