Climatic and topographic changes since the Miocene influenced the radiation and biogeography of tent tortoises (Psammobates tentorius) species complex in southern Africa as inferred from mitochondrial and nuclear genes, microsatellite markers as well as ecological niche modelling


 Background

Climatic and topographic changes function as key drivers in shaping genetic structure and cladogenic radiation in many organisms. Southern Africa has an exceptionally diverse tortoise fauna, harbouring one-third of the world’s tortoise genera. The distribution of Psammobates tentorius covers two of the 25 biodiversity hotspots of southern Africa, the Succulent Karoo and Cape Floristic Region. The highly diverged P. tentorius represents an excellent model species for exploring biogeographic and radiation patterns of reptiles in southern Africa.
Results

We investigated genetic structure, population dynamics and radiation patterns against temporal and spatial dimensions since the Miocene in the Psammobates tentorius species complex, using multiple types of DNA markers and niche modelling analyses. Cladogenesis in P. tentorius started in the late Miocene when populations dispersed from north to south to form two geographically isolated groups. The northern group diverged into a clade north of the Orange River (OR), followed by the splitting of the group south of the OR into a western and an interior clade. The latter divergence corresponded to the intensification of the cold Benguela current, which caused western aridification and rainfall seasonality. In the south, tectonic uplift and subsequent exhumation, together with climatic fluctuations seemed responsible for radiations among the four southern clades since the late Miocene. We found that each clade occurred in a habitat shaped by different climatic parameters, and that the niches differed substantially among the clades of the northern group but were similar among clades of the southern group.
Conclusion

Climatic shifts, and biome and geographic changes possibly functioned as major driving forces in shaping cladogenesis, population dynamics and genetic structure in southern African tortoise species. Our results revealed that the cryptic cladogenic radiation of the P. tentorius species complex was probably shaped by environmental cooling, biome shifts and topographic uplift in southern Africa since the late Miocene. The Last Glacial Maximum may have impacted the population dynamics of P. tentorius substantially. We found the taxonomic diversify of the P. tentorius species complex to be highest in the Greater Cape Floristic Region. All seven clades discovered warrant conservation attention, particularly C1- C3, C5 and C6.

To help clarify the relationship between ecological processes and speciation, various ecological niche modelling (ENM) approaches have risen to prominence in the past decade that enable the answering of deeper evolutionary questions about niche evolution, speciation and the accumulation of ecological diversity within clades [39,40,41,42]. These ENM methods also became essential for identifying suitable habitat and for predicting habitat suitability temporal-spatially, making them critically important to investigations in evolutionary and conservation biology [43].
In this study, we investigated whether the microsatellite DNA based genetic structure derived for the P. tentorius species complex corroborated the phylogeny retrieved using mtDNA and nDNA sequence data. We hypothesized that climate uctuations and associated habitat changes over the Neogene and Quaternary in uenced cladogenesis and population dynamics in P. tentorius, and that physical barriers, in the form of either unsuitable climate or topography, possibly obstructed dispersal routes, leading to isolation and allopatric diversi cation. We also postulated that the diversi cation rate and taxon diversify would be higher in the Fynbos and Succulent Karoo vegetation of the GCFR than in the other biomes, as was found in plants. From an ecological perspective, we hypothesized that the different candidate species occurred in distinct ecological niches, and that the suitability of the niches of each one varied with time. Finally, we postulated that environmental and climatic variables drove range changes against temporal and spatial dimensions. To test these hypotheses, we calibrated a dated phylogeny to evaluate if the timing of cladogenic events corresponded to speci c climatic and/or landscape evolutionary events. This was followed by a character dependency analysis to assess if barriers such as the GE, CFMs and OR in uenced radiation in the P. tentorius species complex, and a habitat reconstruction analysis to gauge the role of vicariance, dispersal and extinction in its radiation. We also used ENM techniques to evaluate ecological niches across the 7 candidate species against temporal and spatial dimensions. In addition, we did a diversi cation rate analysis to better understand evolutionary patterns over time in this group. Lastly, diversi cation rate and ENM analyses were used to evaluate its conservation status and to predict future trends which could be critically important in that regard.

Genetic structure and phylogeny
The STRUCTURE results were generally congruent with the phylogenetic ndings of [35], showing 5 clusters that corresponded to the ve putative species of the optimal clustering scheme (at K = 6 in all cases) for both datasets partitioned by clades and subclades. The difference between C3 and C6 was, however, not clear (Fig. 2). Furthermore, the STRUCTURE results failed to distinguish between C1 and C4, and C5 and C7. The differences between C3 and C6 were nonetheless signi cant for scenarios K = 7 to K = 11 (results not shown). The DAPC membership cluster analysis, however, retrieved 7 signi cant clusters as optimal membership clustering scheme (Fig. 2), which perfectly corroborated the seven mtDNA clades retrieved in the previous study.
Both ML and BI analysis results for the combined dataset (mtDNA + nDNA) generated a tree topology similar to the mtDNA based phylogenetic tree of [35], except for some differences in branch lengths (Fig. 2). The different nodes were generally well supported, (BP > 70, PP > 0. 95), except the node between (C2 + C3) and C6. The mtDNA + nDNA tree topology also revealed low support values for node C6+(C2 + C3) (BP < 70, PP < 0.95). These ndings were congruent with the results reported in [35].
The taxon diversity in the Fynbos and Succulent Karoo vegetation of the GCFR (six clades in total: C1-C5 and C7) was higher than in the Nama Karoo (three clades in total: C1, C2 and C6) and the marginally occurring Albany Thicket (only C1).

Divergence dating analyses
When comparing the dating results between the gene tree ( Fig. 3) and species tree (Supplementary Figure  S1) for both the mtDNA dataset and the combined dataset, we found that the dating results were nearly identical. Also, the tree topologies of the mtDNA dataset and the combined dataset showed no con icts for both the gene trees and species trees. For both the gene tree and species tree, the combined dataset exhibited older age estimates, possibly due to the much slower evolution rate in the nDNA loci.
The total evidence chronogram (based on the mtDNA + nDNA dataset) revealed that the rst divergence that split Stigmochelys pardalis from the genus Psammobates occurred about 33.28 Ma (at node 1, see Fig. 3 and Supplementary Table S1) in the Eocene-Oligocene. Psammobates tentorius branched off from the sister group comprising P. geometricus + P. oculifer about 23.1 Ma (node 2) at the Oligocene-Miocene boundary. The divergence between P. geometricus and P. oculifer dated back to the early Miocene at about 20.12 Ma (node 3). The radiation of the P. tentorius species complex (node 4) started in the late Miocene, about 10 Ma, before the late Miocene cooling, and after the mid-Miocene Climatic Optimum. Soon after the start of the radiation, at about 9.14 Ma, C6 split off from C2 + C3 (node 5). At about 5. 7 Ma, C2 and C3 diverged (node 6) and more or less simultaneously, C5 segregated from (C1 + C4) + C7 about 5 Ma (node 7). The splitting of C7 from (C1 + C4) occurred at about 3.5 Ma (node 8) whereas C1 and C4 diverged at about 2.5 Ma (node 9). The rest of the radiations within clades occurred within the Pleistocene epoch.

SAMOVA results and habitat reconstruction
The SAMOVA analysis (when K = 7) retrieved seven geographic groups that corroborated the seven clades, with strong support ( xation indices F SC = 0.11, F ST = 0.82 and F CT = 0.8). In total, 80.04% of the total variance occurred among groups (Va: df = 6, sum of squares = 14555.33, variance = 114.24), 17.73% of the variance came from within populations (Vc: df = 112, sum of squares = 2834.23, variance = 25.31) and only 2.23% of the variance came from within groups (Vb: df = 46, sum of squares = 1583.72, variance = 3.18). Signi cance tests with 1023 permutations con rm this strong genetic structure against spatial dimensions (Vc and F ST , p < 0.0001; Vb and F SC , p < 0.0001 and Va and F CT , p < 0.0001).
The BioGeoBEARS model test suggested "DIVALIKE + J" as the best model with the highest AICc weight score (see Supplementary Table S2) and a signi cant p-value in the LRT test (p < 0.0001) for all three Page 7/51 analyses. The analysis based on the seven geographic regions suggested that all cladogenic events were caused by vicariance, but that the cladogenic events at nodes 4, 5 and 7 were also in uenced by dispersal (Table 1 and Fig. 4). Our results revealed that both dispersal and vicariance events caused bifurcation between the two major branches, the northern (C2 + C3) + C6) and southern (C1 + C4) + C7) + C5). Clade 6 and its habitat F was isolated from that of C2 + C3, and C5 separated from (C1 + C4) + C7). The separation between C2 and C3, the branching off of C7 from C1 + C4, and the split between C1 and C4 were all driven by vicariance events. The radiation originated in region BEF, which lay north of the GE and the Little Karoo area. Our biome-based reconstruction analysis suggested the Nama Karoo and Succulent Karoo as the ancestral biomes, in which cladogenic events of P. tentorius took place before those in the Fynbos (see Fig. 4 and Table 1). Our BioGeoBEARs analysis suggested that in the case of node 4, the separation between the Nama Karoo and Succulent Karoo represented a signi cant vicariance event which in uenced the cladogenic radiation of the P. tentorius species complex. It seems as though the northern population was derived from the Nama Karoo biome. The shifting of the latter biome in turn represented a vicariance event that split C3 from C2 and led to the isolation of the Succulent biome. In the case of the southern population, it possibly originated in the Succulent Karoo, and the shift from Succulent Karoo to Fynbos may have been the vicariance event that separated C4 from C1.
Our geographic barrier based analysis generally corroborated the ndings of the geographic region based reconstruction, which suggests that the cladogenic radiation at all nodes was strongly in uenced by geographic barriers, by way of vicariance as well as dispersal (see Fig. 4 and Table 1).
The phylorate plot showed two clear evolutionary rate shifts, one at C1 and the other at C2 (Fig. 5a). The RTT plot for the entire P. tentorius species complex showed a trend of increasing diversi cation rates since its radiation started about 10 Ma, levelling off only in recent time (about 1 Ma, Fig. 5b). The RTT plots at clade level ( Fig. 5c-i) revealed a trend of increasing diversi cation rates at C3, C4, C5 and C6, but especially at C6 (Fig. 5i). A clear decreasing trend was found at C1 and C2.
The BAMM results (Table 2) revealed C1 and C2 as having substantially higher net diversi cation rates (diversi cation rate -extinction rate) and lower turn-over rates than the rest of the clades, though C2 exhibits the lowest diversi cation rate. All clades and groups had higher speciation rates than extinction rates. The overall lineage turn-over rate was relatively high in clades C3, C5-C7 and the two major branches "(C1 + C4) + C7) + C5" and "(C2 + C3) + C6". The entire complex generally had a low overall net diversi cation rate and a relatively high lineage turn-over rate The plots of lambda retrieved from the BAMM analysis results (Fig. 5) revealed a clear increasing diversi cation trend at C3, C4, C5, and C7, but especially at C6 (increasing lambda in all cases), since the slope of linear relationship between lambda and the time line in each case was signi cantly steeper than that in the other clades. By contrast, C1 and C2 showed a remarkable decline in their lambda plots, indicating a decreasing trend in net diversi cation.

Dependent character analysis
The MuSSE analysis suggested that the GE and SM geographic barriers may have substantially in uenced the historical radiation of the P. tentorius species complex, and identi ed the "full model" (lambda, mu and q, which all differed across the three regions), with strong statistical support, as the best scenario (model likelihood scores and criteria provided as Supplementary Table S3). The biome based MuSSE analysis statistically con rmed that biomes affected radiation signi cantly, since the "full model" was favoured according to the LRT test. By contrast, results of the BiSSE analysis indicated that the OR did not in uence radiation in P. tentorius, because results for the best model, "free lambda", were not statistically signi cant.
The BiSSE analysis results revealed a much lower net diversi cation rate, and higher lineage turn-over rate, in populations north of the OR than in ones south of it ( Fig. 6 and Table 2). MuSSE analysis of regions indicated low net diversi cation rates and high lineage turn-over rates in populations south of the SM, and in populations occupying a mixture of Fynbos and Succulent Karoo vegetation ( Fig. 6 & Table 2).

Environmental niche modelling
Our ENM analyses generally provided meaningful predictions of suitable ranges for the different groups of the P. tentorius species complex (in all cases, the average test AUC (under ROC) > 0.90 with a standard deviation < 0.06; for details see Supplementary Table S4), suggesting that the modelling process was reliable. The optimized parameter settings for each of the ENM analyses are given in Supplementary  Table S5. Nevertheless, the mean omission rates of C5 and C7 showed relatively large deviations from the predicted omission rates (plot not shown), suggesting that the ENM range predictions for C5 and C7 should be interpreted with caution. Possible reasons for this may be their smaller ranges and smaller sample sizes.
The majority of the groups received relatively meaningful range predictions, which aligned well with their currently recognized distribution ranges (see Fig. 7 and Fig. 1), but not with the niche boundaries between C1 and C4, C4 and C5, C1 and C7, and C5 and C7, which are still fuzzy. The pairwise niche signi cance analysis results based on I statistical criterion found no signi cant differences between C1 and C5, and between C1 and C7, while all the other pairwise comparisons between clades were signi cant (see Supplementary Table S6).
Overall, when comparing the ENM analysis results across timelines, many groups exhibited substantial shrinking of suitable habitats in the Last Glacial Maximum (LGM) period (see Fig. 7 and Supplementary  Table S7) compared to the Last Interglacial (LIG) period. When comparing the size of currently suitable habitats to that in the Middle Holocene (MIDH), an overall declining trend was observed for the P. tentorius species complex, however, the decline was only signi cant in case of C2, C4 and C5 (Supplementary Table S7). As far as predicting future population dynamics scenarios is concerned, C1, C2, C3, C5 and C6 revealed substantial declines in suitable habitats, with only C4 and C7 exhibiting increasing trends in suitable areas. It is noteworthy that C2 and C5 showed continuous declines in suitable areas since the MIDH period.
In terms of climatic variables that made an impact in the ENM analyses (for details, see Supplementary  Table S8 and Fig. 8), substantial differences were found among clades and across different periods. In most cases, the LGM period appeared to be more greatly impacted by the changing climatic variables.

Discussion
Genetic structure and phylogeny Overall, our microsatellite DNA based DAPC clustering analysis results corroborated the sequence DNA (mtDNA + nDNA) based phylogeny, giving a similar genetic structure pattern as the seven clades.
However, our STRUCTURE clustering analysis did not reveal diagnosable differences between C1 and C4, and between C5 and C7 (Fig. 2). According to (Zhao et al., [35]), some of the species delimitation analyses also revealed that the species boundary between C1 and C4 was not clear, thus treating C1 and C4 as the same candidate species or Operational Taxonomic Unit (OTU) is certainly reasonable. However, treating C5 and C7 as a single candidate species or OTU, will con ict with our phylogenetic structure, since the DNA sequence based phylogeny did not support C5 and C7 as a sister group. Furthermore, the assigning of OTU or species status should not purely rely on genetic markers, but on multiple evidence which includes comparative morphology and ecology [44,45]. Our ENM analyses, for example, showed that the niche overlaps among clades of the southern group (C1, C4, C5 and C7) was substantial, and that their ecological niches did not differ signi cantly. Therefore, when deciding on the most parsimonious strategy, based on genetics and ENM, for assigning OTU status that would be crucial to its conservation management, we recommend that four OTUs or candidate species be recognized in the taxonomic revision of the P. tentorius species complex, namely, C2, C3, C6 and (C1 + C4 + C5 + C7). Furthermore, as a morphological character based study is in progress (Z. Zhao et al. in prep), we suggest that the nal taxonomic decisions also take into account those ndings.

Radiation patterns and biogeography
Cladogenesis in P. tentorius started during the early to late Miocene, a period initially characterised by warm temperatures [6], which may have facilitated dispersal of ancestral P. tentorius over a relatively wide area of southern Namibia and northern South Africa. Ancestral area analyses showed that region BEF (Nama Karoo and Little Karoo) was the most likely regions of origin. The absence of fossil records in the biometric analysis may, however, have compromised the accuracy of estimates, consequently, the divergence dates should be treated as ranges rather than means.
Our calibration dating analyses suggested that the rst divergence that split Stigmochelys pardalis from genus Psammobates matched a global cooling event at the Eocene-Oligocene boundary (about 33.9 Ma), which is believed to have led to aridi cation and major extinctions in many organisms, but without causing changes in biome composition [46]. This event may have led to range contractions, which facilitated the divergence between Psammobates and Stigmochelys. It also suggests that the split between Psammobates and Stigmochelys may not have been caused by biome shifts, but rather by climate change. The timespan when the P. tentorius ancestor diverged from other Psammobates species presumably covered the early Miocene, when temperatures dropped and the Antarctic ice-sheet formed [6]. In addition, the opening of the Drake Passage caused cold water (Benguela Current) to ow northwards along western South Africa [47] so that the southern part of Africa became progressively more arid [7]. This cooling and gradual aridi cation possibly contributed to diversi cation and the separation of ancestral species into distinct geographical regions. The fact that P. tentorius lls the geographic space between P. geometricus in the southwest and P. oculifer in the north, suggests that the ancestor of the three species possibly had a wide distribution across southern Africa during the early Miocene, and that global cooling in the early Oligocene likely restricted their ranges and isolated populations. Evidence shows a signi cant global temperature decline during the 25 − 23 Ma interval [6,48]. This cooling event possibly led to the contracting of ranges and refugia, resulting in the separation between P. tentorius and other Psammobates species. The divergence between P. geometricus and P. oculifer dates back to the early Miocene, coinciding with a warm and humid period (23 − 16 Ma, [48,49]). The warm and humid climate may have resulted in rapid changes in vegetation and habitats, which may have driven the divergence between P. geometricus and P. oculifer. Since the time of formation of the biomes perfectly matched the beginning of cladogenic radiation in P. tentorius [4,50,51], it is possible to speculate about the geographic region in which each taxon developed.
Divergence at node 4 divided the ancestral species into two lineages, a northernmost clade (C2 + C3 + C6) and a southernmost clade (C1 + C4 + C5 + C7). The "DIVALIKE + J" RASP route analysis indicated that the range of the ancestral group rst became restricted to region BEF, from where a southwards dispersal occurred to region BFAEG. The nal step was a divergence between BF and AEG, indicating a geographic split between the northern and southern lineages. From this scenario it would appear that ancestral P. tentorius crossed the OR despite it representing a potential geographic barrier. The analysis suggested that both vicariance and dispersal played a role in this rst divergence. The biome-based reconstruction did not show signi cant geographic con icts, and also revealed a signi cant vicariance event caused by biome segregation between the Nama Karoo and Succulent Karoo, suggesting that biome shifting also drove the rst divergence. Soon after the P. tentorius ancestor divided into two lineages, C6 diverged from C2 + C3 at node 5. The DIVALIKE + J RASP route analysis indicated a second crossing of the OR from region F, followed by dispersal into regions BC on the southern side of the OR in South Africa, and nally a divergence between populations north and south of the OR. This divergence apparently involved both dispersal and vicariance. The biome changes does not seem to have in uenced cladogenesis by vicariance at this point, but rather facilitated dispersal.
Within the P. tentorius species complex, radiation was possibly initiated in the late Miocene (about 9.98 Ma), so it seems as though the mid Miocene Climatic Optimum [49] was not responsible for the radiation of P. tentorius. Instead, the cladogenic events of P. tentorius seem to have coincided with the late Miocene cooling, intensifying of aridi cation of southern Africa [52], the late Miocene global diversi cation of succulent plants [52] and also the Miocene-Pliocene stepwise intensi cation of the Benguela upwelling [53]. The late Miocene cooling and intensifying of aridi cation may thus have caused habitat range shrinking and the retracting of refugia, which possibly initiated the cladogenic radiation of the P. tentorius species complex. As succulent plants are the major food source of P. tentorius [22,54], their diversi cation may have in uenced its habitat range substantially. The fragmentation of habitats with suitable succulent plants may in turn have resulted in isolation of subpopulations, providing the essential driving force for diversi cation. Furthermore, the intensi cation of the Benguela upwelling (at about 10 Ma) may have resulted in increasing aridity on the western side of South Africa [53], where the majority of populations of P. tentorius occurred, particularly along the west costal area. As discussed above, the increasing aridity may have resulted in the retracting of ranges and refugia, which possibly led to its speciation. Several studies using different methodologies indicated that major uplift events in southern Africa occurred during the Cretaceous but that tectonic uplift during the Cenozoic was at a relatively small scale [55,56,57]. Yet, uplifts of up to 250 m during the Miocene [58,59] may have in uenced the distribution routes of P. tentorius. The character dependency analyses suggested that the OR did not play a signi cant role in cladogenesis events, and that vicariance resulted from climate changes rather than from the development of the OR. Equally for node 5, we propose that dispersal related to P. tentorius crossing the OR for a second time and dispersing into the northern regions of South Africa, while vicariance resulted from populations retracting into refugial areas in response to a steep reduction in temperature after the late Miocene cooling [6]. The western GE in Namibia and South Africa may have served as both refugia and corridors for diversi cation in P. tentorius, as it had for other taxa [60]. More research is necessary to see if members of the clades still cross the OR on occasion and perhaps hybridize.
There are few reptile taxa with similar distribution patterns as P. tentorius in southern Africa, for which divergence dates have been estimated. Nevertheless, Tolley et al. [61] proposed that the deep divergences among major clades of the genus Bradypodion date back to the late Miocene, and implicated vicariance associated with both geological and climatic events. The latter included uplift of the western GE [59], formation of the Benguela upwelling system and late Miocene cooling [62].
The next two divergences, at nodes 6 and 7, occurred more or less at the same time in the late Miocene, with both divergences being ascribed to vicariance, and in the case of node 7 also to dispersal. The divergence at node 6 involved the separation of C3 on the western coastal plains (region C) from C2 above the GE (region BD). The most likely explanation for this cladogenic event is increasing western aridity due to the intensi cation of the Benguela upwelling in the period from 6.2-5.5 Ma [9,53,63], which perfectly matches the divergence time between C2 and C3. In parallel with western aridi cation, the initiation of a winter rainfall regime and the development of Succulent Karoo vegetation on the west during the late Miocene, came to full development in the Pliocene [7,8,21]. Our biome-based reconstruction results also corroborate this nding, of a vicariance event caused by the separation of the Nama Karoo and Succulent Karoo biomes which may also have facilitated the divergence between C2 and C3. Vicariance due to climate thus seems more likely than the western GE developing into a physical barrier between the coastal plains and interior plateau. The GE undoubtedly would have restricted gene exchange, but there are low-lying regions, particularly in the northwest, which could have served as dispersal routes between the coastal plains and plateau, and thus between C3 and C2. Nevertheless, there is evidence for a minor uplift in the west in the early Pliocene [59], which possibly strengthened isolation between the two clades.
Psammobates tentorius dispersed to the southern region of South Africa (ADEG) by the early Pliocene.
The rst divergence at node 7 isolated C5 in region E (the western Little Karoo) from the rest of P. t. tentorius. Thereafter, C7 in the eastern Little Karoo region G (Oudtshoorn basin) diverged from P. t. tentorius clades 1 and 4, in regions A and D. In both instances, vicariance seemed to have been of primary importance, although dispersal is also suggested. Both geographic region and geographic barrier based reconstruction analyses supported these assumptions. These divergences may have been caused by either climate change or Cenozoic exhumation [14,64]. In particular, uplift intensi ed greatly in the last 5 Ma for the GE alone, its highest elevation point rising from 300 m to 900 m. Therefore, the intensi cation of uplift may have facilitated geographical isolation. Since climatic and vegetational changes were not restricted to the Little Karoo, it seems unlikely that these changes could have played a major role in vicariance. Our biome reconstruction analysis also implied that biome shifting did not contribute signi cantly to cladogenesis at node 7 and node 8, despite a minor contribution from dispersal at node 8.
Evidence for the effects of geomorphological changes on vicariance in the southern clades seemed strong. Thermochronology results indicated that regions west and east of Worcester in the southwestern Cape underwent burial (up to 1.2 km) prior to exhumation (starting 30 − 20 Ma), which created the current relief with surrounding peaks reaching elevations of 1500 m [14]. It seems likely that over time this exhumation isolated P. t. tentorius populations in the Little Karoo from their western and northern counterparts. In addition, thermochronology results indicated differential uplift followed by exhumation (starting 30 − 20 Ma) in mountain ranges of the Oudtshoorn basin [14]. These events could have isolated C7 in region E. Despite being on a more recent time scale, divergence within a Little Karoo endemic plant species, Berkheya cuneata, shows two distinct lineages associated with the western and eastern Little Karoo, respectively [65], similar to C5 and C7 of P. t. tentorius. Pygmy geckos (Goggia) also show a similar differentiation pattern in the southern Cape [24]. On a broader scale, Cowling et al. [20] proposed that increased topographic heterogeneity in response to moderate uplift in the Miocene played a major role, together with climatic deterioration, in the rapid diversi cation of the Cape plant clades.
The nal divergence between C4 and C1, respectively, on the western GE (region D) and below the southern GE (region A), occurred during the late Pliocene with vicariance as the main driver, possibly due to both tectonic and climatic events. The late Miocene and early Pliocene (7.7-4.0 Ma) was relatively humid [66], which may have facilitated movements and gene exchange between populations below and above the southern escarpment. Also, as discussed above, the uplift intensi ed substantially in the last 5 Ma (for the GE alone) which may have facilitated geographic isolation. A possible contact zone may have been the relatively low slopes of the GE in the northern Tankwa Karoo, southwest of Calvinia. The Tankwa Karoo is currently one of the most arid parts of the Karoo with < 100 mm rain per annum [67].
Aridi cation may have closed this contact zone and isolated C1 and C4 from each other. Green et al. [14] provide evidence for a Cenozoic uplift of around 700 m and subsequent denudation in the Beaufort West GE region, but there is no evidence suggesting that this tectonic event was localised or affected the remaining southern escarpment. Nevertheless, if such an uplift in uenced dispersal routes across the GE, it should have contributed to this diversi cation event. Noteworthy about this is that our ENM results also predicted a substantial contact zone in the Tankwa Karoo between C1 and C4, which supports our assumption. More importantly, this Tankwa Karoo contact zone between C1 and C4 appeared to show a substantial declining trend from LIG to the recent and future, which may imply that the further diversi cation between these two clades will be reinforced.
The phylogeographic analyses detected several Pleistocene radiations within C1, C2, C4 and C6, while the lack of ne structure in C3, C5 and C7 may be due to small sample sizes. Periodic climatic cycles in response to Milankovitch forces caused glaciers to wax and wane over the Pliocene and Pleistocene [68], which in uenced species distribution and diversi cation worldwide [2]. Radiations within C2 and C6 were possibly driven by the southward expansion of Kalahari sands in response to Pleistocene glaciation [16], as has been shown for the Namaqua rock rat (Micaelamys namaquensis) [69] and Chacma baboons (Papio ursinus) sensu lato [70].
Generally, our reconstruction analyses revealed that each clade occurred in a unique geographic region, suggesting that both biome shifting and geographic barriers played important roles in the cladogenesis of P. tentorius. Furthermore, the biome-based reconstruction analyses imply that biome formation and distribution may not be a once-off event, and that biome shifting did not necessarily correlate with allopatric speciation in P. tentorius.
Diversi cation rates, potential cladogenic drivers and regional patterns The radiation in P. tentorius was apparently not in uenced by the OR but differed among regions and among biomes. The latter results, however, become more complex when considering the association between clades and regions or biomes. The lineage turnover rate was low for the region above the GE, which is paralleled by a low turnover rate for C2, but not for C4 or C6. Lineage turnover rate was also low for the Nama Karoo, as well as for C1 and C2 associated with this biome, but the turnover rate for C6 in the Nama Karoo north of the OR was high. The region and biome analyses were not very informative possibly because the division of characters oversimpli ed the great climatic and topographic heterogeneity of the southern African landscape [38,71].
Net diversi cation rate for the entire P. tentorius species complex was low but lineage turnover rate was high. Low net diversi cation rates and high lineage turnover rates applied to all clades, except C1 and C2, occurring in the Nama Karoo in central South Africa. The high net diversi cation rates of C1 and C2 suggested that these clades had accumulated much diversity over a relatively short period of time, possibly due to a lower extinction rate rather than a higher diversi cation rate than in the other clades.
Yet, the evolutionary dynamics of C1 and C2 differed, with signi cant rate shifts in both clades associated with an increase in rate for C1 and a decrease in rate for C2. Species distribution models for Bitis arietans indicated that the interior of South Africa became inhospitable for this species during glaciation periods, and that populations continuously retracted and expanded to and from refugia in the northern and southern regions of South Africa [72]. We propose that C2 retracted northward during glaciation periods and that the relatively low landscape heterogeneity of their habitat (Bushmanland and Upper Karoo bioregions) reduced rate shifts in it. Our ENM results also supported the increasing range suitability from LIG to the recent in the northern population. In contrast, C1 possibly retracted to southern regions near the CFMs in the south, which include several vegetation types such as the Lower Karoo bioregion, Albany Thicket and Fynbos. This also corroborates our ENM results (from LGM to recent). The higher topographic and vegetation heterogeneity of the southern region may explain the increase in rate shift in C1 but a decreasing rate shift in C2.
The evolutionary dynamics were similar for C3, C5, and C7 below the GE, and for C4 and C6 above the GE.
The two biodiversity hotspots of the GCFR cover the full distribution ranges of four P. tentorius clades (C3, C4, C5 and C7), and part of the ranges of C1 and C6 (see Fig. 1 in [73]). Topographic diversity in the habitats of C4, C5, C6, and C7 is high but low for C3, except in the Richtersveld in the north, which includes the arid western escarpment [60]. Associated with this topographic diversity is high climatic diversity, with a strong aridity gradient from north (very arid) to south (more mesic) along the western, winter-rainfall region of South Africa [13]. The southern region is arid south and north of the SM, but becomes more mesic towards the east, and topographic heterogeneity creates localized climates [71] over the range of the southern clades. North of the OR, the habitat of C6 is also diverse not only with regard to topography but also climate. The southwest is highly arid and receives winter rains, whereas rainfall increases towards the east is received in summer rain [74]. The topographic and climatic heterogeneity of the abovementioned regions probably explain their similarities in evolutionary dynamics, characterised by high diversi cation rates, countered by high extinction rates, culminating in low net diversi cation and high lineage turnover rates.
The GCFR is recognised for its great oristic richness and endemism [37,38], with much of its present diversity ascribed to radiation events during the Pliocene and Pleistocene [5]. The Pliocene-Pleistocene glacial cycles involved uctuations between cool-arid conditions during glacial periods and warm-humid conditions during interglacial phases [1,6], causing the contraction and expansion of taxa in and out of refugial areas. Several studies suggested the existence of western and eastern refugia in the southern GCFR [72,75], as well as refugial areas in mountains of the western GCFR, such as the Cederberg and Kamiesberg [60,76]. It appears that large patches of Fynbos and Succulent Karoo biomes persisted as refugia over glacial/inter-glacial cycles, which was not the case with the Nama Karoo; recent expansion of the Nama Karoo may thus explain present contact zones of many taxa [5] and references therein. Such contact zones for P. tentorius include intergradation zones between C1 and C2, and between C2 and C4, which in both instances may have involved recent expansions of Nama Karoo clades.
Perspectives on potential impacts of climate and environmental change Our ENM results revealed that clades C2, C3 and C6 occurred in unique niches. The niche boundaries between C1, C4, C5 and C7 are, however, not clear. These results imply that C2, C3 and C6 occur in unique habitats, in uenced by different climatic and environmental factors compared to clades C1, C4, C5 and C7. The latter four clades occur in habitats shaped by similar climatic and environmental factors. The Tankwa Karoo region (southwestern Karoo) was regarded as a "three-way" contact zone, where all three currently recognized subspecies were believed to occur, given the high polymorphism in morphology found in one area [77]. Zhao et al. [35], however, found that all these morphs belonged to the same clade, C1 (supported by both mtDNA and nDNA sequence data). Our ENM analyses revealed substantial overlapping zones of suitable areas in the Tankwa Karoo for C1, C4 and C5 (corresponding to P. t. tentorius), C2 and C6 (corresponding to P. t. verroxii). This interesting nding implies that the high phenotypic variation in the Tankwa Karoo population was possibly driven by convergent evolution due to high microclimatic heterogeneity, making the Tankwa Karoo suitable for all these clades.
Clearly, the recent Pleistocene climatic oscillations (from the LIG to LGM period, but particularly in the LGM) functioned as important evolutionary drivers that impacted population dynamics. Areas with suitable habitats and potential ranges for the P. tentorius species complex were primarily based on our ENM analyses, with range shifts becoming particularly pronounced in the LGM period. Although there was an overall decline in the potentially suitable range area from the LIG to LGM for the P. tentorius species complex overall, possibly due to cooling which led to the retracting of ranges and refugia, the decline was signi cant only in the case of C3, C4, C6 and C7. The reason for the increasing ranges of C1, C2 and C5 is unclear, but may be correlated with variations in localized climates during the LGM. Studies in tropical and subtropical areas found that different population dynamics scenarios could shape glacial refugia [1,78]. Noteworthy about our ENM results was that they revealed a diagnosable subdivision within C6 (between western and eastern populations), which seems to have been reinforced since the LGM cooling. This was also observed in the microsatellite based cluster analysis and DNA sequence based phylogeny. Besides this, our haplotype based diversi cation rate analysis revealed an increasing trend of diversi cation that corroborated our ENM results for C6. The predicted potential range exhibited an overall decreasing trend. This is probably related to environmental gradients that limited homogenizing through gene ow or it could be a refugial area that has persisted over time, resulting in a build-up of diversity. The range shrinking and refugia retracting can facilitate isolation and speciation [1,79]. Collectively, it implies that further divergence and cladogenic radiation within C6 is very likely. This was also advocated by [35], who found that genetic diversity within C6 was comparatively high.
When comparing the results of diversi cation rate analyses (Fig. 5) and ENM analyses in the LGM (Supplementary Table S7), an interesting trend was observed, namely, that a decrease in suitable range generally correlated with an increase in diversi cation rates, and vice versa. This trend is in-line with the theory that cooling results in the retracting of ranges and refugia, which facilitate diversi cation and speciation.
Our ENM analyses revealed that range suitability of each clade was shaped by different combinations of climatic variables. The impact of climatic variables at different periods was also different for each clade.
The LGM was the period in which the climatic variables had their most signi cant impact on most clades (see Fig. 8). It re ects the high level of ecological differences among clades, particularly between the northern (C2, C3 and C6) and southern groups (C1, C4, C5 and C7).
Overall, potentially suitable range for the P. tentorius species complex appears to be strongly in uenced by annual precipitation, temperature and temperature seasonality, with annual precipitation as the most important one. Rainfall patterns determine primary productivity (plants), and is therefore strongly correlated with the distribution of vegetation and biomes. This is especially crucial for P. tentorius, as a herbivore [22]. These ndings also corroborated those of our biogeographic analyses that cladogenesis was substantially in uenced by biome and vegetation characteristics. Some climatic variables were, however, ltered out after the correlation analysis, which may also in uence range suitability, making further ecological studies crucial to unmasking the full ecological picture.
Prospective conservation management plan Since P. t. tentorius (C1, C4, C5 and C7) and P. t. verroxi (C2 and C6) have been categorized as "Near Threatened", and P. t. trimeni (C3) downgraded to "Endangered" by the IUCN in 2018 [80], they all deserve conservation attention. Our ENM analyses showed a substantial decline in the current suitable range of C2, C4 and C5, and also predicted a substantial declining trend for C1, C2, C3, C5 and C6 by 2070. Both C1 and C2 also showed consistent decreases in diversi cation rates. These ndings highlight the need for conservation attention at the very least for C1 and C2. Apart from the in uence of varying climate, there are also anthropogenic factors such as urbanisation which will likely reinforce future range declines, particularly in the southern group, where such activity is high.
Although most of the range of C3 falls in a protected area (Namaqua National Park), the ndings support the prediction of Branch [22], that it will be impacted by ongoing climate change. The distribution ranges of C1, C2 and C4 overlap with those of other species listed by the IUCN in 2018 [80], namely, the Karoo dwarf tortoise (Chersobius boulengeri, "Endangered"), speckled dwarf tortoise (Chersobius signatus, "Endangered") and the armadillo girdled lizard (Ouroborus cataphractus, "Vulnerable"), so they could therefore indirectly bene t from conservation attention given to those taxa. An area, for example, warranting preservation in this regard is the biological corridor between the Nuweveldberge and Sneeuberge, where C1 and C2 coexist along with the sympatrically occurring Karoo dwarf tortoise. Unfortunately this area is not protected yet.
The home range of C6 certainly deserves conservation attention, given that it has the highest genetic potential for further cladogenesis. Its distribution range overlaps with that of the Nama dwarf tortoise (Chersobius solos) categorized as "Vulnerable" [80]. Fortunately, the region already has protected areas like the Ai-Ais Transfrontier Park and Tirasberg Conservancy.
Lastly, clades occurring in the Fynbos and Succulent Karoo (C3-C5 and C7) generally exhibited low net diversi cation rates and high turn-over rates, particularly those in the Fynbos which that showed negative net diversi cation rates. These ndings emphasize the importance of conserving the Fynbos and Succulent Karoo biomes.

Conclusions
The study provided strong evidence that cladogenesis in P. tentorius can be linked to climatic uctuations and topographic changes in southern Africa since the Miocene, thus supporting our hypothesis. It appears if climate was of greater importance than topographic changes in earlier diversi cation events, but that uplift events together with climate change played a signi cant role in later divergences. The climatic and topographic changes linked here to early divergences in P. tentorius were also proposed for vicariance in other reptiles at the generic level. Ecologically, each clade in the northern group (C2, C3 and C6) occurs in a unique niche shaped by different climatic factors. By contrast, clades in the southern group (C1, C4, C7 and C5), showed no signi cant niche differences.
The results also correspond to other studies showing high taxon diversity (in terms of the number of clades) in the GCFR, not only for plants but also for animals, including reptiles. Diversi cation patterns of P. tentorius in the late Miocene and Pliocene thus seem to have paralleled those of other organisms, supporting the hypothesis of higher diversity in the GCFR than elsewhere over its distribution range. The strong association of P. tentorius clades with particular regions and vegetation types suggests that the clades evolved allopatrically and that contact in restricted areas is recent, following range expansions of some clades. However, although the clades abut, they do not necessarily overlap because vegetation in the regions regarded as possible intergradation zones forms a mosaic, which may still keep clades distinct.
Nevertheless, more research is necessary to establish if the clades hybridize in the so-called intergradation zones. Conservation awareness is warranted for all clades in the P. tentorius species complex, particularly for C1, C2, C3, C5 and C6. This study together with the ndings in [35] provides strong evidence that P. tentorius requires a taxonomic revision, which would impact the Red List Assessment of the species. As a species, the IUCN currently lists P. tentorius as "Near Threatened" [80] and P. t. trimeni as "Endangered" [80].

Taxon sampling
We sampled 404 specimens of P. tentorius from 76 localities throughout its distribution range in southern Africa (see Supplementary Table S9

DNA extraction, PCR ampli cation and sequencing
Tortoises are slow evolving compared to most other reptiles [86]. Thus, fast evolving mtDNA markers rather than slow evolving nDNA markers should be used to reconstruct their phylogenetic histories, particularly those of small-scale species complexes [87]. Studies showed that slow evolving nuclear markers are only suitable for studying the phylogeny of organisms with a wide timespan, whilst, faster evolving mtDNA markers become more useful when focusing on recent events with relatively short timespans [87,88]. Despite this, Townsend et al. [89] found that nDNA PRLR is one the fastest evolving nDNA maker for inferring phylogenies in the order Squamata.
We collected fresh tissue from tail tips (15-25 mg), tissue from dead carcasses or blood (0.1 ml) from the subclavian vein. Fresh tissue was preserved in 96% ethanol and blood was dissolved in 1% of a 10M Sodium-EDTA anticoagulant solution. Connective tissue and bone were kept dry in vials, whereas preserved muscle was stored in 96% ethanol. All samples were stored at -80 ºC in the laboratory until we did the DNA extraction. Genomic DNA was extracted using QIAGEN DNeasy Blood and Tissue kits (QIAGEN, Germany) following the manufacturer's instructions. For the ancient DNA extraction from museum specimens or severely degraded shells, we followed the protocol used in [35].
Polymerase Chain Reactions (PCR) were performed in a BIO-RAD T 100™ Thermal Cycler (Singapore) under the following parameters: an initial 4 min denaturation step at 94 ºC, followed by 37 cycles ( for the concatenated dataset. We also used jModeltest v.2 [98] to determine the best substitution model and parameter settings for each data partition via AIC criterion (Supplementary Table S11). Nucleotide base biases across different partitions were determined through the homogeneity test implemented in PAUP v.4.0 [99]. Substitution saturation tests were performed using DAMBE v.6.1.9 [100] at both gene and partition (on protein coding gene) levels and were visually plotted for transitions and transversions against GTR-model modi ed genetic distance diagrams. This was done to investigate potential substitution saturation, particularly at the third codon position in the protein coding fast-evolving mtDNA.
The DAMBE v.6.1.9 software was used to read sequence frames in order to determine the codon positions.

Microsatellite DNA genotyping
Extracted genomic DNA was used as template DNA in the PCR. Twenty two microsatellite DNA loci, which were successful in previous testudines population genetics studies [101,102,103,104,105,106,107) were tested. They were tested on the target animal using PCR, with all forward primers labelled uorescently with different dye probes (primer details and optimal annealing temperatures are given in Supplementary  Table S12). All PCR reactions were performed using KAPA2G Robust HotStart Readymix, USA, in a BIO-RAD T 100™ Thermal Cycler (Singapore) under the following parameters: an initial 5 min. denaturation step at 94 ºC, followed by 35 cycles of 30 seconds of denaturation at 94 ºC, 30 seconds of annealing (different optimal annealing temperatures were used depending on the loci, details given in Supplementary Table S12), and a 1 min. extension step at 72 ºC, with a nal 10 min. extension step at 72 ºC. The PCR annealing temperatures for the different microsatellite DNA loci were optimized using temperature gradient tests. Only loci showing a positive amplicon peak were used for further genotyping.
We applied multiple mixed PCR's to enable the simultaneous ampli cation of multiple microsatellite DNA markers. We subdivided the microsatellite markers into six multiple mixed PCR groups according to the type of dye probes, fragment lengths and optimal PCR annealing temperatures, to maximize visualization of the genotyping and to minimize the interference among the different loci during genotyping (details given in Supplementary Table S12). The PCR conditions for these six multiple mix reactions were the same as above for the single locus based PCR. The annealing temperature for each multiple mix reaction is also given in Supplementary Table S12.
In order to verify the microsatellite repeat motif unit in the P. tentorius species complex, and to ensure that the amplicons obtained from the dye labelled primer pairs indeed came from the target region, and not from other regions (caused by unspeci c primer binding sites), we sequenced some individuals of the seven clades with the same primer pairs but without the dye labels (with the same PCR conditions as the above genotyping). The PCR products were electrophoresed in 1% agarose gel, visualized under UV light, and puri ed using a BioFlux PCR Puri cation Kit (Bioer Technology, China). Puri ed PCR products were cycle sequenced using BigDye (ABI PRISM® BigDye Terminator v3.1 Cycle Sequencing Kits, USA) and standard methods. The Big-Dye PCR products were puri ed with Zymo DNA Sequencing clean-up kits (Epigenetics Company, USA), prior to sequencing in an ABI 3500 genetic analyser.
We used PCR amplicon length for microsatellite genotyping and analyses throughout the study. In all cases, CONVERT [108] and PGDSpider [109] was used to convert datasets into different input le formats for different analyses. Sanger sequences were analysed using ABI Prism Sequencing Analysis software v.3.7 (Applied Biosystems), then aligned with MUSCLE v.3.2 and manually checked with MEGA v.7.
Among these tested 22 microsatellite loci, 19 showed positive results and were used for carrying out further analyses. To ensure that the microsatellite dataset did not contain null alleles, stuttering or large allele dropout [111], we checked it using Micro-Checker [111] to avoid misleading genotyping results. For all classes a combined probability of p < 0.05 was considered a signi cant indicator of the presence of null alleles. Consequently, 14 out of the 19 loci were detected as negative in exhibiting null allele signal, and were used in further analyses (for details, see Supplementary Table S13).

Genetic structure and phylogenetic analysis
To determine genetic structure and identify signi cant clusters in the microsatellite dataset, the Discriminant Analysis of Principal components (DAPC, [112]) was performed using the R package 'adegenet' [113] of the program R v.3.5. The optimal group membership cluster scheme was mapped using the " nd, clusters" function. Additionally, STRUCTURE v 2.3 [114] was used to map the best cluster scheme by using Admixture as ancestry model to allow individuals the possibility of having mixed ancestry with the Bayesian algorithms. In order to test genetic structure at different levels, we partitioned the dataset into two testing schemes. In the rst scheme, individuals were divided into seven mtDNA clades, thus giving a total of seven populations. In the second scheme individuals were divided into 10 populations. Clade 1 was subdivided into a western and an eastern population. Clade 6 was also subdivided into a western and eastern population. Both C1 and C6 showed clear subdivisions in [35]. As for C2, the northern population was morphologically distinct from the southern population, we therefore subdivided it into a northern and southern population, although the mtDNA loci did not show them as clearly separate according to Zhao et al. [35]. Analyses were run for 5 million generations using MCMC with the rst 50, 000 discarded as burn-in. STRUCTURE HARVESTER [115] was used to summarize outputs of STRUCTURE and to determine the best K value ((the one with the largest In Pr(X|K)). The program Clumpak [116] was used to visualize population structure across multiple cluster schemes (clustering under different K values).
To compare the genetic structures between the microsatellite and sequence DNA datasets (mtDNA + nDNA), we rst removed all sequences from alignments having outstanding microsatellite genotyping data to ensure that the mtDNA dataset completely matched the microsatellite DNA dataset. We then combined the trimmed mtDNA dataset with the nDNA dataset. We used BEAST v 2.4 [117] and the STACEY package [118], respectively, to carry out a Bayesian MSC phylogenetic analysis and to generate a species tree and a gene tree. We optimized the prior setting based on results of PartitionFinder and jModeltest (for details, see Supplementary Table S11). All individuals were assigned to the seven mtDNA clades. Analyses ran with 55 million generations, sampling every 5000 generations and discarding the rst 25% as burn-in. Results were loaded into Tracer to check sample mixing to ensure the majority of the parameters had ESS > 200. Both gene and species trees generated from BI analyses were imported into the TreeAnnotator package implemented in BEAST after the rst 10% of trees was discarded as burn-in. In addition, we also performed a maximum likelihood analysis (ML) using RAxML v.8 [119]. To evaluate the strength of support for each node, we used the same partitioning scheme as with the BEAST analysis, with 1000 nonparametric bootstrap replications. The ML algorithm model selected was GTRCAT. In all cases, bootstrap support (BP) > 70% for ML [120] and posterior probability (PP) > 0.95 for Bayesian inference [121] were considered as strongly supported.
To test whether Fynbos and Succulent Karoo vegetation of the GCFR had greater diversity than other biomes, we used ArcGIS v. 10.4 (ESRI, 2016) to strictly examine the biome of each individual on the basis of genetic structure. We then counted the taxon diversity in the Fynbos, Succulent Karoo and Nama Karoo.

Divergence time dating
In order to generate an ultrametric tree to meet the requirements of habitat reconstruction analysis, we used DnaSP v.5 [122] to determine the haplotype sequences, and removed all the excess identical sequences in the concatenated mtDNA and nDNA datasets. Since the saturation test did not detect nucleotide saturation, we performed the complete dating analyses with the combined mtDNA + nDNA and mtDNA sequence datasets, respectively. We used 171 haplotype sequences to perform the Bayesian calibration dating analysis with the program BEAST v.2.4.8 with StarBEAST package for both the combined and mtDNA datasets. Each partition was applied with different optimized substitution models and parameters determined by another independent test using PartitionFinders v.2 and jModeltest v.2 (see Supplementary Table S14). Trees and clock models were linked, site models unlinked, and dating priors in the most recent common ancestor (MRCA) set based on ve calibration nodes from a dating study using the complete mtDNA genome [123]. For the mtDNA dating analysis, the log normal MRCA prior setting details are given in Supplementary In order to de ne population groups that were geographically homogeneous and maximally differentiated from each other, as well as to identify the potentially crucial genetic barriers that separate these population groups, we used SAMOVA 2.0 [124] with concatenated mtDNA to delineate groups with temporal and spatial dimensions. SAMOVA analyses with ranges of K = 5 to K = 8 were tested, and the results with minimal subdivision in clades were used for further habitat reconstruction analysis.

Habitat reconstruction analysis
This analysis was done to test how the shifting of geographic regions, biomes and topographical barriers (GE and CFMs) in uenced cladogenic radiation process against a timeline. Three independent habitat reconstruction analyses were performed on the seven geographic regions, the three biomes and the three regions separated by the topographical barriers speci ed below.
First, we partitioned the distribution range of the P. tentorius species complex into seven main geographical areas according to the optimal grouping scheme suggested by results of the SAMOVA 2.0 analysis (see Fig. 1). Second, we split its range into the three biomes, namely, Nama Karoo, Succulent Karoo and Fynbos. Last, we partitioned its range into three regions: north of the GE, between the GE and CFMs, and south of the CFMs. Biogeographic reconstructions were then performed on each of these three partitions with the R package BioGeoBEARS [125,126,127] via the program RASP v.4.0 [128], respectively. To unify the methods of biogeographical inference, such as the dispersal-vicariance analysis (DIVA), the Dispersal-Extinction-Cladogenesis (DEC-Lagrange) analysis and the Statistical BayArea analysis, the BioGeoBEARS analysis rst tests the optimal biogeographic reconstruction model by evaluating the founder effect parameter "J" [129,130]. Thus, for each of the three analyses, we rst compared combinations of the three models with and without "J", and selected the optimal model based on AIC corrected (AICc) and AIC weight criteria with the LRT test. The ultrametric species tree generated by BEAST BI analysis with the mtDNA + nDNA dataset was used for all three biogeographic reconstruction analyses. All outgroups were truncated before initiating analyses to improve the accuracy of the reconstruction results [128].

Diversi cation rate analysis
In order to detect and visualize the diversi cation dynamics, evolution rate heterogeneity and signi cant rate shifts across different clades on the calibrated BEAST chronogram (mtDNA + nDNA dataset), we employed the Bayesian statistical framework approach [131,132] of the BAMM software [133], together with the R packages 'BAMMtools' [133] and Ape [134]. We used the "setBAMMpriors" function to determine a proper prior for the expected number of shifts, initial lambda, lambda shift and initial mu.
Analyses were conducted with four simultaneous MCMC chains for 50 million generations, using the "speciation-extinction" model, deltaT = 0.1, Swap period = 1000 and default settings for the rest of the parameters. Pairwise probability comparisons obtained were then used to reconstruct a matrix with macro-evolutionary cohort analysis [135] to look for differences in the macroevolutionary rate regime among the lineages of the P. tentorius species complex based on haplotypes. A phylorate plot with optimal rate shift scheme was generated through BAMMtools to visualize evolutionary rate variation across the Bayesian phylogenetic tree topology. Rate through time (RTT) plots were constructed using BAMM on all groups to visualize the radiation rate shift against temporal scales in different groups. To compute the diversi cation rate (lambda) and extinction rate (mu) independently for each clade and crucial evolutionary group, R packages 'dplyr' [136], 'readr' [137] and 'ggplot2' [138] were used to calculate lambda and mu for each target group. The net diversi cation (r = lambda -mu) and turn-over rates (t = mu/lambda) were used as indicators to evaluate the radiation of the P. tentorius species complex against temporal dimensions [139], as well as predict future trends.

Character dependency analysis
A character dependency analysis was used to verify whether the cladogenesis and radiation of the P. tentorius species complex was in uenced statistically by geographic barriers and biomes. When considering the OR as a geographic barrier, we used the BiSSE approach [140] instead of the GeoSSE approach [141] because there is no clear evidence of population overlap across the OR. We performed the BiSSE analysis (6 parameters) with R packages 'ape', 'diversitree' [142] and 'phytools' [143] to investigate whether the OR signi cantly in uenced the radiation history of the P. tentorius species complex. To do this, we coded populations north of the OR as the "0" state, and populations south of it as the "1" state.
We then built likelihood functions into different scenarios: a) full model ((through function 'all.different', lambda1 ≠ lambda2, mu1 ≠ mu2, q01 (transition rate from character "0" to "1") ≠ q10 (transition rate from character "1" to "0")), b) lambda different model (through function 'free.lambda', lambda1 ≠ lambda2, mu1 = mu2 and q01 = q10), c) variable q model (through function 'free.q', lambda1 = lambda2, mu1 = mu2 and q01 ≠ q10). We used the ANOVA based likelihood ratio test (LRT) to determine the "best scenario" under the criteria ln likelihood value and AIC criterion. The best-t model was used to run a Bayesian MCMC analysis. First a preliminary run with 100 steps was done to estimate parameter settings. The formal long run was done with 100000 steps, sampling every 1000 steps.
We used the MuSSE approach [144] implemented in R packages 'ape', 'diversitree', 'phytools' and 'geiger' [145] to determine whether barriers like the GE and Swartberg Mountain (SM) signi cantly in uenced the radiation of the P. tentorius species complex. We partitioned the complex's distribution range into three geographic regions as "1" (above the GE), "2" (below the GE, but north of the SM including the western region below the GE) and "3" (south of the SM). In total, 12 parameters were included in the analysis.
The same method was used to model the interactions between cladogenesis and biome. To simplify the biome analysis and interpretation, we excluded Albany Thicket and included only Nama Karoo, Succulent Karoo and Fynbos. Furthermore, we considered only Succulent Karoo as relevant in the Little Karoo because tent tortoises prefer this vegetation type in that region. Our biome assessment included 1) Nama Karoo on the northern and southern sides of the GE, 2) Fynbos on the western side of the GE and 3) Succulent Karoo along the west coast and in the Little Karoo, and we encoded them as characters "1", "2" and "3", respectively. The rest of the steps followed the MuSSE analysis of the geographic regions. For all the above BiSSE and MuSSE analyses, the input tree used was the calibrated BEAST chronogram with the mtDNA + nDNA dataset.

Environment niche modelling analysis
To accurately estimate suitable habitats through timelines (from past, current to future), 455 samples with genetic information con rmed were used. We did not include museum or internet records which lacked genetic information, since many of the museum collections have problems of incomplete information. Datasets with detailed coordinate information were partitioned into seven separate datasets, based on the seven clades. In total, eight datasets were independently analysed by ENM, namely, datasets C 1-7 and the combined dataset (including all 455 specimens of P. tentorius used in this study). Each dataset was written in WGS1984 pseudo-projection format. To reduce the sampling bias and spatial auto-correlation, we used the R package 'spThin' [146] to thin each dataset to no more than one sample per kilometre.  [148].
To decrease the correlation between the bio-climatic variables, we rst performed a preliminary ENM analysis with all 19 variables. One hundred bootstrap replications were run in the program Maxent v. 3.4.1 [149], with 10,000 maximum iterations. A jack-knife method was applied to measure variable importance.
The R packages 'maptools' [150], 'rgeos' [151], 'raster', 'rgdal' [152] and 'dismo' [153] were used to calculate the pairwise correlation of the 19 variables. If the correlation coe cient between two variables was greater than 0.7, we omitted the one with the lower contribution based on the average contribution from the preliminary ENM run. After the omission procedure of the correlation test, we retained seven variables: Bio1, Bio3, Bio4, Bio8, Bio9, Bio12 and Bio19 (Supplementary Table S17) to estimate the potentially suitable habitats for each clade and the entire P. tentorius species complex across the timeline.
To optimize the parameter settings, the R package 'ENMeval' [154]       (a) The phylorate plot shows the haplotype diversi cation rate across the BEAST chronogram topology.
The "warm" colours represent fast rates, whilst "cool" colours represent the branches with slow rates (units per event or lineages per Ma). The red spots indicate corresponding nodes with signi cant rate shifts. (b) The rate though time (RTT) plot for the overall Psammobates tentorius complex. (c)-(i): The RTT plots for different clades (C1-C7). The red dot indicates a signi cant rate shift detected by BAMM. Figure 6 The estimated lambda, mu and q (A-C) from the BiSSE analysis: "0"-populations from north of the Orange River, "1" -populations from south of the Orange River; D-F, MuSSE analysis results of three geographic regions: "1" -region above the GE, "2" -region south of the GE, and north of the SM and "3" -region south of the SM. Plots show lambda, mu and net diversi cation rates (r) across the three different regions; G-I, MuSSE analysis modelling the lambda, mu and r of three biomes: "1" -Nama Karoo, "2" -Fynbos mixed with Succulent Karoo, "3" -Succulent Karoo. Units were "events/Myr" in all cases.

Figure 7
The potential ranges of suitable habitat for each of the seven clades of the P. tentorius species complex, combined in southern Africa across different periods ranging from the LIG to the future (up to 2070).
Environmental niche modelling was done with present bioclimatic variables on the basis of extant points of occurrence (black triangles) of the species, using the Maxent program. Information about clades and periods are given at the top-left corner of each analysis.

Figure 8
Heatmaps visualizing the changing impact contributed by each critical climatic variable (Jack-knife AUC > 0.75) against the timeline (ranging from the LIG to the future) across all groups (C1 -C7 independently and combined). The warm colours represent climatic variables with high in uence on niche modelling, whilst, cold colours represent low impact; the numbers in the legend are percentages.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.