Research Article
Print
Research Article
The missing piece of the puzzle: A new and widespread species of the genus Rhynchocalamus Günther, 1864 (Squamata, Colubridae) from the Arabian Peninsula
expand article infoFulvio Licata§, Lukáš Pola|, Jiří Šmíd|, Adel A. Ibrahim#, André Vicente Liz§, Bárbara Santos§, László Patkó¤, Ayman Abdulkareem¤, Duarte V. Gonçalves«, Ahmed Mohajja AlShammari», Salem Busais˄, Damien M. Egan˅, Ricardo M. O. Ramalho¦, Josh Smithsonˀ, José Carlos Brito§
‡ Universidade do Porto, Vairão, Portugal
§ CIBIO, Vairão, Portugal
| Charles University, Prague, Czech Republic
¶ National Museum, Prague, Czech Republic
# Suez University, Suez, Egypt
¤ The Royal Commission for AlUla, Oud Dunes - Amr Aldamri, Riyadh, Saudi Arabia
« University of Porto, Matosinhos, Portugal
» University of Ha'il, Ha'il, Saudi Arabia
˄ University of Aden, Aden, Yemen
˅ Natural History Collective Ventures, Wilayah Persekutuan Kuala Lumpur, Malaysia
¦ King Abdullah University of Science and Technology, Thuwal, Saudi Arabia
ˀ Prince Mohammed bin Salman Royal Reserve Development Authority, Al-Wajh, Saudi Arabia
Open Access

Abstract

Discovery rates of new species are uneven across taxonomic groups and regions, with distinctive and widely distributed species being more readily described than species with secretive habits. The genus Rhynchocalamus includes five species of secretive snakes distributed from Egypt eastwards to Iran, including the Arabian Peninsula. A wide biogeographic gap exists within the genus, which separates R. dayanae found in south Israel from R. arabicus, which occurs in the coastal areas of south Yemen and Oman. We describe Rhynchocalamus hejazicus sp. nov., a small, secretive snake, with a distinctive colouration and a melanistic morph. The new species occurs in the northwestern Hejaz region of the Kingdom of Saudi Arabia (KSA) and fills a large part of the existing distribution gap of the genus in the Arabian Peninsula. Molecular analyses of mitochondrial (12S, 16S, cytb) and nuclear genes (cmos, MC1R, NT3, RAG1) indicate that R. hejazicus sp. nov. is closely related to R. dayanae and R. arabicus, but uncertainty on the deep relationship within the genus remains. The new species has a large distribution range which potentially includes other regions in Jordan and KSA, and is associated with mountainous areas with cold wet seasons. Furthermore, it inhabits sandy and stony soils with varying vegetation cover and can be found in anthropogenically disturbed habitats, suggesting that the species should not be categorised as threatened according to IUCN criteria. The discovery of such a distinctive species highlights the existing gap in the description of rare and secretive species, and the need to enhance sampling efforts and monitoring strategies to fully capture species diversity in unexplored areas.

Key Words

Biogeography, Colubrinae, Middle East, secretive species, Serpentes, species distribution model

Introduction

Closing the gap between extant and described species is a daunting task hampered by the intrinsically slow pace of the taxonomic process and the paucity of resources therein deployed (Engel et al. 2021), and is ultimately determined by the rate at which taxonomists discover species (May 2004). Recent advancements in molecular genetic techniques have boosted new species descriptions, however, discovery rates remain uneven across taxonomic groups and regions, being determined by species’ morphology, life–history, and distribution range. For instance, small–bodied, less obvious, and small–ranged species are generally described later than bigger, more brightly coloured or widely distributed ones. Furthermore, geographical remoteness, along with behaviours and life–histories, may play a major role in determining species’ discovery, and secretive species or species from remote areas are more unlikely to be discovered (Scheffers et al. 2012).

Snakes are a diverse group of reptiles (Uetz et al. 2024) but are notoriously hard to detect due to their secretive behaviour and rarity (Rodda 1993; Kéry 2002). This is reflected in wide taxonomic and knowledge gaps, especially in small, secretive species (Vilela et al. 2014) and undersampled regions of the world, such as the Middle East, where systematic and biogeographic data for several snake taxa is still lacking.

A remarkable example and one of the most obscure Arabian snakes is a small and rather secretly living colubrid genus Rhynchocalamus Günther, 1864. The genus currently comprises five recognized species (Uetz et al. 2024): (i) R. arabicus Schmidt, 1933, from Aden, Yemen and Dhofar in south Oman; (ii) R. dayanae Tamar, Šmíd, Göçmen, Meiri & Carranza, 2016, from the Negev area, south Israel with potential presence also in the Sinai Peninsula, Egypt; (iii) R. levitoni (Torki, 2017), from the Zagros Mountains of western Iran with potential presence also in neighbouring Iraq; (iv) R. melanocephalus (Jan, 1862), from the Sinai Peninsula, Egypt, Israel, Jordan, Lebanon, Syria, Turkey, and recently discovered on Cyprus and in Medina Province, Saudi Arabia; and (v) R. satunini (Nikolsky, 1899), from Turkey, Armenia, Iraq, and Iran (Avcı et al. 2015; Šmíd et al. 2015; Tamar et al. 2016, 2020; Fathinia et al. 2017; Torki 2017; Aloufi et al. 2021).

The phylogenetic relationships within the Rhynchocalamus genus showed a relevant biogeographic gap in the species distribution, with two sister species R. dayanae and R. arabicus separated by more than 2,500 km (Tamar et al. 2016). The existing gap extends along the Kingdom of Saudi Arabia, where sampling effort is still scarce or non-existent for certain regions (e.g., Alatawi et al. 2020; Šmíd et al. 2021).

In this study, we contribute to filling taxonomic and biogeographic gaps in the reptile diversity of the Arabian Peninsula and, in particular, in the understudied genus Rhynchocalamus, by providing (i) a formal description of a new species of Rhynchocalamus snake corroborated by morphological and molecular analyses, and (ii) a distribution model of this new species and its environmental requirements, including a proposal of its conservation status. Lastly, (iii) we report the exceptional finding of a melanistic morphotype of the new species.

Materials and methods

Field sampling

Field sampling was conducted in the Kingdom of Saudi Arabia (KSA) in 2017, 2021, 2022, and 2023. Field surveys in Hail Province, KSA were conducted by AAI, AMS, and SB in May 2017 and by AMS in July 2021. Field surveys carried out in the Prince Mohammad bin Salman Royal Reserve (PMBSRR), Tabuk Province, KSA were conducted by DME and LP in March and April 2022. Field surveys in AlUla County (Medina Province) were carried out by FL, JCB, AVL, BS, and DGV in May, June, and November 2023. GPS coordinates (datum WGS84) and high–resolution photographs were taken for each individual encountered.

We collected and stored tissue samples in 96% ethanol, whereas all vouchers were first fixed in 96% ethanol and then stored in 70% ethanol. We deposited the vouchers in the herpetology collections of National History Museum Prague, Czech Republic (NMP), Muséum National d’Histoire Naturelle, Paris, France (MNHN), and AlUla Museum (RCU-URN). The complete list of tissue samples and specimens analysed in this study is listed in Suppl. material 1: table S1.

Morphological analysis

We selected the morphological characters based on previous taxonomic studies of the genus Rhynchocalamus (Tamar et al. 2016). We measured the following characters on the right side of the specimens using the image processing software ImageJ (vers. 1.53; Schneider et al. 2012) on macroscopic images: snout-vent length, measured from the tip of snout to vent (SVL); tail length measured from vent to tip of tail (TL). In addition, we collected the following pholidotic variables: number of preoculars (PreO); number of postoculars (PostO); number of temporal scales (TS); number of post-temporal scales (PTS); number of loreal scales (LS); number of ventrals (VS); number of subcaudal scales (SCS); number of upper labial scales (UL); number of lower labial scales (LL); number of lower labials in contact with anterior inframaxillar (InfLC). Lastly, we also report the shape of internasal scales (IntN; whether it was triangular or trapezoidal), and whether they were separated by the rostral scale (IntNsep).

Original photographs of all specimens in high resolution have been uploaded to the Morphobank database (https://morphobank.org/; Project no. 5111) where they are publicly available for download.

Sex determination

We used molecular markers located on the gametologous genes (i.e., homologous genes located in non–recombining regions of sex chromosomes) to determine the sex of the specimens, following the method of Laopichienpong et al. (2017). Colubrids have a female heterogametic ZZ/ZW sex-determination system (Bull 1980), with males being identified by the presence of a single band on gel electrophoresis (i.e., the alleles have the same length on the Z chromosomes), whereas females have two bands (i.e., the allele on the W sex chromosome is shorter than that on the Z chromosome; Laopichienpong et al. 2017). We amplified the gametologous gene CTNNB1 using the primers Eq–CTNNB1–11–F1 and Eq–CTNNB1–13–R (Matsubara et al. 2016; Laopichienpong et al. 2017) and following the PCR conditions reported by Laopichienpong et al. (2017).

Molecular analyses

DNA extraction, PCR, and sequence analysis

Genomic DNA was extracted from the ethanol–preserved tissue samples using the Tissue Genomic DNA Mini Kit (Geneaid, Taiwan) or the DNeasy® Blood & Tissue Kit (Qiagen, Germany). We PCR–amplified up to seven genetic markers, three mitochondrial (mtDNA): ribosomal 12S rRNA (12S) and 16S rRNA (16S), cytochrome b (cytb); and four nuclear (nDNA), the oocyte maturation factor MOS (cmos), the melanocortin 1 receptor (MC1R), the neurotrophin–3 (NT3), and the recombination activation gene 1 (RAG1). We used the same primers and PCR conditions as described in detail in Šmíd et al. (2015, 2021), Tamar et al. (2016), and Jablonski et al. (2019). The PCR products were then visualised on electrophoresis in 2% agarose gel, subsequently purified using EXOSAP–IT® PCR Product Cleanup Reagent (Thermo Fisher Scientific, USA) and bidirectionally Sanger–sequenced in Macrogen (the Netherlands) and at the Centre for Molecular Analyses (BIOPOLIS/CIBIO).

Apart from generating sequences for the newly obtained material, we generated an additional 27 new sequences for 11 samples from the studies by Šmíd et al. (2015) and Tamar et al. (2016). We retrieved all available sequences for Rhynchocalamus species from GenBank originating from the two above studies and Avcı et al. (2015), Fathinia et al. (2017), and Tamar et al. (2020). Lytorhynchus diadema was used as an outgroup. Geneious R11 (Kearse et al. 2012) was used to inspect raw sequence files, assemble contigs, generate consensus sequences, and concatenate alignments. Heterozygous positions in the nuclear markers were identified by the Heterozygous Plugin, checked by eye, and coded according to the IUPAC ambiguity codes.

Sequences of each genetic marker were aligned separately using MAFFT online service (Katoh et al. 2019). The Q–INS–i strategy that considers the secondary structure of the RNA was applied for the 12S and 16S gene fragments, while the sequences of the remaining genetic markers were aligned using default settings. Finally, the sequences of protein-coding genes were translated into amino acids using appropriate genetic codes and no stop codons were detected, indicating that no pseudogenes were amplified.

Phylogenetic analyses and nuclear networks reconstruction

We conducted Maximum Likelihood (ML) and Bayesian inference (BI) phylogenetic analyses applying a partition scheme by gene. The ML analysis was carried out in IQ–TREE (Nguyen et al. 2015) using its online web interface W–IQ–TREE (Trifinopoulos et al. 2016). The best substitution model for each gene was selected automatically by ModelFinder (Kalyaanamoorthy et al. 2017) as implemented in IQ–TREE. Branch support was assessed by the Shimodaira–Hasegawa–like approximate likelihood ratio test (SH–aLRT; Guindon et al. 2010) and the Ultrafast bootstrap approximation algorithm (UFBoot; Minh et al. 2013), both with 1000 replicates.

We also carried out a Bayesian inference using MrBayes v.3.2.1 (Ronquist et al. 2012) with the same partitioning strategy as for the ML analysis. The best-fit substitution models were selected using PartitionFinder v.2.1.1 (Lanfear et al. 2017) with the following parameters: branch lengths linked; MrBayes models; AICc model selection; and each gene representing a separate partition specified as a user scheme search. The BI was set to run through the CIPRES Science Gateway (Miller et al. 2010) in three independent runs, each for 10 million generations with parameters and trees sampled every 10,000 generations. All parameters (statefreq, revmat, shape, and pinvar) were unlinked for the partitions. After inspecting that the runs had converged, we discarded as burn-in 25% of posterior trees from each run. A 50% majority-rule consensus tree was then produced from all post-burnin posterior trees. Stationarity was determined by examining the standard deviation of the split frequencies between the three runs (being lower than 0.01), the Potential Scale Reduction Factor diagnostic (PSRF approaching 1.0), and by confirming that all parameters had reached stationarity and had sufficient effective sample sizes (ESS >200) using Tracer v.1.7.2 (Rambaut et al. 2018).

To inspect genealogical relationships and the level of nuclear allele sharing among the Rhynchocalamus species, we reconstructed haplotype networks for the four nuclear loci. To resolve the heterozygous single nucleotide polymorphisms, the alignments of the four nuclear loci were phased separately using the PHASE algorithm (Stephens et al. 2001) as implemented in DnaSP v.6 (Rozas et al. 2017) with the probability threshold set to 0.7. The outgroup was excluded from the alignment for this analysis. Haplotype networks were constructed from the phased alignments using the TCS algorithm (Templeton et al. 1992; Clement et al. 2000) implemented in PopART (Leigh and Bryant 2015). Inter- and intraspecific uncorrected p-distances with pairwise deletion were calculated for 12S, 16S, and cytb in MEGA v11 (Tamura et al. 2021).

Species distribution models

The extent of the study area ranged from the northwesternmost tip of the Sinai Peninsula to the northernmost tip of Israel, Jordan, and Kuwait, including the whole of the Arabian Peninsula. To describe the climatic conditions of the study area, we obtained 19 bioclimatic variables from Worldclim2 (Fick and Hijmans 2017), with a ~1 km2 resolution. Furthermore, we included the soil type (FAO; http://www.fao.org/), and the following topographic variables: altitude (meters), slope (degrees), Topographic Position Index (TPI; i.e., the difference between the elevation of a cell and the eight neighbouring cells), and degree of exposure to the east and the north (hereafter referred to as eastness and northness), with values ranging from –1 to 1 (i.e, complete exposure). Eastness and northness were obtained by computing cosine and sine transformations of the angular direction of the raster cell to the geographical east and north, respectively. We used the function “get_elev_raster” (package elevatr; Hollister et al. 2023) to extract the Digital Elevation Model (approx. 150 m resolution) of the study area, and the function “terrain” (package raster; Hijmans et al. 2021) to compute the other topographical variables. We assessed multicollinearity among variables by performing a Variance Inflation Factor analysis using the function “vif” (package usdm; Naimi et al. 2014) and retained the variables with vif values < 10 (Curto and Pinto 2011; Dormann et al. 2013). As a result, the bioclimatic variables included in the models were the mean diurnal range of temperatures (BIO2), isothermality (BIO3), mean temperature of the wettest and driest quarter (BIO8 and BIO9, respectively), precipitation of the wettest and driest month (BIO13 and BIO14, respectively), precipitation seasonality (coefficient of variation; BIO15), and precipitation of warmest and coldest quarter (BIO18 and BIO19, respectively). All variables were scaled to ~1 km2 resolution.

We built and evaluated our Species Distribution Models (SDMs) using the sdm package (Naimi and Araújo 2016) in R environment (R vers. 4.3.2; R Core Team 2024), using three different modelling techniques: Generalized Linear Models (GLMs), Random Forests (RFs), and Maximum Entropy Models (Maxent), as they are considered the best-performing ones (Kaky et al. 2020; Ancillotto and Labadessa 2023). To run the presence–background SDMs, we randomly generated 10,000 background points selected within the study area. For each technique, we performed 10 replicates. We performed an ensemble forecasting approach, weighting models based on their AUC scores, which allows reducing the uncertainty of predictions by single model algorithms (Watling et al. 2015). We used 30% of randomly selected occurrence data for model performance testing and 70% for model training. Model performance was assessed by inspecting the values of the area under the receiver operating characteristic curve (AUC), True Skill Statistics (TSS) (Araújo and New 2007), and the continuous Boyce Index (Boyce et al. 2002; Hirzel et al. 2006), which was computed using the “ecospat.boyce” function of the R package ecospat (Broennimann et al. 2023). We assessed the effect of each environmental predictor on the probability of occurrence of the new Rhynchocalamus species by inspecting the response curves and estimated each variable’s relative importance by using the “getVarImp” function in the sdm package (Naimi and Araújo 2016). Lastly, we binarised the model predictions by using as thresholds the lowest predicted suitability value among the occurrence points of the species.

Evaluation of conservation status

We used GeoCAT (Bachman et al. 2011; at https://geocat.kew.org) and all the distributional records of the new species (see Suppl. material 1) to compute the Extent of Occurrence (EOO) and the Area of Occupancy (AOO; 2–km grid cell width), two metrics related to the geographic range used in the process of evaluating the threat category of a species by the International Union for Conservation of Nature (IUCN). Lastly, we propose an evaluation of the conservation status of the new taxon following the IUCN Red List guidelines (IUCN SSC Standards and Petitions Committee 2022).

Results

Phylogenetic analyses

The final concatenated alignment of the three mtDNA and four nDNA genetic markers included 48 samples and 172 sequences. The total length was 4,823 base pairs (bp; 625 bp of 12S, 512 bp of 16S, 1092 bp of cytb, 408 bp of cmos, 665 bp of MC1R, 486 bp of NT3, and 1035 bp of RAG1).

Both ML and BI phylogenetic analyses resulted in almost identical topologies (Fig. 1). The clade of all six Rhynchocalamus species was strongly supported (SH-aLRT = 100/UFBoot = 100/pp = 1, support values are given in this order hereafter), as well as all the species within the genus. Rhynchocalamus melanocephalus was recovered to be sister to the remaining species of the genus but with poor support (46/66/0.72). Rhynchocalamus levitoni and R. satunini formed a clade (100/100/1). Rhynchocalamus arabicus was recovered to be a sister species to R. dayanae and the new species described herein (89.8/93/1), which were supported as sister species (99.9/100/1). The deeper relationships of the Rhynchocalamus phylogeny, however, remained unresolved (Figs 1, 2).

The reconstructed haplotype networks (Fig. 2) showed a certain degree of allele sharing in three out of four nuclear loci. The new species of Rhynchocalamus had private alleles (i.e., not shared with other species) only in RAG1, which is also the only gene free of allele sharing among all the Rhynchocalamus species. The new species shared alleles with its sister species R. dayanae in MC1R and NT3. Additionally, the new species shared alleles with R. arabicus, R. melanocephalus, and R. satunini in cmos, indicating signs of incomplete lineage sorting rather than ongoing gene flow among, to our current knowledge, largely allopatric species.

Inter- and intraspecific uncorrected p-distances for all three mtDNA genes are summarised in Table 2.

Figure 1. 

Maximum Likelihood phylogenetic tree reconstructed from the concatenated dataset of 12S, 16S, cytb, cmos, MC1R, NT3, and RAG1 genes (4,823 bp). The tree was rooted using Lytorhynchus diadema (not shown in the figure). Support values are indicated by the black squares (SH-aLRT/UFBoot/pp ≥ 80/95/0.95) or by exact values near nodes. Samples for which new genetic data were generated are highlighted in bold. Complete trees with original ML and BI support values are provided as Suppl. material 1: figs S1, S2, respectively. The specimen depicted is the individual NMP 76815 (sample code LP760) from the PMBSRR, Tabuk / Medina Province, Saudi Arabia.

Figure 2. 

Haplotype networks of the four nDNA genes (cmos, MC1R, NT3, and RAG1). Circle sizes are proportional to the number of samples that share that allele. Short transverse bars on the connecting lines indicate the number of mutational steps between alleles.

Table 1.

Measurements (in mm) of the type series of Rhynchocalamus hejazicus sp. nov. For abbreviations see the Methods section.

Voucher code RCU-URN-93850 MNHN-RA-2023.0013 NMP 76815 RCU-URN-94065 RCU-URN-94064
Sample Code FLI447 JIR544 LP760 JCB222 FLI330
Type status Holotype Paratype Paratype Paratype Paratype
Sex Male Female Male Male Female
SVL 321.11 339.51 214 166.18 226.48
TL 63.60 64.24 45.50 38.32
PreO 1/1 1/1 1/1 1/1 1/1
PostO 2/2 1/1 1/1 2/2 1/1
TS 1/1 1/1 1/1 1/1 1/1
PTS 1/0 1/1 2/1 1/1 1/1
LS 1/1 1/1 1/1 1/1 1/1
VS 247 250 228 227 240
SCS 67 70 68 70
UL 6/6 6/6 6/6 6/6 6/6
LL 8/8 8/8 8/8 8/8 8/8
InfLC 4 4 4 4 3
IntN Triangle Trapezoid Trapezoid Trapezoid Triangle
intNsep Yes No No No No
Table 2.

Mean genetic distances (uncorrected p-distances) between the Rhynchocalamus species based on the 12S and 16S (below the diagonal), and cytb (above the diagonal). Intraspecific distances are shown on the diagonal in bold for 12S, 16S, and cytb, respectively.

R. arabicus R. dayanae R. melanocephalus R. satunini R. hejazicus sp. nov. R. levitoni
R. arabicus NA 10% 10.5% 12.9% 9.5% 9.9%
R. dayanae 5.9 / 4.2% 0.1 / 0.3 / 0.5% 10.1% 11.9% 9.5% 9.8%
R. melanocephalus 7.3 / 3.9% 6.7 / 4% 0.6 / 0.5 / 1.6% 12.9% 10.3% 11.6%
R. satunini 8.9 / 4.6% 7.4 / 5.3% 9 / 3.7% 1.3 / 0.4 / 1.5% 12.3% 7.6%
R. hejazicus sp. nov. 6.1 / 3.7% 3.8 / 3.1% 5.7 / 3.8% 5.5 / 5.1% 0.2 / 0.2 / 0.7% 11.1%
R. levitoni NA NA NA NA NA NA / NA / 2.7%

Taxonomic implications

The results of the phylogenetic analyses, as well as the morphological comparison, indicate that the unknown species of Rhynchocalamus from Saudi Arabia represents a new distinct species in need of formal taxonomic recognition provided here:

Rhynchocalamus hejazicus sp. nov.

Type material

Holotype. RCU-URN-93850 (sample code FLI447; Fig. 3; Table 1). Adult male; Kingdom of Saudi Arabia, Medina Province, Shaaran Nature Reserve (26.85895°N, 38.30387°E), 1035 m asl.; collected by Fulvio Licata on May 24, 2023. MorphoBank accessions: M902248–M902255.

Paratypes. RCU-URN-94064 (sample code FLI330; Fig. 3; Table 1), melanistic adult female collected by Fulvio Licata and Leili Khalatbari on May 11, 2023, KSA, Medina Province, Harrat Khaybar (25.60676°N, 39.96734°E), elevation 1616 m asl, MorphoBank accessions: M902243–M902247; RCU-URN-94065 (sample code JCB222), juvenile male collected by José Carlos Brito on November 2, 2023, KSA, Medina Province, Jabal Al Ward (26.36263°N, 37.38780°E), elevation 1277 m asl, MorphoBank accessions: M902256–M902261; NMP 76815 (sample code LP760; Fig. 4; Table 1), juvenile male collected by Damien M. Egan and Lukáš Pola on April 22, 2022, KSA, Tabuk / Medina Provinces, Prince Mohammed bin Salman Royal Reserve, 26 km east of Al Qurr village (26.6028142°N, 37.1225145°E), elevation 781 m asl, MorphoBank accessions: M902293–M902338; NMP 76816 (sample code JIR551), adult, sex undetermined, collected by Ahmed Mohajja AlShammari and Mubarak Al–Aslami in July 2021, KSA, Hail Province, Jabal Salma, south of An Niayy village (27.165278°N, 42.294444°E), elevation 1052 m asl, MorphoBank accessions: M902290–M902292; MNHN–RA–2023.0013 (sample code JIR544; Fig. 4; Table 1) adult female collected by Adel A. Ibrahim and Ahmed Mohajja AlShammari on May 11, 2017, KSA, Hail Province, Jabal Salma, Wadi Al–Azraq, east of Zekheen village (26.948889°N, 41.973889°E), elevation 1099 m asl, MorphoBank accessions: M902262–M902289.

Other material

Additional unvouchered individuals, assigned to this species based on phenotypical resemblance and unique head pattern, were encountered in May, June, and November 2023 by Vidak Lakušić, Gholam Hosein Yusefi, and Fulvio Licata in Harrat Uwayrid, Al-Gharameel, Wadi Nakhlah, and Sharaan (conservation areas located in the AlUla County, Medina Province). Other individuals were observed in July and August 2023 by Neil Rowntree and Euan Ferguson in NEOM, the northwestern region of Tabuk Province, KSA. One record was found in the public database iNaturalist (https://www.inaturalist.org/observations/20014261). Based on the provided information, the snake was recorded on February 1, 2019 at As Salam, Medina Province (24.464058°N, 39.537669°E). An additional individual was reported to us by Muteb Masad Al–Malki from the vicinity of Adham city (20.4486525°N, 40.8792636°E) from Mecca Province.

Etymology

The species name is a latinized noun in masculine gender derived from the word "Hejaz–" = Hejaz Mountains, a mountain range located in the Hejaz region (an important region located in western Saudi Arabia, where the two holy cities of Islam, Mecca, and Medina are located) where most individuals were observed, and the Latin suffix "–icus" = “belonging to”. We suggest the common name “Hejaz black–collared snake” in English and أبو حناء [Abu Henna] in Arabic for the new species.

Diagnosis

The new species of Rhynchocalamus from the Hejaz Mountain range in western Saudi Arabia is characterised by the following morphological characters: (1) SVL 209.2–339.5 mm in adults; (2) tail length 38.3–64.2 mm in adults; (3) loreal scale present; (4) large 3rd and 4th upper labial scales in contact with the eye; (5) one preocular scale; (6) 1–2 postocular scales; (7) one temporal scale; (8) 0–2 post–temporal scales; (9) six upper labial scales; (10) eight lower labial scales; (11) usually four lower labial scales in contact with the anterior inframaxillars; (12) usually one gular scale in contact with anterior inframaxillars, situated between the posterior inframaxillars; (13) 15 smooth dorsal scales at mid-body; (14) 11–12 dorsal and temporal scales surrounding the margin of parietals; (15) 227–250 ventrals; (16) anal and subcaudal scales divided; (17) 67–70 subcaudal scales; (18) dorsal colouration in life deep reddish with a distinctive black collar extending behind the parietal scales and abruptly stopping or tapering backward in the middle, and a pale reddish band passing behind the eyes, through the middle of the supraoculars and the frontal scale, encompassing the temporal and parietal scales.

Colouration in life (adults)

The upper surface of the head is shiny black from the middle of the supraoculars and the frontal scale to the tip of the snout, which is whitish; a wide band, pale reddish dorsally and fading to whitish ventrally, passes behind the eyes, through the middle of the supraoculars and the frontal scale, encompassing the temporal and parietal scales; a black collar around the neck reaches the ventrals, and abruptly stops or tapers dorsally towards the centre; dorsal surface of the body and tail uniformly deep reddish from the end of the collar to the tail tip; ventral surface of the body deep reddish fading to whitish in the upper part of the body, narrowing in correspondence of the black collar; colour pattern paler in alcohol preserved specimen. A melanistic morph, uniformly black, also occurs.

Comparison

Rhynchocalamus hejazicus sp. nov. is morphologically similar to the other Rhynchocalamus species, and it can be distinguished by slight differences in size, colouration, and head and body scalation.

In comparison with R. arabicus, R. hejazicus sp. nov. has a lower number of subcaudal scales (67–70 vs. 71–81 in R. arabicus). We show that a melanistic morph of R. hejazicus sp. nov. also occurs, therefore the new species can be easily misidentified as R. arabicus.

Rhynchocalamus hejazicus sp. nov. differs from R. dayanae by smaller maximum body size (339.5 vs. 432.1 in R. dayanae), a shorter tail (38.3–64.2 vs. 59.2–94.1 in R. dayanae), by a higher number of ventrals (227–250 vs. 198–229 in R. dayanae) and subcaudals (67–70 vs. 54–62 in R. dayanae), Lastly, R. hejazicus sp. nov. can be further differentiated from R. dayanae by the presence of a pale reddish band passing between the eyes and the neck.

Rhynchocalamus hejazicus sp. nov. differs from R. melanocephalus (both Southern population from the Negev region in Israel and Northern population from the Mediterranean ecoregion; sensu Tamar et al. 2020) by a smaller maximum body size (339.5 vs. 499.2 in R. melanocephalus), a lower number of post-temporal scales (1 vs 2 in R. melanocephalus), a higher number of lower labial scales (8 vs. 7 in R. melanocephalus), a higher number of ventrals (227–250 vs. 164–235 in R. melanocephalus), and subcaudals (67–70 vs. 29–69 in R. melanocephalus). Lastly, R. hejazicus sp. nov. can be further differentiated from R. melanocephalus by a pale reddish band between the eyes and the neck.

Rhynchocalamus hejazicus sp. nov. differs from R. satunini in having a longer tail (64.2 vs. 54 mm in R. satunini), a higher number of ventrals (227–250 vs. 201–226 in R. satunini) and subcaudals (67–70 vs. 53–64 in R. satunini), and a lower number of upper labials (6 vs. 7 in R. satunini). Lastly, R. satunini is characterised by two black patches on a pale reddish/whitish background on the prefrontals and the parietals, and a black band around the neck that does not reach the ventrals.

In comparison with R. levitoni, R. hejazicus sp. nov. has a lower number of upper labials (6 vs. 7 in R. levitoni), and it can be distinguished by the overall different colouration (deep reddish vs. lemon yellow in R. levitoni), and the absence of a V-shaped band on the neck, and dark markings on the parietals.

Description of the holotype

Adult male (voucher code RCU-URN-93850) (Fig. 3). Snout-vent length 321.1 mm, tail length 63.6 mm. Body cylindrical and slender. The head is small and relatively narrow, slightly distinct from the neck. Head scales not keeled. Rostral scale elongated, extending backwards and wedged between the internasals, which are separated, and touching the bottom of the suture between the prefrontals. The rostral is bordered by two upper labials, two nasals, two internasals, and two prefrontals. Nostrils are situated on undivided nasal scales internasals of triangular shape, separated by the rostral scale. A trapezoid loreal scale at either side in contact with the 2nd and 3rd upper labials. The eyes are small, with circular pupils. One preocular scale on both sides. Two postocular scales on both sides. Bell-shaped frontal scale, located between the supraoculars and wedged between the parietals. Six squarish upper labials, the 3rd and 4th in contact with the eye. two large parietal scales; one temporal and one post–temporal scale on the left side, while on the right side, only one enlarged temporal is present. Eight lower labials. Four lower labials are in contact with anterior inframaxillars on each side. Two gular scales are positioned between the posterior inframaxillars and in connection with the anterior inframaxillars. Anterior inframaxillars are almost double in size than posterior inframaxillars. Eleven dorsal and temporal scales surround the posterior margin of the parietals. Fifteen dorsal scale rows at mid-body. 247 ventrals and 67 subcaudals, including a conical scale at the tail tip.

Known and potential distribution and habitat

The known extent of occurrence (EOO) of R. hejazicus sp. nov. is more than 274,674 km2 and the area of occurrence (AOO) is 56 km2. The individual reported from the vicinity of Adham city (20.4486525°N, 40.8792636°E) from Mecca Province represents the southernmost record to our knowledge.

Available material and photographic records suggest that the species is scattered across Tabuk, Medina, Hail and Mecca Provinces, and is likely endemic to KSA. Nevertheless, it should be noted that additional range extensions of many overlooked species are now being reported with each herpetofaunal survey carried out in the unexplored regions of northwestern and western KSA (e.g., Aloufi et al. 2020, 2021, 2022; Pola et al. 2024), suggesting that the range of this rather secretive species could be even wider. Currently, the northwesternmost individual is reported from the area west of Jabal Al Lawz, approx. 70 km south of the Jordanian border. Given the presence of similar habitats harbouring many shared species, we cannot rule out the species occurrence in that area, especially in the mountains of Aqaba.

Rhynchocalamus hejazicus sp. nov. has been observed between 456 and 1610 m a.s.l., in the following habitats: (i) sandy flatland with sparse vegetation (Acacia sp., bushes, and tussock grasses; Fig. 5); (ii) large, stony wadis in lava fields (=harrat) with dense patches of woody vegetation (Acacia sp.; Fig. 5); (iii) sparsely vegetated rocky creek with temporary pools.

The species distribution models achieved good performance levels (AUC: 0.9–0.94; TSS: 0.78–0.87; Boyce Index: 0.42–0.99), resulting in an overall good predictive accuracy of the final ensemble models (Boyce index: 0.71). The potential distribution of R. hejazicus sp. nov. mostly includes the known range of the species, but other extralimital suitable areas are identified to the northwest (i.e., the Sinai Mountains in South Sinai, Aqaba mountains, and Wadi Rum in Jordan), and to the north (i.e., Jabal Tubaiq hills at the border with Jordan and Upper Galilee mountains in North Israel; Fig. 6). The main drivers of environmental suitability in our models were the altitude (26.7% of explained variance) and the mean temperature of the wettest quarter (BIO8; 15.6%), suggesting that R. hejazicus sp. nov. is associated with altitudes above 1,000 m a.s.l., in regions characterised by cold wet months (<20 °C). Furthermore, the precipitation of the warmest quarter (BIO18; 13.5%) and the precipitation seasonality (BIO15; 13%) were also related to the probability of occurrence of the species, indicating that the species is likely to occur in areas with relatively dry hot months (precipitations < 50 mm) and stable rainfall regimes.

Natural history

Little is known about the species’ natural history and behaviour. Although the current number of observations is limited, it appears that R. hejazicus sp. nov. has mainly nocturnal activity as all individuals were encountered active at night. We assume that in many biological aspects, it will be concordant with its congeners (Amr and Disi 2011; Tamar et al. 2016; Torki 2017; Fathinia et al. 2017).

Other reptile species observed in syntopy of the NMP 76815 (LP760) were Ptyodactylus cf. hasselquistii (Phyllodactylidae), Pristurus guweirensis (Sphaerodactylidae), and Echis coloratus (Viperidae). In Hail Province, at Jabal Salma range, Echis coloratus, Cerastes gasperettii (Viperidae) and Walterinnesia aegyptia (Elapidae) were found in sympatry (Alshammari and Ibrahim 2015; Alshammari et al. 2017). In AlUla County, it was found in sympatry with Ptyodactylus guttatus (in prep.), Ptyodactylus cf. hasselquistii, Bunopus tuberculatus, Hemidactylus granosus, Stenodactylus doriae, Stenodactylus slevini, Tropiocolotes yomtovi (the latter five all Gekkonidae), and Echis coloratus.

Conservation status

We propose to classify R. hejazicus sp. nov. as Least Concern, based on the wide EOO, and the diverse and continuous habitats where the species is encountered, including habitats with heavy anthropogenic disturbance (e.g., overgrazing) and urban habitats (see https://www.inaturalist.org/observations/20014261). Although R. hejazicus sp. nov. is at the moment known from 14 localities and its AOO is only 56 km2, we consider it premature to classify it as Near Threatened as the limited number of localities could reflect the generally undersampled distribution and secretive habits of the species. Furthermore, six of the known localities occur within conservation areas located in AlUla County (Medina Province), where protection measures have been implemented.

Figure 3. 

Holotype (RCU-URN-93850, sample code FLI447, bottom left) and paratype (RCU-URN-94064, sample code FLI330, bottom right) specimens in life. Lateral and ventral views of their heads (above each life picture). Photo credit: FL.

Figure 4. 

Colour variation within R. hejazicus sp. nov. (or its lack, thereof). Top row: two unvouchered specimens from NEOM, Tabuk Province, KSA (photo credit: Euan Ferguson and Neil Rowntree); bottom left: paratype NMP 76815 (sample code LP760, photo credit: DME); bottom right: paratype MNHN–RA–2023.0013 (sample code JIR544, photo credit: AAI).

Figure 5. 

Habitats of the holotype and two paratype specimens of R. hejazicus sp. nov. Top: locality of the holotype RCU-URN-93850 at the Shaaran NR, AlUla County, Medina Province, KSA (photo credit: FL); bottom left: habitat of the paratype MNHN–RA–2023.0013 at Wadi Al-Azraq, Jabal Salma, Hail Province, KSA (photo credit: AAI); bottom right: habitat of the paratype RCU-URN-94065 the Harrat Khaybar, Hail Province, KSA (photo credit: FL).

Figure 6. 

Distribution of Rhynchocalamus hejazicus sp. nov. showing the location of the material examined in this study. Circles indicate samples used for the genetic analyses, squares indicate examined voucher specimens, and diamonds indicate photographic records. The black arrow denotes the type locality (Shaaran NR in AlUla County, Medina Province). Geographical sampling is overlaid by the predicted species distribution model categorised into binary predictions of the suitable areas. The threshold used for the categorisation was the lowest predicted suitability value among the occurrence points (= 0.05).

Discussion

We describe a new species of the genus Rhynchocalamus based on phylogenetic and morphological analyses. The meristic and morphometric traits of R. hejazicus sp. nov. are largely similar to the other species of Rhynchocalamus, from which, however, the new species can be distinguished for its unique colouration. The phylogenetic analyses show that the new Rhynchocalamus species is sister to R. dayanae from south Israel and closely related to R. arabicus from Yemen and Oman. The presence of this species in western Saudi Arabia fills a large biogeographic gap between northerly and southerly distributed Rhynchocalamus sister taxa, highlighting the potential relevance of mountainous areas as biogeographic corridors for the spread of this genus in the Arabian Peninsula. Previous phylogenetic studies were unable to resolve deep relationships within the genus due to low support of the deepest nodes (Šmíd et al. 2015; Tamar et al. 2016; Fathinia et al. 2017). Despite the additional three nuclear genes included here in the phylogenetic dataset, the deep relationships within the genus remained unresolved in our analyses, preventing clear conclusions on the relationships between the clades within the genus as well as its biogeographic history. This may be caused by a high degree of allele sharing in most of the nuclear markers analysed, and only genome-wide data (e.g., SNPs) will likely be able to fully resolve the phylogeny of Rhynchocalamus.

The distribution modelling identified potentially suitable extralimital areas north and northwest of the known distribution range of the new Rhynchocalamus species. It is unlikely that the species occurs in the well-sampled Upper Galilee mountains in north Israel, where the congeneric R. melanocephalus instead occurs (Tamar et al. 2016; Werner 2016). The distribution models may overpredict species’ distributions, especially if there is a lack of a priori knowledge of species interactions and dispersal abilities, which are difficult to take into account in the modelling procedure (Briscoe et al. 2019; Velazco et al. 2020). Conversely, our models might also have underpredicted suitable areas for the species, as they did not include fine-scale habitat and soil variables which could represent important determinants of distribution in fossorial snakes (e.g., Wagner et al. 2014).

We report a melanistic morph of the new species (Fig. 3). Melanism may have a significant adaptive role in ectotherms, having, among others, thermoregulatory and antipredatory functions (Trullas et al. 2007). In desert animals, melanism is found to have aposematic effects, except for habitats with dark substrates such as black sands or lava fields, where it has cryptic functions (Cloudsley-Thompson 1979). The melanistic individual of R. hejazicus sp. nov. was found in Harrat Khaybar, a lava field located north of Medina (Fig. 5), where melanism might have been positively selected for antipredatory, cryptic functions. This represents an exceptional finding, but it might not be uncommon, taking into account the limited number of observations of this species and the elusiveness of these snakes. Melanism could have a deeper role in the evolution of the genus, as suggested by the southerly distributed, melanistic R. arabicus.

Although a recent boost in herpetological surveys is rapidly contributing to filling up knowledge gaps in reptile diversity in KSA (e.g., Aloufi et al. 2020, 2021, 2022), still vast areas remain unexplored and secretive, potentially widespread species can be expected to occur. Rhynchocalamus hejazicus sp. nov. is the sixth species of Rhynchocalamus snakes and yet is not the lesser-known, having more occurrence records and a larger distribution range than most congeneric species.

The description of R. hejazicus sp. nov. benefited from previous taxonomic and biogeographical insights into the genus (Tamar et al. 2016). However, the detection of secretive species in underexplored areas, such as the Hejaz region, primarily stems from rigorous sampling efforts. This emphasises the significance of developing efficient monitoring strategies aimed at comprehensively documenting species diversity in uncharted territories. Such efforts are crucial for advancing our knowledge of ecosystems and community functions and thus bolstering global conservation efforts.

Author contributions

FL, LP, JŠ, JCB, and DGV conceived the ideas and designed the methodology; FL, LP, DME, and JCB collected the data; AMA, SB, AAI, DGV, AVL, BS, PL, and JCB contributed to the field sampling; FL, LP, AVL, and BS analysed the data; FL and LP led the writing of the manuscript. JŠ, JCB, PL, and AA contributed critically to the draft. All authors have read critically the draft and gave final approval for publication.

Funding

LP was supported by Charles University grant no. SVV260685/2023. JŠ was supported by the Czech Science Foundation (GACR, project number 22-12757S) and by the Ministry of Culture of the Czech Republic (DKRVO 2024–2028/6.I.a, 00023272). JCB and DVG were supported by FCT (CEECINST/00014/2018/CP1512/CT0001; 2020.03848.CEECIND). JCB, DVG, FL, BS, AVL were supported by the project Inventory of AlUla Fauna (PR6869) from RCU – The Royal Commission for AlUla.

Acknowledgements

We would like to thank Nicolas Vidal (MNHN, Paris, France) for providing the photos of the specimen under his care. LP and DME would like to thank Talal Assaf Al–Atawi for his company, and technical support during the fieldwork. The fieldwork of LP and DME was part of a project focused on biological inventory of the PMBSRR developed and funded by Prince Mohammed bin Salman Royal Reserve Development Authority. AMA would like to thank Mubark Alshammari for the company and technical support during the fieldwork. We would like to express our thanks to Neil Rowntree and Euan Ferguson (United Kingdom) for sharing their records from NEOM, the new region in northwestern Saudi Arabia, with us along with the photographs and the accompanying information, as well as to Muteb Masad Al–Malki (Saudi Arabia) for sharing the record from Mecca Province. We would like to thank Leili Khalatbari, Gholam Hosein Yusefi, and Vidak Lakušić for contributing to the finding of the animals and for sharing their records from the AlUla region, and Nina Séren for contributing to the genetic analyses. Acknowledgements extended to Alaaeldin I Soultan, Ingrid Stirnemann, and the Royal Commission for AlUla for their support in developing fieldwork in the AlUla region.

References

  • Alatawi AS, Gilbert F, Reader T (2020) Modelling terrestrial reptile species richness, distributions and habitat suitability in Saudi Arabia. Journal of Arid Environments 178: 104153. https://doi.org/10.1016/j.jaridenv.2020.104153
  • Aloufi A, Baker MAA, Amr ZS (2020) First record of the Variable Racer, Platyceps variabilis (Boulenger, 1905), for Saudi Arabia. Herpetology Notes 13: 973–975.
  • Aloufi A, Amr Z, Baker MA (2022) Reptiles from ’Uruq Bani Ma’arid and Harat al Harrah protected areas in Saudi Arabia: Reptiles from two protected areas in Saudi Arabia. Herpetology Notes 15: 483–491.
  • Alshammari AM, Ibrahim AA (2015) Lizards and snakes in the historical Faid protected area (Faid Hema), Ha’il region, Saudi Arabia. Herpetological Conservation and Biology 10: 1021–1029.
  • Alshammari AM, Busais SM, Ibrahim AA (2017) Snakes in the Province of Ha’il, Kingdom of Saudi Arabia, including two new records. Herpetozoa (Wien) 30: 59–63.
  • Ancillotto L, Labadessa R (2023) Can protected areas and habitats preserve the vulnerable predatory bush cricket Saga pedo? Journal of Insect Conservation 27(4): 615–624. https://doi.org/10.1007/s10841-023-00484-w
  • Avcı A, Ilgaz Ç, Mahdi Rajabizadeh, Yılmaz C, Üzüm N, Adriaens D, Kumlutaş Y, Olgun K (2015) Molecular phylogeny and micro ct-scanning revealed extreme cryptic biodiversity in Kukri snake, Muhtarophis gen. nov., a new genus for Rhynchocalamus barani (Serpentes: Colubridae). Russian Journal of Herpetology 22: 159–174.
  • Bachman S, Moat J, Hill AW, de la Torre J, Scott B (2011) Supporting Red List threat assessments with GeoCAT: Geospatial conservation assessment tool. ZooKeys 150: 117–126. https://doi.org/10.3897/zookeys.150.2109
  • Briscoe NJ, Elith J, Salguero-Gómez R, Lahoz-Monfort JJ, Camac JS, Giljohann KM, Holden MH, Hradsky BA, Kearney MR, McMahon SM, Phillips BL, Regan TJ, Rhodes JR, Vesk PA, Wintle BA, Yen JDL, Guillera-Arroita G (2019) Forecasting species range dynamics with process-explicit models: Matching methods to applications. Ecology Letters 22(11): 1940–1956. https://doi.org/10.1111/ele.13348
  • Broennimann O, Cola VD, Petitpierre B, Breiner F, Scherrer D, D’Amen M, Randin C, Engler R, Hordijk W, Mod H, Pottier J, Febbraro MD, Pellissier L, Pio D, Mateo RG, Dubuis A, Maiorano L, Psomas A, Ndiribe C, Salamin N, Zimmermann N, Collart F, Smith T, Guisan A (2023) ecospat: Spatial Ecology Miscellaneous Methods. https://cran.r-project.org/web/packages/ecospat/index.html [April 2, 2024]
  • Dormann CF, Elith J, Bacher S, Buchmann C, Carl G, Carré G, Marquéz JRG, Gruber B, Lafourcade B, Leitão PJ, Münkemüller T, McClean C, Osborne PE, Reineking B, Schröder B, Skidmore AK, Zurell D, Lautenbach S (2013) Collinearity: A review of methods to deal with it and a simulation study evaluating their performance. Ecography 36(1): 27–46. https://doi.org/10.1111/j.1600-0587.2012.07348.x
  • Engel MS, Ceríaco LMP, Daniel GM, Dellapé PM, Löbl I, Marinov M, Reis RE, Young MT, Dubois A, Agarwal I, Lehmann AP, Alvarado M, Alvarez N, Andreone F, Araujo-Vieira K, Ascher JS, Baêta D, Baldo D, Bandeira SA, Barden P, Barrasso DA, Bendifallah L, Bockmann FA, Böhme W, Borkent A, Brandão CRF, Busack SD, Bybee SM, Channing A, Chatzimanolis S, Christenhusz MJM, Crisci JV, D’elía G, Da Costa LM, Davis SR, De Lucena CAS, Deuve T, Fernandes Elizalde S, Faivovich J, Farooq H, Ferguson AW, Gippoliti S, Gonçalves FMP, Gonzalez VH, Greenbaum E, Hinojosa-Díaz IA, Ineich I, Jiang J, Kahono S, Kury AB, Lucinda PHF, Lynch JD, Malécot V, Marques MP, Marris JWM, Mckellar RC, Mendes LF, Nihei SS, Nishikawa K, Ohler A, Orrico VGD, Ota H, Paiva J, Parrinha D, Pauwels OSG, Pereyra MO, Pestana LB, Pinheiro PDP, Prendini L, Prokop J, Rasmussen C, Rödel M-O, Rodrigues MT, Rodríguez SM, Salatnaya H, Sampaio Í, Sánchez-García A, Shebl MA, Santos BS, Solórzano-Kraemer MM, Sousa ACA, Stoev P, Teta P, Trape J-F, Dos Santos CV-D, Vasudevan K, Vink CJ, Vogel G, Wagner P, Wappler T, Ware JL, Wedmann S, Zacharie CK (2021) The taxonomic impediment: A shortage of taxonomists, not the lack of technical approaches. Zoological Journal of the Linnean Society 193(2): 381–387. https://doi.org/10.1093/zoolinnean/zlab072
  • Fathinia B, Rastegar-Pouyani E, Rastegar-Pouyani N, Darvishnia H (2017) A new species of the genus Rhynchocalamus Günther, 1864 (Reptilia: Squamata: Colubridae) from Ilam province in western Iran. Zootaxa 4282(3): 473–486. https://doi.org/10.11646/zootaxa.4282.3.3
  • Fick SE, Hijmans RJ (2017) WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37(12): 4302–4315. https://doi.org/10.1002/joc.5086
  • Guindon S, Dufayard J-F, Lefort V, Anisimova M, Hordijk W, Gascuel O (2010) New Algorithms and Methods to Estimate Maximum-Likelihood Phylogenies: Assessing the Performance of PhyML 3.0. Systematic Biology 59(3): 307–321. https://doi.org/10.1093/sysbio/syq010
  • Hijmans RJ, van Etten J, Sumner M, Cheng J, Baston D, Bevan A, Bivand R, Busetto L, Canty M, Fasoli B, Forrest D, Ghosh A, Golicher D, Gray J, Greenberg JA, Hiemstra P, Hingee K, Ilich A Geosciences I for MA, Karney C, Mattiuzzi M, Mosher S, Naimi B, Nowosad J, Pebesma E, Lamigueiro OP, Racine EB, Rowlingson B, Shortridge A, Venables B, Wueest R (2021) raster: Geographic data analysis and modeling. https://CRAN.R-project.org/package=raster [January 12, 2022]
  • Jablonski D, Kukushkin OV, Avcı A, Bunyatova S, Kumlutaş Y, Ilgaz Ç, Polyakova E, Shiryaev K, Tuniyev B, Jandzik D (2019) The biogeography of Elaphe sauromates (Pallas, 1814), with a description of a new rat snake species. PeerJ 7: e6944. https://doi.org/10.7717/peerj.6944
  • Kaky E, Nolan V, Alatawi A, Gilbert F (2020) A comparison between Ensemble and MaxEnt species distribution modelling approaches for conservation: A case study with Egyptian medicinal plants. Ecological Informatics 60: 101150. https://doi.org/10.1016/j.ecoinf.2020.101150
  • Kalyaanamoorthy S, Minh BQ, Wong TKF, von Haeseler A, Jermiin LS (2017) ModelFinder: Fast model selection for accurate phylogenetic estimates. Nature Methods 14(6): 587–589. https://doi.org/10.1038/nmeth.4285
  • Katoh K, Rozewicki J, Yamada KD (2019) MAFFT online service: Multiple sequence alignment, interactive sequence choice and visualization. Briefings in Bioinformatics 20(4): 1160–1166. https://doi.org/10.1093/bib/bbx108
  • Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, Buxton S, Cooper A, Markowitz S, Duran C, Thierer T, Ashton B, Meintjes P, Drummond A (2012) Geneious Basic: An integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics (Oxford, England) 28(12): 1647–1649. https://doi.org/10.1093/bioinformatics/bts199
  • Kéry M (2002) Inferring the absence of a species: A case study of snakes. The Journal of Wildlife Management 66(2): 330–338. https://doi.org/10.2307/3803165
  • Lanfear R, Frandsen PB, Wright AM, Senfeld T, Calcott B (2017) PartitionFinder 2: New methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. Molecular Biology and Evolution 34: 772–773. https://doi.org/10.1093/molbev/msw260
  • Laopichienpong N, Tawichasri P, Chanhome L, Phatcharakullawarawat R, Singchat W, Kantachumpoo A, Muangmai N, Suntrarachun S, Matsubara K, Peyachoknagul S, Srikulnath K (2017) A novel method of caenophidian snake sex identification using molecular markers based on two gametologous genes. Ecology and Evolution 7(13): 4661–4669. https://doi.org/10.1002/ece3.3057
  • Leigh JW, Bryant D (2015) popart: full‐feature software for haplotype network construction. Nakagawa S (Ed.). Methods in Ecology and Evolution 6: 1110–1116. https://doi.org/10.1111/2041-210X.12410
  • Matsubara K, Nishida C, Matsuda Y, Kumazawa Y (2016) Sex chromosome evolution in snakes inferred from divergence patterns of two gametologous genes and chromosome distribution of sex chromosome-linked repetitive sequences. Zoological Letters 2(1): 19. https://doi.org/10.1186/s40851-016-0056-1
  • May RM (2004) Tomorrow’s taxonomy: Collecting new species in the field will remain the rate–limiting step. Philosophical Transactions of the Royal Society of London, Series B, Biological Sciences 359(1444): 733–734. https://doi.org/10.1098/rstb.2003.1455
  • Miller MA, Pfeiffer W, Schwartz T (2010) Creating the CIPRES Science Gateway for inference of large phylogenetic trees. In: 2010 Gateway Computing Environments Workshop (GCE), 1–8. https://doi.org/10.1109/GCE.2010.5676129
  • Minh BQ, Nguyen MAT, von Haeseler A (2013) Ultrafast Approximation for Phylogenetic Bootstrap. Molecular Biology and Evolution 30(5): 1188–1195. https://doi.org/10.1093/molbev/mst024
  • Naimi B, Araújo MB (2016) sdm: A reproducible and extensible R platform for species distribution modelling. Ecography 39(4): 368–375. https://doi.org/10.1111/ecog.01881
  • Nguyen L-T, Schmidt HA, von Haeseler A, Minh BQ (2015) IQ-TREE: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Molecular Biology and Evolution 32(1): 268–274. https://doi.org/10.1093/molbev/msu300
  • Pola L, Crochet PA, Geniez P, Shobrak M, Busais S, Jablonski D, Masroor R, Abduraupov T, Carranza S, Šmíd J (2024) Some like it hot: Past and present phylogeography of a desert dwelling gecko across the Arabian Peninsula. Journal of Biogeography 00: 1–15. https://doi.org/10.1111/jbi.14823
  • Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA (2018) Posterior Summarization in Bayesian Phylogenetics Using Tracer 1.7. Systematic Biology 67(5): 901–904. https://doi.org/10.1093/sysbio/syy032
  • Rodda GH (1993) Where’s Waldo (and the snakes)? Herpetological Review 24: 44–45.
  • Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Höhna S, Larget B, Liu L, Suchard MA, Huelsenbeck JP (2012) MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space. Systematic Biology 61(3): 539–542. https://doi.org/10.1093/sysbio/sys029
  • Rozas J, Ferrer-Mata A, Sánchez-DelBarrio JC, Guirao-Rico S, Librado P, Ramos-Onsins SE, Sánchez-Gracia A (2017) DnaSP 6: DNA sequence polymorphism analysis of large data sets. Molecular Biology and Evolution 34(12): 3299–3302. https://doi.org/10.1093/molbev/msx248
  • Scheffers BR, Joppa LN, Pimm SL, Laurance WF (2012) What we know and don’t know about Earth’s missing biodiversity. Trends in Ecology & Evolution 27(9): 501–510. https://doi.org/10.1016/j.tree.2012.05.008
  • Šmíd J, Martínez G, Gebhart J, Aznar J, Gallego J, Göçmen B, de Pous P, Tamar K, Carranza S (2015) Phylogeny of the genus Rhynchocalamus (Reptilia; Colubridae) with a first record from the Sultanate of Oman. Zootaxa 4033(3): 380–392. https://doi.org/10.11646/zootaxa.4033.3.4
  • Šmíd J, Sindaco R, Shobrak M, Busais S, Tamar K, Aghová T, Simó-Riudalbas M, Tarroso P, Geniez P, Crochet P-A, Els J, Burriel-Carranza B, Tejero-Cicuéndez H, Carranza S (2021) Diversity patterns and evolutionary history of Arabian squamates. Journal of Biogeography 48(5): 1183–1199. https://doi.org/10.1111/jbi.14070
  • Stephens M, Smith NJ, Donnelly P (2001) A new statistical method for haplotype reconstruction from population data. American Journal of Human Genetics 68(4): 978–989. https://doi.org/10.1086/319501
  • Tamar K, Šmíd J, Göçmen B, Meiri S, Carranza S (2016) An integrative systematic revision and biogeography of Rhynchocalamus snakes (Reptilia, Colubridae) with a description of a new species from Israel. PeerJ 4: e2769. https://doi.org/10.7717/peerj.2769
  • Tamar K, Wiedl HJ, Maza E, Jablonski D, Meiri S (2020) Discovery of the Black-headed Ground Snake Rhynchocalamus melanocephalus (Jan, 1862) in Cyprus (Reptilia: Colubridae). Zoology in the Middle East 66(2): 118–123. https://doi.org/10.1080/09397140.2020.1757914
  • Templeton AR, Crandall KA, Sing CF (1992) A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics 132(2): 619–633. https://doi.org/10.1093/genetics/132.2.619
  • Trifinopoulos J, Nguyen L-T, von Haeseler A, Minh BQ (2016) W-IQ-TREE: A fast online phylogenetic tool for maximum likelihood analysis. Nucleic Acids Research 44(W1): W232–W235. https://doi.org/10.1093/nar/gkw256
  • Velazco SJE, Ribeiro BR, Laureto LMO, De Marco Júnior P (2020) Overprediction of species distribution models in conservation planning: A still neglected issue with strong effects. Biological Conservation 252: 108822. https://doi.org/10.1016/j.biocon.2020.108822
  • Wagner RO, Pierce JB, Rudolph DC, Schaefer RR, Hightower DA (2014) Modeling Louisiana pine snake (Pituophis ruthveni) habitat use in relation to soils. Southeastern Naturalist (Steuben, ME) 13(sp5): 146–158. https://doi.org/10.1656/058.013.s514
  • Watling JI, Brandt LA, Bucklin DN, Fujisaki I, Mazzotti FJ, Romañach SS, Speroterra C (2015) Performance metrics and variance partitioning reveal sources of uncertainty in species distribution models. Ecological Modelling 309–310: 48–59. https://doi.org/10.1016/j.ecolmodel.2015.03.017
  • Werner YI (2016) Reptile Life in the Land of Israel with Comments on Adjacent Regions. Frankfurt am Main.

1Fulvio Licata and Lukáš Pola contributed equally to this work

Supplementary material

Supplementary material 1 

Supporting information

Fulvio Licata, Lukáš Pola, Jiří Šmíd, Adel A. Ibrahim, André Vicente Liz, Bárbara Santos, László Patkó, Ayman Abdulkareem, Duarte V. Gonçalves, Ahmed Mohajja AlShammari, Salem Busais, Damien M. Egan, Ricardo M. O. Ramalho, Josh Smithson, José Carlos Brito

Data type: docx

Explanation note: table S1. List of material analysed in this study including information on voucher and tissue codes, country, province and locality of origin, GPS coordinates (datum WGS84), and GenBank accession numbers. Accession numbers of sequences generated for this study are highlighted in bold. Rows with neither voucher nor tissue codes are observations used only for the species distribution modelling. figure S1. Phylogenetic tree resulting from the Maximum likelihood analysis of three mitochondrial and four nuclear markers concatenated. Values by branches indicate SH-aLRT/UFBoot (see Materials and methods for abbreviations). figure S2. Phylogenetic tree resulting from the Bayesian analysis of three mitochondrial and four nuclear markers concatenated. Values by branches indicate posterior probabilities.

Download file (42.33 kb)
login to comment