The possible region of the Late Miocene split of the sandfly subgenus Transphlebotomus Artemiev and the early late Neogene to late Quaternary dispersal of the ancestor of Phlebotomus mascittii Grassi

The distribution of the Mediterranean Transphlebotomus species shows a marked zoogeographical dichotomy in the sense that Phlebotomus mascittii has a wide range in Europe, and the other species are restricted to the East Mediterranean region. The study aimed to investigate how the Neogene to late Quaternary climatic-geographical alterations could influence the split of the sandfly subgenus Transphlebotomus and the speciation of Phlebotomus mascitti. For this purpose, the climatic suitability patterns of the species were modelled for seven Neogene and Quaternary periods and the divergence times of Transphlebotomus clades were estimated. The model results suggest that the common ancestor of the extant Mediterranean-European Transphlebotomus species could be adapted to the Late Miocene climate of Western and Central Europe. Phylogenetic results suggest that the speciation of Ph. mascittii started in the Tortonian period, plausibly related to the rise of the Dinaric land bridge. The Central and Eastern Paratethys Seas could have played an important role in the split of the ancestral Phlebotomus mascittii populations and other Transphlebotomus populations. These other species can be the descendants of ancient Transphlebotomus populations adapted to the hotter and drier climate of the areas south of the Central and Eastern Paratethys. Their divergence could be strongly linked to the formation of the Aegean trench and, later, the Messinian salinity crisis. The Pliocene climatic fluctuations could result in habitat loss of Transphlebotomus populations in Europe which was particularly significant during glacial maxima such as the Last Glacial Maximum.


Introduction
Phlebotomine sandflies (Diptera: Psychodidae: Phlebotominae) are vectors of arboviruses and are the principal, if not the sole vectors of Bartonella bacilliformis (Strong et al. 1913) causing human bartonellosis in Latin America (Herrer and Christensen 1975). Most importantly, however, sandflies are the principal vectors of the protozoan parasites Leishmania spp. (Kinetoplastida: Trypanosomatidae), the causative agents of leishmaniasis (Akhoundi et al. 2016). The disease is endemic in 98 countries in Europe, Africa, Asia, and the Americas. Worldwide, around 12 million people are infected with Leishmania spp. and 50,000 to 90,000 cases of visceral leishmaniasis and 700,000 to 1 million new cases of cutaneous leishmaniasis are estimated to occur annually, with more than 1 billion people living in endemic areas. Leishmaniasis is particularly prevalent in the world's wet and semi-arid tropical regions (WHO 2021). Globally, the most affected regions are the tropical areas of Brazil, North and Northeast Africa, the Middle East, and South and East Asia (Burza et al. 2019). The New World and Old World-dwelling sandfly species have different climatic and breeding site preferences (Azar and Nel 2003).
In the Old World, several sandflies occur in arid, semiarid, (sub-) tropical and temperate areas in which they breed in a large variety of domestic, peri-domestic, and sylvatic 1 3 sites such wet cracks (Orshan et al. 2016), burrows of mammals or termite hills, animal barns close to human dwellings, and even tree holes or leaf litter (Feliciangeli 2004). It is worth mentioning that Phlebotomus species also occur in humid environments in Afro-Eurasia. For example, in Ethiopia, Phlebotomus (Synphlebotomus) martini Parrot, 1936 also inhabits the border zone of savanna and subtropical highland climate regions , where the annual mean precipitation reaches the 800-1000 mm value (Asefa et al. 2020). In Southwest Hungary, Phlebotomus (Larroussius) major subsp. neglectus Tonnoir, 1921 occurs in Nagyharsány (Tánczos 2012;Trájer et al. 2018a), where the annual rainfall is 572 mm (Trájer 2017). In the Americas, Lutzomyia França, 1924 (genus as proposed by Young and Duncan 1994) species occur mainly in tropical forests and savannas and breed predominantly in leaf litter (Dutari and Loaiza 2014) or within the tabular roots and bases of trees (Vivero et al. 2015). Female sandflies take blood on various vertebrate species, including mammals (e.g. lagomorphs; González et al. 2021) and birds (e.g. chickens; Sant'Anna et al. 2008), but the extant sandfly fauna also includes herpetophilic species, which predominantly feed on lizards (Pombi et al. 2020). In the case of mammals, burrowing species such as rodents (Yaghoobi-Ershadi and Javadian 1996), rabbits (Sáez et al. 2018) and hyraxes (Bsrat et al. 2015) play the most important role in sandfly ecology and the maintenance of local transmission cycles of leishmaniasis in the Old World.
Evolutionary hypotheses suggest a relationship between leishmaniasis' landscape epidemiology and sandflies' coevolution with Leishmania species of tropical and subtropical mammals and lizards (Ready 2013). Sandflies are a subfamily of the order Diptera L, whose members appear first in the geological record with ancient Nematoceran species in the early Middle Triassic epoch (Lukashevich et al. 2010). They emerged about the same time as the first lepidosauromorph reptiles, the ancestors of the extant lizards, appeared (Renesto and Posenato 2003). The first appearance of Psychodidae Newman, 1834, the family of the sandfly-related insect taxa is dated to the Upper Triassic period (Blagoderov et al. 2007), which coincides with the emergence of the modern ecosystems (Andrade Filho and Brazil 2003). The Late Triassic was also when mammals evolved from small, advanced cynodonts, the mammal-like reptile ancestors of the Mammaliaformes clade (Wallace et al. 2019). The parallel success of dinosaurs and other archosaur clades resulted in the extinction of most of the non-mammalian mammal-like reptiles at the end of the period (Benton 1983). Possibly, the earliest mammals hid from agile carnivorous reptiles in burrows since it is known that burrowing could have been common already among the end-Permian cynodonts (McLoughlin et al. 2020) and Late Triassic cynodonts created complex burrow systems (Benton 2021). These facts can be related to the evolution of sandflies, and Leishmania parasites, as several Old-World sandfly species are known to breed in burrows of rodents and hyraxes and feed on mammalian blood. The Jurassic, parallel to the break-up of the supercontinent of Pangea, could have been the time of the speciation of Psychodidae, including the differentiation of the tribes Hertigiini and Phlebotomini (Andrade Filho and Brazil 2003). Furthermore, based on fossil records, the early members of the Phlebotomus genus appeared in the middle part of the Cretaceous period (Kaddumi 2007;Stebner et al. 2015). The geological and climatic changes of the Tertiary period led to the rise of the ancestors of the present-day clades. The radiation of the oriental phlebotomine sandflies could have started about 50 mya in the Eocene (Ilango 2011). The rapid Oligo-Miocene climatic changes and the collision of Eurasia and the terranes of the former Neo-Tethys Ocean could have played a crucial role in determining the climatic requirements and the range of the emerging Mediterranean sandfly taxa. Among them, the late Neogene aridification of the Mediterranean and peri-Mediterranean region (Esseghir et al. 2000;Cruaud et al. 2021;Trájer et al. 2022a), the permanently changing coastlines of the former Paratethys Sea (Depaquit et al. 2002;Trájer 2022b), as well as the Messinian Salinity Crisis (Kasap et al. 2015;Trájer et al. 2021) are worth mentioning.
The members of the sandfly genus Transphlebotomus occupy large regions of the Mediterranean area in Europe, Northern Africa, and Western Asia. Among them, Phlebotomus (Transphlebotomus) mascittii Grassi, 1908, occurs in oceanic and continental climate regions of Central Europe such as Belgium (Depaquit et al. 2005), Germany (Naucke and Pesson 2000), Austria (Kniha et al. 2020) and Hungary (Farkas et al. 2011), but also the subtropical and Mediterranean climate regions of France, Italy, or Slovenia (Veronesi et al. 2007;Prudhomme et al. 2015;Praprotnik et al. 2019). Therefore, it is plausible that Ph. mascittii is the northernmost sandfly species in Europe and it is an assumed but unproven vector of Leishmania infantum Nicolle, 1908(Zanet et al. 2014, Obwaller et al. 2016. Surprisingly, the range of the other thermophilic Transphlebotomus species, Phlebotomus (Transphlebotomus) anatolicus Erisoz Kasap, Depaquit and Alten, 2015, Phlebotomus (Transphebotomus) canaaniticus Adler and Theodor, 1931, Phlebotomus (Transphlebotomus) economidesi Léger, Depaquit, and Ferté, 2000a, b, Phlebotomus killicki Dvořák, Votýpka and Volf, 2015 simonahalepae Cazan, Erisoz Kasap and Mihalca, 2021, are restricted to the relatively mild winter and hot summer climate regions of Southeast Europe (Ivović, et al. 2007), Asia Minor (e.g. Arserim et al. 2022), and the Levantine (e.g. Maroli et al. 2009) and do not occur in the oceanic and continental regions of Europe. Kasap et al. (2015) found that the first split of Transphlebotomus occurred around 10 million years ago. The authors suggested that the present-day occurrence and the climatic needs of species of the subgenus can be explained by the post-Miocene geographical climatic and geographical changes of the Aegean region. However, the main drivers of the diversification of the subgenus were not investigated in model environments.
The study aimed to test the Eastern Mediterranean speciation theory of the Mediterranean sandfly species using modelling tools. A special focus was given to the possible role of the Central Paratethys in the early split of these subgenera and the late Neogene-early Pleistocene zoogeographical changes of the ancestor of Phlebotomus mascittii and the other, non-Ph. mascittii Transphlebotomus species in Europe, Asia Minor, Levant, and North Africa.

Outline of the study
To investigate the speciation drivers of the sandfly subgenus Transphlebotomus in the Mediterranean region, the following tasks were performed: 1) The known occurrences of Ph. mascittii and the united distribution area of the other Transphlebotomus taxa of the eastern Mediterranean Basin were georeferenced. 2) Based on the georeferenced occurrence data and climate models, the climatic limits of the species were determined. The determined number and type of climatic constraints depended on the nature of the available palaeoclimatic data. Fort the heterogenic nature of the climatic data used in the study, the products of all model environments were presented in different subsections in the Results. 3) Tortonian and Messinian climatic models were created based on existing, site-like palaeoclimatic reconstructions. 4) Three model environments were established: two different Late Miocene -Tortonian and Messinian -model environments and a mid-Pliocene interglacial (MIS19), Last Interglacial Period and Last Glacial Maximum agerelated model environment. 5) Divergence time estimates of Transphlebotomus species were calculated to compare climatic models with potential diversification events of the member species.
The QGIS free and open-source geographic information system was used for data processing and modelling purposes (Lacaze et al. 2018).

Sandfly data sources
A total of 184 Transphlebotomus occurrences were involved in the study in Europe, Asia Minor, and the Levant. The occurrence data of Ph. mascittii (n sites =157) and the other non-Ph. masctitii Transphlebotomus occurrences (Σn sites =29) were based on the available data in the literature ( Fig. 1; red points). The other, non-Ph mascittii species were Ph. anatolicus (n site =2), Ph. canaaniticus (n site =6), Ph. economidesi (n site =5), Ph. killicki (n site =8), and Ph. simonahalepae (n site =1). Furthermore, five additional non-Ph. mascittii Transphlebotomus occurrences from Asia Minor were also georeferenced ( Fig. 1; yellow points). Table 1 shows the occurrence data sources of the studied Transphlebotomus species.
In 72% of the occurrences, the exact coordinates were known. In the remaining cases, heterogeneous methods were used to identify the location of the original trappings as accurately as possible. In some cases, the exact coordinates of the catch sites could be reconstructed. For example, the coordinate determination of the occurrence of Ph. economidesi in Mandria cave, Cyprus, was based on the coordinates of the entrance of the cave (Kasap et al. 2015); in Monteggio and Sessa, southern Switzerland, Grimm et al. (1993) provided occurrence data on a small-scale schematic map with the pattern of respective streets and even houses. By identifying the geographical position of the streets and individual houses within the small settlements, the coordinates extraction was possible within a resolution accuracy of ten meters. In the case of other small settlements, where a schematic map was not provided, the centre coordinates of the villages were considered in georeferencing. For example, Mazeris et al. (2010) collected Ph. economidesi in Dora, Cyprus, a village smaller than 1x1 km.
Sometimes, the resolution was between 1x1< to ≤10x10 km. For example, Léger et al. (2000a) found Ph. economidesi in Armenochori, Cyprus. It is also a very small village, with a maximum extension of 5x2 km. Therefore, it can be concluded that the used coordinates were accurate or sufficiently accurate. For example, in the case of Ph. mascittii occurrence sites, it means that of all occurrence data points, 102 were exact point coordinates, 44 were at a scale of 1x1 km or lower, and eleven were at a scale of 10x10 km or lower. From the point of view of distribution modelling, even the accuracy of ≤10x10 km resolution is generally acceptable because these model experiments aim to predict large-scale occurrence patterns and not local habitats. Furthermore, the cut of 2-2 percentiles from the sampled climatic data was utilised to filter the possibly appearing non-relevant climatic values. Finally, it is worth mentioning that the resolution of the used current period's climate models, which were utilised for climatic data extraction, is 2.5 arc-minutes. This resolution in the mid-latitudes is equal to ~5x5 km. Therefore, it means that the accuracy of the site coordinates fits the resolution of the reference period's climatic models used in the study.

Climatic and topographic data sources
The Tortonian and the Messinian climate data were based on the palaeoclimatic reconstructions of Bruch et al. (2006Bruch et al. ( , 2011. The authors applied the Coexistence Approach of Mosbrugger and Utescher (1997) using the fossil plant assemblages to estimate the former palaeoclimatic conditions of the former West and Central Paratethys, as well as the Aegean area and Asia Minor. The Coexistence Approach method or concept is based on the current distribution arealimiting values of the closest extant relatives of the fossil plant species. This method can be used mainly for the Cenozoic period, but due to the uncertainties of the similarity between the living and the extinct floral elements, the method requires robust statistical tools (e.g. Mosbrugger and Utescher 1997;Erdei et al. 2007). Bruch et al. (2006) produced the values of four bioclimatic variables. Table 2 shows the modelled periods and the related data sources.
To create continuous, georeferenced climatic data of the Mediterranean region in the wider sense, the method of Trájer (2022a, b, c) was used who utilised the following steps to produce the palaeoclimatic maps: 1) The sites of the palaeoclimatic reconstructions were georeferenced in QGIS software, and the site-related values were added to the sites in the attribute table. Not the absolute values, but the differences between the palaeoclimatic and the modern values were applied. 2) The differences between the palaeoclimatic and the modern values were interpolated using the IDW interpolation function of the QGIS software.
3) The modern climatic values were modified by the interpolated difference values. 4) Finally, the former sea level was reconstructed. It should be noted that this study did not aim to create exact geological or palaeogeographical reconstructions. Figure 2 shows the interpolated difference values compared to the reference periods of the Tortonian annual precipitation sum (bio12), annual mean temperature (bio1), January T m 01), and July mean temperature (T m 07) values. The Tortonian values were produced by modifying the present-day values by the IDW interpolated differences. Figure 3 shows the interpolated difference values compared to the reference periods of the Messinian temperature annual range (bio5-bio6; bio7), annual precipitation (bio12) and July precipitation sum (P07) values. The Messinian values were produced by modifying the present-day values by the IDW interpolated differences. As the potential Messinian range of the non-Ph. mascittii Transphlebotomus species were observed to cover former coastal regions, and the palaeoclimatic conditions of the formerly desiccated Mediterranean Basin's floor were also reconstructed following the method of Trájer et al. (2021).
In this process, the steps were as follows in the case of July and the Annual Precipitation Sum values: 1) Several coastal points were designated along the entire Mediterranean basin.  The sea surface temperature on the formerly desiccated abyssal plain of the Mediterranean Sea was calculated  . 2 The interpolated difference values compared to the reference periods of the Tortonian annual precipitation sum, annual mean temperature, January, and July mean temperature values according to the dry adiabatic lapse rate. The global sea level was 25 m higher than today (Dwyer and Chandler 2009). The thermal patterns of the Mediterranean Basin can be calculated according to the following formula during the Messinian Salinity Crisis and the dry adiabatic lapse rate (Ritter 2006): where T agl is the abyssal ground-level temperature (°C) of the Messinian seafloor, Δh is the altitude below sea level (m), and T rsl is the reconstructed Messinian temperature at sea level. Figure 4. shows the thermal surplus of the desiccated Mediterranean seafloor in the Messinian period.
As the topographic model, the ETOPO1 1 Arc-Minute Global Relief Model of the Earth's surface was utilised, integrating land topography and ocean bathymetry (Amante and Eakins 2009). In the case of the mid-Pliocene and the MIS19 Quaternary period's model, existing palaeoclimatic models were used.
Model environments and the determination of the extrema As mentioned, the available data sources made different numbers and types of climatic extrema possible. Model 1 environment was attributed to the creation of the Tortonian models; the mid-Pliocene models and three different, MIS19, Last Interglacial Period and Last Glacial Maximum's Quaternary period-related models were performed in model environment 2, and model environment 3 was the source of the Messinian period's models. The climatic values used in the different model environments and determined climatic extrema can be seen in Table 3.
For acquiring the range-limiting extrema of Ph. mascittii and the united group of the other East Mediterranean Transphlebotomus species, the 2.5 arcmin resolution model of Karger et al. (2017) related to the reference period of 1979-2013 was used. As was previously mentioned, 2-2 percentiles were cut from the absolute maximum and minimum values of the factors to avoid the involvement of unrealistic climatic contracts according to the generally applied considerations of environmental modelling (Trájer et al. 2013). Table 4 shows the determined climatic distribution-limiting extrema.

Modelling the former sandfly ranges
The modelling of the former distribution areas followed the logic of Boolean algebra. The area of all range areas was modelled according to the following general equations: Where v n represents the n th climatic constraint of the distribution area of a species, v n_limit_min and v n_limit_max are the T agl = 6.5 × 10 3 × (Δh + 25) + T srl The interpolated difference values of the climatic factors, compared to the reference periods of the Messinian temperature annual range, annual precipitation and July precipitation sum values lower and upper distribution-limiting values related to the climatic constraint. The potential area-based suitability patterns were determined according to the following mathematical formalism: where A(v 1 ; v2…v n ) shows the potential distribution area of the given species, which contains the remaining areas after taking into consideration the factor-related limitations.
Then, the modelled values were transformed into percentage (%) values and were colourized.

The potential mid-Pliocene distribution areas
The potential mid-Pliocene cold period's range of a Ph. mascittii-like species based on the former climatic suitability values is focused predominantly on Western Europe, including, e.g. the territory of the present-day British Islands, the Dogger Bank region, and South France. The potential mid-Pliocene cold period's range of a non-Ph. mascittii-like Transphlebotomus sandfly taxa could have shown the highest affinity to the area of present-day Eastern England, France, the Iberian Peninsula, Asia Minor, the South Apennine Peninsula, and certain regions of North Africa (Fig. 8).
In the mid-Pliocene warm period, the potential range of a Ph. mascittii-like species could have covered large regions of West and Central Europe, including, e.g. the northern and central regions of present-day France, Central Germany, the southwest part of Scandinavia, the Transdanubian part of Hungary and Transylvania. In contrast, the potential mid-Pliocene warm period's climatic suitability values of a non-Ph. mascittii-like Transphlebotomus sandfly taxa show a markedly different picture, being restricted to certain Mediterranean regions of Europe, including certain parts of the Iberian, the Apennine and the Balkan Peninsulas, Asia Minor, the Mediterranean islands like Corsica, Sardinia and Sicily and certain smaller Northwest African and Levantine regions (Fig. 9).

The potential Pleistocene distribution ranges
The modelled MIS19 period's potential range of a Ph. mascittii-like species is somewhat like the present-day occurrence of the species, which covers large regions of West and Central Europe (Fig. 10a). However, the model suggests that in this period, the climatic suitability of the Mediterranean peninsulas could be lower for the species than today. In the Last Interglacial Period, a Ph. mascittiilike species could occur in Western Europe, including the wider coastal regions of the Bay of Biscay, North France and the southern area of the Benelux states, the milder regions of the British Isles and maybe some mountainous regions of the Apennine Peninsula (Fig. 10b). During the Last Glacial Maximum, the climatically suitable areas for a Ph. mascittii-like species could be restricted to the narrow, currently partly submerged shelf and coastal regions of the Mediterranean Sea including the territories around  (Fig. 10c).
As of today, in the MIS19 period, the potential range of a non-Ph. mascittii-like Transphlebotomus sandfly taxa could cover the Mediterranean territories, including large regions of the South Balkan, the Apennine, and the Iberian Peninsulas, as well as some North African and Levantine territories and the western and central part of present-day France (Fig. 10d). In the Last Interglacial Period, a non-Ph. mascittii-like Transphlebotomus species could inhabit large regions of present-day southern Europe and Western Europe (Fig. 10e). For the Last Glacial Maximum, the potential area of a non-Ph. mascittii-like Transphlebotomus species could be restricted to South Europe (Fig. 10f).

Divergence time estimates
Based on the rejection of the null hypothesis of the equal evolutionary rate at a 5% significance level (p<0.0001), BEAST analysis was conducted under a lognormal, uncorrelated relaxed clock model. The determined tree topology divides the studied Transphlebotomus species into two major clades. The first clade consists of Ph. anatolicus and Ph. canaaniticus. Phlebotomus mascittii, Ph. killicki, Ph. economidesi and Ph. simonahalepae form the second clade. In the second clade, Ph. mascittii is the sister taxon of the other sandfly species. Divergence time estimates of the six Transphlebotomus species ranged from 5.4 to 10.9 mya. The divergence between clade I and clade II could happen ca. 10.9 mya (95% HPD interval = 7-16). Within clade I, the estimated split of Ph. anatolicus and Ph. canaaniticus happened 5.6 mya (95% HPD interval=3.2-8.8). Within clade II, the first split between Ph. mascittii and the other lineages (Ph. killicki, Ph. economidesi, and Ph. simonahalepae) was estimated at approximately 9.5 mya (95% HPD interval = 5.9-14.3), the second split between Ph. killicki and Ph. economidesi + Ph. simonahalepae at around 7.4 mya (95% HPD interval = 4.4-11.2), and the last split between Ph. economidesi and Ph. simonahalepae at approximately 5.4 mya (95% HPD interval = 3.1-8.3) (Fig. 11).

Discussion
The model results suggest that if the ancient Transphlebotomus populations had similar climatic tolerance to their extant offspring, their Late Miocene range could have covered both the North and the South sides of the Central and east Paratethys Seas. Several authors claimed that the Paratethys realm played a very important role in the Neogene evolution of the Mediterranean sandfly species and the speciation of the present-day clades (Steininger and Rogl 1984;Léger and Pesson 1987;Marchais 1992;Esseghir et al. 1997;Depaquit et al. 1998Depaquit et al. , 2002Trájer 2022a, b). The Paratethys Sea became a semi-locked sea approximately 30-35 Ma during the lower Oligocene epoch when it had several narrow connections to the Atlantic and the Mediterranean Seas via, e.g.
the Rhine Graben, the Rhone Strait, and the Slovenian Strait (Palcu 2018). Later, due to the rise of the Alpine and Dinaric Mountain systems, the Paratethys lost its direct connections to the epicontinental seas. After the close of the Slovenian Strait in the early Serravallian, it became an intermittently isolated or semi-isolated sea (Palcu et al. 2015). The connection between the Paratethys Sea and the global seas was reestablished due to the formation of the Mid-Aegean trench ca. 12 mya (Popov et al. 2006), which also formed a significant barrier against the migration of invertebrate species (Simaiakis and Mylonas 2008;Papadopoulou et al. 2010). From this era, the Central and the East Paratethys partly divided the Eastern European and Southwest Asian fauna.
However, at the turn of the late Middle -early Late Miocene epochs, the isolation only became partial in the westernmost part of the Central Paratethys area since a narrow land bridge was formed between the Balkans and the Alps due to the rise of the mountain systems (Kováč et al. 2017). Based on the phylogenetic results of this study, it is very plausible that this geographical change led to the north-westward spread of ancient Transphlebotomus populations and the speciation of Ph. mascittii. With the formation of the Aegean trench, these large-scale geographical processes could also explain the appearance of the Ph. killicki-Ph. economidesi-Ph. simonahalepae clade in the Aegean area. It should be noted that the phylogenetic results of the study are following the ML analysis-supported previously published cytb tree topologies (Kasap et al. 2015;Cazan Fig. 11 Phylogenetic analysis of cytb from Ph. anatolicus, Ph. canaaniticus, Ph. economidesi, Ph. killicki, Ph. mascittii, Ph. simonahalepae, and Ph. chinensis  The calculated appearance times suggest that both the divergences of the Ph. anatolicus-Ph. canaaniticus (5.59 mya) and the Ph. economidesi-Ph. simonahalepae (5.36 mya) approximately coincides with the Messinian Salinity Crisis (~5.96 to 5.33 mya; Gautier et al. 1994). The possible role of the Messinian Salinity Crisis on the speciation of sandfly species was suggested by many authors (Kasap et al. 2015;Cruaud et al. 2021;Trájer et al. 2021;Pavlou et al. 2022).
The distribution of terrestrial mammals strongly supports zoogeographical and land connections in the Mediterranean between the mainland areas and the Mediterranean islands (Azzaroli and Guazzone 1979). It is known that after backstripping the Pliocene-Quaternary and the evaporite layer, the Pre-Messinian seafloor was generally 1800-1900 m deep (Netzeband et al. 2006) However, notably deeper areas also could have existed, e.g. in the Hellenic Trench (the current deepest point: Calypso Deep, 5,110 m u.s.l.). Considering the physiological-thermal tolerance of sandflies, it seems problematic because it is hard to explain how the relatively sensitive sandflies survived and spread in the bottom of the extremely hot Mediterranean basin . For comparison, on the coast of the present-day Dead Sea, the mean annual temperature is 24°C, the average temperature of the warmest quarter is 32°C, and the annual precipitation is about 300 mm (Harris et al. 2014). On the Gaza coasts, around 100 km to the west of the Dead Sea coasts, at sea level, these values are about 21°C, 27°C and 480 mm. Differences in the climatic conditions exist at the approximately -435 meters elevation difference between the sea-level elevation sites of the Eastern Mediterranean Sea at Gaza and the coastal regions of the salt lake.
The desiccated Mediterranean Basin generally could have been more than four times deeper than the present-day Dead Sea or even deeper. Since -as was mentioned -the adiabatic temperature change is 0.65°C per 100 m (Ritter 2006), it can be calculated that the annual mean temperature of the Mediterranean seafloor could have been at least 12°C higher than the former 0 m level territories which were 25 m higher than today (Dwyer and Chandler 2009). These harsh conditions could have posed unsolvable problems for ancient sandfly species even if we know that they can live in certain hot, semi-arid regions like the Negev in Israel (Schlein et al. 1984) or Central Tunisia (Barhoumi et al. 2016). On the other hand, large areas of the seafloor -including the Algerian, the Balearic, the Ionian, and the Levantine subbasins -could have been covered by hypersaline residual lakes according to the occurrence of large covered evaporitic formations (Rouchy et al. 2006). The most probable free migration route could have existed between present-day Tunisia, Sicily, and Calabria.
Based on the ideas, and results discussed above, the following hypotheses can be set related to the Late Miocene Late Pleistocene speciation and migration events of the Mediterranean Transphlebotomus species: It must be acknowledged that apart from only one record, Ph. simonahalepae might be cavernicolous (Cazan et al. 2021), which means that it might have special requirements (climatically and ecologically in general). It also explains that the climate of the mid-Pliocene warm period seems unsuitable for Transphlebotomus in Romania and Bulgaria; however, Ph. simonahalepae could have survived in caves, a more stable environment. Caves can provide adequate shelter for extrazonal sandfly populations to survive in continental areas. For example, Phlebotomus (Larroussius) perfiliewi Parrot, 1930 from Esztergom, North Central Transdanubia, Hungary and the Gellért Hill cave, Budapest, were found in Hungary in the 1960s (Tánczos 2012;Trájer 2017) and Ph. neglectus uses the fissures of an abandoned quarry in Southern Transdanubia, Hungary, as shelters (Trájer et al. 2018a) plausibly to survive the cool winters of the area. As another example, Ph. mascittii is active year-round in Corsica as it lives in a cave-like habitat (tunnel) with stable climatic conditions all year (Naucke et al. 2008). As was mentioned, the presence of Ph. mascittii in Northern Africa raises questions. At first, it can be mentioned that the species is only known from one localization from the area (Berdjane-Brouk et al. 2011).
It is not surprising that Mediterranean sandflies, which occur in Southern Europe, can also be found in Northern Africa. For example, Phlebotomus (Paraphlebotomus) alexandri Sinton, 1928, Phlebotomus (Larroussius) perniciosus Newstead, 1911, Phlebotomus (Phlebotomus) papatasi (Scopoli, 1786, Phlebotomus (Larroussius) ariasi Tonnoir, 1921, Ph. perfiliewi and Phlebotomus (Paraphlebotomus) sergenti Parrot, 1917 can be found on both, the European and the African sides of the Mediterranean Sea (ECDC 2022). Among them, Ph. perniciosus and Ph. ariasi have a somewhat similar distribution in the West Mediterranean region in that sense that they are absent in the East Mediterranean Basin but occur on the West Mediterranean coastlines both in Africa and Europe. This can be a consequence of a similar evolutional and palaeozoogeographical history. It is plausible that before the Messinian epoch, Ph. ariasi and Ph. perniciosus, which together form a so-called 'West-Mediterranean group' were adapted to the more humid and balanced climate of Southwestern Europe (Trájer et al. 2018b).
It is plausible that Ph. mascittii were also adapted to somewhat similar conditions in the Late Miocene era in the West Mediterranean Basin. The presence of the species in Corsica (Naucke et al. 2008), Sardinia (Biocca et al. 1977), Sicily (Lisi et al. 2014), and Montecristo (Zanet et al. 2014) raises the question of when and how this species dispersed to these islands. The cooler or humid climatic events of the Messinian stage might have formed a dried-up corridor for migration events; however, this was not shown in the models for average climatic conditions. In addition, the close genetic distances (cytb) of Ph. mascittii from Corsica and mainland Europe suggest a more recent dispersal to the island rather than during the Messinian (Kniha et al. 2020). Model results show that during the humid and tempered climate of the Late Pliocene Glacial Event (Marine Isotope Stage M2) in Europe, Ph. ariasi and Ph. perniciosus-like ancestral sandfly species could have covered large regions of Europe ca. 3.3 mya, the large potential range area of which collapsed suddenly at the start of the warm mid-Pliocene warm period 3.205 mya (Trájer 2020). The mid-Pliocene climatic fluctuations could also have caused notable stress for ancient populations of Ph. mascittii which can explain why isolated populations remained in North Africa. However, it is questionable whether the later glacial-interglacial climatic fluctuations could support the survival of the North African populations. It should be noted that in the case of island populations, late anthropogenic transport also cannot be excluded from the Neolithic period. However, natural transport media like trees floating on the water surface could also have promoted the dispersal of Mediterranean sandflies (Trájer 2021).
It is striking that since Ph. mascittii diverged from the Ph. killicki-Ph. economidesi-Ph. simonahalepae clade as early as ca. 9.5 mya in the Tortonian stage of the Late Miocene, and the model results suggest that the climatic conditions of the warmer periods of the Neogene theoretically could make it possible for the eastward migration of the species to the Balkans, it has never been recorded in continental Greece. Phlebotomus mascittii plausibly formed from the division of the ancestral Transphlebotomus population due to the separation effect of geographic barriers. The barrier originally could be the Central Paratethys, and later, it was formed by the emerging Dinarids and the Alpine Mountain ranges. The investigation of the Late Miocene mammal faunas (e.g. Ataabadi et al. 2016) and palaeovegetation models (e.g. Pound et al. 2011) revealed that in the Tortonian, two major biomes existed in Western Eurasia: the zone of dense, humid, subtropical temperate forests which flourished predominantly in Central Europe and the Apennine Peninsula, and a woodland-shrubland biome zone with relatively open vegetation, which was characteristic to the South Balkan and Anatolia. At that time, temperate forests were present above 60°N, and warm-temperate mixed forests covered much of Europe and South-East Asia (Pound et al. 2011). As described above, the ancestors of Ph. mascittii plausibly were adapted to this environment.
In the latest Miocene, the increasing aridity caused the decline of subtropical forests thorough the mid-latitude regions of Eurasia (Nelson 2021;Habinger et al. 2022). This event could play a central role in the Late Miocene diversification of Transphlebotomus species, establishing the differences between the climatic requirements of the ancestors of Ph. mascittii and the other Transphlebotomus species. Although in the first part of the Pliocene period, the climate of such Central European regions like the Carpathian Basin was still mild and humid, supporting mesophytic forest (Erdei et al. 2007), which could be a suitable habitat for a Ph. mascitti-like sandfly species (see Fig. 9a). However, the Late and mid-Pliocene cooling events created relatively cool and dry climatic conditions in Central Europe, which narrowed the occurrence of the ancestral Ph. mascittii populations to South-Eastern Europe and Anatolia; and the non-Ph. mascittii populations to the oceanic climate-influenced regions of Western Europe (Fig. 8b). For the MIS19 interglacial, a wide range gap could exist between the ancestors of Ph. mascittii and the other Transphlebotomus species, which became wider for the Late Pleistocene (Fig. 10). This distribution gap and the different climatic adaptation of Ph. mascitii and its relatives can probably explain the distribution differences seen today, including the distribution patterns of Ph. mascitii and other Transphlebotomus taxa in the Balkans.
The late Quaternary models suggest that even in such interglacial periods as the Last Interglacial, a Ph. mascittii-like sandfly species could not reach its current, wide European range, which covers Central Europe. Instead, the potential area of the species covered the present-day milder oceanic climate regions of Western Europe. The glacial maximums could lead primarily to the serious habitat loss of Ph. mascittii, the other Transphlebotomus populations could find relatively large refugia in South Europe and Anatolia. The more expressed vulnerability of a Ph. mascittii-like to the period can be explained by the fact that the extant species dominantly occupy temperate seasonal forest biome-covered regions in Europe (compare the Ph. mascittii occurrences in Fig. 1 with the present-day potential natural vegetation map of Hickler et al. 2012). However, the temperate mixed broad-leaved forests and thermophilous mixed broad-leaved forests were significantly shrunk and retreated to Southern Europe during the Last Glacial Maximum (Svenning et al. 2008;Trájer 2022a). The glacial retreat of the humid, mild-climate forest habitats should have caused a significant genetic bottleneck effect in the case of the ancestral Ph. mascittii and other Transphlebotomus populations. It is worth mentioning that similar glacial refugia, as shown in Fig. 10c and f can also be hypothesised in the case of other European sandfly taxa in Southern Europe, including the Iberian Peninsula, territories north of the Pyrenees, the coastal regions of the Apennine Peninsula, the South Balkan, and the Aegean islands and Western Anatolia (Aransay et al. 2003;Mahamdallie et al. 2011;Depaquit et al. 2015;Trájer and Sebestyén 2019). Finally, it can be added that further phylogenetic studies are needed to elucidate where glacial refugia existed in southern Europe and how postglacial warming impacted the northward migration of the ancestral populations of the extant Transphlebotomus species.

Conclusions
Transphlebotomus species could have emerged in the Late Miocene in the Aegean-Asia Minor area. In the Tortonian stage, certain Transphlebotomus populations may have crossed Central Paratethys over the Dinaric land bridge, and this population became the founder of the later Ph. mascittii species. In parallel, the formation of the Aegean trench could have led to the emergence of other Transphlebotomus species. These geological processes dissected the ancestral Transphlebotomus populations into a European, Balkan and Asia Minor-Levantine group. Due to the desiccation of the Mediterranean Sea in the Messinian stage, the ancestors of Ph. mascittii perhaps could have migrated to Northern Africa during the cooler and more humid periods. However, the model results, which were based on the general climatic conditions, do not support this hypothesis. The Messinian Salinity Crisis triggered the speciation of several East Mediterranean Transphlebotomus taxa.
Acknowledgements The research presented in this article was carried out with the support of the NKFIH-471-3/2021 project. We also acknowledge the critical notices and useful suggestions to the study reviewers, namely Jérôme Depaquit and the second reviewer, who wishes to remain anonymous.
Funding Open access funding provided by University of Pannonia.

Data availability statement
The datasets generated during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.