Using species range data, phylogenies, and fossil records of thousands of extinct and extant terrestrial vertebrates, He et al. quantify temporal changes in beta dissimilarities among zoogeographical regions surrounding the Tibetan Plateau during the Cenozoic era. This study shows how geological events underpin the evolutionary history of zoogeographical regions.
The Tibetan Plateau (TP) and surrounding regions have one of the most complex biotas on Earth. However, the evolutionary history of these regions in deep time is poorly understood. Here, we quantify the temporal changes in beta dissimilarities among zoogeographical regions during the Cenozoic using 4,966 extant terrestrial vertebrates and 1,278 extinct mammal genera. We identify ten present-day zoogeographical regions and find that they underwent a striking change over time. Specifically, the fauna on the TP was close to the Oriental realm in deep time but became more similar to the Palearctic realms more recently. The present-day zoogeographical regions generally emerged during the Miocene/Pliocene boundary (ca. 5 Ma). These results indicate that geological events such as the Indo-Asian Collision, the TP uplift, and the aridification of the Asian interior underpinned the evolutionary history of the zoogeographical regions surrounding the TP over different time periods.
IntroductionThe Tibetan Plateau (TP) uplift was one of the most important geological events in the Cenozoic era (~65 Ma–present1,2). It substantially modified the topography3 and atmospheric circulation4 of Asia (Fig. 1) and resulted in one of the most complex biotas on Earth5. Eight major zoogeographical regions were recently identified surrounding the TP5,6, namely, the Mongolian Plateau, Central Asia, North Asia, West Asia, South Asia, Southeast Asia, South China and North China (Fig. 1). However, it is unknown how these present-day zoogeographical regions evolved over geological time, even though this information is crucial for understanding the origin and evolution of life in Asia.Fig. 1Palaeogeographical processes in the Tibetan Plateau and surrounding regions during the Cenozoic.Palaeotopography was derived from a palaeo-digital elevation model (palaeoDEM, 1° × 1° resolution) developed by Scotese & Wright52. Changes in the atmosphere–ocean climate system were compiled from the data of Sun & Wang50. Present-day zoogeographical regions were adapted from Kreft & Jetz6 and Holt et al.5. CA Central Asia, MP Mongolian Plateau, NA North Asia, NC North China, SA South Asia, SC South China, SEA Southeast Asia, TP Tibetan Plateau, WA West Asia.There is growing evidence that the present-day zoogeographical regions surrounding the TP are the products of geological processes and past climatic changes7–9. A common hypothesis is that the TP uplift created species-dispersal barriers during the Cenozoic10,11, and subsequent climatic changes in Asia increased environmental heterogeneity (Fig. 1, see refs. 12,13); also both events caused the geographical isolation of resident lineages and facilitated the differentiation of zoogeographical regions14. Other biogeographical analyses, however, revealed that historical events, such as the Eocene Indo-Asian Collision15,16, the intercontinental biotic exchange between Eurasia and North America17 and the Pleistocene glaciation cycle18, might have expanded species’ ranges, promoted dispersal and attenuated the faunistic dissimilarities among regions. These processes have been proved to facilitate dispersal and vicariance for many lineages that might increase or decrease the number of taxa common to different regions19, and ultimately alter their pairwise faunistic relationships over time. However, it remains uncertain how these processes shaped the evolution and emergence of the present-day zoogeographical regions surrounding the TP.Recent phylogeographical analyses on the TP have associated biogeographical and evolutionary lineage relationships with specific geological events and periods14. Unfortunately, most available empirical studies have relied upon the interpretation of single-taxon analyses. They have inferred the influences of geological processes and climatic shifts on genus- or species-level distributions of specific taxa10,11,13. However, responses to common geological events might greatly vary among lineages owing to their biological and ecological differences20. These differences would result in incongruent biogeographical patterns across different taxonomic lineages over space or time7,14. Furthermore, present-day zoogeographical regions were structured by a combination of multiple speciations, extinctions and dispersal processes at several time periods21,22. Thus, biogeographical meta-analysis16,22 and community-level analyses5,23, which integrate individual taxon histories into shared biotic area histories, were more promising to clarify the processes shaping biogeographical regions over time19,24.To date, two primary analyses of community-level data have been used to reconstruct the evolutionary history of zoogeographical regions. One tracked temporal changes in beta diversity between extant communities over a phylogenetic timescale23,24, and another compared compositional dissimilarities among fossil assemblages over geological time9,25. However, both methods have their pros and cons26. The former method provides a finer resolution regarding the spatial and temporal changes in communities23, but always fails to deal with past extinctions and distribution changes27 and, therefore, provides only indirect evidence. Although the inclusion of ancestral range reconstruction in quantifying phylogenetic dissimilarity can improve estimates of evolutionary history, it is difficult to incorporate extinct lineages into the analysis (ref. 24, but see refs. 28,29). In contrast, palaeontological materials can provide a direct record of past changes in communities, but they always suffer from incomplete preservation30, which possibly conceals some important signs of biogeographical events31. Nevertheless, despite the limitations of the respective methods, phylogenetic and palaeontological analyses can usefully complement each other in biogeographical studies26, and their correspondence has begun to emerge in large-scale biogeographical contexts32,33.In this study, we reconstructed the evolutionary history of the zoogeographical regions surrounding the TP using 4966 extant terrestrial vertebrates along a phylogenetic timescale and 1278 extinct mammal genera over geological time. By comparing analyses implemented over phylogenetic and geological timescales, we aimed to explore the timeframe within which the present-day zoogeographical regions evolved during the Cenozoic era and estimate the time when the present-day spatial structure of the zoogeographical regions emerged. To reconstruct historical changes in the zoogeographical regions, we quantified the phylogenetic beta dissimilarity using extant species along the phylogenetic timescale. For extinct lineages, we calculated beta dissimilarity based on mammal fossil assemblages over geological time. We assessed the changes in assignments and topologies of hierarchical clustering dendrograms and relative positions in ordinations based on beta dissimilarity at different phylogenetic depths and geological periods. Finally, we explored the relationships between the evolutionary history of the zoogeographical regions within the context of geological and climatic events. Our study reveals that the zoogeographical regions underwent a striking change during the Cenozoic era, and broadly emerged in the Miocene/Pliocene boundary (ca. 5 Ma) owing to a series of geological events such as the Indo-Asian Collision, the TP uplift and the aridification of the Asian interior.ResultsZoogeographical regions over phylogenetic timeTen present-day zoogeographical regions were delineated by unweighted pair-group method using arithmetic average (UPGMA) clustering based on pβsim matrix, namely the Tibetan Plateau, Mongolian Plateau, Central Asia, North Asia, West Asia, South Asia, Indochina, Malay Peninsula, South China and North China & Korea (Fig. 2). At a phylogenetic depth of 60 Ma, nine clustered zoogeographical regions were identified (Fig. 2). They roughly corresponded to the Palearctic realm, North China & Korea, West Asia, South China, Hengduan Mountains, the northern part of Indochina, the southern part of Indochina, Malay Peninsula and South Asia (Fig. 2a). At the phylogenetic depth of 40 Ma, the northern and southern parts of Indochina and the Malay Peninsula were merged into a united region. The southern part of the TP was separated from South Asia. The boundary between South China and North China & Korea moved from ca. 30 °N to 40 °N. At the phylogenetic depth of 20 Ma, the most striking change was that Central Asia combined with the Mongolian Plateau and emerged as an independent region. Then, Central Asia was separated from the Mongolian Plateau at a phylogenetic depth of 10 Ma. The spatial structures of the present-day zoogeographical regions are broadly similar to those at a phylogenetic depth of 5 Ma, when the whole TP was identified as an independent region. These findings are broadly consistent with the analyses performed on the whole-region species list for all terrestrial vertebrates (Supplementary Fig. 1). Notably, these biogeographical processes varied among taxonomic groups, as reflected by the spatial patterns (Supplementary Fig. 2) and the Mantel correlation test (Supplementary Fig. 3). The pβsim structure between mammals and all terrestrial vertebrates showed the highest correlation in the present day (Supplementary Fig. 2 and 3). In contrast, their correlations were gradually weaker than those between ectotherms (i.e., reptiles and amphibians) and their combined counterparts in deeper phylogenetic time bins (Supplementary Fig. 3).Fig. 2Temporal changes in the zoogeographical regions surrounding the Tibetan Plateau at successive phylogenetic depths during the Cenozoic.a Map showing zoogeographical regions based on pβsim dissimilarity between pairs of grid-based terrestrial vertebrate communities at different phylogenetic depths. The width of the boundary was scaled to the βsim dissimilarity, with thinner lines showing lower βsim dissimilarities. CA Central Asia, IC Indochina, ML Malay Peninsula, MP Mongolian Plateau, NA North Asia, NCK North China & Korea, SA South Asia, SC South China, TP Tibetan Plateau, WA West Asia. b Dendrograms plotted by the unweighted pair-group method using arithmetic average clustering. c Coordinates for non-metric multidimensional scaling ordination based on the pβsim dissimilarity matrix for grid-based terrestrial vertebrate communities.Temporal changes in the spatial structures of the zoogeographical regions were reflected in the UPGMA dendrograms and non-metric multidimensional scaling (NMDS) ordinations (Fig. 2b, c). Interestingly, the relationships between the zoogeographical regions over phylogenetic timescales were well illustrated by the Procrustes analysis (Fig. 3). For example, in the deep branches, the grid cells of North China & Korea and North Asia, and those of the Mongolian Plateau and Central Asia largely overlapped. In contrast, their differences were clearer in the shallow branches (Fig. 3). However, the relationships between West Asia and other regions underwent less change over phylogenetic timescales (Fig. 3). Interestingly, the species assemblage on the TP was similar to that of the Oriental realm in the past, but it became closer to the Palearctic realm towards the present day (arrow of the TP in Fig. 3 pointing to the Palearctic realm). This shift was also illustrated by the spatial patterns showing that the boundary between the TP and Mongolian Plateau became shallower in the present day, while the boundary between the TP and South Asia gradually strengthened (Fig. 2).Fig. 3Temporal changes in terrestrial vertebrate communities in zoogeographical regions on a phylogenetic timescale.Each gridded terrestrial vertebrate assemblage is represented by a small arrow linking the six coordinates through time (0 Ma, 5 Ma, 10 Ma, 20 Ma, 40 Ma and 60 Ma, respectively). Large arrows represent average gridded species assemblages across zoogeographical regions. Arrow colours correspond to the zoogeographical regions in the inset map. CA Central Asia, IC Indochina, ML Malay Peninsula, MP Mongolian Plateau, NA North Asia, NCK North China & Korea, SA South Asia, SC South China, TP Tibetan Plateau, WA West Asia.Zoogeographical regions over geological timeThe UPGMA clustering based on the fossil records showed that relationships between zoogeographical regions underwent a striking change through geological time (Fig. 4). In the Eocene (56–33.9 Ma), the TP was merged with South China and North China & Korea and then grouped with Central Asia and the Mongolian Plateau. South Asia emerged as the most distinct region, showing strong dissimilarities to the other regions (Fig. 4). During the Oligocene (33.9–23.0 Ma), after the initial collision between India and Eurasia, the TP was first grouped with South Asia. During the Early Miocene (23–15.9 Ma) and Mid–Late Miocene (15.9–5.3 Ma), however, the TP was first merged with the Mongolian Plateau. Then, the combination of the TP and Mongolian Plateau was grouped with the Oriental realm (i.e., South China and South Asia) in the Early Miocene, whereas they merged with the Palearctic realm (i.e., North China & Korea, North Asia and West Asia) in the Mid–Late Miocene. During the Pliocene–Pleistocene period (5.3 Ma–11.8 Ka), the division between the Palearctic and Oriental realms emerged. When we quantified βsim dissimilarity based on the extant mammal lists for the whole region, four groups of zoogeographical regions were identified, namely group 1: Central Asia + Mongolian Plateau + North Asia, group 2: South Asia + West Asia, group 3: Indochina + Malay Peninsula and group 4: North China + South China + Tibetan Plateau (Fig. 4).Fig. 4Temporal changes in mammal fossil assemblages between zoogeographical regions during the Cenozoic.a Beta dissimilarities between zoogeographical regions inferred from UPGMA clustering are indicated by lines, with thinner lines indicating higher βsim dissimilarities. Palaeomaps were derived from a palaeo-digital elevation model developed by Scotese & Wright52. b Dendrograms from the unweighted pair‐group method using arithmetic average hierarchical clustering of regional mammal lists during five time intervals and the present. c Non‐metric multidimensional scaling ordinations based on the βsim matrices of regional mammal assemblages. CA Central Asia, IC Indochina, ML Malay Peninsula, MP Mongolian Plateau, NA North Asia, NCK North China & Korea, SA South Asia, SC South China, TP Tibetan Plateau, WA West Asia.Comparison between phylogenetic and palaeontological inferencesThe faunistic relationships between zoogeographical regions based on phylogenetic information and fossil data yielded considerable differences. For instance, Central Asia combined with the Mongolian Plateau emerged as an independent region at a phylogenetic depth of 20 Ma, while this pattern was not detected in the fossil data. The present-day zoogeographical regions broadly emerged at a phylogenetic depth of 5 Ma (Fig. 2; Supplementary Fig. 1), whereas the spatial structures of the zoogeographical regions between the Pliocene–Pleistocene period and present day displayed some differences based on fossil data (Fig. 4). Nevertheless, we found some consensuses among the phylogenetic and palaeontological inferences. For example, South Asia emerged as a distinct region in the early Cenozoic (compare Figs. 2 and 4). The faunistic similarity between the TP and Oriental realms was close in the early stages after the Indo-Asian Collision, whereas the TP has become more similar to the Palearctic realms since the Early Miocene (ca. 23–15.9 Ma, Figs. 2c–4).DiscussionOur findings identified ten present-day zoogeographical regions surrounding the TP, which are broadly consistent with the regions identified by previous global regionalisation studies5,6, despite different taxonomic groups being used as inputs. This indicates that different lineages share the similar ecological and historical drivers that underpin their co-occurrence distributions8,34. Notably, our results did not recognise the Sino-Japanese realm at a higher classification level as proposed by Holt et al.5, as the present-day TP was first grouped with the Mongolian Plateau and then merged into other regions within the Palearctic realm (Fig. 2a, b). This pattern corroborated the assumption that the distinctiveness of the Sino-Japanese realm is rather weak and may be easily altered by slight changes in the methods or taxonomic groups used35. Nevertheless, as this study focused only on the TP and its surrounding regions, further studies are required to assess the validity of the Sino-Japanese realm.When we tracked the evolutionary history of the zoogeographical regions back to the early Cenozoic, South Asia (i.e., the Indian plate) was identified as an independent region based on both the phylogenetic and palaeontological results (Figs. 2 and 4). Although the phylogenetic distinctiveness of South Asia was lower than the palaeontological estimate, the independence of South Asia is still valid (Fig. 2). Tectonic studies have suggested that the Indian plate was part of Gondwana and was completely isolated from Eurasia before the Eocene continental collision (ca. 55–50 Ma36). As expected, it harbours several relict taxa closely related to the African and Madagascar lineages37,38 that represent obviously distinct fauna from that of Eurasia. Notably, our results showed that the faunistic relationships between South Asia and the other zoogeographical regions within the Eurasian plate became closer towards the present day (Figs. 2 and 4). This shift is likely to reflect the imprint of the extensive biotic interchanges between South Asia and Eurasia since the Eocene15,16. Consequently, present-day South Asia shares a large proportion of extant lineages with Eurasia, and is placed in the Oriental realm instead of the African realm at a higher classification level5,6.It has long been suggested that the uplift of the TP reconfigured the spatial structures of the biota in Asia7,14, but the faunistic relationships between the TP and its surrounding regions have not been adequately assessed. Our results revealed that the fauna on the TP was closer to the Oriental realm in deep time, but became more similar to the Palearctic realms towards the present time. This finding was independently confirmed by the pβsim dissimilarity of extant species (Figs. 2 and 3) and βsim dissimilarity of mammal fossils (Fig. 4). This trend is probably because, during the Eocene to Oligocene (ca. 56–23 Ma), the southern part of the TP began to emerge above sea level owing to the Indo-Asian Collision, whereas its northern part was still part of the Tethys Ocean1. Accordingly, some Palaeogene-aged lineages could have been expected to co-occur within both the Indian plate and TP39. However, as the uplift of the TP proceeded, the high, steep mountains in the southern TP began to act as a barrier to biotic exchanges between the TP and the Oriental realm since the Middle Miocene (ca. 15.9–11.6 Ma40). In contrast, multiple species could disperse from the Palearctic realm into the TP via the northern routes due to the moderate topographic gradient41,42. Interestingly, only the southern part of the TP was identified at the phylogenetic depth of 40–10 Ma (Fig. 2), indicating that its species assemblage was more phylogenetically distinct in deep time. Both of these patterns provided biogeographical evidence that the TP underwent a south-to-north uplift process7,43, and supported the recent palaeobotanical findings that the elevation of the southern part of the TP reached its present height at the Eocene–Oligocene boundary (ca. 34 Ma44).We also found a profound change in that Central Asia combined with the Mongolian Plateau became an independent region at the phylogenetic depth of 20 Ma (Fig. 2). This timing is approximately contemporaneous with the aridification of the Asian interior during the early Miocene (ca. 24–22 Ma45,46), which was characterised by the contraction of the wet–humid biome and expansion of the dry–cool biome47. The aridification has been proposed to have created new habitats for speciation or dispersal barriers, and resulted in extensive vicariance events according to multiple taxon-specific studies13,48. Other evidence from mammalian fossils49, palaeobotanical records50 and sedimentary data45 support such an important regional climatic change during this period and corroborate our phylogenetic results. However, the emergence of Central Asia and the Mongolian Plateau region was not observed in our fossil data (Fig. 4), which likely reflects the limitations of our palaeontological inferences (see details below).Present-day spatial structures of the zoogeographical regions broadly emerged at the phylogenetic depth of 5 Ma (Fig. 2), as evidenced by the phylogenetic analyses based on both the gridded assemblages (Fig. 2) and whole-region species lists (Supplementary Fig. 1). According to the fossil data, although the spatial patterns during the Pliocene–Pleistocene period (5.3 Ma–11.8 Ka) were rather incongruent with those of the present day (Fig. 4), the correlation between their βsim dissimilarity matrices remained significant (Mantel test, P < 0.05; Supplementary Table 1). Previous palaeontological analyses in China yielded a similar timing and showed that the spatial structure of mammalian communities originated during the Pliocene9,51. This probably happened because the emergence of modern orographic patterns52 and monsoon systems4 has formed dispersal limitations and accelerated lineage diversification. Many in situ radiation events within zoogeographical regions have been reported, including the pika Ochotona on the TP53, the gibbon Hylobates in Indochina54 and the palm squirrel Funambulus in South Asia55, as well as several examples reported for other taxa (e.g., plants56, fishes57 and birds58). Although the subsequent Pleistocene glaciation cycle was expected to influence the geographical ranges of species and biodiversity distributions59, our results detected negligible influences of this event on the spatial structure of the zoogeographical regions surrounding the TP (Fig. 2).Notably, some discrepancies in the historical changes in zoogeographical regions emerged between taxonomic groups and analytical methods. For example, the ectotherms were more important than the endotherms in structuring the zoogeographical regions at deeper phylogenetic time bins (Supplementary Figs. 2 and 3). This probably results from differences in their life history strategies and evolutionary processes. The ectotherms are generally confined to fewer zoogeographical regions due to higher environmental sensitivity and weaker dispersal capacity60. Niche shifts in endotherms are faster than those in ectotherms61, which have likely altered more pβsim patterns of endotherms towards the present62. Furthermore, the discrepancies in the phylogenetic and palaeontological inferences highlight the different outcomes of these methods. For example, the fauna of South Asia was clearly distinct from that of Eurasia during the early Cenozoic, but its phylogenetic dissimilarity was much lower than the palaeontological dissimilarity. This pattern illustrates that the past extinction of ancient endemism and geographical range shift caused by biotic interchanges might conceal considerable phylogenetic signals in deep time63, and thus potentially biases the dissimilarity estimates19,24. Since the extant species represent only a subset of lineages in the phylogenetic tree, a comprehensive picture of the early history of a biota and its temporal changes cannot be resolved by the extant species alone63. However, analyses based on the fossil record fail to detect some biogeographical changes, such as the emergence of Central Asia and the Mongolian Plateau during the Early Miocene (ca. 23–15.9 Ma). One possible reason is that we had to merge the fossil records to coarse zoogeographical regions owing to the low number of fossils, which, however, only reflect the temporal changes in pairwise faunistic relationships between regions. Alternatively, the available fossil collection inevitably suffered from uneven sampling (Supplementary Table 2; see ref. 30) and time averaging31, making it impossible to clarify the finer-scale position of the zoogeographical boundaries and to reconstruct their successive temporal changes as the phylogenetic methods do23. Nonetheless, although there are some differences in the methods and inferred patterns between phylogenetic and palaeontological estimates, the comparisons between these two methods are still informative. For example, the changes in the relationships between the TP and its surrounding regions were effectively resolved by both methods. Overall, using community-level data to reconstruct the temporal changes in biogeographical regions is still a challenging and ongoing mission. Further studies are needed to integrate the molecular phylogenies and fossil data into a combined dataset64, together with ancestral area estimates65, to enable a more comprehensive understanding of evolutionary histories of present-day biogeographical patterns.In conclusion, based on the long-term changes in beta dissimilarity inferred from the palaeontological data and phylogenetic information, this study reconstructed the evolutionary history of the zoogeographical regions surrounding the TP during the Cenozoic Era. Our study demonstrated that the faunistic relationships among these regions underwent a substantial reconfiguration during the Cenozoic as a consequence of several biogeographical events during different periods. These events included the Indo-Asian Collision, the TP uplift and the aridification of the Asian interior. The present-day zoogeographical regions surrounding the TP originated during the Miocene/Pliocene boundary (ca. 5 Ma) when the modern geographical pattern and climatic systems were established. The present study highlights the importance of comparing phylogenetic and palaeontological inferences to reconstruct the history of biogeographical regions. In this way, we may enhance our comprehension of the origin and evolution of life driven by various eco-evolutionary processes over space and time.MethodsSpecies dataWe obtained extant species distribution maps from the IUCN Red List website (http://www.iucnredlist.org) for mammals and amphibians, Birdlife International and NatureServe (http://www.birdlife.org) for birds and Roll et al.66 for reptiles. We excluded introduced, marine and domestic species. Species geographical ranges were transformed into presence and absence data in a matrix of 110-km × 110-km grid cells with the Behrmann projection. We removed grid cells with a land area <50% and species richness <5 to minimise the negative influences of the unequal sampling area and statistical uncertainty. We obtained the most comprehensive dated phylogenies available online (http://vertlife.org/phylosubsets) for each vertebrate group. For mammals, we used a phylogenetic tree (5911 species) from Upham et al.67 that used two levels of Bayesian inference (backbone relationships and species-level phylogenies) to constrain the age and topological uncertainty. For birds, we used a phylogenetic tree (9993 species) from Jetz et al.68 based on the Hackett family-level backbone. For reptiles and amphibians, we used phylogenetic trees from Tonini et al.69 and Jetz & Pyron70, comprising 9574 squamate species and 7238 amphibian species, respectively. In these phylogenies, the topology of species with molecular data was fixed, and the remaining species unsampled for DNA-sequence data were assigned randomly within their genus or higher-level groups based on morphology69,70, resulting in a distribution of 10,000 trees. We downloaded a set of posterior distributions of trees (n = 1000) online using complete lists of all available species, and obtained the maximum clade-credibility phylogenies using the ‘maxCladeCred’ function from the ‘phangorn’ package71 in R version 3.6.072. After combining the distributional and phylogenetic data, our dataset comprised a total of 4966 extant terrestrial vertebrates, including 1022 mammals, 1741 birds, 1453 reptiles and 750 amphibians (Supplementary Data 1).Fossil dataWe obtained fossil records from four databases: Institute of Vertebrate Paleontology and Paleoanthropology, Beijing (http://www.ivpp.ac.cn/), the Paleobiology Database (https://www.paleobiodb.org/), the New and Old Worlds database (https://www.helsinki.fi/science/now/) and the Fossilworks database (https://fossilworks.org/), accessed in April 2018. We focused only on the mammal fossils owing to their relatively good preservation and samples73. In addition, we used genus rather than species as the analytical unit because the fossil records at the genus level included more complete sampling and reliable identification64. We standardised the taxonomy according to the Paleobiology Database and excluded taxa unidentifiable at the genus level. We removed the duplicated records and merged spatially closest collection localities by combining those within 0.1 latitude and longitude9. Our final fossil dataset consisted of 5880 fossil occurrences of 170 families and 1278 genera (Supplementary Data 2). We reconstructed the fossil records from present-day coordinates back to their palaeo-position based on the mean age of the fossil in a temporal range using the ‘reconstruct’ function in the ‘chronosphere’ package74.Delineation of present-day zoogeographical regionsTo delineate present-day zoogeographical regions, we used Simpson’s phylogenetic beta diversity (pβsim) to generate pairwise dissimilarities between all pairs of grid cells using R package ‘betapart’75. We calculated four pβsim matrices for individual taxonomic groups (mammals, birds, reptiles and amphibians) and generated combined pβsim matrices for all terrestrial vertebrates by taking the mean pβsim values5. We compared eight hierarchical clustering methods on the pβsim matrices and assessed the performance of different algorithms in transferring the dissimilarity matrices into dendrograms using cophenetic correlation coefficients6. As the UPGMA achieved significantly better performance than the other clustering algorithms (Supplementary Fig. 4), we only used UPGMA clustering for further analyses. We selected suitable cut-off points in the dendrograms using the ‘recluster.region’ function in the R package ‘recluster’76 based on the explained dissimilarity and mean silhouette width5 considering the number of regions ranging from 2 to 15 (Supplementary Data 3). We defined the zoogeographical regions as the grid cells were geographically coherent and could be clearly delineated in space. We also ran NMDS ordination to investigate the relationships between zoogeographical regions based on the community compositions in two-dimensional space.Zoogeographical regions over phylogenetic timeTo assess the changes in zoogeographical regions surrounding the TP over phylogenetic timescales, we quantified pβsim between gridded species assemblages at different phylogenetic depths19,23. This method cuts a phylogenetic tree at a specified depth and collapses all descendent leaves into ancestral branches23. When the geographical distributions of the descendent leaves were merged into their ancestral branches, a branch × site matrix emerged for a predefined depth. We cut the phylogenetic trees into different time bins from 60 Ma to 0 Ma and generated four pβsim matrices for four individual taxonomic groups. We employed UPGMA clustering and NMDS ordinations to investigate the relationships among the gridded species assemblages based on the combined pβsim matrices for four taxonomic groups in each time slice. Again, we used the explained dissimilarity and mean silhouette width to determine suitable cut-off points in the dendrograms. We investigated the evolutionary history of the zoogeographical regions based on the topological and assignment changes in the UPGMA clustering dendrogram and NMDS ordinations. In addition, we assessed the strength of the relationship between the present-day pβsim matrices and those at different phylogenetic depths using the Mantel correlation test. To visualise the relationships between zoogeographical regions over phylogenetic time, we ran the NMDS for various time periods and maximised the correspondence between ordination pairs using Procrustes analysis19 via the ‘procrustes’ function in the R package ‘vegan’77. To assess the cross-taxon congruence in biogeographical processes, we performed these analyses separately for the four individual taxonomic groups.Zoogeographical regions over geological timeTo explore changes in the zoogeographical regions over geological time, we used Simpson’s beta diversity (βsim) to generate pairwise dissimilarities between fossil assemblages. Because sample completeness of the fossil records varies considerably in space and time (Supplementary Fig. 5), we assigned the fossil records to coarse spatial and temporal scales to strengthen the sampling intensity for each assemblage. To maximise the comparisons in the analyses based on the phylogenetic information, we assigned each fossil record to one of five time intervals: Eocene (56.0–33.9 Ma), Oligocene (33.9–23.0 Ma), Early Miocene (23.0–15.9 Ma), Mid–Late Miocene (15.9–5.3 Ma) and Pliocene–Pleistocene (5.33–11.8 Ka), and to one of the coarse‐grained zoogeographical regions identified by the present-day phylogenetic dissimilarity. We performed UPGMA clustering analyses and NMDS ordinations based on the βsim matrices in different time intervals to explore changes in the zoogeographical regions over geological time. For present-day structures of the zoogeographical regions, we ran analyses based on the extant mammal lists for the whole region to maximise the comparisons of the fossil data.Statistics and reproducibilityWe used Wilcoxon signed-rank tests to compare eight hierarchical clustering methods on the βsim matrices. We used UPGMA clustering analyses and NMDS ordinations based on the βsim matrices in different time bins to explore changes in the zoogeographical regions over time. We used Mantel correlation tests to calculate the correlation coefficients of the βsim between each taxonomic group and all terrestrial vertebrates in different time slices. Statistical significance was calculated with a permutation test. A P value of <0.05 was considered statistically significant. All statistical analyses were performed in R version 3.6.072. All raw data and custom R codes are available from the Dryad Digital Repository (10.5061/dryad.5x69p8d1078).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article.Supplementary informationSupplementary InformationSupplementary Data 1Supplementary Data 2Supplementary Data 3Description of Additional Supplementary FilesReporting SummaryPeer Review File
Biogeography Macroecology