The role of climatic cycles and trans-Saharan migration corridors in species diversification: Biogeography of Psammophis schokari group in North Africa

Highlands, hydrographic systems and coastal areas have been hypothesised to form corridors across the hyperarid Sahara desert in North Africa, allowing dispersal and gene flow for non-xeric species. Here we aim to provide a genetic test for the trans-Saharan corridor model, and predict the location and stability of ecological-corridors, by combining phylogeography and palaeoclimatic modelling. The model was the Psammophis schokari (Schokari sand racer) group, fast-moving and widely distributed generalist colubrids occurring mostly in arid and semiarid scrublands. We combined dated phylogenies of mitochondrial and nuclear markers with palaeoclimatic modelling. For the phylogeographic analysis, we used 75 samples of P. schokari and P. aegyptius , and Bayesian and Maximum-Likelihood methods. For the ecological models, we used Maxent over the distribution of P. schokari and West African lineages. Models were projected to past conditions (mid Holocene, Last Glacial Maximum and Last InterGlacial) to infer climatic stable areas. Climatic stability was predicted to be mostly restricted to coastal areas and not spatially continuous. A putative temporary trans-Saharan corridor was identified in Eastern Sahara, with a more stable one along the Atlantic coast. Six parapatric lineages were identified within P. schokari , four occurring in North Africa. These likely diverged during the Pliocene. The Tamanraset River might have been a vicariant agent. African lineages may have experienced further subsequent diversification during the late Pleistocene. The main P.schokari refugia were probably located along the northern margins


Introduction
Numerous geological and climatic events have affected the geographic and biological diversity of North Africa in the last few million years (Le Houérou 1997;Fabre 2005). Geological events include the opening of the Mediterranean to the Atlantic 7-9 million years ago (Ma), the subsequent closure 6 Ma and re-opening 5.3 Ma (MSC, Krijgsman et al., 1999), recurrent episodes of desiccation and refilling in the Red Sea area (Girdler 1991;Bosworth et al., 2005), marine transgressions (Tawadros 2011), or the Atlas mountains uplift (de Jong, 1998). These had climatic repercussions, but the most wide-ranging climatic event was a shift from tropical to arid environments around mid-Miocene (Zachos et al., 2001) that eventually led to the appearance of the Sahara desert between 7 Ma and 2.5 Ma (Schuster et al., 2006;Swezey 2009). Arid and humid conditions then alternated during the last few million years, causing a series of expansions and contractions of climatic zones (Le Houérou 1992;Swezey 2009) that largely determined current biodiversity patterns .
While during humid phases, the hyper-arid regions were reduced and probably isolated, in arid phases the mesic species were pushed towards coastal areas and mountains (Le Houérou 1992;Messerli and Winiger, 1992). This resulted in disjoint distributions and allopatric diversification, currently best observed in mesic species, for example Mediterranean-Sahel separations (e.g. Gonçalves et al., 2012;Guillaumet et al., 2008), isolated populations in highlands (Geniez and Arnold, 2006;Metallinou et al., 2015), rock pools (Brito et al., 2011a;Vale et al., 2015) or desert-border refugia (Dobigny et al., 2013). Arid phases conversely allowed range expansions of xeric taxa Kissling et al., 2016;Leaché et al., 2017;Pook et al., 2009), later broken during humid phases (Metallinou et al., 2015;Pook et al.,  Landscape features modulate gene flow (Brown and Lomolino, 1998). Depending on climaticcycle phase and species' ecological requirements, a geographic feature can constitute a barrier to gene flow, an ecological corridor (connecting two larger similar areas), a filter bridge (or barrier, a more selective connection), or a refugium (where a species survives during unfavourable periods). Lake Chad, for instance, has been identified as a refugium (Granjon and Dobigny, 2003) or a corridor connecting the Sahel to the Tibesti (Drake et al., 2011;Dumont 1982) for mesic taxa and a vicariant barrier for xeric ones (Pook et al., 2009;Metallinou et al., 2015). Mountains, presently working as biodiversity hotspots and refugia for mesic species (Brito et al., , 2011aVale et al., 2015), likely constitute barriers for lowland or xeric species. Lastly, mountains, coastal areas and hydrographic systems can be linked, forming ecological corridors for mesic species. Several areas likely to constitute North-South oriented ecological corridors have been proposed (Dumont 1982); these geographic features are hereby referred to as trans-Saharan corridors (tS-corridors), to avoid confusion with the ecological corridor feature. They include the more transitory river drainages from central Sahara Mountains (Drake et al., 2011), or the more stable (thus possibly refugia for some species) Red Sea, Nile River or Atlantic Sahara . However, these have been proposed based on ecological/geological data and species' distributions, and no genetic-level assessments have been conducted so far.
Our goal is to provide a genetic assessment of the validity of the tS-corridor model. Psammophis schokari (Schokari sand racer) group in North Africa as a model since it is widely distributed in North Africa, occurring mostly in arid and semiarid habitats; it presents a continuous distribution along the Atlantic coast, thus apparently making use of this tS-corridor ( Fig. 1). These snakes are large and mobile animals with likely good dispersal abilities. Thus, we do not expect limited dispersal or narrow distribution to restrict its use of ecologically suitable corridors. Previous work on P. schokari (Rato et al., 2007) identified several lineages but no conclusions were drawn regarding trans-Saharan dispersal.
Here we propose to model the species' potential distribution in different climatic phases to assess the suitability of candidate tS-corridors as ecological corridors for Psammophis, thus allowing us to build clear hypothesis on the persistence of gene flow along tS-corridors for mesic species during the past climatic cycles. The general aim of this study is to assess the role of corridors in trans-Saharan dispersal, with particular focus on the Atlantic Sahara tS-corridor, in an integrative framework joining phylogeography and palaeoclimatic modelling. Using P.
schokari, we aim to answer the questions: 1) where are the areas with higher climatic stability throughout the species range and particularly West Africa and where are the potential dispersal routes across the Sahara?; 2) how is the genetic variability spatially structured?. By combining results from these two sections, we expect to find phylogeographic patterns coherent with refugia close to the Mediterranean coast and in the Saharan mountains, and tScorridors connecting them.

Sampling and study areas
The snake genus Psammophis includes 34 diurnal fast-moving species occurring mostly throughout tropical Africa, with some species reaching the Middle East and South-Central Asia (Sindaco et al., 2013, Uetz andHošek, 2016). Psammophis schokari (FORSKÅL, 1775), commonly known as Schokari sand racer, is a common colubrid occurring from West Africa to 7 India ( Fig. 1), mostly in desert and xeric scrublands, marginally in sandy habitats in dry Mediterranean zones (Kelly et al., 2008;Schleich et al., 1996;Sindaco et al., 2013). Its sister species Psammophis aegyptius MARX, 1958, commonly known as Egyptian sand snake, is a typical Saharan species renowned for inhabiting even the driest areas of the eastern Sahara (Baha El Din, 2006).
We used 68 samples of P. schokari covering a representative part of the species distribution, particularly in West Africa ( For ecological models, a total of 748 observations (Figs. A.3, A.4) were collected from fieldwork (n = 244), museum collections (n = 110) and bibliography (n = 394). These observations were used to create two datasets: (i) Global, with 629 records at five arc-minute resolution (~10x10 km); and (ii) Regional (Northwest Africa), with 379 records at 30 arc-second resolution (~1x1 km). A 150 km buffer around minimum convex polygons including each dataset was used to delimit two corresponding study areas (see 2.6). All spatial analyses were conducted in ESRI ArcGIS 10. In order to reduce bias from uneven sampling, and to geographically and environmentally homogenize datasets (Merow et al., 2013), localities were randomly removed 8 from clusters of species occurrence (e.g. Martínez-Freiría et al., 2015). The nearest-neighbour (NN) index (ArcGIS 10) was used as assessment. Low clustered distributions were obtained for both global (z-score = -19.67; NN-ratio = 0.478) and regional (z-score = -11.57; NN-ratio = 0.599) datasets, keeping 388 and 225 records, respectively (Table 2).

Climatic variables
Nineteen variables for current climatic conditions at 30 arc-second resolution (~1x1 km) were downloaded from WorldClim (www.worldclim.org; Hijmans et al., 2005). Variables were clipped to each study area (Global and Regional) and, for the Global dataset, upscaled to five arc-minutes (~10x10 km). After visual inspection, five variables were excluded due to the presence of spatial artefacts. The remaining 14 variables (Table A.2) were considered for ecological models. Bivariate correlations among the 14 variables were tested within Global and Regional datasets. We retained the same five slightly correlated (R<0.7) variables in both datasets (BIO 4, 10, 12, 14, 19; Table A upscaled to 5 arc-minute resolution for Global models and kept at original pixel size for Regional models.

DNA extraction and amplification
DNA was extracted from ethanol-preserved tissue using DNeasy Blood & Tissue Kit (Qiagen) as per manufacturer's instructions. Amplifications were performed using MyTaq™ Mix. To benefit from sequence data available from GenBank (Kelly et al., 2008;Rato et al., 2007), specimens were bi-directionally sequenced for two mitochondrial (NADH dehydrogenase subunit 4, ND4; cytochrome b, CYTB) and two nuclear (

Phylogenetic analyses
Sequences were aligned using MAFFT v7 (Katoh and Standley, 2013), with Auto option, then proofread by eye. No stop codons were found in coding genes. Each marker was individually analysed inferring independent ML trees in RAXML v8.1.21 (Stamatakis 2014), in order to detect topological incongruences suggesting sample curation errors.
For both sequence datasets, the best-fit partitioning scheme and models of molecular evolution were selected using PARTITIONFINDER v.1.1.1 (Lanfear et al., 2012)  MRBAYES was run for 10 7 generations in two independent runs sampling every 1000 generations. Parameters of sequence evolution (statefreq, revmat, shape, pinvar) were unlinked for all partitions and the overall rate (ratepr) variable among them. BEAST was run in CIPRES gateway (Miller et al., 2010) in three independent runs of 5x10 7 generations, sampling at every 5000, with unlinked substitution and clock models, under an uncorrelated lognormal relaxed clock (Drummond et al., 2006), and considering ambiguities in nuclear sequences (manually editing the xml file to UseAmbiguities=true). A constant population size coalescent tree prior (Kingman 1982) was used for dataset 1, and a Yule speciation tree prior (Yule 1925;Gernhard 2008) for dataset 2. Burn-in was determined using Tracer v1.6 (Rambaut et al., 2014), upon stabilisation of log likelihood, average standard deviation of split frequencies, and ESS for all parameters. For BEAST analysis, runs were combined with LogCombiner and a maximum credibility tree was generated with TREEANNOTATOR (both in the BEAST package).
Burn-in was determined using TRACER v1.6 (Rambaut et al., 2014). ML analyses were performed in RAXML v8.1.21 (Stamatakis 2014) through RAXMLGUI 1.5b1 (Silvestro and Michalak, 2012), with partition schemes as above and GTR+G model of sequence evolution. The program was set to perform 10 ML searches and 1000 thorough bootstrapping replicates. To test the paraphyly against the monophyly of African populations, the best tree in which Algerian/Tunisian populations are sister taxa to the rest of African and Middle East populations was compared with the alternative constrained topology in which Africa was forced monophyletic. Per-site log likelihoods were obtained with RAxML, then used to run the Shimodaira-Hasegawa (SH) (Shimodaira and Hasegawa, 1999) and Approximately-Unbiased (AU) (Shimodaira 2002) topology tests in CONSEL (Shimodaira and Hasegawa, 2001 (Stephens et al., 2001) implemented in DNASP was used to infer haplotypes for nuclear sequences. The software was run for 10 4 iterations with a thinning interval of five and a burn-in value of 1000, and repeated three times. Consistency was checked across runs by analysing haplotype frequency estimates and goodness-of-fit measures. Haplotypes were used to produce haplotype networks using NETWORK 5.0 (Fluxus-engineering.com) with medianjoining algorithm and default parameters (Bandelt et al., 1999).

Time calibration
No calibration constraints are available within the genus Psammophis or the family Psammophiidae. Therefore, to obtain estimates of the divergence times among P. schokari lineages, the dataset was expanded to span the superfamily Colubroidea plus Achrochordus, and thus make use of the calibration scheme used by Wüster et al. (2008) and Pook et al. (2009). Other authors have used a slightly different prior set (Sanders and Lee, 2008;Kelly et al., 2009), so we also dated the phylogeny according to the revised version by Sanders et al. (2010). RAG2 was excluded given that it was not available for most of the dataset.

Palaeoclimate modelling
Palaeoclimate models were computed at a Global scale, which includes the whole species distribution and aims at a species-wide average niche, and at a Regional scale, including just the West African populations and targeting local environmental species preferences. Models were generated using the maximum entropy approach with MAXENT v3.3.k . This algorithm requires only presence data, performs well comparing to other methods (Elith et al., 2006), and has been used successfully in modelling snake species distributions (e.g. Brito et al., 2011b;Martínez-Freiría et al., 2015). Thirty replicates were run with random seed and 70%/30% training/testing partition and using bootstrap with replacement. Models were 12 run with auto-features, and the area-under-the-curve (AUC) of the receiver-operating characteristics (ROC) plots was taken as measure of individual model fit (Fielding and Bell, 1997). The importance of climatic variables in explaining the species' distribution was determined by their mean percentage contribution to the models.
The individual model replicates were used to generate a mean forecast of probability of species occurrence under current conditions (Marmion et al., 2009). Standard deviation was used as indication of prediction uncertainty (e.g. Brito et al., 2011b;Martínez-Freiría et al., 2015). Individual model replicates were projected to past climatic conditions (midHol, LGM and LIG) and subjected to the same procedure. Stable climatic areas, i.e. stable potential areas of occurrence that could serve as refugia in different time periods (Carnaval et al., 2009), were calculated by averaging the probability of occurrence across all time phases.

Niche overlap test
In order to assess the potential role of ecological processes in the evolutionary history of the lineages along the Atlantic corridor, we tested niche overlap, equivalency and similarity between regional lineages (see Warren et al., 2008). Tests were based on a 2D representation of climatic space of lineages (retrieved using 14 climatic variables; Table A.   (Table A.3). Uncorrected p-distances and diversity measures for the mitochondrial markers of P. schokari and P. aegyptius lineages can be found in Table 1.

Phylogenetic relationships
Genetic identity of specimens was coherent with prior morphology-based species assignment: specimens identified as Psammophis schokari and P. aegyptius formed two reciprocally monophyletic groups in the concatenated phylogenetic ML and BI trees (Fig. 2). The nuclear DNA did not exhibit a clear separation in the network analysis, with a shared basal haplotype in c-mos, and some RAG2 schokari haplotypes being closer to aegyptius haplotypes than to the remaining schokari alleles (Fig. 2) Moroccan populations (Fig. 2).

Time calibration
Dataset 2 resulted in an alignment of 2379 positions (Table A.3). Using the alternative calibration scheme by Sanders et al. (2010) produced very similar results (Fig. A.2a).
Divergence events among the major lineages of P. schokari were all placed in the Pliocene.
African P. schokari lineages and P. aegyptius may have suffered roughly contemporaneous

Palaeoclimate models
Predictions of ecological models were robust (Table 2), identifying most of the species' occurrence with high probability (Figs. A.3-A.6). Accuracy was lower in the most arid parts of the Algerian Sahara and in portions of the Asian range. Contribution of variables changed according to modelling approach ( Table 2).
The major areas predicted with climatic stability (potential refugia) in the Global model were mostly restricted to coastal areas and were not spatially continuous (Fig. 4A). In the Regional

Discussion
The history of North Africa is rich in climatic and geological changes with significant impact on local biodiversity, making it hard to link timing of diversification with particular past events. Nevertheless, by combining geological/climatic knowledge and phylogenetic data, it is possible to approximate the processes that shaped Sahara colonization and the diversity of species in and around it. The models here presented include only climatic conditions, as no accurate geological data are available to past conditions for such large area. However, climatic variables should also be a reasonable surrogate for past habitat changes, as the distribution of desert species (particularly ectotherms) is strongly dependent on temperature and rainfall (Ward, 2009). Although geological features (e.g. soil type) also determine habitat and species distributions, we did not include such data in the models, as the large fluctuation in aeolian sand deposition regimes and scarce paleontological data caused by study area remoteness and stratigraphic discontinuities caused by erosion hamper predictions about past land-cover or soil types. Also, our paleoclimatic reconstructions are referent to the Pleistocene. During Pliocene, similar cycles seem to have occurred both globally (Lisiecki and Raymo, 2005) and in North Africa (Rohling et al., 2015), but uncertainties regarding the climatic conditions at the time and out-of-pleistocene-bounds temperature variation (Snyder, 2016) limit the extrapolations of the models beyond the Pleistocene.

Migration corridors across the Sahara
Trans-Saharan corridors (tS-corridors) can be more or less permanent (i.e. active ecological corridors), with the coastal ones expected to be more permanent . The climatically stable areas mostly fit this prediction, with the Atlantic Sahara suggested even as potential refugium area. In a simple scenario, stronger gene flow along the Atlantic than from the coast to the Algerian Mountains would then be expected, but the distribution of genetic variability seems to indicate the opposite.
The same mitochondrial haplotype was found in E-Tunisia and Hoggar (Figs. 2, 4). This  Brito et al., 2014. But see Nicolas et al., 2015). The species' broad niche can partly explain it, and this shows that the predicted level of connection of northern areas with south Algeria Mountains seems therefore sufficient to allow the species' dispersal.
According to the ecological models, the link is more likely to occur under current conditions, which suggests a very recent (re-) colonization or secondary contact. The most probable path to southern Algerian mountains seems to be through Fezzan region in Libya, but a tS-corridor from there to the Sahel seems unlikely. A putative tS-corridor seems more likely in Eastern Sahara, active under LGM-like conditions (CCSM and MIROC, Fig. A.3), linking the Mediterranean to Gilf Kebir, Jebel Uweinat, Ennedi/Borkou, and possibly Tibesti and Marra mountains. However, identification of suitable areas in the Eastern Sahara is affected by the low availability of precise observational data (Fig. 1), resulting in less clear distribution patterns.
The Atlantic coast was the sole permanent tS-corridor for mesic species well supported and it seems to have been the most climatically stable dispersal path across the Sahara. The split of the Moroccan and Mauritanian lineages seems to contradict this pattern, but these lineages might have been separated by the Tamanrasset palaeoriver basin, which opened on the Atlantic coast in the north of present-day Mauritania ( Fig. 4; Skonieczny et al., 2015). This river system was activated in the humid periods of the Sahara, which seem to have occurred cyclically since before the Pliocene (Rohling et al., 2015). With a drainage area comparable to the ones of Niger or Zambezi rivers (Vorosmarty et al., 2000), it could represent a barrier during humid phases, as well as constitute an unsuitable sandy expanse during warm ones. However, it is still possible the observed pattern is a coincidence, thus phylogeographic studies with other taxa are needed to verify this hypothesis. Continued isolation may have subsequently led to some niche divergence among lineages, although the observed pattern could not be dissociated from the differences in local climatic conditions (Fig. A.7). Predicted 17 conditions predict the narrowest suitable area and still allow species persistence along the coast seems to indicate that, at least climatically, some gene flow is always possible. Individual Moroccan-Mauritanian models predicted that lineage contact is most likely under LGM-like conditions (Figs. A.5, A.6), similar to the Eastern Sahara corridor. The Red Sea coast was expected to be clearly represented in the palaeomodels, given that the topography allows altitudinal range shifts according to climate fluctuations (Messerli and Winiger, 1992). Contrary to previously hypothesised, the Red Sea tS-corridor seems not permanent for mesic species, especially during hot and dry phases similar to LIG (Fig. A.3). Nevertheless, pixel size could be constraining model sensibility.
All considered, it seems that the Atlantic coast works as a corridor when considering the climatic cycles altogether, but certain portions (e.g. Tamanrasset river basin) may work as filters depending on the climatic phase. For example, the predicted climatic suitability of northern costal Mauritania was not mirrored by availability of observations/samples, despite considerable local sampling efforts during the last decade (Sow et al., 2014), and is probably

18
The status of sister species is supported for P. schokari and P. aegyptius, in agreement with Rato et al. (2007), although differentiation in the nuclear markers was not complete (Fig. 2).
Although c-mos exhibits a typical pattern of shared ancestral polymorphism where the most common haplotype in the two species is still the ancestral haplotype, RAG2 exhibits clear separation between species, so the nuclear data do not contradict reproductive isolation between these two species. The genetic distance found between the five P. schokari mtDNA lineages, compared in Table 1, was substantial (4.5-7.1%), and while not uncommon in reptiles, similarly deep evolutionary lineages were used to request taxonomy revisions in other species (e.g. Vipera latastei; Velo-Antón et al., 2012). Still, higher resolution nuclear markers or increased genomic coverage (e.g. Velo-Antón et al., 2014) would be needed to support it and to resolve the phylogenetic relations among lineages. Additional sampling would also be necessary to define lineage distributional limits and assess gene flow in contact zones.

Biogeography and diversification
Psammophis schokari and P. aegyptius seem to have diverged around 10 Ma (Fig. 3). The latter exists only in Africa, but P. schokari ranges from the Atlantic coast to India, and the earliestsplitting lineage identified here is found in Oman. If the African diversity is indeed embedded within Asian diversity, an Asian origin of P. schokari, for instance through vicariance from African populations (the ancestor of P. aegyptius), is plausible. However, clarification is dependent on covering sampling gaps, since an Oman lineage migration across the Red Sea cannot be ruled out, as exemplified by other faunal exchanges between Africa and Arabia during the progressive aridification of late Miocene (Metallinou et al., 2012;Smíd et al., 2013;Tamar et al., 2016a). Ecological divergence into new more arid habitats could also have been involved, as suggested for other reptiles (Carranza et al., , 2002Metallinou et al., 2012;Pook et al., 2009;Tamar et al., 2016b).
The successive divergences among major lineages after the Messinian Salinity Crisis (Fig. 3 mirror the same diversification patterns in the Late Miocene and Pliocene found in many vertebrate and reptile groups Geniez and Arnold, 2006;Gonçalves et al., 2012;Wagner et al., 2011), and have been associated with the original onset of the Sahara desert. Lineage distributions are broadly coherent with the predicted climatically stable areas, which could suggest a vicariant effect of climate (Fig. 4). However, our models are based on Pleistocene conditions and do not cover the Pliocene. Still, dust flux records indicate that notwithstanding a slight increase in the last Ma, humid-arid cycles have been occurring similarly since 5 Ma (Trauth et al., 2009).
The roughly contemporaneous and independent divergences of North African sub-lineages after mid-Pleistocene (Fig. 3) may be linked to persistence and subsequent isolation in Quaternary climate refugia during humid-dry cycles. In Morocco for example, the sub-lineages are separated by mountain ranges (Fig. 2), which agrees with a break in climatically stable areas found in Global, Regional, and Moroccan models (Figs. 4;). This suggests a role for climate in the group diversification and matches the refugia patterns identified in the region (Husemann et al., 2014;Martínez-Freiría et al., 2017). Similarly, the decreased suitability predicted for north-eastern Algeria/Tunisia could be linked with the divergence of the two sub-lineages, a pattern also found in other taxa (Guiller and Madec, 2010;Husemann et al., 2014;Nicolas et al., 2015). In Mauritania, sporadic decreases in climatic connectivity or disruptive effects of rivers during humid phases (e.g. Dobigny et al., 2005) could help explain the distribution of the sub-lineages. However, given most of the phylogenetic signal is from mtDNA, effects of coalescence cannot be ruled out without employing multiple nuclear markers. The observed haplogroups in P. aegyptius are probably due to considerable geographic distance (isolation by distance) between sample clusters (Egypt and Niger).

Conclusions and future research
The work here presented contributes to a better understanding of how Saharan mesic landscape in gene-flow dynamics at local scales is also crucial. In addition to vicariance, ecological adaptation to different gradients of aridity (e.g. Tarentola along the Atlantic coast; Carranza et al., 2002) or habitat specialization (e.g. sympatric speciation in Jaculus; Boratyński et al., 2014) can also lead to species diversification. We contributed to a better understanding of the history of P. schokari in North Africa, but still further sampling and additional nuclear markers are needed to understand the history of the species along its full distribution. Of particular importance are the role of the interplay of Arabia Peninsula and Africa, and the possibilities of secondary contacts between species and lineages, as well as their evolutionary and taxonomic status.

Acknowledgements
Authors acknowledge BIODESERTS group members for their assistance during fieldwork. P.

Figure 2
Phylogeographic relationships among the lineages of P. schokari and P. aegyptius, and their position within the genus (RAxML and MRBAYES, dataset 1, all markers). The colours in the tree, map, and haplotype networks represent the lineages identified in the phylogenetic tree analysis. The RAXML tree is represented; posterior probability (pp) and bootstrap values (bss) are indicated as pp/bss; asterisks represent values of 100; values below 95%pp or 70bss are not shown, except for the node in light grey. Haplotype networks were calculated using the Median-Joining algorithm in NETWORK; small black dots represent mutation steps; red diamonds represent predicted missing haplotypes.

Figure 4
Stable climatic areas for P. schokari derived from four modelling frameworks: (A) Global, all the species range; (B) Regional, the NW-African lineages; (C) Moroccan lineage only; and (D) Mauritanian lineage only. Warmer colours depict areas with higher stability. Tamanrasset paleoriver (Skonieczny et al., 2015) is depicted in B-D.