The role of Sahara highlands in the diversification and desert colonization of the Bosc's fringe‐toed lizard

The biogeographic history of the Sahara‐Sahel desert is tightly linked to its extreme and fluctuating palaeoclimate and diverse topography. For the mesic species inhabiting the region, coastal areas and the Nile Valley are perceived as the main pathways to disperse through desert habitats, but past connections may have also occurred throughout currently isolated mountain regions. Herein, we test the trans‐Sahara mountain corridor hypothesis (i.e., mesic connectivity across Central Sahara highlands) and its role in the diversification of a small terrestrial vertebrate.


| INTRODUC TION
A central goal in biogeography is to understand how past climate and geography shaped the distribution and evolution of current species (Lomolino et al., 2006). In this context, deserts have historically attracted little scientific attention, based on the misconception that they are homogeneous landscapes of depauperate biodiversity (Durant et al., 2012). However, desert biotas provide excellent study systems to explore the role of climate variability on biodiversity patterns, owing to the tight links between species distribution and an extreme climate and the quick and wide-ranging environmental shifts that occurred in arid biomes through the past (Ward, 2016). In North Africa, the Sahara-Sahel warm desert ( Figure 1) has experienced strong climatic oscillations since the Late-Miocene aridity onset (Zhang et al., 2014), driven by variations in the African monsoon belt (Skonieczny et al., 2019). These events triggered cyclic land cover shifts between hyperarid and savanna-like conditions which greatly altered the desert range limits (Palchan & Torfstein, 2019), shaping biodiversity distribution and genetic structure through an interplay of range shifts, vicariance and adaptation processes (reviewed by Brito et al., 2014). As a result, the contemporary biodiversity of the Sahara-Sahel is characterized by a mixture of wide-ranging species (e.g., Gonçalves, Martínez-Freiría, et al., 2018) as well as species with narrow fragmented ranges (e.g., Gonçalves, Pereira, et al., 2018), the occurrence of locally adapted ecotypes (e.g., Leite et al., 2015), and high rates of endemism (Brito et al., 2016). This diversity has also been heavily affected by recent population extirpations induced by increasing aridity since the middle Holocene (Drake et al., 2011) and overhunting of large vertebrates during the last century (Durant et al., 2014).
Climate-induced range shifts are dependent on species' tolerance to aridity . For mesic taxa (i.e. occurring within Aridity Index of 0.05-0.50; Ward, 2016), the present-day Sahara Desert is a patchwork of unsuitable hyper-arid regions and isolated patches of suitable habitats, where distribution patterns reflect historical range retractions to mild, stable climate areas. These areas are primarily located in micro-climatic refugia in mountain regions (Martínez-Freiría et al., 2017), around mountain rock pools (Vale et al., 2015), in oases (Shaibi & Moritz, 2010) and along desert fringes . Mountain regions in particular condense more ecological variety per geographic space and retain more humidity than surrounding areas, and thus are key for generating and maintaining diversity (Perrigo et al., 2020;Rahbek et al., 2019), driving adaptive and allopatric diversification in isolated populations (Garcia-Porta et al., 2017). For example, several ectotherms show high rates of intraspecific diversity around the Atlas, at the north-western edge of the Sahara (e.g., Martínez-Freiría et al., 2017). The Central Sahara highlands ( Figure 1) constitute the main mesic patches among massive hyper-arid extensions and potentially play a key role in the persistence refugia, suggesting the existence of alternative trans-Sahara dispersal routes to the putative coastal and Nile corridors.

K E Y W O R D S
arid mountains, biodiversity refugium, climatic cycles, cryptic diversity, dispersal corridor, historical biogeography, North Africa, phylogeography, Pliocene-Pleistocene F I G U R E 1 Study area, species and sampling localities used in this study. The white polygon depicts the current distribution of Acanthodactylus boskianus (modified from Sindaco & Jeremčenko, 2008). White circles and black dotes represent, respectively, observations and DNA samples of A. boskianus (specimen photograph inserted at the bottom left). Blue circles and black dots represent observations and DNA samples for A. schreiberi syriacus, restricted to North Sinai and the Levant. Contemporary and historical geographical features mentioned along the text are also indicated (mountains in bold; water bodies in italics). The map inserted at the top left depicts the location of Sahara highlands (areas >600 m) and the trans-Sahara corridor scheme (arrows). Projection: WGS84 of desert endangered fauna (Brito et al., 2016;Elsen & Tingley, 2015).
However, knowledge on regional biodiversity is generally absent or outdated, as remoteness, armed conflict and socio-political instability have continuously hampered local research . The scarce molecular studies available have uncovered lineages endemic to southern Algerian [e.g., Uromastyx alfredschmidti (Tamar et al., 2017)] and Mauritanian massifs [e.g., Pristurus adrarensis (Geniez & Arnold, 2006)] or mountain meta-population systems in southern Mauritania [e.g., Crocodylus suchus (Velo-Antón et al., 2014)], advancing the presumably role of these mountains as important diversification hotspots. Nevertheless, potential diversity reflecting past evolutionary dynamics remains unexplored across most of the region.
Permanent land connections since the Middle Miocene prompted several faunal exchanges between North Africa and Arabia (Tejero-Cicuéndez et al., 2021), followed by climate-dependent radiations across arid environments (e.g., Šmíd et al., 2020). In North Africa, the main routes allowing for trans-Sahara dispersals of mesic species are restricted to the peripheral coasts, along the Atlantic Sahara  and the Red Sea (Tamar, Scholz, et al., 2016), and to the Nile River basin (Drake et al., 2011). However, the location of Central Sahara highlands makes them potential pivotal hubs of putative inland corridors, active during humid Sahara phases ). This corridor model presumably integrates mountain regions together with expanded hydrological systems and adjacent savanna-like habitats (Figure 1; Drake et al., 2011). The presence of Mediterranean (e.g., Pelophylax saharicus) and Tropical (e.g.,

Hoplobatrachus occipitalis) relict amphibian populations in Central
Sahara supports the existence of recent humid connections with the desert fringes (Drake et al., 2011). The scant molecular data available on mesic taxa recovers a mixed pattern, with reports of high connectivity with the Mediterranean [Pelophylax saharicus (Nicolas et al., 2015)], the Red Sea [Psammophis aegyptius ] and the Sahel [Ptyodactylus togoensis (Metallinou et al., 2015)], but also isolated species [Agama tassiliensis (Gonçalves, Pereira, et al., 2018)] and localized barrier effects (e.g., Tamanrasset palaeoriver ]. However, these are only partial pictures, and obtaining a good framework for trans-Sahara connectivity requires broad scale studies on species (or complexes) with region-wide distribution.
Although it is a lowland species, it occurs in mountain foothills, presumably making use of hypothetical trans-Sahara mountain routes.
The species falls on the more arid side of the mesic spectrum, making it a good candidate to detect more intermittent and less obvious corridors, refugia and diversification centres compared to truly mesic species, which are more easily extirpated. Previous data suggest that this taxon constitutes a species complex (Tamar et al., 2014).
The main aim of this study is to test whether the Sahara highlands acted as diversification nodes and climatic refugia facilitating trans-Sahara connectivity in A. boskianus, by integrating genetic analyses and palaeo-projections of species' contemporary climate-niche models. Specifically, we evaluate: (i) the existence of intraspecific lineages; (ii) how they are spatially structured; (iii) the location of stable and more intermittent areas of suitable climate (refugia and corridors, respectively); and iv) the correlation between genetic diversity, topography and palaeoclimate oscillations. We expect that: (1) refugia will be located along the desert periphery and around mountains; (2) mountains will work as diversification hotspots, hosting high genetic diversity and endemic lineages derived from interplaying topography and long-term climate-driven population persistence in refugia; (3) permanent and temporary North-South corridors will be located, respectively, along the coasts and throughout Central Sahara highlands; (4) some lineages will spread across corridors as a result of gene flow during recent suitable periods. Our integrative analysis should provide novel insights into the historical biogeography of the Sahara-Sahel, particularly the more intermittent dispersal routes used by nonxeric species to overcome the desert arid barriers.

| Study area, organism and data collection
Our study area encompasses the Sahara, Sahel and Arabia arid belts Acanthodactylus boskianus is the most broadly distributed species within a highly specious genus (>40 currently recognized species), expanding throughout most of Sahara, Sahel and Sudanian Africa, the Arabian Peninsula, and the Middle East to the Cizre region, near the Turkish-Syrian border, and the Zagros Mountains, in West Iran ( Figure 1; Sindaco & Jeremčenko, 2008;Uetz et al., 2021). It is one of the most conspicuous inhabitants of arid and semi-arid environments in these regions, being associated to open grounds with scattered shrubby vegetation or vegetated wadis, typically on hard to soft substratum but absent from truly sandy habitats, where it is replaced by congeneric species, and avoids both savanna grasslands and hyper-arid areas (Trape et al., 2012). The species shows extremely high morphological variability, even within a markedly variable genus (Arnold, 1983;Salvador, 1982), and high degree of genetic diversity (Tamar et al., 2014).
Previous studies unveiled intricate intraspecific relationships and geographic clusters (Heidari et al., 2014;Khannoon et al., 2013;Poulakakis et al., 2013;Psonis et al., 2016;Tamar et al., 2014;. However, irregular spatial sampling limited insights into the species' biogeographic history, particularly in North Africa. For phylogenetic analyses (see Section 2.3), we used a set of 86 georeferenced tissue samples mostly from North Africa, which were collected by the authors between 2002 and 2018 or acquired from museum collections and collaborators. We also retrieved GenBank sequences for additional 106 georeferenced samples, which complemented the species' Middle-East range. Our final genetic dataset ( Figure 1; Table S2.1 in Appendix S2) included the four recognized African subspecies of A. boskianus (A. b. asper, A. b. boskianus, A. b. nigeriensis, and A. b. kattensis; no data were available for A. b. euphraticus, from Iraq-Iran) and covers a representative part of the species' range, filling key geographic gaps in North Africa (e.g., Algeria, western Libya, Atlantic Sahara). Previous studies (Tamar et al., 2014) reported A. boskianus as paraphyletic, with A. schreiberi syriacus from Sinai and the Levant embedded within the species' diversity. Thereby, we added samples of A. s. syriacus (n = 13) to the dataset in order to prevent potential implications of unresolved taxonomy (i.e. whether the latter constitutes a separate taxon or an ecotype of the former).
For climate niche models (see Section 2.5), we considered the monophyletic unit composed of A. boskianus and A. s. syriacus, excluding the samples of an evolutionarily independent lineage from Syria-Jordan. Although this lineage has been attributed to A. boskianus, its relationship with the main species' clade is unclear and strong divergence is confirmed by nuclear data, implying that it may constitute a distinct species (Tamar et al., 2014). Moreover, the scarcity of available records (n = 5) impedes robust niche models for this lineage/species. In total, our observations dataset ( Figure 1; Table S2.1) included 507 records gathered from fieldwork (n = 220), museums (n = 75) and bibliography (n = 212).

| DNA extraction, amplification and sequencing
Extractions of genomic DNA were performed from ethanol-preserved tail tips or tissue samples using QIAGEN's EasySpin Kit and stored at −20°C prior to amplification. A first screening of the species' genetic variability and structure was performed by amplifying and sequencing a short fragment of the non-coding, highly conserved mitochondrial marker 12S rRNA (12S; 366 base pairs) for all samples, and combined with 12S sequences retrieved from GenBank. Main 12S haplotype groups (haplogroups) were inferred based on an haplotype network built with Tcs v.1.2.1 (Clement et al., 2000), using 95% parsimony threshold, which resulted in 23 haplogroups separated by at least three mutation steps (Figure 2b). A subset of 44 samples was selected to represent the distribution area and all 12S haplogroups, except for three from the Middle East that were detected only in GenBank sequences (Table S2. Table S2.2 for PCR conditions and primers). PCR products were cleaned with ExoSAP and the sequencing for the forward strand was outsourced to Genewiz (https://www.genew iz.com/). Chromatograms were verified in Geneious Pro v.4.8.5 (Drummond et al., 2010) and sequences were aligned using MaffT v.7 (Katoh & Standley, 2013). Gene alignments were proofread and the absence of stop codons in coding regions was verified. Topological incongruences between markers were primarily checked by inferring Neighbour-Joining trees for each marker in MeGa v.6 (Tamura et al., 2013).

| Phylogenetic analysis and estimation of divergence times
Sequences were concatenated into mtDNA (n = 40; 2169 bp), nuDNA (n = 34; 2169 bp) and mt-nuDNA (n = 34; 4338 bp) datasets, after retaining only samples that were successfully sequenced for at least 75% of the target markers (Table S2.1). Samples of A. opheodurus and A. scutellatus were used as outgroups. The best-fit partitioning schemes and models of molecular evolution were inferred using ParTiTionfinder  (Miller et al., 2011) in two independent runs of 10 8 generations, sampling at every 10 4 , with unlinked substitution and clock models (Drummond et al., 2006), a constant population size coalescent tree prior (Kingman, 1982), and considering ambiguities in nuDNA partitions (manually editing the xml file to Ambiguities = true). Detailed settings of phylogenetic analyses are given in Table S2.3. The convergence of chains and Effective Sample Sizes >200 for all parameters was verified using Tracer v.1.7.1 . Log and tree files of the independent runs were combined using LoGcoMBiner and the subsequent maximum clade credibility summary trees with posterior probabilities for each node, using mean values, were obtained using TreeannoTaTor (both from the BeasT package). The resulting trees were visualized and edited with fiGTree v.1.4.3 (Rambaut, 2016).
Nodes were considered supported if they received a posterior probability ≥0.95. Haplotype networks were constructed for the mtDNA marker with the highest sampling coverage (12S; n = 204; see section 2.2) and for each of the four nuclear genes using Tcs v.1.21 (Clement et al., 2000), with a 95% parsimony threshold, and visualized using TcsBu (Santos et al., 2015). The nuDNA haplotype networks were built using the most likely haplotypes calculated using the PHASE implementation in dnasP v.5.10.01 (Librado & Rozas, 2009), run five times with 10 4 iterations and 10 3 of burn-in.
Main mtDNA lineages were defined as reciprocally monophyletic clades in the mtDNA tree that split before the Pliocene. The 12S haplogroups represented in the mtDNA tree formed distinct, wellsupported phylogenetic groups, which were defined as mtDNA sublineages. The single exception is the haplogroup formed by all A. s. syriacus samples, which was integrated within the sub-lineage Sinai

| Ancestral area reconstruction
To infer the biogeographic history of A. boskianus, we performed a discreet ancestral area reconstruction with the Bayesian Stochastic Search Variable Selection (BSSVS; Lemey et al., 2009)

, implemented in
BeasT v.1.10.4 . We targeted the time-calibrated mt-nuDNA dataset, using previous settings (Table S2.3). Six regions were proposed based on general biogeographic assumptions, species' current distribution and hypothetical ancestral dispersal routes: (1) the Middle East, including Sinai, where the species was likely originated ; (2)

| Modelling of historical climatic niche stability
Initial occurrence records were tested for spatial clustering to account for sampling bias and avoid model overfitting (Boria et al., 2014) using the ecospat.mantel.correlogram function of 'ecospat' R package (Di Cola et al., 2017), and thinned accordingly with 'spThin' R pack-  Table S2.1). For modelling purposes, nineteen climatic variables for current conditions were downloaded from WorldClim at 2.5 arcminute resolution (Fick & Hijmans, 2017) and cropped to the study area. Eight variables were initially excluded due to the presence of spatial artefacts (Bio8, Bio9, Bio15, Bio18, Bio19) or low ecological relevance (Bio2, Bio3, Bio4). The optimal set of low-correlated, highly explanatory variables (bio12, bio11, bio5, bio17) was inferred using the 'SDMtune' R package (Vignali et al., 2020). For palaeo-projections, two time-scales were used: Late Pleistocene-Holocene (last ∼120 kya, where more data on North African palaeogeography is available) and Plio-Pleistocene (last 5.4 Mya, almost since Sahara climatic cycles began; Zhang et al., 2014). In the first time-scale, ensembles of the previously selected WorldClim variables were produced for the mid-  (Brown et al., 2020;Gamisch, 2020), results from both palaeoclimate datasets were compared to verify potential inconsistencies.
The contemporary realized climatic niche (Soberón & Peterson, 2005) of A. boskianus was modelled using the machine-based learning algorithm MaxenT (maximum entropy; Phillips et al., 2006). Despite the choice of modelling technique affecting the results (Thuiller et al., 2019), MaxenT performs well compared to other methods (Elith et al., 2006), requiring only presence data, and has been successfully implemented in other reptiles (e.g., Gonçalves, Martínez-Freiría, et al., 2018). The 'SDMtune' R package was run to optimize model settings (Vignali et al., 2020). One hundred model replicates were run, where presence data for each replicate were selected randomly by bootstrapping, and using 10,000 background sample points. The back-

| Geographic distribution of genetic diversity
To identify environmental drivers of genetic diversity, we assessed the relationship of nucleotide diversity with climatic stability and topography, considering discrete regions throughout the study area. These regions were defined using same-area circular buffers (diameter 550 km), each one including at least five genotyped localities. For the buffers' distribution, we opted to use a criterion of maximum representation rather than random location, after verifying that a random-placing approach resulted in a very low final number.

| Phylogenetic analyses and cyto-nuclear discordances
The BI mtDNA tree revealed four main well-supported, spatially- Nevertheless, it is worth noticing that Egypt 1 and Egypt 2 share alleles in nuclear MC1R and FGB.

| Ancestral area reconstruction
The ancestral biogeography of A. boskianus recovers a most-likely Middle East origin (77% probability; Figure 4)

| Climatic suitability stability
Climate-niche models had an overall high accuracy (AUC = 0.77; Figure S3.6), but lower in the Sahel ( Figure S3.2). Bio12 (annual precipitation) was the variable that mostly explained species' distribution (Table S2.5). Low mean annual precipitation and temperature correlated with range expansions ( Figure S3.7), resulting in higher overall connectivity and full North-South corridors across Central Sahara in LGM projections ( Figure S3.4). Partial trans-Sahara connections, between Central Sahara and the Mediterranean, were recovered for midHol. The projection to LIG conditions generally matched with the present, with complete isolation of suitable climate patches in Central Sahara. A latitudinal corridor along the Sahel was not evidenced for any of these periods, with only occasional connections between central and eastern suitable fragments.
In North Africa, major historical refugia were continuously distributed along the coasts and the southern foothills of the Atlas, but also in isolated Central Sahara mountain patches ( Figure 5). In the Middle East, refugia also matched with areas of higher elevation, being located in the Sinai and a large part of the Arabian Peninsula.
Continuous Mediterranean-Sahel connections were found along coastal refugia (particularly the Atlantic Sahara) and, with a lower frequency, inland across eastern Sahara. Partial trans-Sahara connections, between Central Sahara and the Mediterranean, were inferred across the western side of the desert. For historical stability, the occasional connectivity between Central Sahara and the Sahel was suggested only by Worldclim ( Figure 5).

| Drivers of genetic diversity
The genetic diversity of A. boskianus was significantly affected by absolute altitude and, in lower terms, Plio-Pleistocene climate-niche stability (Table 1)

| DISCUSSION
The integration of genetic and ecological tools has helped to understand which ecological mechanisms modulated the ancient biogeography and diversification of North African desert fauna Martínez-Freiría et al., 2017;Nicolas et al., 2017;Velo-Antón et al., 2018). Nevertheless, several biogeographic questions at this broad geographic scale remain unresolved.
Here, we identify potential mesic dispersal routes, climatic refugia and diversification hotspots across the Plio-Pleistocene climatic cycles in the Sahara-Sahel and Arabian arid regions. Particularly, we discuss the historical role of the barely known Central Sahara mountainous refugia in past trans-Sahara faunal exchanges, which comprise an alternative radiation scheme to the putative coastal corridors, and provide further genetic grounds for the role of arid mountains as diversification centres. Our phylogeographic results largely complement previous studies on A. boskianus (Harris & Arnold, 2000;Khannoon et al., 2013;Psonis et al., 2016;Tamar et al., 2014;. Although our molecular dataset is the most comprehensive one to date, some sampling gaps persist, particularly in Mali, Chad, coastal Libya and south-western Egypt. Nevertheless, the genetic patterns discussed here can certainly be used as baseline for further hypothesis testing in the region, particularly on mesic taxa. Our models assume a conserved ecological niche and are based only on climate, due to the lack of accurate spatially explicit geological and land cover data for past conditions. However, climatic variables should constitute a reasonable surrogate for past land cover variations, due to the feedback mechanisms between rainfall reduction and shifts in vegetation communities (Ward, 2016). The fit of models was reasonably high, considering that the target is a generalist species ranging over a very large area. Populations falling outside average ecological optimums and the inclusion of regions that are not accessible due to barriers or historical effects in the model training (e.g., northwest of Atlas Mountains, Horn of Africa, east of Zagros Mountains) probably had a negative effect on model accuracy , but these caveats are hardly avoidable considering the spatial scale of the analysis.

| Historical biogeography and diversification
Our ancestral biogeography inference locates the origin of A. boskianus in Arabia and suggests two successful independent Pliocene expansions into Africa (Figure 4), in accordance with previous findings . These dispersals probably constitute back-to-Africa events (Tejero-Cicuéndez et al., 2021), since the genus is most likely of African origin . Extant representatives of both dispersal events currently inhabit the Lower Nile region, and are hence harbouring high nuclear genetic diversity ( Figure 3a). The earlier colonization event resulted in an extensive radiation across the Sahara, with diversification at the western side of the desert taking place in the last four million years, as illustrated by local lineages crown age (Figure 2b). In contrast, the lineage involved in the second colonization event remained restricted to the Lower Nile region. Inter-lineage competition likely explains the different outreach of each dispersal event, impeding the latter to expand throughout already filled ecological niches, but potential breeding barriers need to be explored to verify this hypothesis. Alternative explanations include environmental opportunity (different local conditions at each colonization time) or higher adaptive potential in the first colonizing group leading to a higher ability to disperse.
Climatic cycles and sand desert expansions, probably starting at least 7 Mya in eastern Sahara (Zhang et al., 2014), are consistent with our estimated crown age for A. boskianus and its rapid Plio-Pleistocene diversification (Figure 2b), prompting repeated range

TA B L E 1 GLM of the effect of
Plio-Pleistocene climatic stability and topography on the nucleotide diversity of Acanthodactylus boskianus, based on discrete spatial units across the study area (n = 14) contractions/expansions and subsequent episodes of population isolation/re-connectivity ( Figures S3.7). Although aridity presumably induces diversification in desert mesic species , our palaeo-range shifts suggest that increasing mean annual precipitation across the Sahara-Sahel were followed by population fragmentation in A. boskianus. This pattern probably results from the particular habitat requirements of A. boskianus (a desert species absent from both sandy habitats and grasslands) in comparison to other Sahara mesic species, typically less tolerant to aridity [e.g., Agama sp. (Gonçalves, Pereira, et al., 2018)]. Accordingly, the expansion of savanna habitats during humid Sahara phases likely fragmented populations into pockets of semi-arid, rocky and less vegetated places. In turn, increasing aridity and lower temperatures during dry periods allowed expansions throughout suitable habitats while avoiding sand dune areas.

| Dispersal routes
Pelophylax saharicus (Nicolas et al., 2015)], providing further evidence for the existence of historic mesic connections. Secondly, putative corridors linking the area south of the Atlas Mountains to the Hoggar and Tassili n'Ajjer across central Algeria are suggested by the palaeomodels but not well-supported by genetics, as occurred in the above studies. The potential barrier effect of palaeohydrological features (e.g., Tamanrasset palaeoriver; Ahnet-Mouydr palaeolake; Figure 1) and especially vast dune fields (e.g., Algerian Grand Ergs) likely resulted in the lack of genetic connectivity between areas linked by occasional suitable climate. Thirdly, the presence of continuous areas of suitable climate in the past between the Tibesti and coastal Libya was inferred by the models, but gaps in genetic sampling and incomplete knowledge on the species' distribution prevents testing their effects on past population connectivity. In the Oscillayers climate-niche stability inference, all the above-mentioned Central Sahara mountains were predicted to remain isolated from the Sahel, but this pattern was not fully corroborated by the LGM Worldclim projection ( Figure S3.4).
In addition, ancestral biogeography supported the Sahel origin of extant populations around the Aïr Mountains (Figure 4), a similar pattern also found in Agama boueti (Gonçalves, Pereira, et al., 2018). Further sampling from the Malian and Nigerien Sahel are needed to verify this and other putative corridors.
An intermittent inland Mediterranean-Sahel corridor across eastern Sahara, previously suggested by data on Psammophis aegyptius , is here corroborated by intermittent climatic suitability and shallow genetic differentiation between the Nigerien and Egyptian populations (Figures 2b and 5).
Mountain regions in north-eastern Chad, north-western Sudan and south-western Egypt [i.e., Ennedi, Marrah and Gilf Kebir (Figure 1) (Mouline et al., 2008), V. niloticus (Dowell et al., 2016)]. Further sampling in Libya, Sudan and Chad, combined with the use of fastevolving genetic markers, is necessary to verify potential gene flow and dispersal routes across eastern Sahara.
The hypothesized North-South coastal corridors  were supported in A. boskianus by continuous areas of climatic stability ( Figure 5). In addition, a single sub-lineage was found around the Red Sea (Figure 2a), suggesting genetic connectivity and a matching pattern with other mesic taxa [e.g., Pseudotrapelus chlodnickii (Tamar, Scholz, et al., 2016), Ptyodactylus hasselquistii (Metallinou et al., 2015)].
Nevertheless, two allopatric and nearly endemic sub-lineages were found inhabiting the Atlantic Sahara. It is possible that the Tamanrasset palaeoriver (cyclically active during humid phases; Skonieczny et al., 2015) separated populations of a Pleistocene common ancestor, a similar barrier effect detected in P. schokari  and Acanthodactylus aureus (Velo-Antón et al., 2018). We found genetic support for recent connectivity between the Atlantic Sahara and the Sahel but restricted to southern Mauritania ( Figure 4), with no evidence of mtDNA gene flow from central Sahel populations. The lack of historical climatic suitability and the barrier effect of sand dunes between Mauritania and Mali (Metallinou et al., 2015) or the Niger River likely precluded connectivity between these two regions. The apparent barrier effect of the Tamanrasset and Niger water bodies contrasts with the unexpected permeability of the Nile River (Khannoon et al., 2013), since local lineages can be found on its both sides. This does not mean that lizards can currently cross the wide Nile River, but that the barrier was occasionally permeable as a result of changes in river bed or historical desiccations [e.g., reduction of the Nile system to intermittent seasonal flow 15 kya (Lamb et al., 2007)]. The role of the Nile River in the biography of desert ectotherms is still uncertain, with available examples of acting as a major biogeographic barrier (e.g., Machado et al., 2020) and as a permeable one (e.g., Tamar, Scholz, et al., 2016).
The unresolved tree topologies limit insights into the biogeographic history of the north-western side of the Sahara Desert.
Despite this, our ancestral area reconstruction suggests a Pliocene, West-Sahara origin of the extant Mediterranean lineages (Figure 4).
There were two likely expansion routes: one south of the Atlas favoured by climatic suitability, subjected to historical climateinduced population expansion/isolation events ( Figure S3.7); and an intermittent one along the southern margin of the Tamanrasset palaeoriver (Skonieczny et al., 2015), which is supported by genetic connectivity from the Atlantic coast to central Algeria (Figure 2a) but not by the models (Figure 5). Interestingly, assuming that the Mediterranean was colonized later than West Sahara rules out the Maghreb coast, which is climatically stable and an expected corridor, from the basal radiation scheme of A. boskianus from the East to the West side of the desert. Increasing surface water flow [e.g., Sahabi and Wadi Nashu palaeorivers (Drake et al., 2008)] and sand dunes (e.g., Libyan Desert; Figure 1) may have precluded connectivity across northern Libya, restricting dispersals to permeable, inland routes around mountains. Additional sampling in Libya and Chad and resolving phylogenetic relationships are needed to fully understand the original colonization scheme of the Sahara.

| Systematic implications
Deep mtDNA ( Figure 2b) and nuDNA (Figure 3b) divergences and lack of nuclear allele sharing between some of the A. boskianus lineages indicate a species complex in need of systematic revision.
Nine nuDNA units were identified, each potentially reflecting cryptic species. The lack of nuDNA data for three mtDNA sub-lineages precludes a complete overview of the cryptic diversity contained in our dataset. The presence of undescribed diversity in remote areas is frequent (Ficetola et al., 2013), and several recently described taxa are distributed in mountain regions of North Africa (Gonçalves & Brito, 2019;Miralles et al., 2020) and Arabia (Carranza & Arnold, 2012; Simó-Riudalbas et al., 2017). Extensive molecular and morphological analyses, especially targeting contact zones between lineages to assess their level of reproductive isolation, are needed to understand and describe the cryptic diversity within the A. boskianus complex.

| CONCLUSIONS
Our integrative study provides novel genetic and ecological insights into the biogeographic history of North Africa, laying the ground for future hypothesis testing. We highlight the role of the still poorly known Central Sahara highlands as isolated mesic refugia and skyisland diversification centres. We also support coastal areas, particularly the Atlantic Sahara, as recurrent North-South climatic corridors, while eastern Sahara acted as an intermittent connection throughout. Comprehensive population genetic/genomic data are needed to investigate reproductive isolation and resolve unclear lineage relationships. Additional sampling in remote areas of Mali, Libya, Chad and Sudan will provide the baseline to fully understand the radiation scheme of A. boskianus throughout the Sahara Desert.
Lastly, our work underlines the importance of considering Sahara highlands in conservation planning, emphasizing their critical role as centres of endemism and climatic refugia for highly vulnerable desert ectotherms to face global warming.

ACKNOWLEDG EMENTS
We thank BIODESERTS group members for their assistance dur-

CONFLIC T OF INTERE S T
The authors declare that they have no conflict of interest.

DATA AVA I L A B I LIT Y S TATEM ENT
GenBank accession numbers and sampling localities are indicated in