Pre-Quaternary divergence and subsequent radiation explain longitudinal patterns of genetic and morphological variation in the striped skink, Heremites vittatus

Background Many animal and plant species in the Middle East and northern Africa have a predominantly longitudinal distribution, extending from Iran and Turkey along the eastern Mediterranean coast into northern Africa. These species are potentially characterized by longitudinal patterns of biological diversity, but little is known about the underlying biogeographic mechanisms and evolutionary timescales. We examined these questions in the striped skink, Heremites vittatus, one such species with a roughly longitudinal distribution across the Middle East and northern Africa, by analyzing range-wide patterns of mitochondrial DNA (mtDNA) sequence and multi-trait morphological variation. Results The striped skink exhibits a basic longitudinal organization of mtDNA diversity, with three major mitochondrial lineages inhabiting northern Africa, the eastern Mediterranean coast, and Turkey/Iran. Remarkably, these lineages are of pre-Quaternary origin, and are characterized by p-distances of 9–10%. In addition, within each of these lineages a more recent Quaternary genetic diversification was observed, as evidenced by deep subclades and high haplotype diversity especially in the Turkish/Iranian and eastern Mediterranean lineages. Consistent with the genetic variation, our morphological analysis revealed that the majority of morphological traits show significant mean differences between specimens from northern Africa, the eastern Mediterranean coast, and Turkey/Iran, suggesting lineage-specific trait evolution. In addition, a subset of traits exhibits clinal variation along the eastern Mediterranean coast, potentially indicating selection gradients at the geographic transition from northern Africa to Anatolia. The existence of allopatric, morphologically and genetically divergent lineages suggests that Heremites vittatus might represent a complex with several taxa. Conclusions Our work demonstrates that early divergence events in the Pliocene, likely driven by both climatic and geological factors, established the longitudinal patterns and distribution of Heremites vittatus. Subsequent radiation during the Pleistocene generated the genetic and morphological diversity observed today. Our study provides further evidence that longitudinal diversity patterns and species distributions in the Middle East and northern Africa were shaped by complex evolutionary processes, involving the region’s intricate geological history, climatic oscillations, and the presence of the Sahara. Electronic supplementary material The online version of this article (doi:10.1186/s12862-017-0969-0) contains supplementary material, which is available to authorized users.


(Continued from previous page)
Conclusions: Our work demonstrates that early divergence events in the Pliocene, likely driven by both climatic and geological factors, established the longitudinal patterns and distribution of Heremites vittatus. Subsequent radiation during the Pleistocene generated the genetic and morphological diversity observed today. Our study provides further evidence that longitudinal diversity patterns and species distributions in the Middle East and northern Africa were shaped by complex evolutionary processes, involving the region's intricate geological history, climatic oscillations, and the presence of the Sahara.
Keywords: Phylogeography, Western Palearctic, Cline, Latitude, Longitude, Climate oscillation, Heremites vittatus, Intraspecific variation Background Intraspecific patterns of biological diversity are often the result of geological and climatic processes in the past. In the western Palearctic, the best understood example is the widespread presence of latitudinal gradients in genetic diversity in European taxa [1][2][3][4]. During interglacial periods in the Pleistocene, these species rapidly and repeatedly extended from Mediterranean refugia into central and northern Europe. Genetic bottlenecks during the cyclical expansions led to allele loss and decreased heterozygosity, and ultimately a reduction in genetic diversity in northern areas. While latitudinal gradients can be found in many species in Europe, the geographical shape of southern Europe made longitudinal movements more infrequent. Longitudinal patterns of genetic diversity are thus often confined to the different Mediterranean peninsulas, and can be difficult to disentangle from latitudinal gradients [5,6].
Compared to continental Europe, relatively little is known about the structure and evolution of biological diversity in other regions of the Palearctic, notably in the Middle East and northern Africa [7][8][9][10][11][12]. Several species in northern Africa display post-glacial latitudinal range expansions, often into Europe despite the potential barrier of the Mediterranean Sea [13]. Remarkably, many species in this region share a roughly longitudinal distribution extending from Iran and Turkey along the eastern Mediterranean coast into northern Africa [6]. Unlike in Europe, strong sea barriers are absent in this region, which potentially facilitated longitudinal movements. In addition, the presence of the Sahara as a southern barrier for many species may have further supported longitudinal migration. European species with more pronounced longitudinal patterns of genetic diversity often originated in Asia Minor and expanded westwards along the north Mediterranean coast. As such, Asia Minor has become known as a center of genetic diversity in the Middle East [14][15][16]. However, whether species that expanded from Asia Minor into northern Africa instead of Europe are also characterized by longitudinal patterns remains poorly understood. More generally, few studies have addressed the biogeographic mechanisms and underlying evolutionary timescales in species with longitudinal distributions in this region [5,6,17,18]. Here, we examined these questions in the striped skink, Heremites vittatus, a scincid lizard with a predominantly longitudinal distribution across the Middle East and northern Africa.
Lizards of the genus Mabuya sensu lato (s.l.) are some of the most widely distributed skinks with a circumtropical distribution. For a long time, these skinks were collectively allocated to the genus Mabuya s.l., until the Afro-Malagasy Mabuya taxa were assigned to the genus Euprepis [19,20], and later re-assigned to Trachylepis [21]. The Middle Eastern Mabuya s.l. species were preliminarily included into Trachylepis, although they form a separate radiation [22]. To account for this, the genus Heremites was recently revalidated for these species, including H. vittatus, H. auratus, and H. septemtaeniatus [23]. Of these, H. vittatus has the largest distribution range, occurring in Algeria, Tunisia, Libya, Egypt, Israel, Jordan, Lebanon, Cyprus, Syria, Turkey, Iraq and Iran [24][25][26][27][28][29]. Despite its wide distribution, the striped skink has been regarded as a monotypic species [30]. In several parts of this large distribution range, e.g., in Turkey, Iran, and Lebanon, morphometric studies uncovered considerable morphological variation within local populations and differences between populations on a regional scale [31][32][33][34][35][36]. Part of this variation is likely due to adaptation to local environmental conditions, rather than phylogenetic divergence. For example, in a population in southeastern Turkey, uniform skinks are more abundant in habitats with low grass cover, while striped specimens are more abundant in habitats with high grass cover, indicative of disruptive selection by visual predators [36]. A distribution model for H. vittatus demonstrated that the species is predominantly found in areas with high winter precipitation (>300 mm), and rainy winters may be the driving factor behind its distribution [37]. Populations in northern Africa inhabit wetlands and oases [24], while those in Iran and Turkey live in mountainous areas [26]. Due to its poor dispersal skills and its dependence on humid habitats in the arid Middle Eastern environment, the striped skink is a sensitive indicator of geographic processes that have driven the distribution and intraspecific evolution of animals in this region [38].
In this work, we combined investigations of mtDNA (cytochrome b) sequence and multivariate morphological variation across the species' distribution range, and provide the first comprehensive analysis of the biogeography and evolutionary history of the striped skink. We chose to study both genetic and morphological variation because evolutionary processes can affect genetics and morphology in different ways and at different times during the period of divergence. Unlike putatively neutral molecular markers, morphological traits are under varying degrees of selection; thus, it may be expected that morphological traits can show both clinal variation as a result of local adaptation, as well as trait conservation as a result of deep divergence events [39]. Together, insights from genetic and morphological variation can thus provide for a richer and more integrated understanding of the underlying evolutionary processes [39,40].

Morphological analysis
We first performed a non-metric multidimensional scaling (NMS) analysis to assess the morphological variation in an unbiased way. A 3D plot of the first three NMS axes did not reveal any discrete grouping of samples in the morphospace (Fig. 6). To test for an effect of geographical location on the position in the morphospace, we color-coded samples according to their latitude ( These results could be explained by morphological divergence between the three main mitochondrial lineages (Tunisian/Libyan, Turkish/Iranian, eastern Mediterranean) in our phylogeographic analysis, because these lineages are distributed in an allopatric, roughly longitudinal way. To examine this hypothesis in our morphological dataset, we assigned specimens based on their collection localities to three morphological groups (western, central, eastern) coinciding with the putative ranges of the three mitochondrial lineages (Fig. 7a). Their ranges are delimited by prominent barriers for animal dispersal: The Tunisian/Libyan lineage is separated from the eastern Mediterranean lineage by the arid Marmarica region (located between the Cyrenaica and the Nile delta) [41], while the eastern Mediterranean lineage is separated from the Turkish/Iranian lineage by the Amik plain in southern Turkey (Fig. 7a) [42,43]. We then tested traits individually for significant mean differences between these groups. We excluded binary and meristic traits that showed little variation and are unlikely to contribute to the pattern of variation, and focussed on five meristic and 14 mensural traits. Consistent with the hypothesis of morphological divergence between the mitochondrial lineages, we found that the majority of traits showed significant mean differences between the western, central, and eastern group (13/19 traits, 68%; Table 2, Fig. 7f).   Alternatively, the effect of latitude/longitude on the position in the morphospace could be explained by clinal variation within geographical groups. While the western (Algeria, Tunisia, Libya) and eastern (Turkey, Iran) groups have a roughly longitudinal distribution, the central group (Egypt, Israel, Jordan, Lebanon, Cyprus, Syria) is predominantly distributed along the latitudinal dimension (Fig. 7a). To account for these different geographical extensions, we regressed traits against longitude for the western and eastern group, and against latitude for the central group. 10 of 19 traits (53%) showed evidence of clinal variation in one geographical group. No trait showed evidence for clinal variation in more than one group. Of the 10 traits, seven showed clinal variation in the central group. Two traits showed clinal variation in the western group, and one trait in the eastern group (Fig. 7f).
We next investigated how mean differences between groups overlap with clinal variation within groups. Of 19 traits, three (16%) showed evidence of neither mean nor clinal variation (Fig. 7b). Six traits (32%) showed evidence of only mean differences (Fig. 7c), while three traits (16%) showed evidence of only clinal variation. Seven traits (37%) showed evidence of both mean and clinal variation. Of the latter seven traits, five showed clinal variation in the central group (Fig. 7d-e). Taken together, our results suggest that the majority of morphological traits in H. vittatus exhibit mean differences between geographical groups with either no evidence for clinal variation or clinal variation along the eastern Mediterranean coast.

Historical biogeography
Our mitochondrial phylogeography suggests that Heremites vittatus diverged from H. auratus in the late Miocene (13.3 mya; range 7.9-19.7 mya), consistent with a previous estimate that dated this split to 10.54 mya [44]. Heremites vittatus likely emerged as a distinct species in southwest Asia, where the Heremites radiation originated [23]. In the mid Pliocene (3.6 mya; range 2.3-5.2 mya), the ancestors of the eastern Mediterranean and northern African lineages split off from the ancestral Irano-Anatolian lineage and became allopatric. This divergence event coincided with prominent geological events in Anatolia during the late Neogene. The Anatolian mountain ranges, e.g., the Amanus and Taurus mountains, created geographical barriers between the Anatolian and eastern Mediterranean regions in the South. The vast lake system in central Anatolia delimited the distribution into western Anatolia. Other emerging dispersal barriers subsequently supported this divergence, notably the Amik plain in southern Turkey. As a result, vicariant cladogenesis, frequent among other animals in this geologically complex region (e.g. green lizards), may explain the divergence of the Irano-Anatolian lineage and the ancestors of the eastern Mediterranean and northern African lineages [14,43,45,46].
During much of the Pliocene, the climate was relatively warm and humid in the Palearctic [47], thus potentially facilitating the colonization of northern Africa from the eastern Mediterranean. In the late Pliocene, however, the climate suddenly aridified and became colder [41,48,49]. Around this time, the northern African lineages split off (2.5 mya; range 1.6-3.6 mya). Consistent with the preference of the modern species for more humid habitats in northern Africa [24,37], this climatic turn may have restricted the species to the remaining humid regions, and isolated them from populations inhabiting the eastern Mediterranean coast. Similar divergence events coinciding with the transition from hotter/wetter to colder/drier climate in the late Pliocene/ early Pleistocene have previously been observed in other animal species of the Mediterranean (e.g., tree frogs, Typhlops vermicularis) [8,11,50]. Notably, H. vittatus is the only Heremites species in northern Africa today, suggesting that other Heremites species never migrated into northern Africa, or were unable to withstand the drastic climatic changes at the Pliocene/Pleistocene transition.
During the Pleistocene, the Sahara increasingly aridified, which restricted more humid habitats to a few coastal areas (e.g., Cyrenaica) and oases [41,51]. This aridification may have led to the further fragmentation of the northern African lineage, as evidenced by the deep split between the Tunisian and Libyan haplotypes in the Cyrenaica (1.9 mya; range 1.0-2.9 mya). Although these climatic cycles must have influenced distributions during the Pleistocene, the mitochondrial lineages remained separate and apparently did not mix. The extant, allopatric populations in Libya, Tunisia, and Algeria are thus potential remnants of a once more extended distribution in northern Africa. Similar fragmentation processes in the Pleistocene likely affected other animal species in northern Africa currently limited to the Cyrenaica and coastal areas of Tunisia and Algeria, such as lizards of the genus Ophisops [12].
In the late Pleistocene, the Irano-Anatolian and eastern Mediterranean lineages further radiated, establishing  Based on our dataset, most of these lineages occupy nonoverlapping regions, suggesting that these haplotypes are not admixed. One interesting exception is the cooccurrence of the Jordanian-Lebanese-Syrian MHG with the Cypriot-Israeli-Jordanian MHG at one locality in Jordan (Al Himma), raising the possibility of admixture between these lineages in this region. In our dataset, the monophyletic lineage on the island of Cyprus is sister to the lineages in Israel and Jordan and diverged approximately 0.5 mya (range 0.3-0.8 mya). During the Pleistocene, Cyprus was likely never connected to the mainland through a land bridge [25,52], so that the species must have reached Cyprus via overseas dispersal, which has also been suggested for other species in Cyprus, e.g., Hyla savignyi [8,25,53]. However, these hypotheses may not be correct if the true mainland source of the Cypriot population was not included in our dataset. Although our study shows again the usefulness of fastevolving mtDNA genes to dissect more recent intraspecific divergence in Palearctic reptiles (e.g., [54][55][56][57][58]), we cannot rule out the possibility of discordant evolutionary histories of the mitochondrial and nuclear genome [59][60][61]. The mtDNA lineage evolution described in this study should thus be corroborated with nuclear sequence variation before additional inferences (e.g. for taxonomic purposes) are made on the evolutionary history of this complex.

Morphological evolution
Our multidimensional scaling analysis suggests that morphological variation in Heremites vittatus is linked to latitude and longitude. At the level of individual traits, most show mean differences between three morphological groups (western, central, eastern) coinciding with the putative ranges of the three old mitochondrial lineages. We note that while these results are suggestive, proof of concordant genetic and morphological variation would require data collection from identical specimens. Interestingly, we further find that these traits show either (1) no evidence for within-group clines, or (2) clinal latitudinal variation in the central (eastern Mediterranean) group.
The first category, mean differences with no evidence for within-group clines, argues for a close association of trait variation with the three old mtDNA lineages. For example, the western (northern African) group has an average of 18.2 subdigital lamellae, while the central (eastern Mediterranean) and eastern (Turkish/Iranian) group have on average 16.4 and 16.1, respectively, with no evidence for a clinal increase of lamellae in the central group (Fig. 7c). Thus, although they are likely not closely related the eastern and central group share a similar number of subdigital lamellae, which may reflect an ancestral state that was independently maintained in both lineages. By contrast, the higher number of subdigital lamellae in the western group may be an adaptation to the xeric environmental conditions across northern Africa because sand-dwelling lizards often possess pedal specializations that evolved as (convergent) adaptations to sandy habitats [62]. Mean differences between geographical groups can thus be explained by lineage-specific ancestral or derived trait states.  The second category, mean differences between groups with latitudinal clines within the central (eastern Mediterranean) group, suggests that some traits are affected by selection gradients in addition to lineage-specific trait states. One such trait is relative head length. Lizards in the western group have on average longer heads relative to their body length than the eastern group, with no evidence of clinal variation in either group. By contrast, relative head length in the geographically intermediate central group decreases gradually from the western (northern African) to eastern (Turkish/Iranian) trait states: the farther north specimens were collected along the eastern Mediterranean coast, the shorter their heads are relative to body length. The large extension (roughly 1000 km) of this cline suggests that the central group is likely not a hybridization zone of the western and eastern group, but rather a distinct population affected by gradual environmental change. Remarkably, we found little evidence that similar clines exist in the western or eastern group, pointing to unique environmental changes at the transition from northern Africa to Anatolia [45,46].
Taken together, our morphological analysis reveals that morphological variation in Heremites vittatus has a clear geographical organization, and is likely both affected by lineage-specific trait conservation and local trait adaptation.

Conclusions
Many species in the Middle East and northern Africa exhibit predominantly longitudinal distributions, and are potentially characterized by longitudinal patterns of biological diversity, but the underlying biogeographic mechanisms and evolutionary timescales remain largely unclear. In this paper, we provide the first comprehensive assessment of the mitochondrial phylogeography and multivariate morphological variation of the striped skink, Heremites vittatus, a species with a predominantly longitudinal distribution across the Middle East and northern Africa. We suggest that H. vittatus evolved as a separate species in southwest Asia in the late Miocene. We show that the species diverged into three main mitochondrial lineages roughly 2.5-3.6 mya in the eastern, central, and western regions of the distribution range, introducing a longitudinal organization of genetic diversity. This divergence was likely driven by a combination of geological changes in Anatolia and climatic changes in northern Africa. The morphological variation reflects this pattern of mitochondrial divergence: the majority of morphological traits show significant differences between these regions, supporting the idea of independently evolving lineages. We also find that subsequent local radiation in the Pleistocene, likely driven by climatic oscillations and locally emerging geographical barriers, led to the mitochondrial diversity observed today. Morphologically, clinal variation along the eastern Mediterranean coast in a subset of traits suggests that ongoing selection pressures and local adaptation may play an important role in shaping populations at the transition from northern Africa to Anatolia. Our findings of allopatric, morphologically and genetically divergent lineages raise the possibility that Heremites vittatus represents a complex with several undescribed taxa. However, taxonomic decisions should be based on nuclear, mitochondrial, and morphological data collected in the context of future studies on the inter-and intraspecific systematics and taxonomy of Heremites. In conclusion, our study suggests that longitudinal patterns and species distributions in the Middle East and northern Africa may be the result of complex evolutionary processes, driven by the region's geological and climatic history, geographical setup, and the presence of the Sahara.

DNA sequence analysis
The DNA sampling consists of 46 H. vittatus specimens (Additional file 3), including one sequence from GenBank, covering most of the distribution range (Fig. 1a). We also sequenced five H. auratus specimens. We included Trachylepis quinquetaeniata from the African lineage of Mabuya s.l. as an outgroup, because this lineage is one of the closest relatives of Middle Eastern Heremites [23,[63][64][65]. We did not include in the analyses sequences of H. vittatus from GenBank that only partially overlap the alignment generated in this study (but see Additional file 2).
Small pieces of tissue were sampled either from museum specimens (toe or tongue) of various ages and storage conditions, or were provided by colleagues based on field collections (part of tail or liver). Samples were stored in 70% ethanol or EDTA buffer, and genomic DNA was extracted with standard protocols [66]. A portion of the mitochondrial cytochrome b gene was amplified by polymerase chain reaction (PCR) with the primers mt-c-emys (5′-CCG GAT CAA ACA AYC CAA CAG G-3′) [67] and mt-fs-h (5′-CCA GTA GAA CAC CCA TTC ATC ATC ATT GGC CAA CTA-3′) [68]. This resulted in an amplicon with a length of 394 bp, corresponding to positions 14,762-15,155 in the mitochondrial genome of Eumeces egregius (positions 633-1026 in the cytochrome b gene of E. egregius) [69].
DNA sequences were aligned with the ClustalW algorithm implemented in BioEdit [70], and SNPs were verified by checking the sequence chromatograms. No insertions or deletions were found in the alignment. We deposited new sequences in GenBank under accession numbers MF101923-MF101972 (see Additional file 3). We used MEGA5.2.2 [71] to determine nucleotide diversity and to calculate uncorrected p-distances based on pairwise deletion of ambiguous sites, and DnaSP [72] to identify haplotypes and determine haplotype diversity.
We then divided the dataset by codon position, and ran PartitionFinder 1.1.1 [73] to find the best partition scheme and associated substitution models. The best partition scheme included the three codon positions each as a separate partition (lnL = −1930.14357, BIC = 4553.24905498), with the following models: K80 + G (position 1), HKY + I (position 2), TrN + G (position 3).
We used BEAST v1.8.1 and associated software tools [74] to obtain time-calibrated estimates of phylogenetic divergence. We used the partitioned dataset, and set the substitution models available in BEAUTi equivalent to the PartitionFinder model selection for the partitions separately (partition 1: HKY with gamma site heterogeneity; partition 2: HKY with invariant site heterogeneity; partition 3: TN93 with gamma site heterogeneity). In accordance with previous recommendations [75], we used an uncorrelated lognormal relaxed clock model with a constant coalescent tree prior and a random starting tree for all three partitions. Since dated fossils are not available for any Heremites species or related genera (Trachylepis, Chioninia, Eutropis, Mabuya sensu stricto), we instead generated time-calibrated divergence times through estimates of the substitution rate. Substitution rates for cytochrome b in other skinks (Eumeces, Chalcides, Scincus), geckos (Hemidactylus), and lacertid lizards (Lacertini) range from 1.15-1.35% per million years [64,[76][77][78]. We used a normally distributed prior of the substitution rate (ucld) with an initial mean of 1.25% and a standard deviation of 0.5% for all partitions in BEAUTi, and then optimized these values by checking for convergence and high ESS scores in Tracer v1.6. We then ran five separate rounds of each 10 7 MCMC iterations and logged parameters every 1000 steps, which generated 50,000 trees. We combined these runs in LogCombiner v1.8.1 with a burn-in of 10% of each run, generated a consensus tree with TreeAnnotator v1.8.1, and produced the final tree with FigTree v1.4.2.
In addition, we calculated a maximum likelihood (ML) tree with the program RAxML 7.0.4 [79] using the rapid hill climbing algorithm [80]. The dataset was partitioned into the three codon positions corresponding to the PartitionFinder partitioning strategy, and run under the suggested [79] GTR + G substitution model in RAxML.
Sequence evolution within a species does not necessarily follow dichotomous patterns, because co-existence of ancestral and derived haplotypes may introduce reticulate patterns of evolutionary relationships. To account for this phenomenon, we computed a haplotype network with the software NETWORK v4.611 [81] (http:// www.fluxus-engineering.com), using the median-joining (MJ) algorithm with default settings (epsilon = 0). We assigned haplotypes to major haplogroups (MHGs) with the Bayesian implementation of the Poisson tree processes (PTP) model [82], a method to delimit putative species on a phylogenetic tree, using the BEAST tree (Fig. 4) and the PTP default settings (100,000 MCMC generations, thinning 100, burn-in 0.1, seed 123).

Morphological analysis
We examined 146 H. vittatus specimens from across the distribution range of the species (Additional file 4, Fig. 1b). Specimens were only included in the dataset if the collection locality could be reliably geo-referenced. We did not distinguish among sexes of the examined specimens. For bilateral traits, the left character state was recorded. We measured 60 traits and excluded traits that showed no or only sporadic variation, leaving two binary, eleven meristic and 14 mensural traits for further analyses (Additional file 4). All mensural characters, except for snout-vent length, were first normalized for size by dividing through the snout-vent length (TL, HT, VT, ALL, HL, TH, LH, LF, LT, LA) or the head length (SL, HB, INA), respectively. Data were then ln(x + 1)-transformed. We first conducted a non-metric multidimensional scaling analysis (NMS), which is a statistical approach to reduce the dimensionality of the data set by estimating the similarity of every pairwise comparison of samples based on a similarity coefficient. We used the Gower index (= rangenormalized Manhattan index) for meristic and mensural traits, and the Jaccard index for binary characters as defaulted in the program PAST [83]. We retrieved the first three NMS axes and performed a multivariate multiple regression to test for the effect of geographical location on the 3D position in morphospace. We summarized the results with a type 2 MANOVA using Pillai's test to obtain separate η 2 point estimates for the effect of latitude and longitude, and conducted a Bootstrap analysis with 10 5 repeats to obtain 95% confidence intervals for the η 2 point estimates. We restricted subsequent analyses to five meristic and 14 mensural traits that showed substantial variation, and assigned specimens to three geographical groups (western, central, eastern). We conducted pairwise multiple comparison tests for significant mean differences between these groups using Dunnett's modified Tukey-Kramer test for uneven sample sizes and heterogeneous variance as implemented in the "DTK" R package. We also regressed individual traits separately for each group against the longitude (western and eastern groups) and latitude (central group). All statistical analyses were done in R software v3.2.4 [84].

Additional files
Additional file 1: PTP tree used to delimit major haplogroups. Values on nodes are Bayesian support (BS) values. Higher BS values indicate that descendants from this node are more likely to be from one species. Nodes and their descendant branches that were assigned to the same species by PTP are red; singleton species are left in blue. (PDF 248 kb) Additional file 2: Maximum Likelihood (ML) tree based on combined sequences of this study and cytochrome b sequences of Heremites vittatus from Turkey (in green, retrieved from GenBank [44]) that partially overlap (positions 1-187) with the cytochrome b alignment (394bp) presented here. The tree was calculated with RAxML 7.0.4 using the climbing hill algorithm. The dataset was partitioned into the three codon positions, and run with a GTR+G substitution model in RAxML. (PDF 127 kb)