Population genetics and phylogeography of Tabanus bromius (Diptera: Tabanidae)

Background Tabanus bromius (Diptera: Tabanidae) is one of the most notable Tabanidae species of veterinary and medical importance distributed throughout the Palearctic region. In this study, we investigate the genetic diversity and the phylogeographic structure of T. bromius sampled from Turkey, Croatia, and Iran. Methods For this purpose, a 686-base-pair (bp) fragment of mitochondrial DNA cytochrome oxidase I gene (COI) and 1339 bp of the nuclear DNA internal transcribed spacer (ITS) were sequenced from 247 individuals representing 15 populations. Results The sequences generated 169 COI haplotypes and 90 ITS alleles. A higher haplotype/allele diversity (h = 0.9909 for the COI gene and Ad = 0.8193 for the ITS region) compared to a low nucleotide diversity (π = 0.020605 for COI gene and π = 0.013667 for the ITS region), present for a high number of singleton and private haplotypes/alleles imply population expansion in the past. The results of phylogenetic analysis led to the uncovering of geographically significant groupings of lineages with regard to the entrance of the species into Anatolia and the location of major geographic barriers. According to current data, the species appears to have entered Turkey from Caucasia and Iran. A molecular clock applied to the COI data suggests that T. bromius diverged from the outgroup species nearly 8.83 million years ago, around the end of the Miocene era. Conclusions The results of this study indicate remarkable genetic diversity across the studied range of the species. High haplotype/allele versus low nucleotide diversity and demographic analyses implied that the T. bromius populations have undergone a series of expansions and retreats in the past. Our current findings suggest that T. bromius split from outgroups around the Late Miocene. Subsequent diversification events during the climatic and environmental fluctuation times of the Late Pliocene and Early Pleistocene periods also significantly influenced the species, resulting in the formation of some major genetic lineages. The phylogenetic analyses indicate that T. bromius most likely entered Turkey from the Caucasus region and Iran. Graphical Abstract Supplementary Information The online version contains supplementary material available at 10.1186/s13071-021-04970-5.

Investigating the geographic distribution of widespread pest and vector species is crucial for determining infectious diseases, vector control, and pest management [6]. To date, 171 horsefly species, classified under nine genera from three subfamilies, have been recorded from Turkey [7,8]. Tabanus bromius is an important common pest species in the Palearctic region, including Turkey. The species was first recorded from Turkey in 1957 by Moucha and Chvála [9]. Its importance is associated with both the transmission of diseases and economic significance with its appearance in large numbers, as well as its persistence. Females suck blood from animals using their stout mouthparts and can inflict painful bites during their feeding, which can irritate grazing animals considerably, resulting in allergic responses, loss of weight, and reduction in growth rates and milk production [10,11]. Despite its medical and economic significance, apart from a few faunistic studies [7,8], there has been no detailed study on T. bromius from Turkey.
Located in the Western Palearctic region, Turkey is one of the most significant places in the world due to its high biodiversity and number of endemic species [12]. The high diversity is partly associated with its continental features changing surprisingly in short geographic distances in terms of climate and habitats [13]. In this context, both species and genetic diversity are high in Turkey, which is in line with expectations [14]. Anatolia (the Asian part of Turkey) is considered a region where much hybridization is observed due to the topography [13]. It is also the homeland of many species and has been a refuge area for those species especially affected by geological and climatic changes in the past, and it has much more biological importance than many other places in the world considering its role as a natural bridge between Europe, Asia, and Africa [15]. Moreover, the presence of areas used as a refuge by terrestrial organisms in the west and east of Anatolia during the Pleistocene period makes Turkey even more important in terms of biogeography [16,17]. Climatic oscillations throughout the Quaternary period have particularly affected the Northern Hemisphere in different ways, and the response of species in different geographies has been dissimilar [18]. Nonetheless, studies over the last few decades suggest that Anatolia might have acted as a shelter for many species during the glaciation periods, and that these species may have spread to Europe or the Caucasus again (perhaps for the first time) during the post-glacial periods [16,19,20].
In recent years, there has been an increase in research involving Turkish fauna to explore the phylogenetic and phylogeographic structure of various animal groups [19][20][21][22][23][24][25]. Although there have been certain faunistic studies to determine Tabanidae diversity [7,8,[26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41], there has been only one molecular-based study on horse flies in Turkey [42]. Because it has been implicated as a vector of various pathogens, and due to the lack of any molecular-based study, here we aim to investigate the population genetic structure, phylogeography, and demographic history of T. bromius to gain greater insight into this species. Therefore, in the present study, we use cytochrome oxidase I (COI) of the mitochondrial genome and the internal transcribed spacer (ITS) region of the nuclear genome as the combination to (i) reveal the genetic diversity of T. bromius; (ii) determine the distribution of genetic variation across the sampled area; (iii) reconstruct phylogenetic trees to reveal evolutionary relationships of haplotypes/alleles; (iv) estimate divergence of lineages; and (v) reveal possible factors that generate the contemporary phylogeographic pattern in T. bromius. Our results contribute to both the fauna of Tabanidae in Turkey and the vector epidemiology.

Sample collection and DNA extraction
Tabanus bromius females were collected from the 13 locations (Antalya, Artvin, Bitlis, Çanakkale, Elazığ, Eskişehir, Giresun, Hakkari, Hatay, Kayseri, Muğla, Sinop and Thrace) covering most of the distribution range of the species in Turkey. Fresh specimens of T. bromius were caught using a malaise trap or collected manually while sucking blood from cattle. We also used museum materials registered in the inventory of the Zoology Museum of Eskişehir Technical University, Faculty of Science Department of Biology. Specimens from Iran and Croatia were also included in this study (Fig. 1). Despite all our efforts, we were not able to obtain specimens from other parts of the T. bromius distribution range from other countries. Therefore, in this study, a total of 247 samples (83 specimens were museum materials and 164 were fresh samples) representing 15 populations (234 samples from Turkey, 10 specimens from Croatia, and three individuals from Iran) were used to isolate the total DNA. All sampling locations, their abbreviations, and the sample size of the studied populations are given in Table 1. The old museum materials were preserved in ethyl acetate, which is known to cause problems in DNA extraction [43,44]. Therefore, we modified [45] the total genomic DNA isolation method. After extraction, DNA samples were run at 120 V and 50 mA in 1% agarose gel, using 1× TAE (40 mM Tris pH 7.6, 20 mM acetic acid, 1 mM EDTA) buffer. The extracted DNA samples were stored at −20 °C until DNA amplification.
Chromatograms sent by the company were first visually checked, and subsequently transferred to Geneious 11.2.1 software [49] for editing, alignment, and the collapse of similar sequences into haplotypes/alleles. To ensure that the obtained COI sequences were genuine mitogenomes in their origin, internal stop codons and nonsense mutations were carefully checked after translating the sequences into amino acids using DnaSP 5.10.1 software [50]. All haplotypes/alleles were deposited in Gen-Bank (MK941663-MK941831 for the COI haplotypes and MK936073-MK936162 for the ITS1-5.8S-ITS2 alleles). Despite all our efforts, we were not able to amplify the COI gene of 14 of the Turkish museum specimens, and we, therefore, included only 233 COI sequences and 247 ITS sequences in our data analysis. Sequences from five other tabanid species were used as outgroups: Tabanus atratus (GenBank: KM243515), T. rufofrater (GenBank: DQ631993), and T. bifarius [42] for the COI haplotypes, T. bifarius-1 and T. bifarius-2 [42] for the ITS alleles.

Estimating genetic diversity and population structure
Molecular diversity indices, including the number of polymorphic sites (S), nucleotide (π) and haplotype (h) diversity [51], the number of substitutions, and the pairwise nucleotide differences (k) [52], were estimated using the DnaSP 5.10.1 [50] and Arlequin 3.5.2.2 programs [53]. Pairwise comparisons to determine genetic differentiation between populations were estimated using the fixation index, F ST , and their significance was tested by 10,000 permutations. The hierarchical distribution of genetic diversity was calculated as an analysis of molecular variance (AMOVA) with 20,000 permutations, with statistical significance of P ≤ 0.01, using the Arlequin 3.5.2.2 program [53]. The analysis was performed at three levels: among groups, among populations within groups, and within each T. bromius population. Population demographic analyses and any deviations from neutrality were characterized using mismatch distribution analysis as implemented in DnaSP 5.10.1 [50], the sum of squared deviations (SSD), raggedness index (Hri) [54], Tajima's D [55], and Fu's F S [56]. Statistical significance was generated using 1000 simulations, and the significance level was set to P ≤ 0.05.

Genetic diversity of T. bromius
All of the sequences of the COI gene of 233 T. bromius specimens from Turkey, Iran, and Croatia generated 169 unique sequences. A total of 220 T. bromius samples from the Turkish localities yielded 157 haplotypes, 10 specimens from Croatia generated nine haplotypes, and three samples from Iran produced three distinct sequences. The haplotypes and their frequencies are shown in Additional file 2: Table S2. In 169 haplotypes, there were 117 polymorphic characters (17.0%), of which 94 (13.7%) were parsimony informative. The nucleotide frequencies were 38.0% (T), 29.5% (A), 16.7% (G), and 15.5% (C), respectively. The transition/transversion rate ratios were k1 = 5.1 (purines) and k2 = 6.0 (pyrimidines). When the COI sequences were translated into 228 amino acids, only 12 amino acids showed polymorphism. Pairwise differences between the haplotypes ranged between 1 bp (0.1%) and 39 bp (5.9%). Haplotype (h) diversity in T. bromius was h = 0.9909 ± 0.1091, and the nucleotide diversity was π = 0.020605 ± 0.012850. Among all of the COI haplotypes, seven haplotypes were shared between populations, and there were 139 singleton and 23 private haplotypes. In the Turkish populations, there were 128 singleton and 22 private haplotypes. However, in the Croatian population, there were eight singleton haplotypes, and there was only one private haplotype. Of the Iranian specimens, all three produced three unique haplotypes. No haplotypes were shared among Turkey, Croatia, and Iran.
The amplified size of the nuclear genome of T. bromius, including the ITS1 + 5.8S rRNA + ITS2 region, was 1339 bp in length. Since the 5.8S region (123 bp) showed no polymorphism, we eliminated it from our data set; therefore, only the ITS1 and ITS2 regions were combined and used to determine T. bromius alleles. Therefore, the concatenated sequence of the ITS1-ITS2 region was 1216 nucleotides in size without any length differences among 247 samples, which were collapsed into 90 unique alleles (Additional file 3: Table S3). A total of 234 Turkish T. bromius individuals yielded 84 alleles, while the 10 Croatian samples generated six and the three Iranian specimens produced two distinct alleles. In the ITS alleles, 538 characters (44.0%) showed polymorphism, of which 390 sites (32.0%) were parsimony informative. The average base composition was 43.0% (A), 35.0% (T), 12.4% (G) and 9.4% (C), and the transition to transversion ratio was R = 1.38. The pairwise difference between alleles ranged from 1 bp (0.1%) to 326 bp (3.3%). The average allele diversity was Ad = 0.8193 ± 0.1417 and the nucleotide diversity was π = 0.013667 ± 0.007385. The numbers of singleton and private alleles were high, where 63 alleles were singleton, 23 were private alleles, and only four alleles were shared among the localities. The number of singleton alleles was 59, while the number of private alleles was 21 in the Turkish populations. In the Croatian population, four alleles were determined as singleton alleles, with one private and one shared allele being found. All of the Iranian specimens produced one private and one shared allele. There was no singleton allele in the Iranian population.

Demographic history of T. bromius populations
Deviations from neutrality and population size changes were investigated through several analyses. Analyses of all of the COI haplotypes produced significantly nega-  (Table 2). Similarly, Tajima's D and Fu's F S values of individually analyzed Croatian and Iranian localities were not statistically significant. On the other hand, a mismatch distribution analysis of all of the COI haplotypes produced a bimodal profile (Fig. 2a). Individually analyzed Turkish populations showed that the Eskişehir and Iran populations generated bimodal curves; however, the remaining populations produced multimodal curves (Additional file 4: Figure S1).
For the ITS1-ITS2 data sets for all of the populations, Tajima's D value was significantly negative (Tajima's D = −1.9563, P ≤ 0.05), implying deviations from neutrality, while Fu's F S (F S = −5.5854, P ≥ 0.05) value was non-significant. Hri (0.0045) and SSD (0.0040) values calculated for all of the alleles were low. When each population was analyzed separately, Fu's F S values were non-significant for all of the populations, except for Antalya. The Antalya population also suggested population growth by a significantly positive value  (Table 2). Fu's F S and Tajima's D values calculated for the Croatian and Iranian populations were not statistically significant. Other than this, all of the ITS1-ITS2 datasets produced a multimodal curve (Fig. 2b). While the Antalya, Bitlis, Çanakkale, Muğla, and Iran populations showed bimodal curves, the remaining populations generated multimodal curves, and all these populations suggested either stable or diminished populations (Additional file 5: Figure S2).

Phylogenetic and phylogeographic analyses
All three analyses conducted for both data sets produced similar tree topologies with different bootstrap/posterior probability values. Therefore, only Bayesian consensus trees were provided here for both the COI haplotypes and ITS1-ITS2 alleles (Figs. 3,5). The resulting BEAST tree is also compatible with other trees. For the COI data set, all of the trees produced two main clades; the first clade (Clade A) comprised an Iranian haplotype (H113) and Turkish haplotypes from the eastern part of Turkey (Hakkari, Elazığ, and Bitlis). The second main clade (Clade B) splits into two subclades (B1 + B2); the subclade B1 with a northern basal haplotype (H21) sampled from the Artvin population, and most of the remaining haplotypes were from the Giresun and Artvin populations. Moreover, the subclade BI is composed of haplotypes representing eastern, central, northern, and western Turkey. The subclade B2 has a basal haplotype (H111) sampled from the Iranian population, with the remaining haplotypes from Croatia, and it is mainly dominated by the western and central Turkish haplotypes.
Application of a molecular clock showed that the ingroup haplotypes of T. bromius diverged from the outgroup haplotypes around 8.83 million years ago (MYA). The separation of the most basal haplotype (H44) from the main clade was around 6.00 MYA (Clade A). A small clade formed by Sinop haplotypes (H160, H158,   The minimum spanning network analysis supports the general pattern revealed by the phylogenetic trees of the COI gene (Fig. 4). Three main haplogroups are apparent in the network. The first shown by A in the network is coherent with Clade A in the phylogenetic trees. The haplogroup in this cluster is composed of both Iranian (H113) and predominantly eastern Anatolian haplotypes. The other haplotypes that formed subclade B1 and B2 in Fig. 4 The results of the minimum spanning network analysis of the COI haplotypes of T. bromius using HapStar 0.7. Size of the circles is proportional to the frequency of the haplotype. The branch lengths are proportional to the number of hypothetical ancestors. Haplotypes are shown in Additional file 2: Table S2 the BEAST tree are also obvious as B1 and B2 in the network. The B haplogroup is dominated by haplotypes from the central, northern, and western Anatolian populations and a small star phylogeny has H21 as the central haplotype detected from Artvin, the locality that is located in the north eastern part of Turkey. The B2 haplogroup comprises haplotypes, mostly from western and central Turkey. H48, a shared haplotype among Elazığ, Eskişehir, and Muğla with the highest frequency, is connected to other haplotypes found mainly in the western Anatolian populations.
All three phylogenetic analyses of the concatenated ITS1-ITS2 region of T. bromius generate similar trees with different bootstrap/posterior probability values; therefore, only a single tree is provided here (Fig. 5). Analysis of the ITS1-ITS2 region yields predominantly polytomous structure with relatively low resolution. The A46 allele (from Giresun), at the most basal of the tree, and the A11 allele (from Artvin) form a small clade. There is an A90 allele (from the Thrace region of Turkey) at the most basal of the clade showing polytomy. The clade is divided into two sister subclades; the first subclade is composed of Eskişehir and one Hatay (A68) allele, while the other main clade structure consists of many polytomic alleles or allele groups. Furthermore, in the main clade, the Croatian (A63-A64) and the Turkish alleles from Hakkari, Antalya, Sinop, and Eskişehir (A40-A41) are grouped together. The Iranian and Hatay alleles form a monophyletic group, and all of the remaining alleles are from the eastern, central, northern, and western parts of Turkey.
A minimum spanning network of the ITS region is given in Fig. 6. A15 (n = 40 individuals), being the most frequent and most common allele shared among the Artvin, Çanakkale, Eskişehir, Kayseri, and Sinop populations of T. bromius, produces a star-shaped structure. A22 is also a shared allele among Turkey (Çanakkale, Eskişehir, Sinop, and Thrace), Croatia, and Iran. It is remarkable that the A43 allele (n = 1 from Eskişehir) is connected to eight alleles, including the southern, northern, eastern, and western Turkish populations, and Croatia.

Partitioning of genetic variation and population differentiation
The pairwise F ST values show significant genetic differentiation for the COI data set (Additional file 6: Table S4). The highest F ST value is between Iran and Giresun (from northern Turkey F ST = 0.759, P ≤ 0.05). Iran also produces the next highest statistically significant F ST value with another northern Turkish population (from Artvin) (F ST = 0.751, P ≤ 0.05). In addition, the Iranian population is genetically different from other Turkish populations (Antalya, Çanakkale, Eskişehir, Kayseri, Hatay). On the other hand, Iran-Hakkari and Iran-Bitlis populations, which are geographically relatively close to each other, produce low F ST values (F ST = 0.343 and F ST = 0.312, P ≤ 0.05, respectively). For the Croatian population, the highest F ST is observed in comparisons between Croatia and Artvin (F ST = 0.68, P ≤ 0.05), and Giresun (F ST = 0.674, P ≤ 0.05), respectively. Pairwise comparisons also show that the highest differentiation between the Turkish populations is between Çanakkale and Artvin (F ST = 0.684, P ≤ 0.05), and between Çanakkale and Giresun (F ST = 0.670, P ≤ 0.05), respectively.
Genetic differentiation between populations, using the ITS data sets, was also employed through pairwise F ST analysis (Additional file 7: Table S5). The highest An analysis of the hierarchical distribution of genetic diversity shows that the highest genetic variation is at the "within population" level (46.0%, P ≤ 0.01), followed by "among populations" (38.3%, P ≤ 0.01) and "among groups" (15.7%, P ≤ 0.01) for the COI data set when all populations are grouped under three clusters (Group1: HRV Group2: IRN, Group3: the Turkish populations of ANT, ART, BIT, CAN, ELA, ESK, GIR, HKR, HTY, KAY, MUG, SNP, TRA). The genetic differentiation between the groups is F ST = 0.54, P ≤ 0.001. On the other hand, the ITS data set shows that 65.2% of the genetic variation is within populations and 34.7% is among populations when each locality is considered as a separate group. The genetic differentiation between the groups is F ST = 0.34, P ≤ 0.001 (Table 3).  Table 1 are considered, two northern Anatolian populations, Sinop and Giresun, show the highest diversity estimates. Moreover, the Sinop population shows the highest genetic diversity for both the COI gene and the combined ITS1-ITS2 region. Therefore, it seems that Sinop is an important location for the T. bromius genetic variation. Similar results have been observed in another Tabanus species studied from Turkey [42]. In fact, the presence of high genetic diversity has been reported from other animal groups, including insects from Turkey studied [21,[67][68][69][70], and explained by its topographic structure, the presence of microhabitats, and climatic characteristics [71]. Despite the low number of samples included from Iran and Croatia, both localities also display high diversity. However, more specimens covering a larger sampling area across the distribution range of T. bromius are necessary.

Genetic diversity and population demography of T. bromius
Among 169 T. bromius haplotypes, we detected 139 as singleton and 23 private haplotypes. There are a number of sequences that were found either only in the Iranian Fig. 6 The results of the minimum spanning network analysis of the combined ITS1-ITS2 region of T. bromius using HapStar 0.7. Size of the circles is proportional to the frequency of the allele. The branch lengths are proportional to the number of hypothetical ancestors. Allele are shown in Additional file 3: Table S3 or the Croatian populations. Although the number of samples is not equal among Iran, Croatia, and Turkey, there were 128 singleton and 22 private haplotypes in the Turkish populations. Despite the low number of specimens included from Iran and Croatia, it may be correct to suggest that Turkey has the conspicuous level of variation across the sampled area, and that the species has been in Turkey for a prolonged time period, so that it possessed the high number of variants. The presence of both high singleton and high private haplotype numbers imply that these haplotypes have been derived relatively recently. In addition, the high rate of derived haplotypes may be associated with the presence of young sequences in the data set of expanding populations [72]. In T. bromius the nucleotide diversity is low compared to the haplotype diversity for both COI and ITS data for the Turkish populations. It is known that the high haplotype, but low nucleotide diversity, and the private-singleton haplotypes suggest that the population has shown fluctuations in the past, and that it generally points to population growth in the more recent past following population bottlenecks [18].
The frequency and geographic distribution of T. bromius haplotypes/alleles imply the possible effects of certain biogeographic barriers in Turkey. For instance, the allele with the highest frequency and the most common allele (the A15 allele shared among Artvin, Çanakkale, Eskişehir, Kayseri, and Sinop in Turkey) is distributed among the Black Sea region, western, and central Anatolia. The presence of the A15 allele only in the western part of a major mountain range known as the Anatolian Diagonal dividing Anatolia into East and West [73] suggests that T. bromius might have also been affected by this barrier. Moreover, the absence of shared haplotypes/ alleles between the eastern and western sides of the mountain line supports the importance of the Anatolian Diagonal on T. bromius. In fact, it has been shown that the Anatolian Diagonal plays a key role in the distribution of species and genetic lineages, not only in plants [74,75], but also in animals [23,25,68,[76][77][78]. A geographic structuring across the distribution range of T. bromius in Turkey is further supported by the pairwise F ST analysis predicting the levels of genetic differentiation of populations and reveals significant differences from each other in many populations. Iran is one of the entrance doors into Turkey of T. bromius and the population representing this region is quite different from the Black Sea populations (Giresun and Artvin), while they differ slightly from the eastern populations (Hakkari and Bitlis). Another entrance gate, the Black Sea region (Artvin), has an extremely high differentiation rate with its population of Çanakkale and Croatia (representing the European region) proposing the presence of geographic separation of lineages after entering into Turkey [22,68,71].
The effects of historical factors can be traced by employing several population demographic analyses [18]. For T. bromius, the mismatch analysis of the COI data, where all of the populations included in the analysis generated a bimodal profile with Hri = 0.0011 and SSD = 0.0038, but were not statistically significant, indicates population expansion [79]. Fu's F S (−23.7386) is significantly negative, indicating an excess of rare haplotypes over what would be expected under neutrality in an expanding population is also implied by the mismatch distribution. Negative and non-significant Tajima's D imply that populations of T. bromius did not deviate from neutrality. In our results, although Tajima's D and Fu's F S value seem to conflict, the reason is different statistical power, because Fu's F S is more powerful and more sensitive than Tajima's D to population size changes [80]. On the other hand, the combined ITS1-ITS2 region Table 3 AMOVA analysis of Tabanus bromius. Analysis of the COI data includes three groups: Group 1 (HRV), Group 2 (IRN) and Group 3 (ANT, ART, BIT, CAN, ELA, ESK, GIR, HKR, HTY, KAY, MUG, SNP, TRA). In the analysis of the ITS1-ITS2 data, each population was treated as a distinct group. Population abbreviations are shown in Table 1 Source generated a multimodal profile with Hri = 0.0045 and SSD = 0.0040, but it was non-significant, indicating population expansion. For the ITS1-ITS2 region, it can be seen that the species has a stable population distribution. This is also supported by the negative significant Tajima's D = −1.9563 value implying population expansion. Overall, the mismatch distribution analysis and other population demographic analyses, in addition to the combination of high haplotype and low nucleotide diversity of T. bromius populations for both the COI gene and the combined ITS1-ITS2 regions, suggest rapid demographic expansion from a small population in the past [81].

Estimation of divergence times of T. bromius lineages
Horseflies are a cosmopolitan group spread across all continents and are reported to be of Pangea origin. It is thought that they emerged in the late Triassic period just before the drifting of the Laurasia and Gondwana continents [82]. Based on morphological and molecular evidence, the oldest species of a true tabanid is thought to have derived during the late Jura [83], and the emergence of the tabanid species coincided with the early Cretaceous period [84]. The Tabanidae species is rich in the equatorial belt, both in terms of species diversity and population density, but it gradually decreases in number towards the poles [10] because the adult activity period of the species increases depending on climate temperature. The majority of the terrestrial fauna of the Palearctic region including Turkey is thought to have originated from Laurasian fauna [13]. Based on the distribution of the Tabanidae species in the Palearctic region, the Turkish Tabanidae species is reportedly composed of Mediterranean, European, Asian, Eurasian, Ethiopian and endemic species [85][86][87][88][89][90]. This view is based only on the geographic distribution of the T. bromius species in the Palearctic region, and is not focused on any assessment of phylogeography. Our study species, T. bromius, shows a cosmopolitan distribution in the Palearctic region. The species has been recorded from all European countries [91] (except Ireland), North Africa (Algeria and Morocco), the Near East and a number of Middle Eastern countries, Central Asian countries (Kazakhstan, Afghanistan) up to the Northern Urals [10], India [92] and Turkey [7,8,91,93]. When the distribution of both the T. bromius and other Tabanus species [1,10] in the Palearctic region is taken into consideration, it has been estimated that the origin of the group may have been Central or West Asia, and that they might have entered Europe from north of the Black Sea (the Caucasus region), and from Turkey in the south.
Our current findings indicate that T. bromius diverged from outgroups approximately 8.83 MYA around the Tortonian period of the Late Miocene. During the Miocene, more than a quarter of the world was covered with pasture, which may have given the tabanids an advantage in terms of their feeding and egg laying. It is estimated that while this increased the number of warm-blooded animals feeding on this type of pasture, it also increased the density of female tabanids feeding on the blood of these animals. The Miocene was a period of intense climatic and environmental change; therefore, it is possible that the great fluctuations that occurred during the Late Miocene period may have resulted in the separation of T. bromius from the outgroup species. Indeed, concurrent findings have also been suggested in a number of other studies for the separation of species throughout the Miocene period [2]. It seems that after the splitting of the T. bromius lineage from the ancestral stock species, it has undergone certain major diversification events. Our analysis implies that the estimated age of the most basal haplotype of T. bromius (H44 from Çanakkale) is around 6.00 MYA, coinciding with the end of the Miocene period. In the ingroup, an eastern haplotype from Elazığ (H54) split from its sister clade almost 4.30 MYA around the Early Pliocene period. Similarly, the estimated divergence time of the two main clades (Clade A and Clade B) was approximately 3.10 MYA, which indicates the Late Pliocene period, and ongoing diversification events have continued from the beginning of the Pleistocene period until the present day. This coincides with the effects of glacial and interglacial changes in the Quaternary (in Pleistocene) [17,94,95]. Quaternary climatic oscillations caused the expansion or contraction of species that played an important role in shaping the genetic models of populations [96]. Such diminishing/expansion and the isolation following multiple glacial refugia can maintain high intra-genetic diversity [97]. These models can be well illustrated by a phylogeographic analysis with hereditary data in and among the population [98]. It appears that the Quaternary climate oscillations were effective in shaping the populations and lineages of the T. bromius species.
Our phylogenetic analyses of T. bromius suggest that the species most likely entered Anatolia through both the Caucasus region and Iran. The lineage that entered through Iran has shaped the eastern populations (Elazığ, Bitlis, and Hakkari) up to the Anatolian Diagonal, and the lineage that entered through the Caucasus spread westward through the Black Sea region and crossed the Anatolian Diagonal. Thereafter, Central Anatolia (Kayseri, Eskişehir), south of Turkey (Hatay and Antalya), and Western Anatolian populations (Muğla, Çanakkale and Thrace) were relatively separated in the western part of Anatolia. Furthermore, the placement of the haplogroup, consisting of the Thracian and Croatian populations in the same clade, suggests that the species provided a genetic source to the Croatian population. This suggests that the Central European T. bromius populations may have originated from Thracian populations or that a portion of the T. bromius European populations originated from the lineage from Turkey. The Black Sea region populations (Artvin and Giresun) are in subclade B1; it may be correct to assume that the collected samples are related to the distribution in a narrow area north of the North Anatolian mountain range. This is also supported by the co-grouping of Black Sea populations in one of the three-star phylogenies in the haplotype network. The Sinop population haplotypes present an interesting pattern due to phylogenetic trees, suggesting that the Sinop population is interacting with other populations. Similar results have been observed for T. bifarius [42]. Antalya and Hatay haplotypes formed subclade B1 together with the Central Anatolian haplotypes. In particular, all Hatay haplotypes are grouped together with Kayseri (H117) and Sinop (H148) haplotypes and, as a consequence, it has been proposed that the Hatay population might have originated from Turkey. The absence of the species in countries such as Lebanon and Jordan, close to the Hatay border, and the fact that the lineage relationship with the Hakkari population is not seen in the phylogenetic trees, may support this idea. The ITS1-ITS2 analysis produced similar phylogenetic tree topologies. There is also an allele structure with geographically significant groupings among the polytomic groups and Artvin and Giresun alleles that are found at the basal part of all tree topologies. In addition, due to the presence of the T. bromius species in countries such as Russia and Ukraine, it may be assumed that a lineage in the north of the Black Sea (from the Caucasian), which was a freshwater lake at that time, reached Europe. Moreover, these results agree with our results that the Turkish T. bromius populations originated in the Caucasus region (Artvin) and Iran.

Conclusions
In conclusion, the T. bromius populations show quite remarkable genetic diversity. Demographic analyses and high haplotype/allele versus low nucleotide diversity imply that the T. bromius populations might have undergone a series of expansions in the past. The T. bromius haplotypes split from outgroup species around the Miocene period, and current results clearly indicate that the structuring of T. bromius genetic lineages was significantly influenced by climatic and environmental fluctuations during the Late Pliocene and Early Pleistocene periods. The phylogenetic analysis findings show that the species appears to have entered Turkey from Caucasia and Iran, and that the presence of geographic barriers, such as the Anatolian Diagonal, separates geographically important lineage groupings. Larger sampling across the western Palearctic region and a more detailed study are necessary to draw a general conclusion regarding the phylogeographic structure of T. bromius.