Investigation of the possible role of the Central Paratethys as a migration route and speciation area of the ancestors of Mediterranean Larroussius, Paraphlebotomus and Phlebotomus species

The Oligocene and Miocene tectonic and biogeographical alterations of the peri-Mediterranean area could strongly impact the speciation processes and migrations of Mediterranean sandfly species. To understand the possible former role of this palaeobiogeographic factor on ancestral sandfly species, the potential suitability values of five Larroussius, two Paraphlebotomus and one Phlebotomus species were modelled from the Rupelian to the Tortonian stage in the Central Paratethys. The bioclimatic extrema of this sandfly species and the Coexistence Approach-based palaeoclimatic reconstructions made the basis of hypothesis testing. In the case of the Tortonian stage, a georeferenced climatic model was created. The models indicate that the suitability values could notably vary by species and periods. The monsoon-influenced humid subtropical climate of the Oligocene epoch could be less suitable for the ancestors of Mediterranean sandfly taxa than the later, drier humid subtropical climates-characterised Miocene stages. The Central Paratethys area could be less suitable for ancestors of the extant Paraphlebotomus, and Phlebotomus in the Miocene epoch compared to Larroussius species. It may indicate that the Central Paratethys formed a barrier against the east-to-west migration of the ancestors of Phlebotomus and Paraphlebotomus species. In contrast, Larroussius species could prefer the Miocene climate of the Central Paratethys. In the Tortonian stage, mainly the coastal areas of the sea could be colonised by sandflies. These results indicate that the coasts of the Central Paratethys should be not only considered as a potential former migration route but could be a part of the speciation area of Larroussius sandflies.


Introduction
Leishmaniasis is a parasitic disease caused by different Leishmania species and transmitted by phlebotomine sandflies in the Old World. About 2% of the global cases of leishmaniasis occur in the WHO European Region (EUR). However, it is plausible that the real burden on leishmaniasis is underestimated. Based on the reported cases, cutaneous (CL) and visceral leishmaniasis (VL) forms are present in 22-27 of the 53 member states of the EUR (WHO 2018). For example, in 2016, 2862 autochthonous CL and 353 VL cases were reported in the region. About 80% of the CL cases occur in Israel, Turkey, Turkmenistan, and Uzbekistan in the WHO European Region, and VL is the most common in Albania, Georgia, Italy, and Spain (WHO fact sheet leishmaniasis). In 2014-2016, the highest annual case numbers of autochthonous CL cases were reported in Turkey, Uzbekistan, Israel, and Spain. The highest VL incidences were reported in Spain, Georgia, Italy, and Uzbekistan (WHO 2018). Leishmaniasis emerges in Europe and the adjacent regions due to climate change (Dujardin et al. 2008;Ready 2008Ready , 2010. Model projections indicate that since global warming triggers the northward spread of the sandfly vectors, the leishmaniasis-endemic areas in Europe will increase Koch et al. 2017;Chalghaf et al. 2018). Dog travelling (Maia and Cardoso 2015) and HIV coinfection are also well-known factors of the spread and the increase in the prevalence of leishmaniasis in Europe (Monge-Maillo et al. 2014).
The coevolution of sandflies with Leishmania species of mammals and lizards (Ready 2013) and the effects of biogeographical changes on the environmental requirements, as well as the distribution of extant sandfly taxa is still in the focus of sandfly-related ecological and evolutionary investigations (Ready 2000;Cruaud et al. 2021). The Afrotropical-Palearctic-Indomalayan distribution Phlebotomus genus (Diptera: Psychodidae) (Akhoundi et al. 2016), represents a phylogenetically old clade of Diptera. It might appear at the beginning of the Late Cretaceous epoch. The oldest, 94-100 Myr old fossil species, Phlebotomus vetus Stebner et al. 2015, was found in Cenomanian-age amber in Myanmar (Stebner et al. 2015). Based on the fossil record and the present distribution of sandflies, it is likely that the family is native to the humid, tropical palaeoregions (Akhoundi et al. 2016). In contrast to, e.g. mosquitoes and many other members of the family, like Clogmia albipunctata Williston, 1893 (Trájer and Juhász 2017), phlebotomine sandflies generally do not breed directly in waters, but lay their eggs on such wet surfaces as microbial mats or damp decaying materials (Moncaz et al. 2012;Killick-Kendrick 1999). Sandfly larvae develop in moist, windprotected and shaded environments such as the abandoned wells covered with plants, the humid understory environment of tropical herbs-formed bushes, tree holes, the wet surfaces of old stacks of bricks covered with plants, organically rich moist soils of rain forest floors, the wet and faecal-contaminated soils of animal shelters like rodent burrows or cattle huts, the algal math-covered above-ground roots of trees, the wet crevices of caves, rocks and buildings (Feliciangeli 2004;Dinesh et al. 2009;Alencar et al. 2011;Rangel and Shaw 2018;Wijerathna and Gunathilaka 2020).
The optimal temperature for the development of the Mediterranean sandfly larvae is relatively high, maybe referring to their speciation area. The larval optimum temperature is 25°C in the case of Phlebotomus (Larroussius) neglectus Tonnoir, 1921, Phlebotomus (Larroussius) perfiliewi Parrot, 1930 and Phlebotomus (Larroussius) perniciosus Newstead, 1911, 28-34 in the case of Phlebotomus (Phlebotomus) papatasi (Scopoli, 1786) and 31-33°C is the optimum temperature for the larval development of Phlebotomus (Paraphlebotomus) sergenti Parrot, 1917 (Lindgren et al. 2006;Killick-Kendrick 1999). To understand the significance of this data, it is necessary to know that July daily mean temperature in such present-day Mediterranean cities like, e.g. Athens, Málaga, Marseilles, Palermo, and Rome is between 24.2-29.3°C values (KNMI Climate Explorer data, Trouet and Van Oldenborgh 2013). It indicates that nowadays the climate of Southern Europe is sub-optimal for the larval development of the species of subgenera Phlebotomus and Paraphlebotomus, but optimal for Larroussius species. The above-described differences between the optimal larval development temperatures of the different subgenera can be interpreted from a phylogenetic point of view based on the results of Grace-Lema et al. (2015). Based on their results, it can be hypothesised that while the ancestors of Paraphlebotomus and Phlebotomus species might have specialised under hot and semi-arid climatic conditions, the ancestors of Larroussius species could have a subtropical origin. The preferred relative air humidity values of adults also fit this hypothesis, since it is 60-80% in the case of the Larroussius species Ph. neglectus and Ph. perfiliewi, but it is only 36-45% in the case of the Phlebotomus species Ph. papatasi and less than 45% in the case of the Paraphlebotomus species Ph. sergenti (Lindgren et al. 2006;Killick-Kendrick 1999).
Several studies were performed to understand better the Neogene-Quaternary migrations and evolutionary history of the Mediterranean sandfly taxa, mainly on the phylogenetic basis (e.g. Depaquit et al. 1998Depaquit et al. , 2000Depaquit et al. , 2002Esseghir et al. 2000;Alten 2010;Mahamdallie et al. 2011;Aransay et al. 2003;Kasap et al. 2015) and some study was performed also to study the environmental aspect of this problem using modelling tools (e.g. Trájer et al. 2018;Trájer and Sebestyén 2019). All these studies emphasise the importance of the co-evaluation of the climatic and the palaeozoogeographic factors in the investigation of former environmental alterations on the past and present sandfly occurrences. The major problem in the confirmation of the purely phylogenetic-based studies is that while the genetic results provide hard evidence for the relationships and possible divergence times of sandfly species, it is difficult to put the phylogenetic results in an adequate past environmental context due to the sparse fossil record. Unfortunately, the fossil record of Phlebotomus is too sparse to model the former occurrence of the species in Europe. Phlebotomus fossils are known from the Cretaceous amber of Myanmar, as well as subfossil specimens from the Holocene copal of Tanzania (Ilango 2011;Ross 2018).
Although Paraphlebotomus is only one subgenus of phlebotomine sandflies, it plastically depicts the speciation and migration issues that arise in connection with the Miocene epoch's climatic and geographical history of the Mediterranean Basin and the evolutionary history of the related sandfly subgenera. Depaquit et al. (1998Depaquit et al. ( , 2000Depaquit et al. ( , 2002 emphasised the importance of the Paratethys and Tethys Seas as primarily geographical factors and barriers in the speciation of the members of Paraphlebotomus, and we have no rational reason to assume that this would not be the case for other European sandfly subgenera as well. Depaquit et al. (1998Depaquit et al. ( , 2000Depaquit et al. ( , 2002 mainly investigated the role of the Eastern Paratethys Sea; Trájer et al. (2018A) raised the possibility that the Central Paratethys also could play an important factor in the emergence of the subgenus Paraphlebotomus, based on the observations of Depaquit et al. (1998Depaquit et al. ( , 2000Depaquit et al. ( , 2002 because the Central Paratethys could form with the Adriatic Sea an important barrier against the migration of western and (south) eastern European sandfly taxa.
The Adriatic Sea could constitute an important and probably insurmountable east-west and the west-east limit for the migration of Mediterranean sandfly species in the late Neogene period. In the subgenus Larroussius: Phlebotomus (Larroussius) ariasi Tonnoir, 1921 and Ph. neglectus which could be vicariant species now possibly have an overlapping area. It can also be observed in the case of Ph. perniciosus and Phlebotomus (Larroussius) tobbi Adler and Theodor 1930 (VECTORNET 2019). The Adriatic Sea is also the western limit of the Phlebotomus subg. Adlerius Nitzulescu, 1931 and the western limit of Sergentomyia (Sergentomyia) dentata (Sinton 1933) in Europe. However, it would be a mistake to assume that the Adriatic Sea alone was the most important geographical barrier in the past. From the late Eocene to the late Miocene epoch, the actual extension of the Adriatic Sea and the former Central Paratethys Sea together determined the migration pathways between the west and the eastern Mediterranean Basin regions. However, the geographical barrier function of the Central Paratethys always altered. This sea can be characterised as a changing efficiency barrier depending on the period, in contrast to the Adriatic Sea which could form an absolute barrier for sandfly taxa. Around the boundary of the Eocene and Oligocene epochs, the northward tectonic shift of Africa resulted in the subduction and the final disintegration of the former Western Tethys into three major sub-basins and its partial separation from the early Mediterranean Sea and the Indo-Pacific Ocean (Piller et al. 2007). The Central Paratethys Sea was the middle member of the Paratethys basins in the geographical sense. The Central Paratethys Sea which was bordered by the mountain ranges of the Alps, the Carpathians and the Dinarides, covered large areas in the present-day Pannonian Basin, the foredeep basins of the Carpathians, the Transylvanian Basin, and some parts of the Moravian Basin during the Oligocene and Miocene epochs (Kováč et al. 2007). It can be concluded that the Central Paratethys area was the gate for the west to (south) east and from (south) east to and west migrations (Fig. 1).
The Miocene Mediterranean palaeogeography and palaeoclimate of the Paratethys Realm could strongly influence the present-day distribution of Circum-Mediterranean entomofauna (Oosterbroek and Arntzen 1992;Liebherr 1994;Gaudant 2003;Casale et al. 2009;Wappler et al. 2009), including the members of the insect family Psychodidae. The investigation of multiple genetic divergences showed that the former population expansions and contractions of Mediterranean sandflies can be corresponded with past geographical and climatic changes (Mahamdallie et al. 2011). Several authors have claimed that the environmental conditions of the former coastal areas of the Paratethys Sea could play an important role in the Miocene diversification of the Mediterranean sandfly fauna (Depaquit et al. 1998(Depaquit et al. , 2000(Depaquit et al. , 2002Esseghir et al. 1997Esseghir et al. , 2000. The distribution-limiting climatic requirements of European sandfly species also confirm the assumption that their present habitat preference can be a consequence of the Miocene speciation of the ancient sandfly populations around the Paratethys Realm (Trájer et al. 2018). The speciation and dispersion of the Mediterranean Phlebotomine sandflies might be impacted by such serious Miocene regional ecological catastrophes as the Messinian salinity crisis (Alten 2010) or the Sarmatian isolation and sea-level fluctuations of the Central Paratethys (Palcu et al. 2015). These former environmental events can explain the present intraspecific variability of sandflies (Dvořák 2008).
The Central Paratethys not only included the Pannonian Basin System, but the Dinaridic fresh-water basins also belonged to the Central Paratethys Realm (Jamičić 2002) which could explain how the genus Congeriathe enigmatic shell genus of the Pannonian Lakesurvived the final retreat of the sea in the karst systems of the Northwest Dinarids (Stepien et al. 2001). About 12 Ma, due to extensional backarc tectonic deformation, the Pannonian Basin began to sink rapidly (Magyar et al. 1999;Matenco and Radivojević 2012). With the parallel uplift of the bordering mountain chain, the new back-arc basin reached its maximal extent of about 10 Ma. For this time, Central Paratethys became a brackish or almost freshwater lake (Pannonian Lake or Sea; Csató 1993;Matenco et al. 2016). The final retreat of the Central Paratethys Sea occurred in the Messinian (Krijgsman et al. 2010). Although a remnant lake, Slavonian Lake still existed until 4.5 Ma (Kázmér 1990). After the disappearance of the Central Paratethys, the endemic Diptera species were able to move to colonise the lower altitude regions of the Carpathian Range (Keresztes et al. 2012).
Based on the palinspastic reconstructions of the Central Paratethys (Popov et al. 2004;Kováč et al. 2017), during its existence, short-living straits separated, or narrow land bridges connected the Southeast European and the Western Mediterranean drylands across the sea. It indicates that the Central Paratethys could play an important role in the speciation and former migrations of the Mediterranean sandfly fauna. However, due to the lack of sandfly fossils in the Neogene Circum-Mediterranean fossil record, this possible palaeobiogeographic significance of the Central Paratethys directly cannot be investigated.

Aims
Because the investigation of the former environmental alterations on the past migrations and occurrences of sandfly species can help to create better projections for the future spread of this vector group, it was aimed to clarify the potential role of the Central Paratethys in the speciation and/or migration of sandfly species between western and south-eastern Europe by modelling the Oligocene-Miocene potential occurrence of the species. The aim of the study was to test whether the Central Paratethys area itself could be the speciation/migration area of the ancestors of certain sandfly species or not. It can be hypothesised that if the climate of an area, prior to the emergence of an extant sandfly species, would not have allowed the presence of it in an area, then it is very likely that this area was not part of the range of the ancestors of the present-day species during their speciation.

Theoretical considerations
Due to the limitations of the fossil record, we have no exact knowledge about the ecological requirements of ancient insect species. The usage of the distribution-limiting extrema of extant taxa is basically a question of the taxonomic longevity of phlebotomine species. It is known that the average life span of a species is 5-10 My (May et al. 1995). Important information in this regard is that, e.g. Phlebotomus (Paraphlebotomus) similis Perfiliev, 1963 and Ph. sergenti presumably already existed at the time of the Messinian salinity crisis. It is because the continental and island populations of both species were separated at about the same time based on the most parsimonious phylogram of Phlebotomus (Paraphlebotomus) kazeruni Theodor and Mesghali 1964, Phlebotomus (Paraphlebotomus) jacusieli Theodor, 1947, Ph. similis and Ph. sergenti (Depaquit et al. 2002). There is currently no more substantiated reason as to how these relatively poor flying capability species -e.g. the maximal travel distance of Ph. papatasi (Scopoli 1786) is only 1.9-1.5 km according to Orshan et al. 2016 -made the sea voyage between the Middle East and Cyprus or Northwest Africa and Sicily (in the case of Ph. sergenti) or between the Balkan Peninsula and Malta (in the case of Ph. similis) about at the same time unless we hypothesise that they made the trip with "dry feet" in the largely dry Mediterranean basin in the Messinian stage (Trájer 2021a). The role of the Messinian salinity crisis in the evolution of the Mediterranean sandfly taxa was also considered e.g. by Depaquit et al. (2002), Kasap et al. (2015) and Cruaud et al. (2021). Because this event ended in about 5.33 Ma, it can be assumed that the two species may have developed earlier than the late Miocene epoch. According to Cruaud et al. (2021) Even if we accept the ecological conservatism of certain sandfly taxa, naturally, in the case of the earlier periods we cannot be sure that the climatic limitations of the ancestor of a species, e.g. in the Miocene or the Oligocene epochs were the same as in the present times. It is also questionable how former geographic and meso-and microclimatic factors impacted the phenotypic variation of ancient sandflies which could result in differences between sandfly populations. Nevertheless, Ph. ariasi exhibits notable phenotypic plasticity according to altitude for both sexes and to slope and station for females (Prudhomme et al. 2016). Also, significant size variation was detected for both sexes between months in the case of Ph. tobbi (Oguz et al. 2017). Although meteorological factors affect sandfly phenotype, it is not clear whether it is not known what differences exist or existed in the past between different populations in terms of their climatic needs. Based on the above-described uncertainties, the climatic extrema of an extant sandfly species can only be used as a model of the ancestral sandfly taxa. The idea is similar to the basic concept of the Coexistence Approach of Mosbrugger and Utescher (1997) in the sense that it assumes that the climatic needs of (plant) taxa did not change substantially in the last about 30-40 million years. A certain level of niche conservatism can be hypothesised in the case of sandfly species because the Mediterranean sandfly species reach their absolute northern distribution about at the border of the warm temperate (Cfa/Cfb) and humid subtropical (Dfa/Dfb) climates in Europe even though their populations have been repeatedly suffered bottleneck effect during the Pleistocene climate changes due to their inability to adapt to cooler temperate climates (Esseghir et al. 1997;Moo-Llanes et al. 2019).
The logical framework and methodological considerations of the study The following work plan has been identified ( Fig. 2): 1) Eight extant sandfly species as model species of the ancestors of their subgenera were selected for modelling purposes, namely: Ph. perfiliewi, Ph. perniciosus, Ph. tobbi, Ph. sergenti and Ph. similis, Ph. ariasi, Ph. neglectus, Ph. papatasi. Phlebotomus sergenti and P. similis were frequently investigated species of their subgenera in the last two decades in palaeozoogeographical studies. It is possible that the Miocene geographical changes in the Tethys and the Paratethys realm, as well as such events like the Messinian salinity crisis, could play an important role in the speciation and the migration of these species and their relatives (Depaquit et al. 2002;Dvořák et al. 2006;Trájer 2021a). The contours of a similar evolutionary history arose in relation to Ph. papatasi (Depaquit et al. 2000(Depaquit et al. , 2001(Depaquit et al. , 2008. In the case of the members of Larroussius, it was found by Esseghir et al. (2000) that the speciation of members of the subgenus Larroussius also could be strongly impacted by Miocene palaeogeographicalpalaeoclimatic changes and their important evolutionary events could also coincide, e.g. with the late Miocene-Pliocene aridification of the Mediterranean subregion. It can be interesting to compare some   (Hijmans et al. 2005) was modified by the calculated IDW-interpolated values in the case of all factors. As a result of this process, a vector kind of palaeoclimatic model was created. 5) As an additional factor, the so-called Thornthwaite agrometeorological (aridity) index was also calculated both for the pre-Tortonian stages and the Tortonian stage also (the description of the index is seen below). 6) Based on the Climate Envelope Modelling method, the past suitability values of the extant sandfly species were calculated based on the point-like Kiscellian-Sarmatian Central Paratethys stages' climate data and the created georeferenced Pannonian (Tortonian) palaeoclimate model. This model was constructed for the purposes of the present study based on the palaeoclimatic data of Bruch et al. (2006). The relative sandfly activity models of the warmest quarter were based on the same climatic data. 7) The Pannonian (Tortonian) relative abundance values of four sandfly species from the studied eight extant sandfly taxa was modelled related to the former temperatures of the warmest quarter based on the published temperaturecaptured sandfly number data of Tsirigotakis et al. (2018).

The studied periods
In this study, the studied stratigraphic units were denominated both according to the regional and in brackets, the equivalent global units. This is indicated because palaeogeographical and slate climatological information on the area can be searched in the literature primarily based on regional units. The correlation of the regional (Central Paratethys stages) and global chronostratigraphic units was based on the publication of Kováč et al. (2017). The following stages were considered: the Kiscellian (Rupelian; lower Oligocene), the Egerian (middle and upper Chattian and lower and middle Aquitanian; upper Oligocene and the lowermost part of the lower Miocene), the Eggenburgian (upper Aquitanian and the lower part of the Burdigalian, the middle part of the lower Miocene), the Karpatian (an upper, but not the uppermost part of the Burdigalian, the upper part of the lower Miocene), Badenian (the uppermost Burdigalian, the Langhian and the lower part of the Serravallian; upper lower and middle Miocene), the Sarmatian (upper Serravallian; middle Miocene) and the Pannonian (the Tortonian; upper Miocene). Figure 3 visualizes the temporal connection of the regional and global stratigraphic units with the palaeoclimate data-bearing sites and the approximate palaeo-topographic conditions of the periods.

The used sandfly distribution limiting extremes
To create the distribution models, the occurrence limiting climatic variables of Trájer et al. (2013) were used. The Authors performed their analyses based on the sandfly distribution data of the European Centre for Disease Prevention and Control in 2012. The European Centre for Disease Prevention and Control uses several occurrence categories, like 'present' and 'introduced'. The applied reference period climate model of this study was the Regional Climate Model (REMO) regional climate model. The reference period of REMO was 1961REMO was -1990. The somewhat earlier reference period compared to the publication time of the sandfly distribution data was since 1) the several occurrence data related to sandfly species was based on earlier publications and 2) the dispersal rate of sandflies is plausibly low. Due to the accelerating warming of the global climate and the relatively low dispersal rate of sandfly taxa, there is a gap between the observed and the theoretically possible distribution area of the species, which is clearly visible, e.g. in the comparison of the model results of Trájer et al. (2013). Trájer et al. (2013) handled only the 'present' category data as positive areas for permanent occurrence. For the explanation of the employed distribution limiting factors, see the model identification subchapters. Since the original extremes were expressed as monthly mean temperature and Fig. 3 Correlation of the studied periods in regional and global units according to Kováč et al. (2017), the sites of fossil floras which provided the past climatic data and the palinspastic reconstructions (P: period, E: epoch, s-E: sub-epoch, gS: global stage, rS: regional [Central Paratethys] stage) precipitation sum values, those factors were used to determine which values were convertible to bioclimatic values. It was necessary because palaeoclimatic reconstructions use bioclimatic variables.

Palaeogeographical and palaeotopographic reconstructions
The crust of the Central Paratethys is a mosaic of different origin terranes and until the late Miocene epoch, it also contained oceanic crust elements. The European and the Moesian continental crusts are the only firm and stable tectonic units of the area. The AlCaPa, Tisza and Dacia terranes of the Central Paratethys were in constant motions, including also notable rotating movements during the entire Paleogene period to the end of the Middle Miocene epoch. It should be noted that the crustal blocks are still moving today due to the rotation and northward movement of the Adriatic Plate (Handy et al. 2015). The terranes approximatively reached their present relative positions in the early Miocene epoch (Kvaček et al. 2006;Kováč et al. 2017). However, in the middle Miocene epoch, due to an intensive extensional back-arc deformation, the crustal blocks of the Pannonian Basin extended and sunk, and the sea transgressed to the coastal lowlands and marshes of the basin (Matenco and Radivojević 2012). The topographic model of the Central Paratethys from the Kiscellian to the Sarmatian Central Paratethys stages was based on the palinspastic reconstructions of Kováč et al. (2017). Naturally, in a primarily palaeobiogeographic study, the accurate reconstruction of the former geographical conditions cannot be expected, which belongs to the field of geology. From 10 Ma, the geographic patterns of the circum-Paratethys ranges become like the present conditions. Due to this fact, the present-day Carpathian Basin can be held in the topographic sense as the basin of the Pannonian Sea filled with alluvial sediments. These facts indicate that the present topographic conditions can only form the base of correct palaeoclimatic, palaeobiogeographic reconstructions from the Pannonian Central Paratethys stage. Earlier in this period, the relief patterns were not known enough to create topographically adequate models. The shorelines of the Pannonian Central Paratethys, the so-called Pannonian Lake have also been continuously changed. The approximate shoreline patterns of the Pannonian Lake used in this study were based on the reconstruction of Magyar et al. (1999) and represent the approximate shoreline patterns of about 10 Ma. The reconstructed topography of the Pannonian Central Paratethys stage was based on the modification of the present conditions, which method can be justified by the fact that the development of the Carpathian Basin occurred prior to the Pannonian Central Paratethys stage, in the turn of the early and middle Miocene epochs (Ustaszewski et al. 2008).

Palaeoclimatic reconstructions
The Kiscellian to the Sarmatian Central Paratethys stages' climate values used in this study were based on the palaeoclimatic reconstructions of Erdei et al. (2007). The authors adopted the Coexistence Approach of Mosbrugger and Utescher (1997) considering the fossil plant assemblages of the Pannonian domain (the area of the former Central Paratethys) to estimate the former palaeoclimatic conditions of the fossil sites. The Coexistence Approach method or concept is based on the range-limiting values of the closest living relatives of the fossil flora elements. This method can be utilised for the Cenozoic period due to the relative time proximity, if the fossil record is treated with caution (e.g. Erdei et al. 2007) discredited the Ceratozamia and Quercus sect. Cerris taxa from the analyses due to their uncertainty of known climatic needs), many fossil species are used, and the climatic requirements of the extant taxa were processed with appropriate statistical methods (Mosbrugger and Utescher 1997). Erdei et al. (2007) produced the values of the mean annual temperature (Tam,°C), the temperature of the coldest month (Tmcq,°C), the temperature of the warmest month (Tmwq,°C) and the mean annual precipitation (Pas, mm). The Pannonian palaeoclimatic data of the same region was gained from Bruch et al. (2006). As Erdei et al. (2007), Bruch et al. (2006) also used the Coexistence Approach (Mosbrugger and Utescher 1997) to reconstruct the Pannonian climatic conditions. From the produced palaeoclimatic values, the mean annual temperature (Tam,°C), the temperature of the coldest month (Tmcq,°C ), the temperature of the warmest month (Tmwq,°C) and the mean annual precipitation (Pas, mm) were used for modelling purposes. The palaeogeographical reconstruction of the Pannonian Central Paratethys was based on the findings of Kázmér (1990), Magyar et al. (1999) and Matenco and Radivojević (2012). Figure 4 shows both the fossil palaeoflora sites of the Pannonian Central Paratethys stage and the palaeogeographical reconstruction of the Central Paratethys used in the modelling. Table 1 shows the mean climatic values of the fossil flora sites in Hungary by periods.

Model environments
The first model environment (model environment A) was performed to model the suitability and the warmest quarter's relative activity values of the involved sandfly species by sites for the Kiscellian, the Egerian, the Eggenburgian, the Karpatian, Badenian, the Sarmatian Central Paratethys stages. Because the exact topographical conditions of the pre-Pannonian Central Paratethys are only approximately known, only the model results were expressed as site (point) -like data according to the published localities of Erdei et al. (2007). In contrast, the Pannonian (Torontonian) climate data were used to create a continuous climate model for the Central Paratethys area. Since high-resolution climate models were not available for the Pannonian Central Paratethys area, which is identical to the region of the former Pannonian Lake, a climate model was developed based on the modification of the reference climate maps according to the palaeoclimatic data of Bruch et al. (2006).
The second model environment (model environment B) was established to investigate the potential suitability, the warmest quarter's relative activity patterns of sandfly species in the Pannonian stage. Based on the palaeoclimatic reconstructions of Erdei et al. (2007). In the case of the creation of the climatic model of the Pannonian stage, the fossil plant species-inferred palaeoclimate data was thought to be characteristic of the former areas close to sea level. It was hypothesised that plant fossils were deposited in river deltas or in the calderas of maar craters in which the bottom was under or close to the former sea level (the sediments of higher elevation calderas generally were eroded in the last million years). The late Miocene Pannonian age model was based on the palaeoclimatic reconstruction of Bruch et al. (2006). The Pannonian is a Central Paratethys Stage that lasted from 11.6 to 6.15 Ma and is equal to the Tortonian and lower Messinian stages, but the analyses of Bruch et al. (2006) covered the ∼11 to ∼7 Ma section, Pannonian Central Paratethys stage of these regional stratigraphic intervals which is equal to the Tortonian global stratigraphic unit. Because in this case the palaeotopographic conditions of the Central Paratethys could be reconstructed, the arealbased modelling purposes could be performed. Table 2 shows the calibration of the Pannonian climate model.

The climatic factors
The used climatic variables for suitability modelling purposes were as follows: the lower and upper limits of the annual mean of the temperature (Tam;°C; Bioclim1), the mean temperature of warmest quarter (Tmwq limit.min,max , mm; Bioclim10) and the mean temperature of coldest quarter (Tmcq limit.min,max , mm; Bioclim11); as well as the lower and upper limits of the annual precipitation sum (Pas limit.min,max , mm; Bioclim12) and the lower and upper limits of the mean monthly Thornthwaite agrometeorological (aridity) index (TAIm limit.min,max , mm°C -1 ); As can be seen in Table 3, the limitation factor values -the distribution-limiting climatic extrema of the species -are climate model-independent. This means 5x2 bioclimatic factors in the model in the case of both species, which is an acceptable number for distribution modelling purposes. For the numerical modelling of the potential distribution areas of species, these 5x2 extrema of bioclimatic factors were used.
In the case of the Pannonian age model, the palaeoclimatic conditions were coupled with the palaeogeographical conditions shown in Fig. 4. Because the original data of Bruch et al. (2006) were given as site-based climatic values, the climatic values were interpolated to produce areal-based climatic models for the factors described above. At the first step, to create a topography-sensitive model, the reference climate model of WordClim 1.4 climate database ) (Hijmans et al. 2005) was sampled for the same bioclimatic values in the fossil flora sites to calculate the differences in the bioclimatic values between the present and the former climatic conditions. Then, the sites were georeferenced and weighted with the difference values using the Inverse Distance Weighting (IDW) interpolation method. The former palaeoclimatic values were modelled according to the following formula (Keckler 1994): Table 1 The mean climatic values of the fossil flora sites in Hungary by Central Paratethys stages. Tam: annual mean temperature (°C; equal to bioclimatic variable 1), Tmwq: mean temperature of warmest quarter (°C; equal to bioclimatic variable 10), Tmcq: mean temperature of coldest quarter (°C; equal to bioclimatic variable 11), Pas: annual precipitation (mm; equal to bioclimatic variable 12) and TAIm (mm°C -1 ; it is the (annual) mean monthly Thornthwaite agrometeorological (aridity) index) where Z p is an interpolated palaeoclimatic value at the grid node in cm; Z i is a palaeoclimatic value at a site in cm (xi, yj), 1/d p i is a weighting function; n is the number of fossil flora sites that were the basis of palaeoclimatic reconstructions. The d p i is the distance between Z p and Z i sites.
To increase the number of climatic variables, a derived variable, the annual monthly average Thornthwaite agrometeorological (aridity) index (sing. abbr.: TAI) was introduced which can be derived from the temperature and precipitation values (Kemp 1990). The use of an aridity index as additional distribution limiting factor is justified by the fact that sandfly species are sensitive to the combination of strong solar radiation and drought (Trájer et al. 2018;Elnaiem et al. 1998).
The formula of the TAI index is as follows: where TAI is the Thornthwaite agrometeorological (aridity) index in mm°C -1 , P is the monthly precipitation sum in mm, T is the monthly mean temperature in°C.
The (annual) mean monthly TAIm values were calculated according to the following formula: where TAImis the mean monthly Thornthwaite agrometeorological (aridity) index in mm°C -1 , Pas is the annual precipitation sum in mm, Tam is the annual mean temperature in°C. Finally, the reference period's climatic model was modified with the interpolated values by bioclimatic factors and TAIm. The results were projected to the Pannonian topographic models of the Central Paratethys and were represented as heat maps.
It is a common element of model environments A and B that the modelling was based on the binary logic of the Boolean algebra: 1 Tam ð Þ¼ 0 if Tam limit:min > Tam > Tam limit:max 1 if Tam limit:min ≤Tam≤Tam limit:max ð4Þ Where Tam represents the lower and upper limits of the annual mean temperature, the Tmwq is the georeferenced climate model data of the lower and upper limits of the mean temperature of the warmest quarter, and Tmcq is the georeferenced climate model data of the lower and upper limits of the mean temperature of the coldest quarter. Pas represents the georeferenced climate model data of the annual precipitation sum. TAIm is the mean monthly Thornthwaite agrometeorological index.
In the case of model environment A, the suitability values were calculated as point-like values according to the sites of palaeoclimatic reconstructions.
In the case of model environment B, the equations were displayed as areal-based reconstructions because the Pannonian site to site climatic values was interpolated by climatic factors as was described above. The potential area-based suitability patterns were determined according to the following mathematical formalism: A Tam; Tmwq; Tmc; Pas; TAIm ð Þ ¼ Tam⋂Tmwa⋂Tmcq⋂Pas⋂TAIm Where A(Tam;Tmwq;Tmcq;Pas;TAIm) shows the potential distribution area of the given species, which contains the remaining areas after taking into consideration the temperature and precipitation limitations. Summarising the used climatic factors, a total of 14 variables was involved, but it varied by the model environments in which variables were used in the modelling in each geological period.
Calculation of the relative activity values Tsirigotakis et al. (2018) sorted the number of captured sandfly individuals by temperature conditions during their trapping efforts by species. Because Tsirigotakis et al. (2018) reported their trapped sandfly number data in 3-degree intervals, the mean temperature values were paired with the modified values of the collected sandfly numbers. This modification means that in the present study, the function between the temperature and the modified number of the trapped individuals of Ph. neglectus, Ph. perfiliewi, Ph. similis and Ph. tobbi was used to calculate the temperature-related relative (%) numbers, or in other words, as a kind of temperaturerelated sandfly activity patterns. In order to, the number of captured individuals was transformed as follows: where the RA is the relative (%) sandfly activity; N i is the number of the trapped sandflies; N max is the maximum number of sandflies trapped in a case of a temperature value.
Then, the correlation between the physical factors and the relative sandfly activity values were determined. The correlation between the temperature and relative sandfly activity values were approximated by a Gaussian function (r 2 =0.9483, p=0.00493 in the case of Ph. neglectus; r 2 =0.9999, p<0.0001 in the case of Ph. perfiliewi; in the case of Ph. tobbi where RA is the relative activity (%); T mwq is the mean temperature of the warmest quarter.
The RAvalues were calculated according to the mean temperature of the warmest quarter. The applicability of the mean temperatures is since outside the tropics, sandfly species are predominantly active in the evening and early night hours (Kravchenko et al. 2004;Boussaa and Boumezzough 2006;Trájer et al. 2018). In the tropics, the diurnal activity differs from this, because the peak nocturnal activities of sandflies can also occur in the night hours, even after midnight (Aklilu et al. 2017). Because there is a lag between the run of the solar cycle and the run of temperature, the temperature regime of early night hours can be well-approximated by the daily mean temperature (Ephrath et al. 1996), or in a wider sense, by the weekly, monthly, and quaternary mean temperatures -depending on the time frame of the investigated period and the purpose. Since the activity model was only performed in the case of the Tortonian stage and the former climate of the Pannonian age Central Paratethys could be humid subtropical kind of (Cfa) Hably 2002;Utescher et al. 2007;Bruch et al. 2011;Ivanov et al. 2011;Utescher et al. 2017), the use of the mean temperature values for the modelling of the relative activity patterns is an acceptable approximation.
In the case of the older (pre-Tortonian) stages, the climatic suitability values were calculated by the individual sites of palaeoclimatic reconstructions. Model results were displayed using Quantum GIS 3.4.4 software. The Lambert Azimuthal Equal Area (EPSG:3035) was used as a projection system in each case.

Oligocene-middle Miocene suitability patterns
The Kiscellian suitability values were modelled for three locations. In this period, the suitability values could be generally low in each model species compared to the later periods. The suitability values did not reach 100% value except for a sandfly species with Ph. perfiliewi-like environmental needs. The modelled suitability values of a sandfly species whose climatic needs may have been like those of Ph. ariasi do not exceed the 70% value. The suitability values of Ph. sergenti and Ph. similis-like sandfly species could range between 70 and 80%. The maximum suitability values for a sandfly species with Ph. neglectus, Ph. papatasi, Ph. perniciosus and Ph tobbi -like climatic needs reached the 90% value, although it ranges between 70-80 and 90% values. The Kiscellian suitability values were modelled for five locations. In the Egerian Central Paratethys stage, the suitability values showed generally the highest values in the case of all species. Except for a Ph. ariasi-like ancient sandly species, the modelled maximum suitability values by species reached 100% value. The modelled suitability values of Ph. ariasi-like sandfly species could range between 70-90%. In each location, the modelled values of sandfly species whose distribution determining climatic extrema may have been like the limits of Ph. papatasi, Ph. perfiliewi and Ph. tobbi reached the 100% value. In the case of Ph. neglectus-like sadfly species, the modelled values ranged between 90-100% values by location. The modelled values of sandfly species with the climatic limits of Ph. sergenti and Ph. similis were between 80-100 and 70-100% values, although it can be added that in the case of a Ph. similis kind of sandfly species the other values ranged between 70-80% values (Fig. 5).
The Badenian suitability values were modelled for one location. In the case of sandfly species with similar climatic needs which can be found in the case of the extant Ph. sergenti and Ph. similis, the modelled suitability values could be 90%.
The suitability values of the other sandfly species could reach 100% value. The Sarmatian suitability values were modelled for three locations. The maximum suitability values reached 100% value in each species, and except for a Ph. similis-like ancient sandfly species in one location, it was true for each modelled value (Fig. 7).

Relative temperature-dependent activity values of location-based models
The Kiscellian temperature-dependent relative activity values in the warmest quarter of the year ranged between the high 39-43% values. In the Egerian Central Paratethys stage, the values ranged between 34-43% values. The modelled activity values of the Eggenburgian and Karpatian Central Paratethys stages also reached 39-43% values. In contrast, the modelled activity value of the Badenian and Sarmatian Central Paratethys stages were only 24-28 and 19-33% (Fig. 8).

Model environment 2 results
The reconstructed climatic patterns in the Pannonian stage The highest mean annual and the warmest quarter's temperatures occurred along the southern and eastern coastlines of the Pannonian age Central Paratethys. In the coldest quarter of the year, the highest temperatures could be measured along the southern coasts of the sea and in the islands. The wettest areas could be found on the southwest coasts. Figure 9 shows the IDW interpolated maps of the climatic variables.

Pannonian suitability patterns
In the Pannonian stage, the west, north and east shorelines and islands of the Central Paratethys (the Pannonian Lake) showed high suitability values (suitability values were between 90-100%) for ancient sandfly species that had similar environmental needs to Ph. ariasi, Ph. neglectus, Ph. perfiliewi Ph. perniciosus and Ph. tobbi. The climatic conditions could be less suitable for species with Ph. papatasi and Ph. sergenti (suitability values are between 80-90%), and even less for Ph. similis-like climatic requirements (suitability values were between 70-80%). A Phlebotomus ariasi-like species could be present along almost the entire shoreline and the islands of the Central Paratethys, excluding the southwest edge of the sea. Phlebotomus neglectus, Ph. perfiliewi, Ph. perniciosus, and Ph. tobbi-like sandfly might have been lacking from the southwest coasts of the Central Paratethys (Fig. 10).

The relative activity values in the Pannonian stage
The modelled temperature-dependent relative activity values in the warmest quarter of the year show heterogenic patterns depending on the topography and region and species. The  highest activity values can be seen at the south, west, and east shorelines of the former Pannonian Lake and in the islands. Along the shorelines, the relative activity values could reach the 78-100% value during the Pannonian stage in the case of a Ph. neglectus-like ancient sandfly species. In contrast, in the higher elevation regions, the relative activity values in the warmest quarter of the year could be below the 50% value. However, in the case of ancient sandfly species with similar temperature-dependent activity functions like the present-day Ph. perfiliewi, Ph. similis and Ph. tobbi, the highest relative activity values also could be observed in the islands of the Pannonian Sea and along the southern shorelines, the maximums of the warmest quarter's mean temperature-based relative activity values could less than 90% even in most of the south coastal sites, and in the inner areas of the drylands, these values could be less than 50-60%. Because relatively high mountain ranges bordered Central Paratethys, the sandflyinhabited areas could be restricted to the relatively narrow coastal plains and meridional valleys. In accordance with these results, the migration of sandflies could be possible along the shorelines and coastal valleys parallel to mountain ranges. The islands of the Paratethys Sea could also be suitable for sandflies (Figure 11).

Discussion
The tropical or hot subtropical climatic conditions of the early Oligocene epoch, in the Kiscellian Central Paratethys stages, would not have been suitable for an ancient sandfly species  whose climatic needs would have been similar to the applied Mediterranean sandfly species, in general. In the second part of the Oligocene, due to the drier subtropical conditions, the modelled maximum suitability values could reach the highest valueexcept for a species with Ph. ariasi-like climatic needs. In the early Miocene era, the modelled values could be lower for ancient sandfly species with Ph. ariasi, Ph. similis and Ph. sergenti-like environmental needs due to the wetter subtropical climate of the era but could reach the maximum values for such kind of ancient sandfly species which occupied similar niche like the present-day Ph. neglectus, Ph. papatasi, Ph. perfiliewi, Ph. perniciosus and Ph. tobbi. In the middle Miocene period, high (90-100%) suitability values were modelled for all model species, although the climate of the Badenian stage could be less suitable for a Ph. sergenti and Ph. similis-like ancient sandfly species. Because the Tortonian was the last period when the Central Paratethys existedapart from the Messinian age Slavonian Lake which was a relatively small remnant of the Pannonian Sea (Kázmér 1990)it could be a critical period in the aspect of the migration of sandfly species between western and eastern Europe. Surprisingly, the model results showed lower suitability values for the Slovenian land bridge (corridor) than for the Vienna Basin. It is an important result because it could be hypothesised that before the disappearance of the Central Paratethys, the Slavonian land bridge and the North Adriatic area could be a potential migration route for sandfly species. These results contradict the hypothesis that the former Slovenian land bridge could be the main migration route of ancient Paraphlebotomus-like species between the (south)east and the western parts of Europe in Tortonian. As it was hypothesised by Depaquit et al. (2000) that the ancestors of Ph. sergenti and Ph. similis may appear during the late Miocene epoch along the north and south shorelines of the Paratethys. Considering the more novel results of Cruaud et al. (2021), it can be said that the Tortonian stage could be the time of the speciation and radiation of the extant Paraphlebotomus species. It may indicate that it is very improbable that the ancestor of the extant Paraphlebotomus species could migrate between East and West Europe through the Central Paratethys. The opening of the migration routes for sandfly taxa could be the consequence of the Messinian salinity crisis-induced climatic and palaeogeographical changes (Trájer 2021a(Trájer , 2021bTrájer et al. 2021).
The modelled Tortonian relative activity values for Ph. neglectus, Ph. perfiliewi, Ph. similis and Ph. tobbi-like ancient sandfly species indicates that although these values generally could be lower in the northwest part of the Tortonian Central Paratethys, it cannot be excluded that the Vienna Basin and the northern foredeep of the Alpsthe remnant of the former Fig. 8 The modelled temperature-dependent relative activity values of the hypothetical ancient climatic analogues of extant sandfly species in the warmest quarter of the year during the Oligocene and Miocene epochs. a Kiscellian, b Egerian, c Eggenburgian, d Karpatian, e Badenian, f Sarmatian Fig. 9 The IDW interpolated differences of the Pannonian age climatic variables from the reference period's values (the figures of the 1st column) and the reconstructed Pannonian climatic patterns (the figures of the 2nd column) ca. in 10 Ma (a1-a2 Tam: mean annual temperature, b1-b2 Tmwq: mean temperature of the summer months, c1-c2 Tmcq: mean temperature of the winter months and d1-d2 Pas: the annual precipitation sum; e1-e2 TAIm: mean monthly Thornthwaite (aridity) index)  Western Paratethyscould be main migration route between the west and (south)east parts of Europe during the Tortonian and not the northeast Adriatic region. The results indicate that the Tortonian climate of the Central Paratethys could not have been equally favourable for two ancient sandfly species, even if otherwise the climatic conditions would have allowed both species to occur in the study area. This circumstance is clear in the comparison of the modelled relative suitability values of Ph. neglectus and Ph. perfiliewi for this period.
Considering the results, it is hard to explain that the present distribution of Ph. ariasi in Western Europe is limited to the boarding areas of France, Piedmont, and Liguria regions (Maroli et al. 1997). On the other hand, it is also curious that Ph. ariasi and Ph. neglectus, those species that were supposed to share a common origin and have similar ecological requirements (Pesson et al. 1994), have only a small overlapping distribution area in Northwest Italy. One of the possibilities is that the divergence of the lineages of the ancestors of the two species occurred later than the Tortonian stage. It was also found that the climate of the Pannonian Central Paratethys was rather suitable for Ph. ariasi than Ph. neglectus at that time and the climatic suitability of the more remote areas, like the north and south foothills of the Alps or the northwest coasts of the former Mediterranean Sea, were not modelled. Another possibility is that geographical barriers prevented the eastward migration of Ph. ariasi's predecessor. Because the Central Paratethys could not be suitable for such a kind of hot and semi-arid climate preferring ancient Paraphlebotomus species like Ph. sergenti and Ph. similis today, this circumstance could explain why these species lack from the Apennine Peninsula and most parts of western Europe. Depaquit et al. (1998Depaquit et al. ( , 2000Depaquit et al. ( , 2002 showed that Ph. sergenti originated in the Middle East and migrated along the northern shorelines of Africa reaching the Iberian Peninsula through Gibraltar. The most distant area where Ph. sergenti has reached is South France. Because even in the coldest periods of the Quaternary period such as the Last Glacial Maximum ca. 21 ka, refugia existed for sandfly species at the eastern coastal foothills of the Pyrenees (Trájer and Sebestyén 2019), these populations could have survived to the present day. Another problem is the presence of Ph. sergenti in Sicily. Because large glacial refugia existed along the shorelines of the Ligurian and Adriatic coasts of the Apennine Peninsula, it is very unlikely that Ph. sergenti or Ph. similis died out in the entire Apennine Peninsula in the last million years. It is more plausible that this species never reached the Apennine Peninsula and Ph. sergenti might pass through the Pantelleria Graben, for example, during the Messinian (Depaquit et al. 2002). It could happen because it is plausible that the most ancient population of Ph. sergenti appeared at the earliest around in the Messinian period. However, in the case of Ph. similis, the impact of the Messinian salinity crisis on the speciation and migrations of this sandfly species is questionable because it is likely that it could appear later than the late Miocene epoch (Cruaud et al. 2021). Following this logical line, similar migration and subsequent speciation events were hypothesised about the ancestors of the members of subgenus Larroussius by Esseghir et al. (2000) and the ancestors of the extant species of subgenus Transphlebotomus by Kasap et al. (2015) during the Messinian salinity crisis. However, it should be noted that from a physiological point of view, it has never been examined whether the environment-sensitive sandfly species would have been able to migrate along with the extremely hot (greater than 35°C; Murphy et al. 2009), dry plain (Cita, 2001) of the desiccated Mediterranean Basin or not.
However, it should not be forgotten that the biogeographic history of sandflies which determine the present-day occurrence of the species could be both strongly modified natural factors and human activity. It can be hypothesised that these effects could differently influence the survival of the eastern and western Mediterranean populations, including the occurrence of the relictual localities. For example, the large, near-supervolcanic eruptions of the Mediterranean volcanoes, e.g. the Pozzolane Rosse Tephritic Ignimbritecreating VEI 7, Ultra-Plinian eruption in 456-439 ka in Central Italy (Giordano and Doronzo 2017) or the also super-colossal eruption of Santorini in 1610 BC in the Aegean Sea could also heavily affect the sandfly populations at the regional level (Trájer 2021a). It is also plausible that the effects of the glacial phases on the South European sandfly populations could vary by area in the Mediterranean region. A model study showed that while during the Last Glacial Maximum (~21 ka), sandfly species could disappear from South France and from the western regions of the Balkan Peninsula, stabile populations could persist in the Mediterranean islands, the Aegean Archipelago and along the eastern and southern coasts of the Iberian Peninsula (Trájer and Sebestyén 2019). The former nautical activity of the ancient merchant states of Greece and Phoenicia, the Roman Empire or medieval Venice and Genoa could also modify the presently known occurrence of sandfly species since maritime trade still plays a role in the distribution of arthropod vectors (e.g. Hofhuis et al. 2009).
It should be noted that the number of fossil sites involved in the study was different. One site represents the different strati-graphic units of the Langhian and the Burdigalian, three sites the Serravallian, four sites the Rupelian, five sites the Chattian+Aquitanian era and the Pannonian climate model was based on the IDW interpolation of the climatic data of nine sites. This condition is the consequence of the nature of the original climate models of Bruch et al. (2006) and Erdei et al. (2007). Seeing those periods in which at least two sites were involved, certain fluctuation exists between the suitability and other modelled values of the sites. Although, because all the palaeoflora-bearing strata which became the basis of the palaeoclimatic reconstructions sedimented originally along the coasts of Central Paratethys. It means that although the model results of one timeperiod can be handled with a certain caution, these models have an indicative value. To evaluate the possible evolutionary importance of the model results, it is inevitable to involve the results of palinspastic reconstructions. Based on the palinspastic reconstructions of Kováč et al. (2017), the Slovenian corridor was open in the Egerian, Karpatian and Badenian Central Paratethys stages. The status of the corridor is questionable during the Kiscellian and Sarmatian Central Paratethys stages. These facts indicate that the free migration between central-western Europe and the Balkan Peninsula could be given during the Eggenburgian, and the Pannonian and maybe, not permanently, in the Kiscellian and Sarmatian Central Paratethys stages. In addition, in the Kiscellian and the Eggenburgian Central Paratethys stages, the Western Paratethys has separated the Eastern Alps from the continent, and during the Egerian and the Sarmatian Central Paratethys stages, the Anatolian-Balkan landmass was separated from Eastern Europe (Kováč et al. 2017 perniciosus, and Ph. tobbi are the extant analogues of ancient sandflies in the sense of their possible climatic distribution-limiting extrema, but the wide Slovenian strait formed an insurmountable obstacle for sandfly species. 5) In the Badenian Central Paratethys stage, the climate, maybe except for a Ph. sergenti or a Ph. similis-like ancient sandfly species, could be suitable for the studied species, although the wide Slovenian strait still formed an insurmountable barrier for the migration of sandfly species between the Balkan and the north Adriatic area. 6) In the Sarmatian Central Paratethys stage, the climate could be suitable for any ancient sandfly species with climatic needs like the applied eight extant sandfly taxa. Because the barrier role of the Slovenian corridor is questionable, it cannot be excluded that sandfly species could migrate between the Balkan and the North Adriatic region. Although at that time, the connections between the Central and East Paratethys could be a limit against the communication between the sandfly fauna of East Europe and the Balkan and Asia Minor. 7) In the Pannonian era, the climate of the shorelines of the central Paratethys could be suitable for most of the used model sandfly species, except for ancient sandfly species whose climatic needs would have been the same as those of Ph. papatasi, Ph. sergenti and Ph. similis, although, the Slovenian land bridge does not seem to be a suitable area for the studied species during the Tortonian.
To decide which speciation and migration possibilities could occur in the past, the exact determination of the speciation of the Mediterranean sandfly species could be needed. Although the former literature neglected the Central Paratethys only as a speciation area of sandfly species, the results show that the central part of the Paratethys could also be an adequate habitat for the ancestors of the present European sandfly taxa. It is less plausible for ancient Paraphlebotomus species, but the Central Paratethyan origin of certain Phlebotomus and Larroussius species cannot be excluded. It is also an important argument in favour of the above statements that due to the diverse former palaeogeographical conditions of the Central Paratethys, this area could provide adequate territory for the evolution and long-term survival of sandfly taxa. It should be noted that based on the palaeoclimatic reconstructions of Erdei et al. (2007), while in the late Oligocene climate of the Central Paratethys area could be monsoon-influenced humid subtropical (Köppen-Geiger cwa; mean annual temperature: about 15.6-20.6°C, mean annual precipitation: about 1300-1500 mm), in the early and middle Miocene epoch, it could be a somewhat dryer, possibly not a monsooninfluenced type humid subtropical climate (Köppen-Geiger cfa; mean annual temperature: about 14-16.5°C, mean annual precipitation: about 820-1300 mm). Based on the results of Erdei et al. (2007), the climate of the Central Paratethys in the late Miocene era could become a somewhat drier humid subtropical climate (the transition between the Köppen-Geiger cwa and csa; mean annual temperature: about 14-16.6°C, mean annual precipitation: about 760-1180 mm). Considering the present occurrence of Paraphlebotomus species which predominantly covers the arid, semi-arid and Mediterranean areas of North Africa, South Europe, and the Middle East, Cruaud et al. (2021) found that the most probable ancestral range of Paraphlebotomus is equivocal with such arid and semi-arid regions like Irano-Turania or North Africa. However, the direct comparison of today's vegetation conditions and biomes with the occurrence of today's sandflies has pitfalls. While many extant Mediterranean sandfly populations live in often shruband-bush areas where vegetation may have been secondary and anthropogenic. In the same areas, originally extensive forests and groves could occur in the early Holocene epoch as it was found in the case of the ancient Greek city Ephesus (Stock et al. 2020). This hypothesis is supported by the fact that the climatic needs of Mediterranean sandfly species correspond to the requirements of such forests-forming ligneous plant taxa like Quercus ilex Linneaus 1753 or Pinus brutia Tenore 1811 (Bede-Fazekas and Trájer 2013).
Although Paraphlebotomus is only one subgenus of Phlebotomus, this finding is in good accordance with the finding of the present study which suggests that the Central Paratethys could not be among the speciation area of such Paraphlebotomus species like Ph. sergenti and Ph. similis. Grace-Lema et al. (2015) also found in their phylogeneticzoogeographic nature study that Paraphlebotomus species and their relatives have a typical arid North and East Africa, Mediterranean and Middle East distribution and such Phlebotomus species, like the studied Ph. papatasi also belong to this ecological clade of sandflies. It also fits well with the presented model results since the modelled climatic suitability values of this species exhibited very similar changes between the Kiscellian and the Pannonian Central Paratethys stages due to the plausibly similar ancient speciation environment of these species. It should be noted that Ilango, (2004) also confirmed the close sister group relationship between Phlebotomus and Paraphlebotomus species -and as a third member of this group, Synphlebotomus, also. Larroussius is the sister group of Transphlebotomus (Ilango 2004) with a clear Mediterranean and East African distribution area. According to Cruaud et al. 2021) the divergence between the subgenus Larroussius and the Adlerius-Paraphlebotomus-Phlebotomus clade could occur at the turn of the Oligocene and Miocene epochs. Because Larroussius might have adapted to the climate of the Tethys and the Paratethys realms, it could explain why the climatic requirements of the members of this subgenus fit better to the palaeoclimate of the central Paratethys in the studied period.
Naturally, other Mediterranean phlebotomine sandfly species, which can also be found in the narrower and wider area of the former Central Paratethys, could also be used in this study like Phlebotomus (Paraphlebotomus) (Cazan et al. 2021;Dvořák et al. 2020;Hubenov 2021;Kniha et al. 2021). If we extend our field of study to the whole of the Mediterranean or to its wider environment, we should deal with many species of even a subgenus. Due to a large number of the possible modellable taxa that can be used as indicator species related to former climatic conditions, the used eight sandfly species were used bearing in mind that, of course, these species do not cover the full ecological spectrum of Mediterranean sandfly species.
It should also be discussed whether in which habitats the ancient sandfly species could persist along the shores of the Central Paratethys Sea. Of course, this is not to say that in the past, these species existed in the study area, as the used species were handled as indicators or types of sandfly species occupying more or less similar or different niches. It is known that such sandfly species like Ph. neglectus prefer the shaded cracks of karst rocks (Trájer et al. 2018). In the Slovenian coastal-karst area, five Phlebotomine sandfly species were identified by Vladimir et al. (2015) and Ivović et al. (2015): Ph. neglectus, Ph. perniciosus, Ph. mascittii, Ph. papatasi, and Sergentomyia (Sergentomyia) minuta Rondani, 1843. Also, a relatively rich sandfly fauna was found in the coastal district of Bar, Montenegro with the presence of Ph. papatasi, Ph. perfiliewi, Ph. tobbi, Ph. neglectus and S. minuta, which area has a typical Mediterranean coastal karstic landscape (Ivović et al. 2003). In Crimea, sandfly species feed on lizards that reside in caves on the southern coast of the peninsula (Turbanov et al. 2019). These observations confirm that several Mediterranean sandfly species can colonise the relatively dry environment. Geological evidence suggests that in the Pannonian era, the intensive abrasion of the limestone and dolomite ranges created similar coastal karst environments along the shorelines of the Central Paratethys which is characteristic of the present-day Mediterranean area (Popov et al. 2006;Pezelj et al. 2013). Fossilised abrasion coasts and abrasion gravel sediments are known in several sites of the Hungarian Transdanubian Mountains, for example in the Velence Mountains in Pázmánd (Gyalog and Ódor 1983), along the southern slopes of the Balaton Highland (Csillag et al. 2010) or on the Strázsa Hill in Zsámbék (Kercsmár et al. 2020).
Such Mediterranean sandfly species, like Ph. papatasi, use the burrows of small mammals under summer hot and dry Mediterranean climatic conditions as shelters and breeding sites (Wasserberg et al. 2003). The Pannonian coasts of the Central Paratethys were colonised by various small mammal assemblages (Mészáros 2000). Rich fossil small mammal record is known from the karst filling deposits of the coasts of the Pannonian Sea, e.g. from the about 9.5 Mys old palaeokarstfilling deposits of Sümeg, near to the preserved rocks of the former Pannonian abrasion shoreline and other sites in Hungary which could belong to islands and peninsulas in the Pannonian era (Mészáros 2000). A different kind of site was the subtropical swamp forest of Rudabánya, where also rich mammal fauna was found, including rodents (Daxner-Höck 2003). Although under humid subtropical and tropical conditions, sandflies can also be found in the estuarine wetlands and marshes (Brockmeyer et al. 1996;Barata 2017), in Europe, it is not common. For example, in the low alluvial plains and coastal marshes of the Camargue, South France, sandflies are uncommon members of the entomofauna. In this area, Ph. ariasi colonises only the hilltop environments of the delta, reaching its maximum frequency in such drier environments as mixed evergreen oak forests (Rioux et al. 1967). Considering the ecological potential of Ph. ariasi, it should be noted that this species does not prefer the low altitude regions (Aransay et al. 2004). It can be summarised that predominantly the karstic coasts of the Central Paratethys could provide adequate habitats for several sandfly species during the Pannonian stage, although adequate habitats could also exist in the earlier periods of the Paleogene and Neogene periods.

Conclusions
It can be concluded that Central Paratethys could be an important region both in the aspect of the speciation of the Mediterranean sandfly fauna and as a migration route between the west and east, southeast European sandfly faunas in the late Paleogene and the early and middle Neogene periods. The model results indicate that the climate of the former coasts of the Paratethys became more suitable for the ancient Mediterranean taxa in the late Miocene era. The continuously changing geographic and climatic conditions resulted in a very complex distribution history of sandfly species in the area.
Acknowledgements I acknowledge the financial support of Széchenyi 2020 under the GINOP-2.3.2-15-2016-00016 and the support of the NKFIH-471-3/2021 project. I would like to thank both reviewers for their dedicated work and suggestions, which have contributed greatly to improving the publication.
Funding Open access funding provided by University of Pannonia.
Data availability The datasets generated during the current study are available from the corresponding author on reasonable request.

Declarations
Conflict of interest The author declares that he has no conflict of interest.
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://creativecommons.org/licenses/by/4.0/.