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

André Vicente Liz, Dennis Rödder, Duarte Vasconcelos Gonçalves, Guillermo

Velo-antón, Miguel Fonseca, Philippe Geniez, Pierre-andré Crochet, José

Carlos Brito

► To cite this version:

André Vicente Liz, Dennis Rödder, Duarte Vasconcelos Gonçalves, Guillermo Velo-antón, Miguel Fonseca, et al.. The role of Sahara highlands in the diversification and desert colonization of the Bosc's fringe-toed lizard. Journal of Biogeography, 2021, 48 (11), pp.2891-2906. 10.1111/jbi.14250 .
hal-03449129

HAL Id: hal-03449129

https://hal.science/hal-03449129

Submitted on 14 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The role of Sahara highlands in the diversification and desert colonization of the Bosc's fringe-toed lizard

André Vicente Liz<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID

Dennis Rödder<fnref n="3" />

Duarte Vasconcelos Gonçalves<fnref n="1" />,<fnref n="4" /> ID

Guillermo Velo-Antón<fnref n="1" />,<fnref n="5" />

Miguel M. Fonseca<fnref n="4" />

Philippe Geniez<fnref n="6" />

Pierre-André Crochet<fnref n="7" />

José Carlos Brito<fnref n="1" />,<fnref n="2" />

### Correspondence

André Vicente Liz, CIBIO/InBIO, Centro de Investigação em Biodiversidade e Recursos Genéticos da Universidade do Porto, Instituto de Ciências Agrárias de Vairão, Rua Padre Armando Quintas 7, 4485-661 Vairão, Portugal.

Email: andre.vicente.liz@cibio.up.pt

### Funding information

National Geographic Society; Mohammed bin Zayed Species Conservation Fund; Fundação para a Ciência e a Tecnologia; European Regional Development Fund

Handling Editor: Katharine Marske

### Abstract

**Aim:** 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.

**Location:** North Africa and Arabia.

**Taxon:** *Acanthodactylus boskianus* (Squamata: Lacertidae).

**Methods:** We integrated multi-locus mitochondrial and nuclear phylogenies with species' climate-niche modelling, including palaeo-projections. Genetic analyses aimed to assess the species' genetic structure, identify its main mitochondrial lineages and nuclear diversity, and reconstruct its ancestral biogeography. Species' climate-niche stability was modelled independently for the Late Pleistocene-Holocene and the Plio-Pleistocene, to infer historical climatic refugia and dispersal corridors.

**Results:** Four spatially structured mitochondrial lineages, integrating several parapatric sub-lineages, originated during the Plio-Pleistocene. Nuclear data revealed nine potential candidate species. Climatic refugia were located in mountains and desert fringes, remaining consistent for the Late Pleistocene-Holocene and the Plio-Pleistocene. Recurrent North-South climatic corridors were located along the desert periphery, while others less frequent were found across Central Sahara. Ancestral biogeography analyses recovered a recent Pleistocene colonization of the Sahel throughout eastern Sahara and either Sahara or Sahel origin for Central Sahara populations.

**Main conclusions:** Species' diversification was triggered by a combination of Plio-Pleistocene climatic cycles across a complex topographic region, where mountains acted as the main diversification hotspots. The historical role of Central Sahara highlands as main non-peripheral mesic refugia was corroborated. In addition, intermittent climatic connections linked Mediterranean and Sahel ecoregions with Central Sahara.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2021 The Authors. Journal of Biogeography published by John Wiley & Sons Ltd.

<footnote n="1">¹CIBIO/InBIO, Centro de Investigação em Biodiversidade e Recursos Genéticos da Universidade do Porto, Instituto de Ciências Agrárias de Vairão, Vairão, Portugal</footnote>

<footnote n="2">²Departamento de Biologia, Facultade de Ciências, Universidade do Porto, Porto, Portugal</footnote>

<footnote n="3">³ZFMK, Zoologisches Forschungsmuseum Alexander Koenig, Bonn, Germany</footnote>

<footnote n="4">⁴CIIMAR, Centro Interdisciplinar de Investigação Marinha e Ambiental, Terminal de Cruzeiros de Leixões, Matosinhos, Portugal</footnote>

<footnote n="5">⁵Departamento de Ecoloxía e Bioloxía Animal, Universidade de Vigo, Grupo GEA, Vigo, Spain</footnote>

<footnote n="6">⁶CEFE, Univ Montpellier, CNRS, EPHE-PSL University, IRD, Biogéographie et Ecologie des Vertébrés, Montpellier, France</footnote>

<footnote n="7">⁷CEFE, CNRS, Univ Montpellier, EPHE, IRD, Montpellier, France</footnote>
<figure><img src="image_2.png" /></figure>

refugia, suggesting the existence of alternative trans-Sahara dispersal routes to the putative coastal and Nile corridors.

KEYWORDS

arid mountains, biodiversity refugium, climatic cycles, cryptic diversity, dispersal corridor, historical biogeography, North Africa, phylogeography, Pliocene-Pleistocene

# 1 | INTRODUCTION

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 hyper-arid 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 (Brito et al., 2014). 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 (Velot-Antón et al., 2018). 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

<figure><img src="image_3.png" /><figcaption>FIGURE 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 dots 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</figcaption></figure>
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 (Brito et al., 2018). The scarce molecular studies available have uncovered lineages endemic to southern Algerian [e.g., *Uromastyx alfredschmidtii* (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., Šmid 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 (Velo-Antón et al., 2018) 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 (Brito et al., 2014). 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* (Gonçalves, Martínez-Freiría, et al., 2018)] 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 (Gonçalves, Martínez-Freiría, et al., 2018)]. 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.

The Bosc's fringe-toed lizard (*Acanthodactylus boskianus* Daudin, 1802) is a mesic ectotherm exhibiting a wide trans-Sahara range (Sindaco & Jeremčenko, 2008; Figure 1), thus being a promising model to assess biogeographic questions at the North African scale. 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 non-xeric species to overcome the desert arid barriers.

## 2 | MATERIALS AND METHODS

### 2.1 | Study area, organism and data collection

Our study area encompasses the Sahara, Sahel and Arabia arid belts (Figure 1). It covers the ecoregions of Sahara Desert, Sahel acacia savanna, and parts of the Sudanian savanna and Mediterranean steppes in North Africa, and the Arabian Desert and part of the adjacent steppes and xeric shrublands in the Middle East (Dinerstein et al., 2017).

*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; Tamar, Carranza, et al., 2016). 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
<figure><img src="image_4.png" /></figure>

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).

## 2.2 | 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 se-
quences (Table S2.1). These 44 samples were amplified and sequenced
for three additional mitochondrial DNA (mtDNA) fragments: dehy-
drogenase subunit 4 and partial tRNAs (ND4; 803 bp), cytochrome-b
(Cytb; 387 bp), and cytochrome oxidase I (COI; 612 bp); and four nu-
clear DNA (nuDNA) genes: melanocortin-1 receptor (MC1R; 589 bp),
beta-fibrinogen intron 7 (FGB; 1054 bp), ornithine decarboxylase
introns 1-2 and partial coding region (ODC; 467 bp), and phospho-
gluconate dehydrogenase intron 7 (PGD; 526 bp). Polymerase chain
reactions (PCRs) were performed in 10 µl reaction volumes comprising
5 µl of MyTaq (MyTaqTM Mix; Bioline), 3 µl of ultrapure water, 0.5 µl of
forward and reverse primers and 1 µl of template DNA (see 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.genewiz.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).

## 2.3 | 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 v1.1.1 (Lanfear et al., 2012) with the following settings: linked branch lengths, mrbayes models, BIC model selection criterion and greedy search. Partitions schemes by gene and codon were applied. Bayesian Inference (BI) analyses were performed using BEAST v1.10.4 (Suchard et al., 2018) on each dataset. BEAST was run on the CIPRES gateway (Miller et al., 2011) in two independent runs of 10⁸ generations, sampling at every 10⁴, 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 (Rambaut et al., 2018). 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⁴ iterations and 10³ of burn-in.

Phylogenetic trees were time-calibrated using gene substitution
rates. Mean substitution rates of 12S (0.00553 ± 0.00128), ND4
(0.0174 ± 0.004), and Cytb (0.0164 ± 0.00317) were retrieved from
calibrated phylogenies of closely related lacertids (Barata et al., 2012;
Carranza & Arnold, 2012). These rates have been used to date splits
within other *Acanthodactylus* species (Beddek et al., 2018; Velo-
Antón et al., 2018). The substitution rates were incorporated to the
mt-nuDNA and mtDNA trees. Both trees revealed equivalent tree
<figure><img src="image_5.png" /><figcaption>FIGURE 2 (a) Distribution map of the mtDNA lineages and sub-lineages of *Acanthodactylus boskianus* (n = 205; projection WGS84). (b) Left: 12S mtDNA haplotype network (n = 204). White dots represent mutation steps. Right: time-calibrated Bayesian mtDNA phylogeny (n = 40), obtained from the concatenation of four mtDNA fragments (12S, ND4, Cytb and COI; 2169 bp). Black and grey dots represent, respectively, node posterior probabilities (pp) >0.99 and >0.95; lower values are indicated right to the nodes if pp >0.5. Split times between lineages/sub-lineages are given in bold for supported nodes (95% HPD in brackets). Coloured polygons and symbols represent mtDNA lineages/sub-lineages, with name assignation based on geography. Samples of *A. schreiberi syriacus* are indicated with black dots. Haplotype numbers correspond to samples indicated in Table S2.1</figcaption></figure>

topologies, node divergence times and number of resolved nodes (see
Section 3.1; Figure S3.1 in Appendix S3).

Main mtDNA lineages were defined as reciprocally monophy-
letic clades in the mtDNA tree that split before the Pliocene. The 12S
haplogroups represented in the mtDNA tree formed distinct, well-
supported phylogenetic groups, which were defined as mtDNA sub-
lineages. The single exception is the haplogroup formed by all A. s.
syriacus samples, which was integrated within the sub-lineage *Sinai*

due to the lack of distinctiveness in the tree. Moreover, the hap-
logroup formed by the Syrian samples was not represented in the tree
due to failed amplification of additional markers beyond Cytb. The
nuDNA lineages were defined as distinct, well-supported groups in
the nuDNA tree. Mean genetic distances (uncorrected p-distances)
among mtDNA lineages and sub-lineages were calculated in MEGA
v.6 for the marker with the highest sampling coverage (12S; n = 204;
Tamura et al., 2013).
<figure><img src="image_6.png" /></figure>

## 2.4 | Ancestral area reconstruction

To infer the biogeographic history of *A. boskianus*, we performed a discrete ancestral area reconstruction with the Bayesian Stochastic Search Variable Selection (BSSVS; Lemey et al., 2009), implemented in BEAST v.1.10.4 (Suchard et al., 2018). 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 (Tamar, Carranza, et al., 2016); (2) the Lower Nile region, as the most probable gateway into Africa (Tejero-Cicuéndez et al., 2021); (3) the Mediterranean Basin, extending across North Sahara eastwards from the Moulouya biogeographic break (Beddek et al., 2018); (4) West Sahara, integrating the Atlantic coast biogeographic corridor (Velo-Antón et al., 2018); (5) the Central Sahara highlands; and (6) the Sahel ecoregion.

## 2.5 | 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 package (Aiello-Lammens et al., 2015). The final set was composed of 157 uncorrelated observations with minimum distance of ~80 km (Figure S3.2; Table S2.1). For modelling purposes, nineteen climatic variables for current conditions were downloaded from WorldClim at 2.5 arc-minute 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 middle Holocene (midHol; 6 kya) and the Last Glacial Maximum (LGM; 21 kya). Ensembles of projections are widely used to account for the variability between Global Circulation Model scenarios (Figure S3.3; Araujo & New, 2007; Thuiller et al., 2019). Ensembles were produced by averaging the downscaled outputs of all available scenarios (delta method; CCSM4, CNRM-CM5, FGOALS-g2, IPSL-CM5A-LR, MIROC-ESM, MPI-ESM-P and MRI-CGCM3; Fick & Hijmans, 2017) prior to the palaeoprojections (e.g. Eberle et al., 2017). Palaeoclimatic variables were downloaded at the same 2.5 arc-minute resolution and cropped to the study area to decrease computation time. In the second time-scale, the same set of variables was downloaded from Oscillayers dataset for 539 Plio-Pleistocene periods (10 kya time-steps; Gamisch, 2019). Given recent criticism regarding the robustness of Oscillayers

(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 background area was limited by a radial buffer of 500 km from presence records to avoid excessive biases towards regions with no occurrence (Figure S3.2; Anderson & Raza, 2010). For model training/testing, an 80/20% partition of presence records was used. The area-under-the-curve (AUC) of the receiver-operating characteristics (ROC) plots was assumed as measure of individual model fit (Fielding & Bell, 1997). Response curves of climatic variables were produced and verified.

A consensus forecast of the probability of species' occurrence under current climatic conditions was generated by averaging individual model replicates (Marmion et al., 2009). Standard deviation among replicates was verified to account for prediction uncertainty (e.g., Velo-Antón et al., 2018). Mean and standard deviation of the model forecast were projected to past climatic conditions to predict the species' spatial probability of occurrence for each time period. This approach, based on the formula estimated for each model replicate, provides identical results to those obtained by averaging the projections of each model replicate (e.g. Gonçalves, Martínez-Freiría, et al., 2018), but implies lower computation time. Multivariate Environmental Similarity Surfaces (MESS) were inferred for each projection to identify and minimize errors in spatial extrapolations beyond the training range of the models (Elith et al., 2010). The consensus model and projections to past climatic conditions were converted to binary (i.e., suitable/unsuitable) using a 20% training cut-off threshold. Different thresholds were tested in the model projection for the Present, selecting the more concordant one with the known species' range (Figure S3.4). Binary layers were averaged across time periods to obtain a rate map of historical climate-niche stability and identify areas of constant (i.e., refugia; Carnaval et al., 2009) and temporary suitable climate (i.e. corridors). In order to verify the coherence of climate-niche stability between Late Pleistocene-Holocene and Plio-Pleistocene time-scales, projections were averaged separately for WorldClim and Oscillayers.

## 2.6 | 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. In total, 14 units were defined, which hosted ~90% of the genotyped samples (Figure S3.5). For each unit, we calculated: averages of absolute altitude, index of terrain roughness (2.5 arc-minute resolution; Title & Bemmels, 2018) and rate of Plio-Pleistocene climate-niche stability, using the Zonal Statistics tool of ArcGIS v.10.5 (ESRI, 2006); and species' 12S nucleotide diversity (Pi; Nei, 1987), using DNASP 5.10.01 (Librado & Rozas, 2009). A Generalized Linear Model (GLM) with normal errors was fitted in R, using nucleotide diversity as response variable and altitude, climatic stability, terrain roughness, and sample count as low-correlated (Pearson's r < 0.6) predictors (Matthiopoulos, 2011).

## 3 | RESULTS

### 3.1 | Phylogenetic analyses and cyto-nuclear discordances

The BI mtDNA tree revealed four main well-supported, spatially-structured lineages within *A. boskianus*, integrating 18 sub-lineages (Figure 2): (i) an early splitting lineage restricted to inland areas of Syria-Jordan (*Levant*); (ii) a highly divergent lineage extending alongside the Euphrates Basin (*Persian Gulf*); iii) a clade incorporating the remaining ~80% samples from the Middle East, including *A. s. syriacus*, plus two African samples from the adjacent Nile Delta area, which is structured in five sub-lineages (*Arabia*); and iv) a highly diverse clade made exclusively of the African samples and structured in eleven sub-lineages of very different range sizes, with only one distributed across the eastern side of the Sahara-Sahel, from Egypt to Niger, and all the remaining being small-ranged across the western side (*North Africa*). The Afro-Arabian suture zone hosts the largest lineage diversity (Figure 2a). For example, the three Middle East lineages occur in Jordan, and the Nile Delta constitutes a contact zone for Arabian and African lineages. The mtDNA haplotype network (Figure 2b) allocated the four unrepresented haplogroups in the tree within the lineages *Levant* (*Syria*), *Persian Gulf* (*Al Hajar*), and *Arabia* (*Dhofar* and *Northwest Egypt*). Mean 12S genetic distances between lineages ranged from 0.030 to 0.044, while some sub-lineages were highly divergent (e.g., *Euphrates Basin*, *Hoggar*; see Table S2.4). The TMRCA for the first branching within *A. boskianus* (*Levant* vs remaining groups) dates back to the Late Miocene, ~6.91 Mya (95% HPD: 5.0–8.9; Figure 2b). The second split (*Persian Gulf* vs. remaining groups) occurred ~5.78 Mya (4.3–7.4). The split between the two other lineages (*Arabia* vs. *North Africa*), although poorly supported, dates back to the Miocene-Pliocene transition (~5.42 Mya; 4.0–6.8). All sub-lineages of *A. boskianus* diversified throughout the Plio-Pleistocene. Yet, the weak support of various tree nodes implies unclear relationships among some African sub-lineages.

The BI nuDNA tree revealed three main supported branches which largely correspond to the *Persian Gulf*, *Arabia* and *North Africa* mtDNA lineages (Figure 3b). It lacks representatives for the *Levant* mtDNA lineage and two of the *Arabia* sub-lineages, due to failed DNA amplification. These main branches integrate nine supported nuDNA lineages. The position of 13 of the samples and several relationships among the African lineages remained uncertain. The highest diversity in nuDNA lineages is found in Egypt, where three of them occur (Figure 3a). The nuDNA lineages are congruent among the haplotype networks (Figure 3b), with clear differentiation of the three main groups and only occasional allele sharing among the African lineages. Although the assignment of samples to mtDNA and nuDNA lineages/sub-lineages was generally congruent, two major cyto-nuclear discordances are revealed by the nuDNA paraphyly of the mtDNA lineage *Arabia* and the lack of support for the sub-lineage *East Sahara-Sahel*, whose samples are distributed among three nuDNA lineages: *Egypt 1*, *Egypt 2* and *Niger*. Nevertheless, it is worth noticing that *Egypt 1* and *Egypt 2* share alleles in nuclear MC1R and FGB.

### 3.2 | Ancestral area reconstruction

The ancestral biogeography of *A. boskianus* recovers a most-likely Middle East origin (77% probability; Figure 4). The species successfully colonized North Africa twice: the first time during the Pliocene, with a diversification crown age of ~4.7 Mya (HPD 95%: 3.4–6.0) for the last likely Middle East ancestor (50%) and which led to the species' extensive African radiation; and the later taking place within the last ~3.2 Mya (2.2–4.4), with expansion restricted to the Lower Nile region. Sahel populations probably originated at the end of the Pleistocene from two biogeographic regions: the Lower Nile (92%), which led to local diversification in Niger within the last ~1.3 Mya (0.9–1.8); and the West Sahara (85%), with the arrival to southern Mauritania occurring not before ~0.2 Mya (0.1–0.3). Fluxes towards the Central Sahara mountains were at least coming from two different areas: the Nigerien Sahel (94%), dating back to the last ~0.2 Mya (0.1–0.3); and, most likely, areas south of the Atlas Mountains, in West Sahara, in more than one occasion during the Pliocene (60% probability at the crown age of the Sahara groups). Different dispersal events from West Sahara to the Mediterranean may have also occurred since the Late Pliocene (65% probability at the crown age of the Mediterranean groups).

### 3.3 | 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
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

FIGURE 3 (a) Distribution map of *Acanthodactylus boskianus* nuDNA lineages (n = 34; projection WGS84). (b) Left: nuDNA haplotype networks for MC1R (n = 34), ODC (25), FGB (31) and PGD (33). White dots represent mutation steps. Right: Bayesian nuDNA phylogeny (n = 34), obtained from the concatenation of four nuDNA genes (MC1R, ODC, FGB and PGD; 2169 bp). Black and white dots represent, respectively, node posterior probabilities (pp) >0.99 and >0.90; lower values are indicated right to the nodes if pp >0.5. Colours represent nuDNA lineages, with name assignation matching with the corresponding mtDNA sub-lineage if congruent (otherwise, indicated with an asterisk). Main mtDNA lineages are marked using the colour polygons of Figure 2. Haplotype numbers (sorted as MC1R/ODC/FGB/PGD) correspond to samples indicated in Table S2.1

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.
<figure><img src="image_9.png" /></figure>

**FIGURE 4** Time-calibrated Bayesian mt-nuDNA phylogeny (n = 34) implementing the BSSVS method for the ancestral area reconstruction of *Acanthodactylus boskianus*. Black dots represent supported nodes with posterior probabilities (pp) >0.99; lower values are indicated right to the nodes if pp > 0.5. Branch colours represent ancestral range, based on the biogeographic regions indicated in the bottom-left map (regions intersected with the species' range polygon; samples depicted with white circles). Pie charts indicate the probability of each biogeographic region for the main supported nodes. Coloured arrows indicate dispersal events between biogeographic regions, considering supported nodes

<figure><img src="image_10.png" /></figure>

Projection comparisons, only possible for LGM, and historical climate-niche stability maps showed consistent results between Worldclim and Oscillayers datasets. For LGM, the main disparities were restricted to a southern portion of the Arabian Peninsula and the western Sahel, which were inferred as areas of suitable habitat only in Oscillayers (Figure S3.4). For historical stability, the occasional connectivity between Central Sahara and the Sahel was suggested only by Worldclim (Figure 5).

## 3.4 | 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). Climatic stability and altitude positively correlated with nucleotide diversity (Figure S3.8; Table S2.6). Conversely, no significant effect of topographic heterogeneity on nucleotide diversity was found.

# 4 | 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 (Gonçalves, Martínez-Freiría, et al., 2018; 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; Tamar, Carranza, et al., 2016). 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 (Vale et al., 2014), but these caveats are hardly avoidable considering the spatial scale of the analysis.
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 5 Stable climatic areas for *Acanthodactylus boskianus* over time, derived from climate-niche models and respective projections to: (a) two midHol and LGM; and (b) 539 Plio-Pleistocene time steps. Warm colours represent higher climatic stability. Projection: WGS84</figcaption></figure>

<table><thead><tr><th></th><th>Coefficient</th><th>SE</th><th>T</th><th>p</th></tr></thead><tbody><tr><td>Intercept</td><td>-5.536 × 10<sup>-3</sup></td><td>2.89 × 10<sup>-3</sup></td><td>-1.919</td><td>0.087</td></tr><tr><td>Climate-niche stability</td><td>1.357 × 10<sup>-2</sup></td><td>4.47 × 10<sup>-3</sup></td><td>3.037</td><td><strong>0.014</strong></td></tr><tr><td>Absolute altitude</td><td>1.991 × 10<sup>-5</sup></td><td>3.90 × 10<sup>-6</sup></td><td>5.112</td><td><strong>&lt;0.001</strong></td></tr><tr><td>Terrain roughness index</td><td>-2.578 × 10<sup>-4</sup></td><td>1.17 × 10<sup>-4</sup></td><td>-2.205</td><td>0.055</td></tr><tr><td>Sampling size</td><td>7.032 × 10<sup>-5</sup></td><td>7.21 × 10<sup>-5</sup></td><td>0.976</td><td>0.355</td></tr></tbody></table>

Significant results ($p$ < 0.05) are highlighted in bold.

<caption>TABLE 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)</caption>

## 4.1 | 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 (Tamar, Carranza, et al., 2016). These dispersals probably constitute back-to-Africa events (Tejero-Cicuéndez et al., 2021), since the genus is most likely of African origin (Tamar, Carranza, et al., 2016). 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
<figure><img src="image_13.png" /></figure>

contractions/expansions and subsequent episodes of population isolation/re-connectivity (Figures S3.7). Although aridity presumably induces diversification in desert mesic species (Brito et al., 2014), 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.

The positive correlation of genetic diversity with climatic stability and altitude highlights the expected role of mountains as both diversification hotspots and climatic refugia (Table 1; Figure S3.8; Rahbek et al., 2019). Our results revealed absolute altitude as the main topographic factor related to the diversification of *A. boskianus*, contrary with other squamates where genetic diversity is driven by topographic heterogeneity (Šmid et al., 2021). Historical range fragmentation (Figures S3.7) and the presence of endemic lineages/sub-lineages (Figures 2a and 3a) to Hoggar, Tassili n'Ajjer and Fezzan suggest a "sky-island" pattern around Central Sahara highlands (e.g., Gonçalves, Pereira, et al., 2018; Metallinou et al., 2015). These dynamics are also observed south of the Atlas, where the presence of various endemic lineages provide further support for the area constituting a diversification cradle (e.g., Martínez-Freiría et al., 2017). A similar phylogeographic scenario was found around the Arabian mountains, where intense Miocene tectonic activity may have also played an important role in generating diversity (Machado et al., 2020). Our genetic output provides new evidence for the Yemen, Dhofar and Al Hajar mountain diversification centres (Šmid et al., 2021; Figures 2 and 3). The genetic distinctiveness of Levant populations (Table S2.4) also matches with previous studies (e.g., Metallinou et al., 2015).

## 4.2 | Dispersal routes

Genetic patterns (Figure 2a) and climate-niche models (Figure 5) support several historical connections in *A. boskianus* between the Central Sahara highlands and the Mediterranean biome in the western side of the desert. Firstly, the distribution of a mtDNA lineage and the temporary presence of areas with suitable climate from central Tunisia to Fezzan and the nearby Algerian massifs mirrors phylogeographic patterns described in other terrestrial, highly mobile reptile [e.g., *Psammophis schokari* (Gonçalves, Martínez-Freiría, et al., 2018)] and amphibian taxa [e.g., *Bufotes boulengeri* (Nicolas et al., 2017), *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 palaeo-models 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 Air 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* (Gonçalves, Martínez-Freiría, et al., 2018), 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)] could have worked as minor refugia and stepping-stones during range expansion periods (Garcia-Porta et al., 2017). However, the lack of genetic samples across these areas prevents verifying this hypothesis, while the underrepresentation of these mountains in the models might be due to the low number of local observations. An alternative explanation to the recovered genetic pattern is dispersal through the Nile Valley, then along the Sahel, similar to *Varanus niloticus* (Dowell et al., 2016). Our models do not predict the suitability of the Nile Valley for *A. boskianus*, yet the permanent availability of moisturized habitats suggests it as a reasonable scenario. In addition, although a putative Sahel longitudinal connection for *A. boskianus* lacks palaeomodel support, it has been suggested from genetic data of other mesic groups [e.g., *Agama* sp. (Kissling et al., 2016), *Mastomys huberti* (Mouline et al., 2008), *V. niloticus* (Dowell et al., 2016)]. Further sampling in Libya, Sudan and Chad, combined with the use of fast-evolving genetic markers, is necessary to verify potential gene flow and dispersal routes across eastern Sahara.

The hypothesized North-South coastal corridors (Brito et al., 2014) 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* (Gonçalves, Martínez-Freiría, et al., 2018) and *Acanthodactylus aureus* (Velo-Antón et al.,
<figure><img src="image_14.png" /></figure>

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 ecototherms 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 climate-induced 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.

## 4.3 | 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.

# 5 | 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 sky-island 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.

## ACKNOWLEDGEMENTS

We thank BIODESERTS group members for their assistance during fieldwork. We thank H. in den Bosch, E. Didner, M. Geniez, J. Viglione, A. Cluchier, H. Loumassine, M. Beddek, O. Peyre, A. Benlahrech, M. Siol, B. Shacham, Y. Werner and S. Baha el Din for help with the sample collection. Acknowledgements extended to Petros Lymberakis from the Natural History Museum of Crete, Greece and Wolfgang Böhme from the Zoological Research Museum Alexander Koenig, Bonn. This work was funded by Fundação para a Ciência e a Tecnologia – FCT (PTDC/BIA-BEC/099934/2008 and PTDC/BIA-BIC/2903/2012), by FEDER funds through the Operational Programme for Competitiveness Factors - COMPETE (FCOMP-01-0124-FEDER-008917/028276), by AGRIGEN-NORTE-01-0145-FEDER-000007, supported by Norte Portugal Regional Operational Programme (NORTE2020), under the PORTUGAL 2020 Partnership Agreement, through the European Regional Development Fund (ERDF), by National Geographic Society (CRE 7629-04, CRE 8412-08), and by Mohammed bin Zayed Species Conservation Fund (11052709, 11052707, 11052499). AVL, DVG, GVA and JCB were supported by FCT (SFRH/BD/140348/2018, CEECIND/03848/2020, CEECIND/00937/2018, and CEECINST/00014/2018/CP1512/CT0001, respectively). Capture permits were issued by the Haut Commissariat aux Eaux et Forêts (278/2012 and 20/2013/HCEFLCD/DLCDPN/DPRN/CFF) and Ministère de l'Environnement et du Développement Durable de Mauritania (460/MDE/PNBA). Logistic support for fieldwork was given by Pedro Santos Lda (Trimble GPS), Off Road Power Shop, P.N. Banc d'Arguin and P.N. Diawling (Mauritania), Association Nature Initiative (Morocco), and Université des Sciences, de Technologie et de Médecine de Nouakchott.
## CONFLICT OF INTEREST

The authors declare that they have no conflict of interest.

## DATA AVAILABILITY STATEMENT

GenBank accession numbers and sampling localities are indicated in Table S2.1.

## ORCID

André Vicente Liz https://orcid.org/0000-0001-6131-5194

Duarte Vasconcelos Gonçalves https://orcid.org/0000-0003-4299-0375

## REFERENCES

Aiello-Lammens, M. E., Boria, R. A., Radosavljevic, A., Vilela, B., & Anderson, R. P. (2015). spThin: an R package for spatial thinning of species occurrence records for use in ecological niche models. *Ecography*, 38, 541–545. https://doi.org/10.1111/ecog.01132

Anderson, R. P., & Raza, A. (2010). The effect of the extent of the study region on GIS models of species geographic distributions and estimates of niche evolution: preliminary tests with montane rodents (genus Nephelomys) in Venezuela. *Journal of Biogeography*, 37, 1378–1393. https://doi.org/10.1111/j.1365-2699.2010.02290.x

Araujo, M. B., & New, M. (2007). Ensemble forecasting of species distributions. *Trends in Ecology & Evolution*, 22, 42–47. https://doi.org/10.1016/j.tree.2006.09.010

Arnold, E. N. (1983). Osteology, genitalia and the relationships of Acanthodactylus (Reptilia: Lacertidae). *Bulletin of the British Museum (Natural History). Zoology*, 44, 291–339.

Barata, M., Carranza, S., & Harris, D. J. (2012). Extreme genetic diversity in the lizard *Atlantolacerta andreanskyi* (Werner, 1929): A montane cryptic species complex. *BMC Evolutionary Biology*, 12, 1–19. https://doi.org/10.1186/1471-2148-12-167

Beddek, M., Zenboudji-Beddek, S., Geniez, P., Fathalla, R., Sourouille, P., Arnal, V., Dellaoui, B., Koudache, F., Telailia, S., Peyre, O., & Crochet, P.-A. (2018). Comparative phylogeography of amphibians and reptiles in Algeria suggests common causes for the east-west phylogeographic breaks in the Maghreb. *PLoS One*, 13, e0201218. https://doi.org/10.1371/journal.pone.0201218

Boria, R. A., Olson, L. E., Goodman, S. M., & Anderson, R. P. (2014). Spatial filtering to reduce sampling bias can improve the performance of ecological niche models. *Ecological Modelling*, 275, 73–77. https://doi.org/10.1016/j.ecolmodel.2013.12.012

Brito, J. C., Durant, S. M., Pettorelli, N., Newby, J., Canney, S., Algadafi, W., Rabeil, T., Crochet, P.-A., Pleguezuelos, J. M., Wacher, T., de Smet, K., Gonçalves, D. V., da Silva, M. J. F., Martínez-Freiría, F., Abáigar, T., Campos, J. C., Comizzoli, P., Fahd, S., Fellous, A., ... Carvalho, S. B. (2018). Armed conflicts and wildlife decline: Challenges and recommendations for effective conservation policy in the Sahara-Sahel. *Conservation Letters*, 11, e12446. https://doi.org/10.1111/conl.12446

Brito, J. C., Godinho, R., Martínez-Freiría, F., Pleguezuelos, J. M., Rebelo, H., Santos, X., Vale, C. G., Velo-Antón, G., Boratyński, Z., Carvalho, S. B., Ferreira, S., Gonçalves, D. V., Silva, T. L., Tarroso, P., Campos, J. C., Leite, J. V., Nogueira, J., Álvares, F., Sillero, N., ... Carranza, S. (2014). Unravelling biodiversity, evolution and threats to conservation in the Sahara-Sahel. *Biological Reviews*, 89, 215–231. https://doi.org/10.1111/brv.12049

Brito, J. C., Tarroso, P., Vale, C. G., Martínez-Freiría, F., Boratyński, Z., Campos, J. C., Ferreira, S., Godinho, R., Gonçalves, D. V., Leite, J. V., Lima, V. O., Pereira, P., Santos, X., da Silva, M. J. F., Silva, T. L., Velo-Antón, G., Veríssimo, J., Crochet, P.-A., Pleguezuelos, J. M., &

Carvalho, S. B. (2016). Conservation biogeography of the Sahara-Sahel: Additional protected areas are needed to secure unique biodiversity. *Diversity and Distributions*, 22, 371–384. https://doi.org/10.1111/ddi.12416

Brown, J. L., Hill, D. J., & Haywood, A. M. (2020). A critical evaluation of the Oscillayers methods and datasets. *Global Ecology and Biogeography*, 29, 1435–1442. https://doi.org/10.1111/geb.13103

Carnaval, A. C., Hickerson, M. J., Haddad, C. F., Rodrigues, M. T., & Moritz, C. (2009). Stability predicts genetic diversity in the Brazilian Atlantic forest hotspot. *Science*, 323, 785–789. https://doi.org/10.1126/science.1166955

Carranza, S., & Arnold, E. N. (2012). A review of the geckos of the genus Hemidactylus (Squamata: Gekkonidae) from Oman based on morphology, mitochondrial and nuclear data, with descriptions of eight new species. *Zootaxa*, 3378, 1–95. https://doi.org/10.11646/zootaxa.3378.1.1

Clement, M., Posada, D. C. K. A., & Crandall, K. A. (2000). TCS: A computer program to estimate gene genealogies. *Molecular Ecology*, 9, 1657–1659. https://doi.org/10.1046/j.1365-294x.2000.01020.x

Di Cola, V., Broennimann, O., Petitpierre, B., Breiner, F. T., D'Amen, M., Randin, C., Engler, R., Pottier, J., Pio, D., Dubuis, A., Pellissier, L., Mateo, R. G., Hordijk, W., Salamin, N., & Guisan, A. (2017). ecospat: An R package to support spatial analyses and modeling of species niches and distributions. *Ecography*, 40, 774–787. https://doi.org/10.1111/ecog.02671

Dinerstein, E., Olson, D., Joshi, A., Vynne, C., Burgess, N. D., Wikramanayake, E., Hahn, N., Palminteri, S., Hedao, P., Noss, R., Hansen, M., Locke, H., Ellis, E. C., Jones, B., Barber, C. V., Hayes, R., Kormos, C., Martin, V., Crist, E., ... Saleem, M. (2017). An ecoregion-based approach to protecting half the terrestrial realm. *BioScience*, 67, 534–545. https://doi.org/10.1093/biosci/bix014

Dowell, S. A., Portik, D. M., de Buffrénil, V., Ineich, I., Greenbaum, E., Kolokotronis, S. O., & Hekkala, E. R. (2016). Molecular data from contemporary and historical collections reveal a complex story of cryptic diversification in the *Varanus (Polydaedalus) niloticus* Species Group. *Molecular Phylogenetics and Evolution*, 94, 591–604. https://doi.org/10.1016/j.ympev.2015.10.004

Drake, N. A., Blench, R. M., Armitage, S. J., Bristow, C. S., & White, K. H. (2011). Ancient watercourses and biogeography of the Sahara explain the peopling of the desert. *Proceedings of the National Academy of Sciences of the United States of America*, 108, 458–462. https://doi.org/10.1073/pnas.1012231108

Drake, N. A., El-Hawat, A. S., Turner, P., Armitage, S. J., Salem, M. J., White, K. H., & McLaren, S. (2008). Paleohydrology of the Fazzan Basin and surrounding regions: the last 7 million years. *Paleogeography, Paleoclimatology, Paleocology*, 263, 131–145.

Drummond, A. J., Ashton, B., Buxton, S., Cheung, M., Heled, J., & Kearse, M. (2010). Geneious 4.8. 2010. Retrieved from: http://www.geneious.com

Drummond, A. J., Ho, S. Y. W., Phillips, M. J., & Rambaut, A. (2006). Relaxed phylogenetics and dating with confidence. *PLoS Biology*, 4, e88. https://doi.org/10.1371/journal.pbio.0040088

Durant, S. M., Pettorelli, N., Bashir, S., Woodroffe, R., Wacher, T., De Ornellas, P., Ransom, C., Abaigar, T., Abdelgadir, M., El Alqamy, H., Beddiaf, M., Belbachir, F., Belbachir-Bazi, A., Berbash, A. A., Beudels-Jamar, R., Boitani, L., Breitenmoser, C., Cano, M., Chardonnet, P., ... Baillie, J. E. M. (2012). Forgotten biodiversity in desert ecosystems. *Science*, 336, 1379–1380. https://doi.org/10.1126/science.336.6087.1379

Durant, S. M., Wacher, T., Bashir, S., Woodroffe, R., De Ornellas, P., Ransom, C., Newby, J., Abáigar, T., Abdelgadir, M., El Alqamy, H., Baillie, J., Beddiaf, M., Belbachir, F., Belbachir-Bazi, A., Berbash, A. A., Bemadjim, N. E., Beudels-Jamar, R., Boitani, L., Breitenmoser, C., ... Pettorelli, N. (2014). Fiddling in biodiversity hotspots while deserts burn? Collapse of the Sahara's megafauna. *Diversity and Distributions*, 20, 114–122. https://doi.org/10.1111/ddi.12157
<figure><img src="image_15.png" /></figure>

Eberle, J., Rödder, D., Beckett, M., & Ahrens, D. (2017). Landscape genetics indicate recently increased habitat fragmentation in African forest-associated chafers. *Global Change Biology*, 23, 1988–2004.

Elith, J., Graham, C. H., Anderson, R. P., Dudík, M., Ferrier, S., Guisan, A., Hijmans, R. J., Huettmann, F., Leathwick, J. R., Lehmann, A., Li, J., Lohmann, L. G., Loiselle, B. A., Manion, G., Moritz, C., Nakamura, M., Nakazawa, Y., McC. M. Overton, J., Townsend Peterson, A., ... Zimmermann, N. E. (2006). Novel methods improve prediction of species' distributions from occurrence data. *Ecography*, 29, 129–151.
https://doi.org/10.1111/j.2006.0906-7590.04596.x.

Elith, J., Kearney, M., & Phillips, S. (2010). The art of modelling range-shifting species. *Methods in Ecology and Evolution, 1*(4), 330–342.
https://doi.org/10.1111/j.2041-210X.2010.00036.x

Elsen, P. R., & Tingley, M. W. (2015). Global mountain topography and the fate of montane species under climate change. *Nature Climate Change*, 5, 772–776. https://doi.org/10.1038/nclimate2656

ESRI. (2006). ArcMap 10.5. Environmental Systems Research Institute Inc.

Ficetola, G. F., Bonardi, A., Sindaco, R., & Padoa-Schioppa, E. (2013). Estimating patterns of reptile biodiversity in remote regions. *Journal of Biogeography*, 40, 1202–1211. https://doi.org/10.1111/jbi.12060

Fick, S. E., & Hijmans, R. J. (2017). WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. *International Journal of Climatology*, 37, 4302–4315. https://doi.org/10.1002/joc.5086

Fielding, A. H., & Bell, J. F. (1997). A review of methods for the assessment of prediction errors in conservation presence/absence models. *Environmental Conservation*, 24, 38–49. https://doi.org/10.1017/S0376892997000088

Gamisch, A. (2019). Oscillayers: A dataset for the study of climatic oscillations over Plio-Pleistocene time-scales at high spatial-temporal resolution. *Global Ecology and Biogeography*, 28, 1552–1560. https://doi.org/10.1111/geb.12979

Gamisch, A. (2020). A reply to the 'critical evaluation of the Oscillayers methods and dataset'. *Global Ecology and Biogeography*, 29, 1443–1448. https://doi.org/10.1111/geb.13108

Garcia-Porta, J., Simó-Riudalbas, M., Robinson, M., & Carranza, S. (2017). Diversification in arid mountains: Biogeography and cryptic diversity of *Pristurus rupestris rupestris* in Arabia. *Journal of Biogeography*, 44, 1694–1704.

Geniez, P., & Arnold, E. N. (2006). A new species of semaphore gecko *Pristurus* (Squamata: Gekkonidae) from Mauretania, represents a 4700km range extension for genus. *Zootaxa*, 1317, 57–68. https://doi.org/10.11646/zootaxa.13171.1.6

Gonçalves, D. V., & Brito, J. C. (2019). Second Sahelian amphibian endemism suggested by phylogeography of Groove crowned Bullfrog (*Hoplobatrachus occipitalis*) in western Sahel and hints of polyploid species formation. *Journal of Zoological Systematics and Evolutionary Research*, 58, 262–274.

Gonçalves, D. V., Martínez-Freiría, F., Crochet, P. A., Geniez, P., Carranza, S., & Brito, J. C. (2018). The role of climatic cycles and trans-Saharan migration corridors in species diversification: biogeography of *Psammophis schokari* group in North Africa. *Molecular Phylogenetics and Evolution*, 118, 64–74. https://doi.org/10.1016/j.ympev.2017.09.009

Gonçalves, D. V., Pereira, P., Velo-Antón, G., Harris, D. J., Carranza, S., & Brito, J. C. (2018). Assessing the role of aridity-induced vicariance and ecological divergence in species diversification in North-West Africa using *Agama* lizards. *Biological Journal of the Linnean Society*, 124, 363–380. https://doi.org/10.1093/biolinnean/bly055

Harris, D. J., & Arnold, E. N. (2000). Elucidation of the relationships of spiny-footed lizards, *Acanthodactylus* spp. (Reptilia: Lacertidae) using mitochondrial DNA sequence, with comments on their biogeography and evolution. *Journal of Zoology*, 252, 351–362.

Heidari, N., Rastegar-Pouyani, E., Rastegar-Pouyani, N., & Faizi, H. (2014). Molecular phylogeny and biogeography of the genus *Acanthodactylus* Fitzinger, 1834 (Reptilia: Lacertidae) in Iran, inferred from mtDNA Sequences. *Zootaxa*, 3860, 379–395.

Katoh, K., & Standley, D. M. (2013). MAFFT multiple sequence alignment software version 7: Improvements in performance and usability. *Molecular Biology and Evolution*, 30, 772–780. https://doi.org/10.1093/molbev/mst010

Khannoon, E. R., Lunt, D. H., Schulz, S., & Hardege, J. D. (2013). Divergence of scent pheromones in allopatric populations of *Acanthodactylus boskianus* (Squamata: Lacertidae). *Zoological Science*, 30, 380–385.

Kingman, J. F. C. (1982). The coalescent. *Stochastic Processes and their Applications*, 13, 235–248. https://doi.org/10.1016/0304-4149(82)90011-4

Kissling, W. D., Blach-Overgaard, A., Zwaan, R. E., & Wagner, P. (2016). Historical colonization and dispersal limitation supplement climate and topography in shaping species richness of African lizards (Reptilia: Agaminae). *Scientific Reports*, 6, 1–14. https://doi.org/10.1038/srep34014

Lamb, H. F., Bates, C. R., Coombes, P. V., Marshall, M. H., Umer, M., Davies, S. J., & Dejen, E. (2007). Late Pleistocene desiccation of Lake Tana, source of the blue Nile. *Quaternary Science Reviews*, 26, 287–299. https://doi.org/10.1016/j.quascirev.2006.11.020

Lanfear, R., Calcott, B., Ho, S. Y. W., & Guindon, S. (2012). PartitionFinder: Combined selection of partitioning schemes and substitution models for phylogenetic analyses. *Molecular Biology and Evolution*, 29, 1695–1701. https://doi.org/10.1093/molbev/mss020

Leite, J. V., Álvares, F., Velo-Antón, G., Brito, J. C., & Godinho, R. (2015). Differentiation of North African foxes and population genetic dynamics in the desert—insights into the evolutionary history of two sister taxa, Vulpes rueppellii and Vulpes vulpes. *Organisms Diversity and Evolution*, 15, 731–745. https://doi.org/10.1007/s13127-015-0232-8

Lemey, P., Rambaut, A., Drummond, A. J., & Suchard, M. A. (2009). Bayesian phylogeography finds its roots. *PLoS Computational Biology*, 5, e1000520. https://doi.org/10.1371/journal.pcbi.1000520

Librado, P., & Rozas, J. (2009). DnaSP v5: A software for comprehensive analysis of DNA polymorphism data. *Bioinformatics*, 25, 1451–1452. https://doi.org/10.1093/bioinformatics/btp187

Lomolino, M. V., Riddle, B. R., Brown, J. H., & Brown, J. H. (2006). *Biogeography*. Sinauer Associates.

Machado, L., Salvi, D., Harris, D. J., Brito, J. C., Crochet, P. A., Geniez, P., & Carranza, S. (2020). Systematics, biogeography and evolution of the Saharo-Arabian naked-toed geckos genus *Tropiocolotes*. *Molecular Phylogenetics and Evolution*, 155, 106969.

Marmion, M., Parviainen, M., Luoto, M., Heikkinen, R. K., & Thuiller, W. (2009). Evaluation of consensus methods in predictive species distribution modelling. *Diversity and Distributions*, 15, 59–69. https://doi.org/10.1111/j.1472-4642.2008.00491.x

Martínez-Freiría, F., Crochet, P. A., Fahd, S., Geniez, P., Brito, J. C., & Velo-Antón, G. (2017). Integrative phylogeographical and ecological analysis reveals multiple Pleistocene refugia for Mediterranean *Daboia* vipers in north-west Africa. *Biological Journal of the Linnean Society*, 122, 366–384. https://doi.org/10.1093/biolinnean/blx038

Matthiopoulos, J. (2011). *How to be a quantitative ecologist: the* A to R*of green mathematics and statistics*. John Wiley & Sons.

Metallinou, M., Červenka, J., Crochet, P. A., Kratochvíl, L., Wilms, T., Geniez, P., & Carranza, S. (2015). Species on the rocks: Systematics and biogeography of the rock-dwelling *Ptyodactylus* geckos (Squamata: Phyllodactylidae) in North Africa and Arabia. *Molecular Phylogenetics and Evolution*, 85, 208–220.

Miller, M. A., Pfeiffer, W., & Schwartz, T. (2011). The CIPRES science gateway: a community resource for phylogenetic analyses. In: *Proceedings of the 2011 TeraGrid Conference: extreme digital discovery* (pp. 1–8).

Miralles, A., Geniez, P., Beddek, M., Aranda, D. M., Brito, J. C., Leblois, R., & Crochet, P. A. (2020). Morphology and multilocus phylogeny of the Spiny-footed Lizard (*Acanthodactylus erythrurus*) complex reveal two new mountain species from the Moroccan Atlas. *Zootaxa*, 4747, 302–326.

Mouline, K., Granjon, L., Galan, M., Tatard, C., Abdoullaye, D., Atteyine, S. A., Duplantier, J.-M., & Cosson, J.-F. (2008). Phylogeography of a
<figure><img src="image_16.png" /></figure>

Sahelian rodent species *Mastomys huberti*: A Plio-Pleistocene story of emergence and colonization of humid habitats. *Molecular Ecology*, 17, 1036-1053. https://doi.org/10.1111/j.1365-294X.2007.03610.x

Nei, M. (1987). *Molecular evolutionary genetics*. Columbia University Press.

Nicolas, V., Mataame, A., Crochet, P. A., Geniez, P., Fahd, S., & Ohler, A. (2017). Phylogeography and ecological niche modeling unravel the evolutionary history of the African green toad, *Bufotes boulengeri boulengeri* (Amphibia: Bufonidae), through the Quaternary. *Journal of Zoological Systematics and Evolutionary Research*, 56, 102–116.

Nicolas, V., Mataame, A., Crochet, P. A., Geniez, P., & Ohler, A. (2015). Phylogeographic patterns in north African water frog *Pelophylax saharicus* (Anura: Ranidae). *Journal of Zoological Systematics and Evolutionary Research*, 53, 239–248.

Palchan, D., & Torfstein, A. (2019). A drop in Sahara dust fluxes records the northern limits of the African Humid Period. *Nature Communications*, 10, 1–9. https://doi.org/10.1038/s41467-019-11701-z

Perrigo, A., Hoorn, C., & Antonelli, A. (2020). Why mountains matter for biodiversity. *Journal of Biogeography*, 47, 315–325. https://doi.org/10.1111/jbi.13731

Phillips, S. J., Anderson, R. P., & Schapire, R. E. (2006). Maximum entropy modeling of species geographic distributions. *Ecological Modelling*, 190, 231–259. https://doi.org/10.1016/j.ecolmodel.2005.03.026

Poulakakis, N., Kapli, P., Kardamaki, A., Skourtanioti, E., Gökmen, B., Ilgaz, Ç., Kumlutaş, Y., Avci, A., & Lymberakis, P. (2013). Comparative phylogeography of six herpetofauna species in Cyprus: Late Miocene to Pleistocene colonization routes. *Biological Journal of the Linnean Society*, 108, 619–635. https://doi.org/10.1111/j.1095-8312.2012.02039.x

Psonis, N., Lymberakis, P., Poursanidis, D., & Poulakakis, N. (2016). Contribution to the study of *Acanthodactylus* (Sauria: Lacertidae) mtDNA diversity focusing on the *A. boskianus* species group. *Mitochondrion*, 30, 78–94. https://doi.org/10.1016/j.mito.2016.07.001

Rahbek, C., Borregaard, M. K., Colwell, R. K., Dalsgaard, B. O., Holt, B. G., Morueta-Holme, N., Nogues-Bravo, D., Whittaker, R. J., & Fjeldså, J. (2019). Humboldt's enigma: What causes global patterns of mountain biodiversity? *Science*, 365, 1108–1113. https://doi.org/10.1126/science.aax0149

Rambaut, A. (2016). FigTree v1.4.3. Retrieved from: http://tree.bio.ed.ac.uk/software/figtree/

Rambaut, A., Drummond, A. J., Xie, D., Baele, G., & Suchard, M. A. (2018). Posterior summarisation in Bayesian phylogenetics using Tracer 1.7. *Systematic Biology*, 67, 901.

Salvador, A. (1982). A revision of the lizards of the genus *Acanthodactylus* (Sauria: Lacertidae). *Bonner Zoologische Monographien*, 16.

Santos, A. M., Cabezas, M. P., Tavares, A. I., Xavier, R., & Branco, M. (2015). tcsBU: A tool to extend TCS network layout and visualization. *Bioinformatics*, 32, 627–628. https://doi.org/10.1093/bioinformatics/btv636

Shaibi, T., & Moritz, R. F. A. (2010). 10,000 years in isolation? Honeybees (*Apis mellifera*) in Saharan oases. *Conservation Genetics*, 11, 2085–2089.

Simó-Riudalbas, M., Metallinou, M., de Pous, P., Els, J., Jayasinghe, S., Péntek-Zakar, E., Wilms, T., Al-Saadi, S., & Carranza, S. (2017). Cryptic diversity in *Ptyodactylus* (Reptilia: Gekkonidae) from the northern Hajar Mountains of Oman and the United Arab Emirates uncovered by an integrative taxonomic approach. *PLoS One*, 12, e0180397. https://doi.org/10.1371/journal.pone.0180397

Sindaco, R., & Jeremčenko, V. K. (2008). *The Reptiles of the Western Palearctic: Annotated checklist and distributional atlas of the turtles, crocodiles, amphisbaenians and lizards of Europe, North Africa, Middle East and Central Asia*. Edizioni Belvedere.

Skonieczny, C., McGee, D., Winckler, G., Bory, A., Bradtmiller, L. I., Kinsley, C. W., Polissar, P. J., De Pol-Holz, R., Rossignol, L., & Malaizé, B. (2019). Monsoon-driven Saharan dust variability over the past 240,000 years. *Science Advances*, 5, eaav1887. https://doi.org/10.1126/sciadv.aav1887

Skonieczny, C., Paillou, P., Bory, A., Bayon, G., Biscara, L., Crosta, X., Eynaud, F., Malaizé, B., Revel, M., Aleman, N., Barusseau, J.-P., Vernet, R., Lopez, S., & Grousset, F. (2015). African humid periods triggered the reactivation of a large river system in Western Sahara. *Nature Communications*, 6, 1–6. https://doi.org/10.1038/ncomms9751

Šmid, J., Sindaco, R., Shobrak, M., Busais, S., Tamar, K., Aghova, T., & Carranza, S. (2021). Diversity patterns and evolutionary history of Arabian squamates. *Journal of Biogeography*, 1–17.

Šmid, J., Uvizl, M., Shobrak, M., Salim, A. F. A., AlGethami, R. H. M., Algethami, A. R., Alanazi, A. S. K., Alsubaie, S. D., Busais, S., & Carranza, S. (2020). Swimming through the sands of the Sahara and Arabian deserts: Phylogeny of sandfish skinks (Scincidae, *Scincus*) reveals a recent and rapid diversification. *Molecular Phylogenetics and Evolution*, 155, 107012. https://doi.org/10.1016/j.ympev.2020.107012

Soberón, J., & Peterson, A. T. (2005). Interpretation of models of fundamental ecological niches and species' distributional areas. *Biodiversity Informatics*, 2, 1–10. https://doi.org/10.17161/bi.v2i0.4

Suchard, M. A., Lemey, P., Baele, G., Ayres, D. L., Drummond, A. J., & Rambaut, A. (2018). Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10. *Virus Evolution*, 4, vey016. https://doi.org/10.1093/ve/vey016

Tamar, K., Carranza, S., Sindaco, R., Moravec, J., & Meiri, S. (2014). Systematics and phylogeography of *Acanthodactylus schreiberi* and its relationships with *Acanthodactylus boskianus* (Reptilia: Squamata: Lacertidae). *Zoological Journal of the Linnean Society*, 172, 720–739. https://doi.org/10.1111/zoj12170

Tamar, K., Carranza, S., Sindaco, R., Moravec, J., Trape, J. F., & Meiri, S. (2016). Out of Africa: Phylogeny and biogeography of the widespread genus *Acanthodactylus* (Reptilia: Lacertidae). *Molecular Phylogenetics and Evolution*, 103, 6–18. https://doi.org/10.1016/j.ympev.2016.07.003

Tamar, K., Metallinou, M., Wilms, T., Schmitz, A., Crochet, P. A., Geniez, P., & Carranza, S. (2017). Evolutionary history of spiny-tailed lizards (Agamidae: *Uromastyx*) from the Saharo-Arabian region. *Zoologica Scripta*, 47, 159–173.

Tamar, K., Scholz, S., Crochet, P.-A., Geniez, P., Meiri, S., Schmitz, A., Wilms, T., & Carranza, S. (2016). Evolution around the Red Sea: Systematics and biogeography of the agamid genus *Pseudotrapelus* (Squamata: Agamidae) from North Africa and Arabia. *Molecular Phylogenetics and Evolution*, 97, 55–68. https://doi.org/10.1016/j.ympev.2015.12.021

Tamura, K., Stecher, G., Peterson, D., Filipski, A., & Kumar, S. (2013). MEGA6: Molecular evolutionary genetics analysis version 6.0. *Molecular Biology and Evolution*, 30, 2725–2729. https://doi.org/10.1093/molbev/mst197

Tejero-Cicuéndez, H., Patton, A., Caetano, D., Smid, J., Harmon, L., & Carranza, S. (2021). reconstructing squamate biogeography in Afro-Arabia reveals the influence of a complex and dynamic geologic past. *Systematic Biology*, 1–11.

Thuiller, W., Guéguen, M., Renaud, J., Karger, D. N., & Zimmermann, N. E. (2019). Uncertainty in ensembles of global biodiversity scenarios. *Nature Communications*, 10, 1–9. https://doi.org/10.1038/s41467-019-09519-w

Title, P. O., & Bemmels, J. B. (2018). ENVIREM: An expanded set of bioclimatic and topographic variables increases flexibility and improves performance of ecological niche modelling. *Ecography*, 41, 291–307.

Trape, J. F., Chirio, L., & Trape, S. (2012). *Lézards, crocodiles et tortues d'Afrique occidentale et du Sahara*. IRD éditions.

Uetz, P., Freed, P., & Hosek, J. (2021). The reptile database. Retrieved from: http://www.reptile-database.org/

Vale, C. G., Pimm, S. L., & Brito, J. C. (2015). Overlooked mountain rock pools in deserts are critical local hotspots of biodiversity. *PLoS One*, 10, e0118367. https://doi.org/10.1371/journal.pone.0118367

Vale, C. G., Tarroso, P., & Brito, J. C. (2014). Predicting species distribution at range margins: Testing the effects of study area extent, resolution and threshold selection in the Sahara-Sahel transition zone. *Diversity and Distributions*, 20, 20–33. https://doi.org/10.1111/ddi.12115
<figure><img src="image_17.png" /></figure>

Velo-Antón, G., Godinho, R., Campos, J. C., & Brito, J. C. (2014). Should I stay or should I go? Dispersal and population structure in small, isolated desert populations of West African crocodiles. *PLoS One*, 9, e94626. https://doi.org/10.1371/journal.pone.0094626

Velo-Antón, G., Martínez-Freiría, F., Pereira, P., Crochet, P. A., & Brito, J. C. (2018). Living on the edge: Ecological and genetic connectivity of the spiny-footed lizard, *Acanthodactylus aureus*, confirms the Atlantic Sahara desert as a biogeographic corridor and centre of lineage diversification. *Journal of Biogeography*, 45, 1031-1042.

Vignali, S., Barras, A. G., Arlettaz, R., & Braunisch, V. (2020). SDMtune: An R package to tune and evaluate species distribution models. *Ecology and Evolution*, 1–18.

Ward, D. (2016). *The biology of deserts*. Oxford University Press.

Zhang, Z., Ramstein, G., Schuster, M., Li, C., Contoux, C., & Yan, Q. (2014). Aridification of the Sahara desert caused by Tethys Sea shrinkage during the Late Miocene. *Nature*, 513, 401–404. https://doi.org/10.1038/nature13705

on comparative phylogeography. All authors are interested in the assessment of biodiversity patterns and processes in deserts and arid regions.

**Authors' contributions:** AVL, DR, DVG, JCB and PAC conceived the ideas; AVL, DVG, GVA, JCB, PAC and PG collected the data; AVL, DR, DVG, GVA and MF analysed the data; AVL led the writing with inputs from all authors.

## SUPPORTING INFORMATION

Additional supporting information may be found online in the Supporting Information section.

## BIOSKETCH

André Vicente Liz is a Ph.D. candidate at CIBIO-InBIO (University of Porto, Portugal) and Zoological Research Museum Alexander König (Bonn, Germany), and is interested in biogeographic approaches to conservation. His Ph.D. project aims to understand the role of climate-induced dispersal/vicariance as main biodiversity diversification drivers in the Sahara-Sahel desert, based

**How to cite this article:** Liz, A. V., Rödder, D., Gonçalves, D. V., Velo-Antón, G., Fonseca, M. M., Geniez, P., Crochet, P.-A., & Brito, J. C. (2021). The role of Sahara highlands in the diversification and desert colonization of the Bosc's fringe-toed lizard. *Journal of Biogeography*, 48, 2891–2906. https://doi.org/10.1111/jbi.14250