We compare the phylogeographic structure of thirteen codistributed ectotherms including four reptiles (a snake, a legless skink and two tortoise species) and nine invertebrates (six freshwater crabs and three velvet worm species) to test the presence of congruent evolutionary histories. Phylogenies were estimated and dated using maximum likelihood and Bayesian methods with combined mitochondrial and nuclear DNA sequence datasets. All taxa demonstrated a marked east/west phylogeographic division, separated by the Cape Fold Mountain range. Phylogeographic concordance factors were calculated to assess the degree of evolutionary congruence among the study species and generally supported a shared pattern of diversification along the east/west longitudinal axis. Testing simultaneous divergence between the eastern and western phylogeographic regions indicated pseudo-congruent evolutionary histories among the study taxa, with at least three separate divergence events throughout the Mio/Plio/Pleistocene epochs. Climatic refugia were identified for each species using climatic niche modeling, demonstrating taxon-specific responses to climatic fluctuations. Climate and the Cape Fold Mountain barrier explained the highest proportion of genetic diversity in all taxa, while climate was the most significant individual abiotic variable. This study highlights the complex interactions between the Cape Fold Mountains and past climatic oscillations during the Mio/Plio/Pleistocene. The congruent east/west phylogeographic division observed in all taxa lends support to the conclusion that the longitudinal climatic gradient within the Greater Cape Floristic Region, mediated in part by the barrier to dispersal posed by the Cape Fold Mountains, plays a major role in lineage diversification and population differentiation.

In order to determine whether the major phylogenetic patterns uncovered by this study aligned with congruent climatic refugia among the study taxa, a niche modeling approach was used. Maxent 3.4.1 was used to identify persistent Pleistocene refugia based on several contemporary and paleoclimatic scenarios, namely the current (1950-2000), mid-Holocene (6 Ka), last glacial maximum (22 ka) and last interglacial (120-140 ka) time periods. All 19 bioclimatic variables were downloaded from the WorldClim database and cropped to the geographical extent of South Africa using ArcGIS (ESRI, 2012).  Correlated variables (with a Pearson’s correlation coefficient >0.7) for each species were assessed using the r package CorrPlot in RStudio 2021.09.1 and were removed, preferentially retaining variables known to influence the distribution of each taxon and the remaining limiting climatic variables were used to generate the models. A custom R script was used to specify a fixed distance buffer of 50km based on the occurrence records for model building in order to correct the models for overprediction, as well as to produce minimum convex polygons of the distribution for each taxon. ENMeval 2.0 was used to determine the optimum feature class and regularisation multiplier for model construction, after which Maxent 3.4.1was used to generate the SDMs. The models were converted from continuous probability surfaces to binary models representing either suitable or unsuitable climatic conditions based on a 10% training presence cloglog threshold using QGIS 3.24.2. Possible climatic refugia were inferred for each taxon by overlaying the binary models over the various time periods, thereby identifying areas which remained within the habitat requirements of the taxa over time. Finally, the species-specific refugia were intersected in order to determine the shared regions of climatic stability. The final maps were generated using QGIS 3.24.2.

QGIS