Sand fly trapping in the Greek Aegean Islands was conducted during previous field work (2016–2019) (Tsirigotakis et al. 2018; Dvořák et al. 2020) using CDC miniature light traps (John W. Hock Co., Gainesville, FL, USA). The studied Aegean Islands included Crete, Cyclades (Andros, Sifnos, Milos, Folegandros, Santorini & Anafi), Dodecanese (Patmos, Leros, Nisyros & Karpathos), and North Aegean (Ikaria) (Fig. 1) (Tsirigotakis et al. 2018; Dvořák et al. 2020). Due to the presence of specimens with morphological uncertainties (e.g. specimens of P. similis), additional sampling took place during 2018 and 2019 in Cyprus (Agioi Trimithias) and Turkey (Kayseri & Antalya) (Omondi et al. 2020), in order to collect closely related species (e.g. P. sergenti) and enable us to resolve these uncertainties. Overall, 300 traps were used, with each trap being set for one night. Sand fly head, genitalia and wing of all collected specimens were mounted in CMCP-10 mounting medium (Polysciences, Inc., Warrington, PA, USA). At the same time, the rest of the body was stored in 70% ethanol at –20^oC. Species identification was made using the morphological characters of the head, wing and genitalia (Lewis 1982; Killick-Kendrick et al. 1991; Kasap et al. 2015; Dvořák et al. 2020).

All morphologically identified Phlebotomus spp. are shown in Table 1. The collected Phlebotomus fauna from our study areas included 12 species, with the most common being P. neglectus, P. tobbi, P. (Adlerius) simici Nitzulescu and P. similis. Phlebotomus killicki and P. creticus were found only in Crete, P. sergenti only in Cyprus and P. (Adlerius) halepensis Theodor only in Turkey. Also, a male specimen (“C147B”) of the subgenus Adlerius Nitzulescu collected from Andros could not be identified; nevertheless, it had morphological similarities with P. creticus.

The sand fly genomic DNA was extracted from the thorax, legs and the anterior part of the abdomen of each chosen specimen using the Qiagen QIAamp DNA micro kit (Qiagen, Hilden, Germany). Double-stranded PCR was performed to amplify partial sequences of two mitochondrial loci (mtDNA) and four nuclear loci (nDNA). These loci were cytochrome c oxidase subunit 1 (COI), cytochrome b (CytB), internal subscribed spacer 2 (ITS2), domain 9 and 10 of 28S ribosomal RNA (28S), elongation factor 1 alpha (EF1-α) and triose-phosphate isomerase (TPI). All PCR reactions had 25μl volume and were performed in T100 thermal cycler (Bio-Rad Laboratories, California, USA). Primers and PCR conditions are described in Table 2. Single-stranded Sanger sequencing was performed in CeMIA SA (Larisa, Greece) using the same primers as in PCR.

Sequence chromatograms were viewed and edited using CodonCode Aligner v.9.0.1 (CodonCode Corporation, Centerville, USA). Multiple sequence alignments for each locus were performed using MAFFT v.7.475 (Katoh and Standley 2013) with the default settings and the E-INS-i alignment strategy. Our alignments consisted of our samples and some additional sequences retrieved from the GenBank database, belonging to species that are closely related to our samples. The NCBI accession numbers of these sequences are shown in Table S1 of Supporting information. Protein coding sequences were translated into amino acids using Sequence Manipulation Suite (Stothard 2000) to check for the presence of stop codons.

The optimal nucleotide substitution model for each locus was identified using PartitionFinder v.2.1.1 (PF) (Lanfear et al. 2016). For each protein-coding locus, the alignment was partitioned into three blocks containing the 1^st, 2^nd and 3^rd codon positions. We ran PF three times for each locus using the greedy search algorithm with linked branch lengths in calculations of likelihood scores under the Bayesian information criterion (BIC). The difference in those runs was the restriction of the candidate models to those that were available in MrBayes v.3.2.6 (Ronquist et al. 2012), IQ-TREE v.1.6.12 (Nguyen et al. 2015), and BEAST2 v.2.6.2 (Bouckaert et al. 2019). The models which included both gamma distribution and invariable sites were ignored because of the high correlation between the proportion of the invariable sites and the gamma shape, making it extremely difficult to estimate both parameters reliably (Yang 2006).

All retrieved sequences were submitted to GenBank and their NCBI accession numbers are given in Table S1 of Supporting information. The aligned mtDNA, nDNA and concatenated datasets consisted of 1,080, 1,976 and 3,056 base pairs (bp), respectively. Information on conserved, variable and informative sites is given in Table 3. PartitionFinder indicated the best-fit nucleotide substitution models for each locus and each analysis, which are summarised in Table S2 of Supporting information.

Genetic distances were calculated using the Tamura-Nei model (Tamura and Nei 1993) as implemented in MEGA v.7 (Kumar et al. 2016) to better compared with other studies. Phylogenetic trees were constructed using Bayesian Inference (BI) and Maximum Likelihood (ML) on the mtDNA, nDNA gene trees and the concatenated alignments. The BI analysis was performed in MrBayes with four runs and eight chains per run for 10,000,000 generations, with a sampling frequency of 100 using the substitution models selected by PF. To check for convergence and stationarity, we used several convergence diagnostics, such as the average standard deviation of split frequencies (ASDSF-values below 0.01 indicate convergence), the log-likelihood values, the Effective Sample Size (ESS-values below 100 may indicate insufficient sampling) and the average Potential Scale Reduction Factor (PSRF-values should be close to 1 for all parameters). The –lnL stabilised after approximately 10⁶ generations, and the first 25% of the trees were discarded as ‛burn-inʼ phase. A majority-rule consensus tree relied on the remaining trees, and posterior probabilities were calculated as the percentage of samples recovering any particular clade (Huelsenbeck and Ronquist 2001). The ML analysis was performed in IQ-TREE with bio-neighbor joining starting tree under the suggested models selected in PF. Two independent runs were performed with bootstrap (Felsenstein 1985) and ultrafast bootstrap (Minh et al. 2013) values, which were estimated by 1,000 and 10,000 replicates, respectively. Sergentomyia minuta Rondani represented the outgroup in these analyses. FigTree v.1.4.4 (Rambaut 2018) was used to display the phylogenetic trees and Tracer v.1.7.1 (Rambaut et al. 2018) was used to assess the convergence statistics.

We performed a molecular clock test in MEGA before the divergence times estimation analysis was applied to the studied loci. The null hypothesis of equal evolutionary rates throughout the tree was rejected at a 5% significance level (p=<0.05), thus we applied an uncorrelated lognormal relaxed clock for the time estimation (Drummond et al. 2006).

The coalescent species tree and the estimation of divergence times were calculated using Starbeast2 (Ogilvie et al. 2017) as included in the BEAST2 software package (Bouckaert et al. 2019). We used BEAUti v.2.6.2, as implemented in the BEAST2 software package, to create the input file. We used the partitions and substitution models selected by PF (Lanfear et al. 2016). Also, the Yule Model was used for the speciation and the Uncorrelated Lognormal Models for the relaxed molecular clock. The Markov chain Monte Carlo (MCMC) analysis run for 4*10⁸ generations and saved the result every 5,000^th generation. Tracer was used to verify that the analysis had converged and that the satisfactory Effective Sample Size had been obtained. The –lnL stabilised after approximately 4*10⁷ generations, and the first 10% of the trees were discarded. Sergentomyia minuta represented the outgroup, and FigTree was used to display the species tree. To estimate divergence times, we used the available gene-specific substitution rates on CytB from previous studies. Specifically, the substitution rate for CytB was 2.1% (1.57–2,64%), and the clock rate was set at 0.0105 (Esseghir et al. 2000).

Biogeographic analysis was constructed in Reconstruct Ancestral State in Phylogenies (RASP) v.4.2 (Yu et al. 2015) and the ancestral distributions were reconstructed using the Statistical Dispersal-Vicariance Analysis (S-DIVA) (Yu et al. 2010), allowing a maximum of two areas per node. The analysis was conducted using the trees and the dated species tree generated from the BEAST2 software package (Bouckaert et al. 2019). Five geographical regions were defined: A, Middle East; B, Aegean islands; C, East Europe; D, West Europe; E, North Africa.

For conducting species delimitation, we used BP&P v.4.3 (Flouri et al. 2018), which uses the multispecies coalescent model to compare different models of species delimitation and phylogeny in the Bayesian framework. This method accounts for incomplete lineage sorting due to polymorphism and differences between gene trees and species tree (Yang and Rannala 2010). We conducted the A11 analysis for the joint species delimitation and species tree inference of unguided species delimitation (Yang and Rannala 2014), which attempts to group different populations into one species and explore different phylogenetic relationships among the delimited species (Yang and Rannala 2014). Our dataset comprised of five different loci (mtDNA, EF1-α, TPI, 28S & ITS2). The samples were divided into 17 potential distinct species that corresponded to the major clades of the phylogenetic concatenated tree and the morphological identification. The rjMCMC analysis was conducted for 100,000 generations, with a sampling interval of three, ‘burn-in’ phase of 2,500 and each species delimitation model had equal prior probability. We tested three combinations of priors (Leaché and Fujita 2010), where the prior distribution on the ancestral population size (θ) and the prior for the root age (τ₀) were set either to 0.1 [IG(3, 0.2)] or 0.001 [IG(3, 0.002)]. The first combination assumed relatively large ancestral population sizes and large divergences [θ~IG(3, 0.2) & τ₀~IG(3, 0.2)]; the second combination assumed low ancestral population sizes and low divergences [θ~IG(3, 0.002) & τ₀~IG(3, 0.002)] and the third combination assumed relatively large ancestral population sizes and low divergences among species [θ~IG(3, 0.2) & τ₀~IG(3, 0.002)]. Each analysis was run twice with different seed numbers (positive and negative integer) to confirm the consistency between the runs. The starting tree had the topology of the coalescent species tree.

The genetic distances for mtDNA between the local Phlebotomus spp. ranged from 2.78% to 28.35% and for nDNA ranged from 0.4% to 23.22%. Genetic distances among and within species for mtDNA and nDNA (except ITS2) are given in Table 4. Sequence divergence for the COI ranged from 4.32% to 22.04%, for CytB from 3.79% to 32.24%, for EF1-α from 0% to 13.99%, for TPI from 0.93% to 21.81%, for 28S from 0.73% to 9.92% and for ITS2 from 1.94% to 32.8%. The genetic distances for all loci separately are presented in Tables S3–S5 of Supporting information.

Both BI and ML analyses on the concatenated dataset resulted into similar and well-supported phylogenetic trees (BI: lnL = –15,458.60; ML: Bootstrap lnL = –15,314.19 and Ultrafast Bootstrap lnL = –15,314.12) (Fig. 2). According to our results, all local species were monophyletic with very good statistical support. Also, each subgenus representative formed distinct monophyletic clades, with Artemievus being more closely related with Phlebotomus rather than Paraphlebotomus. Moreover, subgenus Transphlebotomus is more closely related to Adlerius than Larroussius. The Adlerius specimen (“C147B”) from Andros Island, which was morphologically similar to P. creticus, was placed in the same clade with P. creticus from Crete; thus, we assigned it to P. creticus, whose closest relative appears to be P. (Adlerius) balcanicus Theodor. As for the species of the P. major complex, P. (Larroussius) syriacus Adler & Theodor was the earliest branching lineage and P. major Annadale was the sister species to P. neglectus. The male specimens “M292” and “M293” (P. cf. major), which were morphologically different from the males of P. neglectus, were grouped with P. major sequences.

Figure 2. 

Bayesian inference concatenated phylogenetic tree of the six loci. Numbers next to the branches represent posterior probabilities (left), bootstrap values (middle) and ultrafast bootstrap values (right).

The mtDNA and the nDNA datasets produced trees with lnL = –7,489.09 and –7,848.33 for BI, respectively, and; lnL = –7,312.67 and –7,679.85 for Bootstrap ML, respectively and lnL = –7,312.06 and –7,679.21 for Ultrafast Bootstrap ML, respectively. The phylogenetic analyses produced less resolved trees but with similar topologies to that of the concatenated gene tree (Fig. S1 in Supporting information). The most important differences were the positions of the P. major complex, P. sergenti and P. similis. More specifically, according to the nDNA phylogenetic tree, P. sergenti and P. similis form a single well-supported clade. In contrast, the topology of these two species in the mtDNA gene tree was identical to the concatenated phylogenetic tree. Regarding the P. major complex, P. neglectus was resolved as paraphyletic in the mtDNA alignment tree since its lineage contained also the specimens of P. major and P. cf. major (Fig. S1 in Supporting information).

The coalescent species tree analysis resulted in good ESS values (lnL = –16,111.83) and the species tree presented very good posterior probabilities (Fig. 3). Phlebotomus alexandri is more closely related to P. papatasi rather than P. sergenti and P. similis, which are closely related. Phlebotomus creticus is more closely related to P. balcanicus rather than P. halepensis or P. simici. Phebotomus perniciosus, which is distributed in the west Mediterranean (Léger and Depaquit 2002), appears to be closely related to P. tobbi (east Mediterranean). The species of the P. major complex appear to be closely related between them, with P. syriacus being the root lineage of the complex.

Figure 3. 

Dated coalescent species tree with posterior probabilities and mean divergence times (above the nodes). The 95% HPD are displayed below the nodes; ancestral geographical distribution is displayed in the up-left corner; letters next to the species correspond to their geographical distribution for the S-DIVA analysis (A: Middle East; B: Aegean Islands; C: East Europe; D: West Europe; E: North Africa).

According to divergence time estimation (Fig. 3), the studied local species of the genus Phlebotomus separated from S. minuta at 34.37 mya (95% HPD 21.79–49.27) and started to diverge in the Oligocene at 29.92 mya (95% HPD 18.89–42.36), with the separation of Artemievus, Phlebotomus and Paraphlebotomus from all the other subgenera. The closely related species P. similis and P. sergenti diverged in the Pliocene at 3.62 mya (95% HPD 1.86–5.55). The estimated age for the split between the representatives of the subgenus Larroussius and those of the subgenera Transphlebotomus and Adlerius was 21.87 mya (95% HPD 14.72–29.73) (early Miocene). Phlebotomus major complex started diverging in the early Pleistocene at 2.45 mya (95% HPD 0.87–4.47), and P. creticus diverged from its closest relative (P. balcanicus) also in the early Pleistocene at 2.11 mya (95% HPD 0.7–3.73).

The results of the biogeographic analysis are shown in Fig. 3 and Fig. S2 of Supporting Information. The S-DIVA analysis proposed that the distributional history for the studied local species of Phlebotomus was dominated by dispersal. More specifically, the results indicated that 33 total dispersal events and only two vicariant events (nodes 5 & 8 in Fig. 3) were necessary to explain the current species distribution. Most dispersal events occurred from A (Middle East) to B (Aegean Islands) (11 events) and from A to C (East Europe) (eight events), while A was also the area with the most speciation events (14). According to the analysis, Middle East was the main contributing geographical area for most Phlebotomus species in this study. The majority of the geographical regions for the ancestral nodes were Middle East.

The results of the species delimitation analyses are given in Table 5. Assuming large ancestral population sizes and large divergences among species, the posterior probabilities of the species in the analysis ranged between 0.89 to 1, with the least supported species being P. major (0.89), P. cf. major (0.89) and P. balcanicus (0.90). The best-supported scheme included all 17 species groups and had a posterior probability equal to 0.76. The next proposed scheme included 16 species groups, with the fusion of P. major-P. cf. major. Assuming low ancestral population sizes and low divergences among species, species delimitation analysis supported a scheme of all 17 species-level clusters with a posterior probability equal to 0.81. The posterior probabilities of the species groups ranged from 0.85 to 1, with the least supported species groups being P. major (0.85) and P. cf. major (0.86). The posterior probabilities of the species for the third combination of priors (large ancestral population sizes and low divergences among species), ranged from 0.84 to 1. The least supported species were P. balcanicus (0.84) and P. major (0.90).

The interspecific genetic distances of the studied local species of Phlebotomus were similar to those found by Esseghir et al. (2000). Our phylogenetic analyses indicated that all studied local species were monophyletic and each subgenus representative formed distinct monophyletic clades. Nevertheless, species complexes such as P. major, P. perniciosus and P. perfiliewi, need further investigation by sampling the remaining taxa, in order to resolve their taxonomic status. Subgenus Artemievus (which is comprised of only P. alexandri) appears to be more closely related with Phlebotomus rather than Paraphlebotomus (where P. alexandri was previously classified). This is in agreement with previously published work, which showed that P. alexandri is more closely related to P. papatasi, P. (Phlebotomus) duboscqi Neveu-Lemarie and P. (Phlebotomus) bergeroti Parrot rather than the species of the subgenus Paraphlebotomus (Krüger et al. 2011; Cruaud et al. 2021). The two species of Paraphlebotomus in our study are closely related to each other, since they formed a single well-supported clade in the nDNA gene tree, while they formed two distinct, well-supported clades in the mtDNA gene tree (Figure S1 of Supporting information).

Subgenus Larroussius is closely related to Transphlebotomus and Adlerius, which agrees with the classification of morphological characters by Rispail and Léger (1998). Larroussius was divided into two major clades, with the first one comprising P. perfiliewi, P. perniciosus Newstead and P. tobbi. Phlebotomus perfiliewi was the earliest branching lineage, and this is congruent with the phylogenetic analysis made by Esseghir et al. (2000). The second clade included the species of the P. major complex with P. syriacus as the earliest branching lineage. Phlebotomus cf. major from Nisyros was placed in the same subclade with P. major, and P. neglectus appears as their sister species. Kasap et al. (2013) observed similar phylogenetic relationships in this species complex. The P. major complex consists of six species with P. neglectus, P. major and P. syriacus present in Turkey (Kasap et al. 2013, 2019). Kasap et al. (2013) examined the distribution of these species in Turkey and found three major lineages: P. syriacus (distributed in southeast Turkey); P. neglectus s.str. (sympatric with P. syriacus and the third lineage); and the third lineage composed of P. major sequences (distributed in central and southwest Turkey). These P. major sequences, which match those in our study, originated from specimens from Iran. Badakhshan et al. (2011) reported that most specimens from Iran were not representatives of P. major s.str. Kasap et al. (2013) suggested that these specimens may belong to a different cryptic species of the P. major complex. We assume that the specimens (seven males) found in Nisyros may belong to the same latter species.

The first lineage of subgenus Adlerius to branch off was P. simici, while the other three species appear to be closely related. According to our analyses, P. creticus formed a single well-supported clade with P. balcanicus as its closest relative. These results were also observed by Dvořák et al. (2020). The Adlerius specimen (“C147B”) from Andros was grouped with the specimens of P. creticus from Crete, and also, the genetic distance from P. creticus was 1.63% for mtDNA and 0% for nDNA. Consequently, we assume it represents P. creticus, which was found in cave entrances across all of Crete, while the Andros specimen was found in an abandoned stone building.

Another important finding of our study is the probable mitochondrial introgression between the species of the P. major complex. In the nDNA gene tree (Figure S1 of Supporting information), three distinct lineages were present, corresponding to the three taxa, and with P. syriacus as the earliest branching lineage. In the mtDNA gene tree (Figure S1 of Supporting information), all taxa were included in a single clade without separating the different species. Pesson et al. (2004) had suggested mtDNA introgression (most likely by hybridization) between P. perniciosus and P. (Larroussius) longicuspis Nitzulescu (Pesson et al. 2004), which are sympatric in southwest Mediterranean (Léger and Depaquit 2002). Esseghir et al. (2000) also suggested introgressive hybridizations between species whose ranges overlap [P. (Larroussius) orientalis Parrot & P. (Larroussius) langeroni Nitzulescu].

According to Bayesian species delimitation analyses, the least supported species are P. balcanicus, P. cf. major and P. major. The taxonomic status of P. balcanicus is unresolved probably because of the only two mitochondrial loci (CytB & COI) that are available in the GenBank database. Further molecular data (nuclear data) of P. balcanicus are needed to resolve its status. The relationship between P. cf. major and P. major is unresolved probably because they may belong to the same species or due to the few molecular data (CytB and EF1-α) available for P. major.

According to our dating analysis, the studied local species of Phlebotomus were separated from S. minuta at 34.37 mya, which is congruent with the suggestion of Tuon et al. (2008) and Akhoundi et al. (2016) that the genus Phlebotomus emerged in the Eocene epoch (33.9 to 55.8 mya). The emergence of Larroussius occurred at 21.87 mya (early Miocene) and started diverging at 13.35 mya (Serravalian age). Léger and Pesson (1987) associated the speciation of Mediterranean Larroussius with vicariance due to tectonic activities during early Miocene, while Esseghir et al. (2000) concluded that the speciation was not dependent on tectonic vicariance, but was affected by palaeoecological changes during late Miocene and that the faunal turnover and dispersion has determined the current species distributions (Léger and Pesson 1987; Esseghir et al. 2000). To extract safer assumptions on the factors that played key roles to the divergence of the subgenus Larroussius, additional sampling for more taxa is needed.

Phlebotomus perfiliewi diverged from P. tobbi and P. perniciosus during the pre-evaporitic stage of the Messinian age (7.25–5.96 mya) (Prista et al. 2015) through dispersal events from an Asiatic origin. However, this hypothesis contradicts the conclusion made by Trájer et al. (2021), who emphasized that during Messinian, the active dispersal of sandflies was very difficult or impossible due to the thermal conditions (Trájer et al. 2021). The only vicariant event in Larroussius lineages occurred at 4.24 mya (Zanclean age), diverging the lineages of P. tobbi and P. perniciosus. This lineage split occurred after the refilling of the Mediterranean Sea and the re-establishment of the marine conditions (Prista et al. 2015). Phlebotomus tobbi is distributed in eastern Mediterranean and Middle East (Léger and Depaquit 2002), while P. perniciosus is distributed in western Mediterranean, from Morocco to Libya and from Portugal to Croatia (Pesson et al. 2004). Therefore, we hypothesize that the ancestor of these two species may have occupied the Mediterranean region sometime before the refilling of the Mediterranean Sea (probably during MSC) and the latter may have acted as the vicariant event that split their lineages.

The speciation within the P. major complex appears to have coincided with the glacial and interglacial periods during the early Pleistocene (2.46 to 2.11 mya) (Poulakakis et al. 2015; Suc et al. 2019), and ancestral reconstruction showed that the complex has an Asiatic origin. During these periods, cycles of alternating hot, dry and cold, wet seasons were present in the Aegean area, with glacial periods outlasting the interglacial (Poulakakis et al. 2015). Consequently, we assume that the glacial and interglacial periods were the key factors for the speciation within the complex and the current distribution of each species was achieved through dispersal. Esseghir et al. (2000) also suggested an Asiatic origin; however, they suggested that the speciation occurred during or shortly after the MSC, which acted as a vicariant event.

The subgenus Transphlebotomus is restricted to the Mediterranean basin (Kasap et al. 2015) and the subgenus Adlerius is distributed mostly in Asia and eastern Mediterranean. Our subgenera representatives of Transphlebotomus and Adlerius, appear to have an Asiatic origin, and their split occurred at 11.84 mya (late Miocene) when the AB began its formation in the Aegean area (Poulakakis et al. 2015). This dating is congruent with that of Kasap et al. (2015). The studied local species of Adlerius have an Asiatic origin and the speciation amongst them was dominated by dispersal. Phlebotomus simici separated from the other species of Adlerius at 8.79 mya (after the formation of the AB). Therefore, we assume that the ancestor of P. simici was able to disperse from the Middle East to the Aegean and East Europe despite the presence of the AB. Due to the apparent crossing of the AB from the ancestor of P. simici, we propose that the AB had a minor influence on this diversification. Species of the woodlouse genus Dolichopoda and the tenebrionid beetle Dichomma dardanum are similar cases where the crossing of the AB was observed (Allegrucci et al. 2009, 2011; Papadopoulou et al. 2009).

The separation of P. halepensis from P. balcanicus and P. creticus coincided with the MSC (5.96–5.33 mya) (Krijgsman et al. 1999), which may have provided the opportunity for the dispersal of the ancestor of P. balcanicus and P. creticus from the Middle East to the Aegean through the formation of steppes and saline deserts. A similar case has been suggested for the separation of bush cricket lineages of the Poecilimon jonicus group (Orthoptera), which they may have taken the opportunity for dispersal during MSC, and were isolated during the refilling of the Mediterranean (Borissov et al. 2020). The ancestor of P. creticus and P. balcanicus had occupied the Middle East and the Aegean, and its diversification was due to a vicariant event during the Pleistocene. This event may have been the alternating glacial and interglacial periods. The current distribution of P. balcanicus consists of the Balkans, Turkey, Armenia and Georgia (Léger and Depaquit 2002; Cazan et al. 2019). Consequently, we assume that these periods may have acted as the vicariant event that separated these lineages, while the ancestor of P. balcanicus was able to disperse in the Balkans.

Finally, the separation of the subgenera Paraphlebotomus, Phlebotomus and Artemievus from each other appears to coincide with the Mid-Miocene Climatic Optimum (~17–15 mya), which represents a geologically warming event (Böhme 2003; Song et al. 2018). As for the closely related P. sergenti and P. similis, they appear to have separated before the establishment of the Mediterranean climate regime as we know it today (3.2 to 2.8 mya) (Blondel et al. 2010) in late Pliocene through dispersal from an Asiatic origin. According to Trájer et al. (2018), the ancestor of these species adapted to the hot and dry summers of the Mediterranean during late Neogene and the lowest minimum tolerable temperature of P. sergenti and P. similis is –7.2°C and –4°C, respectively (Trájer et al. 2018). Subsequently, we assume that this climate regime was the main factor of their divergence through climate adaptation. This assumption is incongruent with the hypothesis made by Depaquit et al. (1998), who suggested that during the Miocene, both species followed different migration routes, one north of the Paratethys and the other south of the Tethys (Depaquit et al. 1998). Our assumption is also incongruent with the alternative proposed hypothesis made by Trájer et al. (2021), in which P. similis, P. sergenti and P. (Paraphlebotomus) jacusieli could have separated after the tectonic subsidence of the Hellenic Orogenic Belt in late Miocene (Trájer 2021). Furthermore, Cruaud et al. (2021) suggested that the separation of the Paraphlebotomus lineages was probably the result of a Messinian vicariant event in which extreme prolonged drought has been instrumental (Cruaud et al. 2021).