28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 The desert vipers of the genus Cerastes are a small clade of medically important venomous snakes within the family Viperidae. According to published morphological and molecular studies, the group is comprised by four species: two morphologically similar and phylogenetically sister taxa, the African horned viper (Cerastes cerastes) and the Arabian horned viper ( Cerastes gasperettii); a more distantly related species, the Saharan sand viper ( Cerastes vipera), and the enigmatic Böhme9s sand viper (Cerastes boehmei), only known from a single specimen in captivity alle gedly captured in Central Tunisia. In this study, we analyzed one mitochondri al marker (COI) as well as genome-wide data (ddRAD sequencing) from 28 and 41 samples, respectively, covering the entire distribution range of the genus to explore the population genomics, phylogenomic relationships and introgression patterns within the genus Cerastes. Additionally, and to provide insights into the mo de of diversification of the group, we carried out niche overlap analyses consid ering climatic and habitat variables. Both nuclear phylogenomic reconstructions and population structu re analyses have unveiled an unexpected evolutionary history for the genus Cerastes, which sharply contradicts the morphological simil arities and previously published mitochondrial approaches. Cerastes cerastes and C. vipera are recovered as sister taxa whilst C. gasperettii is a sister taxon to the clade formed by these two species. We found a relatively high niche overlap (OI > 0.7) in both cl imatic and habitat variables between C. cerastes and C. vipera, contradicting a potential scenario of sympatric s peciation. These results are in line with the introgression found between the northwestern African populations of C. cerastes and C. vipera. Finally, our genomic data confirms the existence of a lineag e of C. cerastes in Arabia. All these results highlight the importance of genome-wide data over few genetic markers to study the evolutionary history of species. 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77

Introduction Species-level phylogenies are crucial for understan ding evolution and speciation (Barraclough & Nee, 2001) . An accurate knowledge of phylogenetic relati onships is critical to study multiple evolutionary processes such as interspecific diversification, tr ait evolution or adaptative responses of species to changing environments, among others (Guo et al., 20 19; Schoch et al., 2009; Zhang et al., 2022) . Traditionally, mitochondrial markers have been wide ly used in phylogenetic and/or phylogeographic approaches to infer the evolutionary history of spe cies (Brown, 2002) . However, mitonuclear discordances have been widely reported in many taxa (e.g., Dinis et al., 2019; Burriel-Carranza et al. , 2023a; Hinojosa et al., 2019; Mochales-Riaño et al. , 2023; Zaidi & Makova, 2019) , mainly due to incomplete lineage sorting, sex-biased asymmetries, introgression, or different selection pressures on mitochondrial and nuclear genomes (Toews & Brelsford, 2012) . Extensive discussion about this topic has led to the general agreement that the use of ge nomic data can efficiently resolve the evolutionary relationships of target species (e.g., Mao et al., 2019) , as well as to infer other mechanisms, like introgression, which may influence the evolutionary trajectory of species (Cahill et al., 2015; Green et al., 2010) . Remarkably, the advent of molecular gen omics has allowed for the identification and characterization of introgression (i.e., the proces s by which genetic material from one species is incorporated into the genome of another species thr ough interbreeding; Harrison & Larson, 2014) with increasing accuracy (Cahill et al., 2015; Leonard e t al., 2013 ; Otten burghs et al., 2017 ), confirming its significant impact on the evolution, adaptation and diversity of many species, including adaptive radiations (Seehausen, 2004) . Adaptive introgression (the transfer of alleles from one species to anot her favored by natural selection) has even been propose d as another potential source of adaptation, aside from new mutations or standing variation (Hedrick, 2013) , and has been reported to be involved in pesticide resistance, color polymorphism or adaptat ion to arid environments, among others (Anderson et al., 2009; Liu et al., 2021; Nanaei et al. 2023;

Rocha et al. 2023; Song et al. 2011).

Arid areas encompass a large proportion of the Eart h9s surface and are important for understanding global biodiversity patterns and especially pattern s of reptile diversity (Brito et al., 2014; Roll etal., 2017; Tejero-Cicuéndez et al., 2022a) . Arid areas a re often inhabited by specialized arid-adapted species, showing high levels of inter- and intra-sp ecific diversity and endemicity (e.g., Burriel-Carr anza 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 et al., 2023 b; Garcia-Porta et al., 2017 ; Metallino u et al., 2012, 2015; Velo-Antón et al 2018 ). Much attention has been focused on arid areas from Austr alia and North America, while North Africa and Arabia have been traditionally neglected despite an increasing number of biodiversity studies carried out in the last decade (e.g., Brito et al., 2014; Carranza et al., 2018; Metallinou et al., 2012, 2015 ; Gonçalves et al., 2018a; Martínez-Freiría et al., 2 017; amíd et al., 2021a; Tejero-Cicuéndez et al., 2022b; Burriel-Carranza et al., 2019, 2023a, 2023b) . These studies have reported high levels of cryptic diversity in these relatively underexplored regions and have emphasized vicariant processes as the primary factor driving speciation/diversification. Additionally, they have provided valuable insights into how the onset and expansion of deserts, in con junction with Plio-Pleistocene climatic cycles, hav e influenced the diversification processes of the region9s biota. However, there has been limited consideration of other evolutionary processes, such as ecological adaptation and lack of geographic isolation (as seen in sympatric speciation), which require further exploration in Palearctic deserts ( Brito et al., 2014) . In these arid regions most studies on reptiles have generally focused on lizards, while other groups, such as snakes, have been relatively neglec ted and under-explored (but see Daboia vipers, Martínez-Freiría et al., 2017 ; Psammophis schokari, Gonçalves et al., 2018 a).

The desert vipers of the genus Cerastes (Family Viperidae, subfamily Viperinae) are a mono phyletic group of medically important venomous snakes adapte d to desert environments (Phelps, 2010; Alencar et al., 2016; amíd & Tolley, 2019) . As other Palear ctic vipers, they are highly dependent on climatic conditions and are susceptible to range shifts and demographic changes driven by climatic oscillations (Brito et al., 2011; Martínez-Freiría et al., 2015, 2017; Lucchini et al., 2023) . The genus Cerastes is widely distributed across the arid regions of North Africa and the Middle East, particularly in the Sa hara Desert and the Arabian Peninsula (Phelps, 2010; Sin daco et al., 2014) . According to morphological (Wagner & Wilms, 2010; Werner et al., 1991) , ecolog ical (Brito et al., 2011) and phylogenetic studies, relying mostly on mtDNA data (Alencar et al 2016; a míd and Tolley, 2019) , the group is comprised by four species: (A) two medium-sized sister species with generalist habitat requirements, the African horned viper (C. cerastes Linnaeus, 1758) and the Arabian horned viper ( C. gasperettii Leviton & Anderson, 1967). Both species exhibit mostly allopa tric distributions, with C. cerastes primarily inhabiting the Saharan desert and having a reported isolated population in southwestern Arabia, while 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

C. gasperettii is exclusively found in arid regions within Arabia. Both species display individuals with and without horns. (B) A more distantly related spe cies, the African sand viper C. vipera (Linnaeus, 1758), specialized in the sandy environments of the Sahara Desert. This species has dorsally positione d eyes and lacks horns entirely. A fourth species, Bö hme9s sand viper (Cerastes boehmei Wagner & Wilms, 2010) is sometimes recognized (e.g., Uetz et al., 2023) but also disputed (e.g., Trape et al., 2023) . It is based on a single specimen with aberra nt horns, allegedly captured in Central Tunisia in 1991 and for which no genetic data was available so far. At the intraspecific level, molecular inferen ces have only been carried out in C. gasperettii (Carné et al., 2020), showing a low genetic variat ion. Morphological data has led to the description of se veral subspecies within the two horned vipers (i.e., C. cerastes and C. gasperettii): C. c. cerastes (widely distributed from Western Sahara to Egypt) and C. c. hoofieni (in southwestern Arabia based only on morphologica l data; Werner et al., 1991) within C. cerastes; and C. g. gasperettii (widely distributed in Arabia) and C. g. mendelssohni (Arava Valley; Werner et al., 1999) within C. gasperettii. A genomic study at the genus level would allow ad dressing the intraspecific variability and systematic of thi s group across North African and Arabia. Moreover, the combination of genomic and ecological data woul d help interpreting the evolutionary history of the group and enhance our understanding of its evolutio n and speciation history.

In this study, we used mitochondrial and double dig est restriction-site associated DNA (ddRAD) sequencing data to infer the evolutionary history w ithin the genus Cerastes and gain new insights into the intraspecific variability of

C. cerastes. Given the remarkable patterns of morphological differentiation and habitat specialization in

Cerastes species, we also implemented niche overlap analyses to explore the potential role of ecologica l adaptation in a putative process of sympatric speciation (e.g., Moutinho et al., 2020) . Finally, following the niche overlap and FineRADstructure results, we tested for ancestral introgression betw een the three species. Our main objectives were: (1) to reassess the phylogenetic relationships within t he genus Cerastes using novel genome-wide data, (2) to discuss possible drivers of speciation within the genus Cerastes across North Africa and Arabia and (3) to clarify the taxonomic status of C. cerastes. 134 135 136 137 138 139 140 141 142 143

Materials and Methods Sample collection A total of 36 samples representing the distribution range of C. cerastes (including C. c. cerastes and C. c. hoofieni), C. gasperettii (C. g. gasperettii) and C. vipera were included in this study (Fig. 1, Fig. S1 and Table S1). No samples were available for C. boehmei and C. g. mendelssohni. Moreover, we included five individuals of Echis omanensis as outgroup (amíd & Tolley, 2019) . Samples consist ed of tail tips or scale clips in live specimens (n = 37) , small fragments of body skin and muscle in roadki ll specimens (n = 1) or stored in museum collections ( n = 3) and were kept in 99% ethanol and stored at -20 ºC upon arrival. For DNA extraction, samples we re extracted using the MagAttract HMW Kit (Qiagen) following the manufacturer9s protocols. ddRADseq library preparation We followed Peterson et al., (2012 ) protocol for do uble digest restriction site-associated DNA (ddRADseq) on a total of 41 samples. We double digested 500 ng of genomic DNA using a combination of a rare and a common cutting restriction enzymes (Sbf1 and Msp1, respectively). Fragments were purified with Agencourt AMPure beads before ligatio n with eight different barcoded Illumina adapters. Samples with unique adapters were then pooled, and each pool of eight samples was size-selected for fragments sized between 415 and 515 base pairs (bp). Illumina multiplexing indices were ligated to individual samples using a Phusion polymerase kit ( high-fidelity Taq polymerase, New England Biolabs). Final pools were single-end sequenced on four lanes to a length of 75 base pairs (bp) Illumi na NextSeq 500, at the UPF Genomics Core Facility (Bar celona, Spain). ddRADseq data processing Raw reads were visually explored with FastQC (Andre ws, 2010) to ensure data quality and absence of adapters. Then, reads for each individual were demu ltiplexed and processed with ipyrad v.0.9.84 (Eaton & Overcast, 2020) , discarding sites with Phred scor e < 33, and reads with more than three missing site s.

Consensus sequences with less than six reads, exces sive heterozygous sites (more than three), or more 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 than two alleles were also discarded. We created th ree different datasets: one to explore population structure analyses and gene flow containing all ind ividuals of Cerastes together with E. omanensis (as outgroup), a second to explore more in-depth population struct ure analyses at the genus level containing only the genus Cerastes samples and a third dataset to explore intraspecifi c relationships with only C. cerastes samples. Then, for each dataset three independent i pyrad runs were performed using different clustering parameters for both filtered and consens us sequences (step 3 and 6 in ipyrad, respectively) : 0.85, 0.89 and 0.92, selecting the clustering that

maximized the number of samples without significant ly reducing the number of loci retained (0.85 with the outgroup and 0.89 for both datasets without the outgroup). Then, following recommendations from O9Leary et al., (2018 ), we used Radiator (Gosselin et al., 2017) , Plink2 (Chang et al., 2015) and VCFr (Knaus & Grünwald, 2017) implemented in a custom script (https://github.com/BernatBurriel/Post_processing_filtering) from Burriel-Carranza et al., (2023a) to filter iteratively and alternatively low -quality samples and loci. For all datasets, values of missing data allowance ranged from 98% to 60% of mi ssing genotype call rate and missing data per individual, decreasing 2% between iterations. Subse quently, we applied a hard threshold of missing genotype call rate of 60%. We removed non-biallelic SNPs, applied a minor allele frequency filter (maf < 0.05) and removed monomorphic sites. For the diff erent analyses, we retrieved loci, haplotype and SNP data (with or without an outgroup, depending on the analyses9 requirements): loci and haplotype datasets were generated with a second round of ipyr ad after removing all individuals that did not pass the previously explained iterative filters and reta ining only loci that were at least in 60% of all individuals. The SNP dataset (with the outgroup) was used for the introg ression analysis, where all previously filtered SNPs were kept. For population structure analyses, we filtered the SNP dataset to only keep unlinked SNPs, selecting the SNP with the highest depth for each locus (from now on uSNP dataset). Finally, the SNP dataset only containing C. cerastes individuals was used for a Mantel test. Population structure analyses We performed a Principal Component Analysis (PCA) w ith the unlinked SNP dataset (a total of 2,856 uSNPS) containing only Cerastes samples using Plink v1.9 (Chang et al., 2015) . Then , we used 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213

ADMIXTURE (Alexander et al., 2009) to analyze the s ame dataset to detect population ancestry from K=2 to K=5. A total of 20 replicates for each K wer e calculated, selecting the best K value after 20 cross validations. We also used fineRADstructure (M alinsky et al., 2018) on the haplotype dataset to explore population structure within the three species, plus the outgroup, setting the minimum number of samples per loci to four. This analysis is optim al to detect different levels of structure within a nd between populations and to explore co-ancestry patt erns for non-model organisms (Malinsky et al., 2018) . We also performed a Mantel test with the

Cerastes cerastes samples (SNP dataset containing only C. cerastes samples). Lastly, we calculated genetic diversity a nd pairwise Fst values with GenoDive v.3.0 (Meirmans, 2020) with 10,000 permutations at the species and subspecies level, using the SNP dataset. Visualization of results from these analys es was performed with R v.3.6.3 (R Core Team, 2021) and the R package ggplot2 (Wickham, 2016) .

COI sequencing and phylogenetic inference We sequenced a COI fragment for a total of 24 individuals for the three species within the Cerastes genus (information regarding the specimens used can be found in Table S1). Primers and PCR amplification conditions were implemented according to Nagy et al., (2012 ). PCR products were purified and Sanger sequenced by Genewiz (UK). Then , we performed a multiple sequence alignment with MAFFT (Katoh & Standley, 2013) for the 24 sam ples together with two previously sequenced samples (13693 and 12953, accession number ON943580 and ON943571, respectively, from VeloAntón et al., 2022 ) as well as two outgroups (two i ndividuals of Echis carinatus, accession numbers MG699965 and MG699966, respectively). We used BEAST 2 v.2.6.4 (Bouckaert et al., 2019) to reconstruct a Bayesian Inference (BI) time-calibrat ed tree. We calibrated the phylogeny dating the deepest node in the phylogeny (between Echis and Cerastes), as suggested by Stange et al., (2018) and also implemented in previous studies focused on reptiles (Burriel-Carranza et al., 2023a, 2023b; Thano u et al. 2023) , with a normal distribution from a mea n age of 28 million years ago (Mya) (amíd & Tolley , 2019) and encompassing a 95% HPD interval (25.65-3 1.55 Mya). We selected a GTR model with four gamma categories (base frequencies and proportion o f invariant sites were estimated), and a relaxed 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 clock LogNormal was used with a coalescent constant population tree prior. We conducted four independent runs of 10 9 generations sampling every 10 4 generations. Convergence was checked with Tracer v.1.7.2 (Rambaut et al., 2018) . Posterior distributions were combined with LogCombiner v.2.6.3, discarding 10% of the posterior trees as burn-in an d a maximum clade credibility tree was obtained calculating median heights in TreeAnnotator v.2.6.3 (BEAST2 v.2.6.4; Bouckaert et al., 2019) . Phylogenomic reconstructions We reconstructed the phylogenomic history of the three species of Cerastes using three different phylogenetic methods. First, with the concatenated loci dataset (531,260 bp, including the outgroup) we reconstructed a Maximum Likelihood (ML) tree. Fo r that, we used RaxML-NG (Kozlov et al., 2019) , with the GTR+I+G model with a total of 100 s tarting trees (50 random and 50 parsimony) and performing 1,000 bootstraps for branch support. Lat er, we filtered the loci dataset to maintain only loci without missing data, obtaining a total of 627 loci and 40,831 bp (including the outgroup). With those, we used BEAST2 v.2.6.4 (Bouckaert et al., 2019) to obtain a time-calibrated tree. We used the same calibration point as described for the mitochondria l phylogeny. We selected a GTR model with four gamma categories (base frequencies and proportion o f invariant sites were estimated), and a relaxed clock LogNormal was used with a coalescent exponent ial population tree prior. We conducted four 8 independent runs of 10 generations sampling every 10 4 generations. Posterior distributions were combined with LogCombiner v.2.6.3. Convergence was assessed with Tracer v.1.7.2 (Rambaut et al., 2018) , applying a 10% burn-in. Finally, and to conf irm the unexpected phylogenomic relationships (see results), we calculated a species tree using SNAPP (Bryant et al., 2012) . This method allows each SNP to have its own history under the multispecies coal escence model, while bypassing the need to sample individual gene trees (Bryant et al., 2012) . Due to computational constraints, we ran an extra ipyrad selecting only two individuals per species (Table S 1) and obtaining a total of 6,668 SNPs (including t he outgroup), allowing a total of 20% of missing data with each SNP present in at least one individual pe r species. We used the same calibration point as desc ribed for the mitochondrial phylogeny. We ran four different SNAPP analyses for 1,000,000 Markov chain Monte Carlo (MCMC) generations, sampling 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 every 1,000 steps, setting mutation rates ( u and v) to 1.0 and using different gamma prior distributi ons for alpha and beta (2, 200; 2, 2,000; 2, 20,000). P osterior distributions were combined with LogCombiner v.2.6.3. Convergence to stationary dist ributions was observed after 10% generations, and hence 10% of posterior trees were removed as burn-i n with Tracer v.1.7.2 (Rambaut et al., 2018) and the posterior distribution was considered adequatel y sampled with an effective sample size higher than 200 for all parameters. Trees were combined with Lo gCombiner v.2.6.3 and visualized in DensiTree v2.2.6 (Bouckaert, 2010) . Alternative prior combina tions produced highly concordant results. We produced a second species tree under the same condi tions but including two specimens of C. c. hoofieni to the same dataset (obtaining 6,379 SNPs) to study possible different evolutionary histories of the two

Niche overlap analyses In order to investigate the potential role of envir onmental adaptation in the diversification of the g enus, we performed pairwise niche overlap tests between t he major clades within the genus, resulting from our ddRADseq inferences, considering two geographic scales and two types of variables. The global scale analyses included a total of 990 records at a spatial resolution of 5 arc minutes (~10x10km; 561 C. cerastes, 213 C. gasperettii and 216 C. vipera) spanning across the whole distribution of the thr ee species (Fig. S2), and two distinct sets of low correlated variables (R < 0.72; Table S2): (1) a clima tic dataset, which included seven precipitation- and te mperature-related variables downloaded from WorldClim v.2.1 (www.worldclim.org), and (2) a land cover dataset, which comprised three percentages of landcover variables, which resulted from calcula ting the number of 10 arc seconds resolution pixels of three landcover categories representative of the Saharan and Arabian deserts (from the ESA GlobCover ver 2.2, Bicheron et al., 2008) in 5 arc minutes resolution pixels, using the function aggregate from ArcGIS v.10.5 (ESRI, 2016). Occurren ce data for the global scale analysis were gathered from distinct sources, including fieldwork conducted by the authors for the period 2001-2021 (n = 241), bibliographic references (n = 180; e.g., Werner et al., 1999; Brito et al., 2011; GarcíaCardenete et al., 2017) , museum collections (n = 10 6; e.g., NHM-London, MNHN-Paris) and GBIF (n 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 = 463; www.gbif.org, 2021). Due to the underrepresentation of the landcover heterogeneity for the Saharan and Arabian deserts as available in the ESA GlobCover dataset, we also approached niche overlap for the Saharan species at a regional scale of the Atlantic Sahara (i.e., Mauritania and Weste rn Sahara), where a detailed landcover map was recentl y derived (Campos et al., 2018) . Consequently, for this regional scale analysis, we considered a subse t of 170 records (118 for C. cerastes and 52 for C. vipera), available at a higher accuracy (i.e., 30 arc sec onds, ~1x1km), since 69% were gathered during fieldwork campaigns (the remaining from bibliograph y, 19% and GBIF, 12%). In addition, we used seven percentages of landcover variables (see Table S2), which resulted from using the function aggregate from ArcGIS over seven landcover categori es or groups of categories with 1 arc second resolution pixel (see Table S2) to new rasters with 30 arc seconds resolution pixel.

For the three approaches (global-climate, global-la ndcover and regional-landcover), occurrence data were used to extract the corresponding variabl e values and then run three distinct PCA. Following a three-dimensional hypervolume approach, the first three components of each PCA were used to describe the environmental niche of each major clad e (i.e., total volume and unique part, intersection and union with the other species9 volume) as obtained from ddRAD data analyses and perform pairwise niche-overlap tests between them (see Martínez-Frei ría et al., 2020 for a similar procedure) . Analyses were performed in R with the 8hypervolume9 package (Blonder et al., 2014) , using a Silverman bandwidth estimator and a set of 999 random points to sample the kernel density to build hypervolumes. To quantify environmental niche overlap, the Sørens on index (K) and the Overlap Index (OI) were used (e.g., Martínez-Freiría et al., 2020) . These statis tics quantify predicted niche similarity, and they both range from 0 (no overlap) to 1 (identical niche). T he Overlap Index (OI) relates the observed and maximum values of K, and is preferred when niches p resent different sizes (Simó-Riudalbas et al., 2018) . In order to test the significance of the ove rlap, a randomization procedure similar to that app lied by Warren et al., (2010 ) was used. By comparing the produced Ks and OIs, we tested whether the observed niche overlap was more different than the overlap that would be obtained by merging occurrences randomly. This procedure was performed 999 times to generate the null distribution, and to evaluate the significance of each overlap. 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322

Introgression analyses To study ancestral introgression, we ran the ABBA-B ABA test implemented in Dsuite (Malinsky et al., 2021) with the SNP dataset (32,254 SNPs), including the outgroup. Todo so, we split our samples per species (and subspecies when possible) and calculat ed the D-statistics using the new phylogenetic topology. Patterson9s D-statistic uses asymmetry in gene tree topologies to quantify introgression between either of two lineages which share a common ancestor (P1 and P2) and a third lineage (P3) that diverged from the common ancestor of P1 and P2 . A fourth individual (P4) is needed to inform the ancestral state for each comparison (Malinsky et al ., 2021) . In here, we set E. omanensis as P4 to set the ancestral state. C. gasperetti or C. vipera were used as P3 while C. c. cerastes, C. c. hoofieni or C. vipera were set as P1 and P2 terminals, depending on the comparison. Finally, we calculated the Fbranch statistics to infer the excess of shared all eles between the groups, using E. omanensis as outgroup.

Results ddRADseq data processing Genome-wide data was sequenced for a total of 41 samples, including five individuals of E. omanensis from regions spanning Western Sahara to Oman (Fig. 1, Fig. S1 and Table S1). However, after the iterative filtering and due to the low DNA quality of some samples, we kept a total of 28 high-quality samples: four C. vipera, seven C. gasperettii, thirteen C. cerastes (seven C. c. cerastes and six C. c. hoofieni) and four Echis omanensis (Fig. S1 and Table S1). Total reads obtained from sequencing added 7 up to 7.7 x 10 . After applying quality filtering, more than 99% o f the raw reads remained, with an 6 average of 2.7 x 10 reads per individual (Table S1).

Population structure analyses Population genomic results of the unlinked SNP dataset (2,856 uSNPs) showed an unexpected result (Fig. 2). PC1 split C. gasperettii from both C. cerastes and C. vipera individuals whilst PC2 separated C. cerastes from C. vipera (Fig. 2A). The C. c. hoofieni individuals clustered together with the rest of 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348

African C. cerastes (Fig. 2A). PC3 sorted C. cerastes following their geographical locations (northwestern Africa, northeastern Africa and Arabi a). The best supported scenario in Admixture recovered the three different species for K=3 (Fig. 2B. For K=4, the samples ofC. c. hoofieni formed a unique cluster, with the two geographically closer African individuals (92 and 12953; Table S1) sharin g a limited amount of coancestry (Fig. 2B). FineRADst ructure analyses showed five different clusters (including the outgroup) representing the different species and with C. c. hoofieni having a different cluster (Fig. 2C). Again, C. gasperettii was the most differentiated within the Cerastes genus and showed high similarity among samples. Cerastes vipera also displayed a high similarity among samples, with only the individual from central Nort h Africa being less similar than the rest. Cerastes cerastes individuals showed three main clusters correspondin g with their geographical locations (northwestern Africa, northeastern Africa and Arabi a). Cerastes vipera samples reported higher shared co-ancestry with the African samples of C. c. cerastes than with the Arabian individuals of C. c. hoofieni (Fig. 2C). Mantel test analysis showed strong isola tion by distance among C. cerastes individuals (r = 0.69, p-value < 0.001) (Fig. S3). The F st comparisons revealed consistently high values in a ll pairwise comparisons, indicating a strong isolation among sp ecies (Table S3), being similar between C. c. cerastes and both C. c. hoofieni and C. vipera. Finally, genetic diversity measures consistently indicated lower observed diversity compared to the expected values (Table S4), withC. c. hoofieni displaying the lowest levels of heterozygosity. 6.48 Mitochondrial phylogenetics The phylogenetic reconstruction of the mitochondrial COI gene fragment split the three different Cerastes species (Fig. 3A). Cerastes vipera represented the deepest split within the genus, dating 12.78 ± 5.62 Mya, while C. cerastes and C. gasperettii were recovered as sister taxa, splitting 9.19 ± 4.25 Mya. Within C. cerastes, there were two lineage splitting 3.76 ± 2.6 Mya: a first one containing individuals from Morocco, Western Sahara, Mauritani a and Algeria, and a second lineage with the individuals from Chad and Egypt together with the A rabian C. c. hoofieni. All individuals within C. 350 351 gasperettii clustered together with an individual from souther n Oman showing a slight separation. Individuals of C. vipera formed a separated lineage.

Phylogenomic reconstruction The ML phylogenomic tree, with 531,260 bp, showed C. vipera and C. cerastes as sister species, whilst C. gasperettii was recovered as the basal species (bootstrap valu es higher than 95%; Fig. S4). Within C. cerastes there were three different lineages: a northwester n African lineage ( C. c. cerastes), a northeastern African lineage ( C. c. cerastes) with two individuals from Chad and Egypt and a th ird lineage with the Arabian samples ( C. c. hoofieni) sister to the northeastern African lineage. For C. vipera, the topology followed a geographical distribution of the four individuals sampled, with closer individuals clustering together. Conversely, in C. gasperettii there were two different lineages, splitting coastal and interior individuals from the Arabian P eninsula (Fig. S4). The time-calibrated Bayesian tree further corroborated the ML topology for the genus and provided support for the intraspecific relationships (Fig. 3B). The split between C. gasperettii and the ancestor of both C. cerastes and C. vipera was estimated around 7.12 ± 3.83 Mya and between C. cerastes and C. vipera around 4.22 ± 2.23 Mya (Fig. 3B). The genomic species tree inferred with SNAPP yielded the same new topology (Fig. S5), as well as the tree including C. c. hoofieni (Fig. S6) and divergence among species were found within the 95% interval of our time-calibrated Baye sian tree. Overall, intraspecific relationships wer e found to be consistent between mitochondrial and nu clear phylogenies (Fig. 3).

Niche overlap analyses The first three components accounted for more than 81% of the environmental variance in the two global PCAs, while this value was 60% in the region al-landcover PCA (Table S5). In the global-climate approach, C. gasperettii showed a much-reduced volume in comparison to the volume exhibited by the clade comprising C. cerastes and C. vipera, and also to the volumes recovered for each specie s. Conversely, in the global-landcover approach, C. gasperettii showed a similar volume to the volumes of the other two species, while it was more reduced when both species were considered as a single clad e 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 (Table S5). In the three approaches, C. cerastes and C. vipera showed rather similar volumes (Table

Niche volume overlap was high in all cases (Fig. 4), with both K and OI indexes ranging from moderate to high (0.49 < K < 0.81; 0.70 < OI < 0.84 ; Table S5). Statistical tests were all significant ( p < 0.05), with the exception of the test of K index for C. cerastes vs C. vipera in the global-landcover dataset (Table S5). Overall, niche overlap analyses show that the three species exhibit non-equivalent but highly overlapping niches.

Introgression analyses We used the filtered SNP dataset containing 32,265 SNPs that were present i n at least 60% of the samples to study introgression between the different clades. Both D-statistics and f4-ratio analyses showed a significant excess of shared alleles betwe en C. vipera and C. c. cerastes individuals (Table 1 and Fig. 5). Moreover, f-branch analyses supported D-statistics and f4-ratio and suggested an introgression of around 10% of the sampled sites be tween C. vipera and C. c. cerastes (Fig. 5). Ancestral introgression was not found between C. gasperettii and C. c. hoofieni, both occurring in Arabia (Table 1 and Fig. 5).

Discussion In this study, we used mitochondrial and genome-wid e data to unveil the phylogenetic relationships of the genus Cerastes, a clade of medically important venomous snakes ad apted to arid environments and mainly distributed across the Saharan and Arabian d eserts. This group has been historically understudied due to their wide distribution and harsh an d remote environments it inhabits, which hampered a systematic field sampling. Our results unravel an unexpected phylogenomic relationship that differs from previous morphological and mitochondrial appro aches, which suggested that C. cerastes and C. gasperettii, two morphologically similar species with shared hab itat requirements, were sister taxa while, C. vipera, smaller and linked to soft sand habitats, was sis ter to them. Moreover, niche overlap analyses did not seem to support a potential scenar io of sympatric speciation between C. cerastes and 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430

C. vipera. Ultimately, our study underscores the significanc e of employing genome-wide data in the examination of speciation processes and the unravel ling of the complete evolutionary history of species. It also highlights the contrasting patterns observe d between genomic data with both morphology and mitochondrial markers.

Unexpected phylogenomic relationships within the genus Cerastes All our phylogenomic analyses support a new phyloge nomic relationship within the genus Cerastes. According to our findings, C. cerastes and C. vipera are sister taxa, while C. gasperettii is more distantly related (Fig. 3B and Fig. S4-S6). This contrasts wtih the previously accepted sister relationship betw een C. cerastes and C. gasperettii, with C. vipera considered as an external group relative to this c lade (Alencar et al., 2016; amíd & Tolley, 2019; Fig. 3A ) . Population structure analyses also showed a similar clustering (Fig. 2). This was an unexpected result, since C. cerastes and C. gasperettii exhibit greater morphological similarities, both with horn and hornless populations, a similar body size, and overlapping habitat requirements, since both specie s are generalists (Werner et al., 1991; Brito et al ., 2011) . Conversely, C. vipera is smaller, strictly hornless and specialized in bu rrowing in sandy substrates, with specific adaptations such as the d orsal position of the eyes (Young & Morain, 2003; Sivan et al., 2013) . This morphological-genomic con trasting pattern has also been reported for other species with complex evolutionary histories present in desert environments as foxes (Leite et al., 201 5; Rocha et al. 2023) .

Prior to our work, all phylogenetic knowledge on these desert-adapted vipers pointed towards C. cerastes and C. gasperettii being sister taxa (Pyron et al. 2013; Alencar et al ., 2016; amíd & Tolley, 2019 and Fig. 3A) . However, these studies predomina ntly relied on a variety of mitochondrial markers, since nuclear data has not been traditionally consi dered within the Viperinae subfamily due to its low variability (e.g., in Eurasian vipers, Freitas et a l., 2020) . One plausible hypothesis to consider is that the introgression observed between C. c. cerastes and C. vipera could have altered the phylogenetic position of C. vipera, potentially bringing it closer to C. c. cerastes. However, if that was the case, an incorporation of C. c. hoofieni into a phylogenetic framework should result in C. cerastes as paraphyletic (with C. c. cerastes and C. vipera clustering together and C. c. hoofieni closer to C. 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 gasperettii), as we did not report introgression between C. vipera and the Arabian subspecies. Nonetheless, this was not the case as coancestry an alyses would show more shared ancestry between C. cerastes and C. gasperettii (Fig. 2C) and the phylogenomic trees incorporating C. c. hoofieni also showed a common evolutionary history within C. cerastes, as both subspecies clustered together (Fig. 3, Fig. S4 and S6).This points towards incomplete lineage sorting of t he mitochondria as the main factor contributing to the observed mitonuclear dis cordance described here. As it has been reported elsewhere, the use of nuclear (and especially genom e-wide data) has been a key factor in unraveling th e phylogenomic relationships among species (e.g., Mao et al., 2019) . Mitochondrial data might not represent the complete evolutionary history of a sp ecies, as it only shows the evolutionary history of the maternally inherited mitochondria. Therefore, i t is recommended to combine NGS techniques together with nuclear data to effectively address t hese questions.

Vicariance, the geographic separation of connected populations eventually leading to different species due to the emergence of a physical barrier, is cons idered one of the key mechanisms driving speciation in low-dispersal reptiles from arid environments (B rito et al., 2014, Tejero-Cicuéndez et al., 2022b) . This type of speciation has been repeatedly reporte d for species from the Arabian Peninsula and North Africa, as these areas present several geographical (e.g., environmental or physiographical) barriers that have had a determinant impact on the distribution a nd divergence of their reptile fauna (Brito et al., 2014; Metallinou et al., 2015; Garcia-Porta et al. 2017; Gonçalves et al., 2018a; 2018b; Liz et al., 2 021; 2023; Tejero-Cicuéndez et al., 2022b; Burriel-Carra nza et al., 2023b) . Among them, deserts are an interesting geographical feature, as they can act b oth as a geographical barrier or a corridor, depend ing on the species9 ecological preference (Brito et al., 2014) . On the one hand, deserts can serve as effective geographical barriers for several reptile species w ith rather mesic requirements, preventing gene flow between the populations on either side and leading to different evolutionary lineages (Rato et al., 20 07; Martínez-Freiría et al., 2017; Gonçalves et al., 20 18a; 2018b; Liz et al. 2021) . On the other hand, de serts have facilitated the dispersal of arid-adapted taxa , expanding the range of xeric species during perio ds of increasing aridity (Amer & Kumazawa, 2005; Brito et al., 2014; Velo-Antón et al., 2018 ; Liz et al. 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 Hemidactylus, Acanthodactylus, Echis, Pseudotrapelus or Uromastyx dispersed from Arabia to Africa, confirming the establishment of a dispersal corrido r connecting Africa and Eurasia through the Arabian Peninsula (Tejero-Cicuéndez et al., 2022b) . These d ispersal events match with a progressive aridification of North Africa and Arabia during the mid-late Miocene, with an hyperaridity period in Arabia, together with the more recent phases of the formation of the Sahara Desert, which began around 7 Mya in northeastern Africa (Schuster et al., 2006 ; Böhme et al., 2021; Crocker et al., 2022) . Assuming an Arabian origin for the genus Cerastes, the rise of aridity and the formation of the Saha ra Desert during the mid-late Miocene could have facil itated the colonization of North Africa by the arid adapted ancestor of C. cerastes and C. vipera, a phenomenon previously observed in other species of reptiles (Amer & Kumazawa, 2005; Metallinou et al., 2012; Liz et al. 2021) . In fact, this period coincides with the split of C. gasperettii and the ancestor of C. cerastes and C. vipera. Although Africa and Arabia started splitting around 30-25 Mya (Boha nnon 1986; Ghebreab 1998; Bosworth et al., 2005) these two geographical entities were still in conta ct through the Sinai Peninsula (Bosworth et al., 20 05; Tejero-Cicuéndez et al., 2022b) . The subsequent for mation of the Red Sea and the Gulf of Aden likely isolated C. gasperettii in Arabia, as these geographical barriers have caus ed strong vicariance speciation in other reptile species (Metallinou et al., 2012; Pook et al. 2009; Tejero-Cicuéndez et al., 2022b) . The dispersal of Cerastes from Arabia to Africa is relatively recent when co mpared to other taxa such as Stenodactylus or Echis, whose estimated colonization events dated around 21.8 and 19.4 Mya, respectively (Metallinou et al., 2012; Pook et al., 2009) . However, there have been other dispersal ev ents of species between Africa and Arabia that suggest t hat dispersion have also occurred at different time s (e.g., 9.8 and 5.9 Mya for Hemidactylus, 8 Mya for Echis pyramidum, 4 Mya forBitis arietans or 1.75 Mya for Naja haje; Pook et al., 2009; amíd et al., 2013) . After the colonization of Africa, t he dry-wet cycles that the Sahara experienced during the Plio- Pleistocene (Le Houérou, 1997; Brito et al., 2014) could have caused an isolation of the widespread Cerastes ancestor and a subsequent speciation into C. cerastes and C. vipera, around 4 Mya (Fig. 3B). Similar processes have been reported for other aridadapted species included in the genera Agama ( Gonçalves et al., 2012 ) or Scincus (Carranza et al., 2008; amíd et al., 2021b) and it is in line with our nich e overlap analyses, as we did not find strong ecolo gical 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 differences between these two species (Fig. 4). Thsi suggests that both species overlap in their habit at use and putatively contradicts a possible scenario of sympatric speciation. This is surprising, as C. vipera has been known to exhibit certain morphological ada ptations, such as the dorsal position of the eyes that allows them to burrow in sandy substrates while ambushing potential preys (Young & Morain, 2003; Sivan et al., 2013) . Conversely, C. cerastes has always been observed in rocky-like deserts (see Brito et al., 2011) . The generalist ecological nich e reported here for C. cerastes may lead to overlapping habitat use between both species (Fig. 4) or indicate that both species exploit the same ecological ni che differently, as the digging habits present in C. vipera. Nonetheless, these unexpected ecological findings could also be attributed to a general absence of sa mpling of these unexplored areas, as recent surveys have found C. vipera in zones with compact sand or C. cerastes in sandy areas without rocks (authors, pers. observation). An alternate interpretation might be that we did not manage to capture the grain of the habitat at the scale we describe it as rocky/ha rd and sandy/soft are often present side by side. The confirmation of the presence of C. cerastes in Arabia using genetic and genomic data (i.e., C. c. hoofieni), with its dispersal occurring around 1.7 Mya (Fig . 3B), could potentially be explained by either northern or southern land bridges. Currently, C. c. hoofieni is exclusively found in the coastal Tihama Desert in southwestern Arabia. In the scenario of a northern land bridge, it would suggest the species9 disappearance from the northwestern Saudi Arabia co astline. Conversely, a southern land bridge would imply its extirpation from Eritrea and Djibouti, a hypothesis that appears more plausible due to competition with other viperid snakes such as membe rs of the genus Echis. To completely unravel the colonization route of C. c. hoofieni, it is essential to conduct a more thorough sampli ng (including individuals from both Sinai and southern Sudan) as well as bioclimatic models encompassing both current and past periods to unveil potential corrid ors for the species dispersal (e.g., Gonçalves et a l., 2018 a; Velo-Antón et al., 2018 ).

Ancient introgression in North Africa We found evidence of ancestral introgression betwee n C. c. cerastes and C. vipera, which aligns with the high niche overlap reported here (Fig. 4 and 5and Table 1). This introgression is supported by th e different ancestry observed between C. c. cerastes and C. c. hoofieni with regards to C. vipera and most 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 likely is due to the allopatric distribution of C. c. hoofieni and C. vipera. FineRADstructure also supported this finding, as C. vipera individuals shared high coancestry with C. c. cerastes samples (Fig. 2C). This discovery opens up avenues for future research to investigate the introgressed DNA regions and how this introgression could have affected the populations of C. c. cerastes compared to C. c. hoofieni, as introgression has been suggested as a source o f adaptive variation (Hedrick, 2013) and has been proven fundamental for the arid adaption of so me species of mammals (Liu et al., 2021; Nanaei et al. 2023; Rocha et al. 2023) . Undoubtedly, addition al individuals of C. vipera from northeastern Africa should be sampled to verify whether this introgression is maintained along the whole species9 range. Insights into the systematics of Cerastes cerastes This study provides novel genetic and genomic information for the subspecies C. c. hoofieni. Moreover, our genome-wide data does not recognize C. c. mutila as an independent genetic unit, as it did not exhibit any genetic differences compared to other i ndividuals of C. c. cerastes (Fig. 2-3 and Fig. S4). This finding further refutes the previous discussion on the validity of C. c. mutila as a subspecies (Martínez del Mármol et al., 2019) . Conversely, our relatively small dataset seems to support the recognition of C. c. hoofieni as a subspecies within C. cerastes (Fig. 2-3 and Fig. S4), withC. c. cerastes being paraphyletic. Nevertheless, in spite of the e xistence of morphological differences (Werner et al ., 1991) , our genomic data so far does not support an elevation to a species status of C. c. hoofieni (Fig. 2-3 and Fig. S3), but future studies should address that with a bigger sampling and using an integrati ve taxonomic approach which combines morphological and climatic data together with genetics (e.g., Martínez-Freiría et al., 2021) . Interestingly, the western populations of C. cerastes show a different genetic and evolutionary history than the rest, pro mpting the need for a taxonomic revision within thi s genetic cluster. Additionally, the taxonomic status of C. g. mendelssohni should be addressed in future research, as we did not manage to obtain samples fo r this study, although Carné et al., (2020) did not find genetic differences between this and the nominotypical subspecies using mitochondrial markers.

Conclusions 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569

In this study, we have substantially increased our understanding on the evolutionary history of Cerastes from the North African and Arabian regions. We prov ide, for the first time, genome-wide information for four medically important venomous snakes from N orth Africa and Arabia (including genome-wide data for Echis omanensis), as well as an unexpected phylogenomic relationsh ip within the genus Cerastes, a crucial step for understanding the speciation p rocesses that have shaped its evolutionary history. Furthermore, we have identified distinct genetic clades within C. cerastes. This finding is of particular significance in the study of venomous sn akes, given the widespread nature of venom variatio n across different taxonomic levels (Chippaux et al., 1991; Rao et al., 2022) . Ultimately, our study emphasizes the importance of employing nuclear and specifically genome-wide data, together with ecological data, to comprehensively infer evolution ary relationships between and within species. Acknowledgements This work was supported by grants PGC2018-098290-B- I00 (MCIU/AEI/FEDER, UE), Spain, PID2021-128901NB-I00 (MCIN/AEI/10.13039/50110001103 3 and by ERDF, A way of making Europe) , Spain, and grant 2021-SGR-00751 from the D epartament de Recerca i Universitats from the Generalitat de Catalunya, Spain to SC. GM-R is supp orted by an FPI grant from the Ministerio de Ciencia, Innovación y Universidades, Spain (PRE2019 -088729), BB-C is supported by FPU grant from Ministerio de Ciencia, Innovación y Universidades, Spain (FPU18/04742), AT is supported by <la Caixa= doctoral fellowship program (LCF/BQ/DR20/117 90007), HT-C is supported by a "Juan de la Cierva

Formación" postdoctoral fellowship (FJC202 1-046832-I) funded by MCIN/AEI/10.13039/501100011033 and by the European Union NextGenerationEU/PRTR, GV-A was supported by the FCT (CEECIND/00937/2018) and r ecently by a Ramón y Cajal research grant (Ref. RYC-2019-026959-I/AEI/10.13039/501100011033) , Ja was supported by the Czech Science Foundation (GA R) under grant number 22-12757S and by the Charles University Research Centre under grant number 204069 and FM-F and JCB are supported by FCT - Fundação para a Ciencia e Tecnologia de Portugal (DL57/2016/CP1440/CT0010, CEECINST/00014/2018/CP1512/CT0001, respectively) and

Conflict of interests: The authors declare no conflict of interest.

Data availability: All new COI sequences can be accessed at NCBI with Accession number SUB13696935 and demultiplex ddRAD sequencing files can be found under BioProject ID 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 Fig. 1: On the top, distribution of the individuals sampled in the stud y, colored per species or subspecies. On the bottom, pictures of the three species considered in this study, from left to right: C. cerastes (Morocco), C. vipera (Mauritania) and C. gasperettii (United Arab Emirates). Authors: F.

Martínez-Freiría and A. Talavera.

Fig. 2: Population structure analyses for the three specie s of Cerastes, with sample colors per species or subspecies. A) Principal Component Analysis (PCA ) with a total of 2,856 uSNPs. B) Admixture analyses at the individual level for clusters with K=3 and K=4. C) Co-ancestry matrix from fineRADstructure, showing pairwise similarity. Dark er colors show higher shared co-ancestry. In bold the putative C. c. mutila individual (16CC008).

Fig. 3: Mitonuclear discordance between mitochondrial and nuclear phylogenies. A) Mitochondrial Bayesian phylogeny for the COI gene f or the species and subspecies included in the study. B) Bayesian time-calibrated tree produced wi th BEAST2 v.2.6.4 based on 627 concatenated loci and covering a total of 40,831 bp. Divergence time between C. gasperettii and the ancestor of both C. cerastes and C. vipera occurred around 7.12 ± 3.83 Mya and between C. cerastes and C. vipera around 4.22 ± 2.23 Mya. Posterior probabilities above 0.95 are indicated with an asterisk. In bold is shown the putative C. c. mutila sample.

Fig. 4: 3D niche overlap plots. Each of the five plots rep resents the three first components of a principal component analyses (PCA) of the climatic/ landcover variability for the occurrence of major clades found within Cerastes, approached at the global and regional scales. Ell ipses represent the 0.95 confidence interval for the scores of each clade in the PCA.

Fig. 5: F-branch analysis indicating an ancestral introgre ssion between the samples of C. vipera with C. c. cerastes over a total dataset of 32,254 SNPs, includingE. omanensis as an outgroup. 995 996 997 results are shown in bold. Abbreviations are as fol lows: Cch, Cerastes cerastes hoofieni; Ccc, Cerastes cerastes cerastes; Cg, Cerastes gasperettii and Cv, Cerastes vipera.

P1

P2

P3

BBAA

ABBA

BABA Cch

Ccc Cch

Ccc

Ccc Cch