Females of T. bifarius were collected from 12 locations across Türkiye, covering most of the species’ known distribution range. Fresh specimens were captured using Malaise traps, and additional samples were obtained from the Zoology Museum of Eskişehir Technical University’s Faculty of Science, Department of Biology (ESTU-ZM). Specimens from Iran and Croatia were also included in the study (Fig. 1).

The specimens used in this study were identified by Prof. Dr. A. Yavuz Kılıç based on morphological keys (Chvála et al., 1972; Kılıç, 1999, 2006). Females are medium-sized horseflies, with body lengths ranging from 12 to 15 mm. Greyish to olive-grey in color, with a relatively narrow frons and well-separated calli; the lower callus is brownish. Eyes haired and with three bands. Antennae are reddish-brown with a dark apex. The thorax exhibits dark brown coloration with distinct yellowish stripes, while the abdomen shows alternating dark and light bands. Notable variability is observed in antennal coloration and abdominal patterns.

In total, genomic DNA was successfully extracted from 187 individuals, including 146 fresh specimens and 41 museum samples, representing 13 populations (182 from Türkiye and 5 from Iran), using a modified version of the extraction protocol by Taylor et al. (1996). Sampling locations, along with geographic coordinates and sample sizes, are provided in Table 1. DNA could not be extracted from the Croatian museum samples, probably due to DNA degradation caused by Ethyl acetate preservation (Dillon et al., 1996; Quicke et al., 1999).

A 650-base pair (bp) fragment of the mitochondrial cytochrome c oxidase subunit I (COI) gene was amplified using the LCO1490–HCO2198 primer pair (Folmer et al., 1994). Additionally, a 1339 bp region spanning the internal transcribed spacer (ITS1–5.8S rRNA–ITS2) of nuclear DNA was amplified using three primer pairs: CS249–FL (Tautz et al., 1988), CAS18sF1–CAS5p8sB1d, and CAS5p8sFt–CAS28sB1d (Ji et al., 2003). PCR amplification, product purification, and sequencing followed protocols described by Sanal-Demirci et al. (2021). Despite extensive attempts, amplification of the ITS region failed in 31 Turkish museum specimens; thus, the final dataset included 187 COI sequences and 156 ITS sequences. All sequences were aligned and curated using Geneious v11.2.1 (Kearse et al., 2012). To avoid including nuclear mitochondrial pseudogenes (numts), COI sequences were translated into amino acids and checked for the presence of stop codons or nonsense mutations using DnaSP v5.10.1 (Librado et al., 2009). All unique haplotypes and alleles have been deposited in GenBank under accession numbers OP934089–OP934169 (COI) and OQ550299–OQ550357 (ITS1–5.8S–ITS2).

The number of polymorphic sites (S), haplotype (h) and nucleotide (π) diversity (Nei 1987), the number of substitutions, and number of pairwise nucleotide differences (k) (Tajima 1983) for both the COI gene and the ITS region were calculated using the DnaSP 5.10.1 and Arlequin 3.5.2.2 programs (Excoffier et al. 2010). DnaSP 5.10.1 was used to do mismatch distribution analysis using pairwise differences between haplotypes/alleles. Expanding populations typically exhibit a unimodal pattern in mismatch distribution analysis, whereas populations in demographic equilibrium display bimodal or multimodal patterns (Rogers and Harpending 1992). The raggedness index (Hri) (Harpending 1994) and the sum of squared deviations (SSD) (Schneider and Excoffier 1999) were also calculated. Besides, Tajima’s D (Tajima 1989) and Fu’s F_S (Fu 1997) tests to check for any deviations from neutrality were performed in Arlequin 3.5.2.2. Statistical significance was generated using 1000 simulations, and the significance level was set to P ≤ 0.05. 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. Genetic diversity was assessed at three hierarchical levels (among groups, among populations within groups, and within populations) using Analysis of Molecular Variance (AMOVA) with 20,000 permutations, with statistical significance of P ≤ 0.01, using the Arlequin 3.5.2.2 program (Excoffier et al. 2010).

Both data sets were analyzed separately to reveal intraspecific phylogenetic relationships among Tabanus bifarius populations. Tabanus maculicornis Zetterstedt, 1842 (GenBank: MZ563331.1), Tabanus shannonellus Kröber, 1936 (GenBank: MZ563357.1) for the COI haplotypes, T. bromius-1 (GenBank: MK936073) and T. bromius-2 (GenBank: MK936162) for the ITS alleles were used as outgroups in all phylogenetic analyses. Maximum parsimony (MP), maximum likelihood (ML) and Bayesian analyses were used to deduce phylogenetic relationships among haplotypes/alleles using PAUP*4.0b10 (Swofford 2002) and MrBayes 3.2.6 (Ronquist and Huelsenbeck 2003), respectively. MP analysis was carried out using the TBR branch-swapping algorithm of random addition of taxa under heuristic search options and was employed with 1000 bootstrap replicates. The JModeltest-2 (Darriba 2012; Guindon and Gascuel 2003) was used to find the best-fit model of substitution to our data sets. The GTR + I + G and HKY + I + G were the best fitting models for the COI gene (AICc = 48,009) and the ITS region (AICc = 65,388), respectively. For the estimation of divergence times of the COI lineages and their confidence intervals, BEAST ver. 1.5.2 was used using a Bayesian Markov chain Monte Carlo (MCMC) method coalescent method by applying GTR + I + G model following the uncorrelated relaxed lognormal clock (Drummond and Rambaut 2007). A 2.3% sequence divergence per lineage per million years mutation rate was used for calibrating molecular clock (Papadopoulou et al. 2010). Both the most recent common ancestors (MRCAs) and most ancient common ancestors (MACAs) calculation and operator optimizations were carried out with BEAUTI v1.8.0 (Hayward and Stone 2006). The BEAST analysis was run for 100 million generations, sampling every 1000 and the convergence to stationary and the effective sample size (ESS) of the model parameters were checked using Tracer v1.6.0. The maximum clade credibility tree was built with TreeAnnotator v1.8.4, discarding the initial 25% samples as burn-in. The results were visualized using FigTree ver. 1.3.1 (Rambout 2009). Since phylogenetic trees may not always be adequate to distinguish genealogical evolutionary relationships because of reticulations (Posada and Crandall 2001), we used network analysis to analyze our data. The pairwise distance matrix between haplotypes/alleles computed by Arlequin 3.5.2.2 was transferred to HapStar version 0.5 (C) (Teacher and Griffiths 2011) to build a minimum spanning tree.

A total of 81 unique COI haplotypes were identified from 187 T. bifarius individuals sampled from Türkiye and Iran. Seventy-seven haplotypes were detected across all Turkish populations, while the five Iranian specimens yielded four distinct haplotypes (Table SS1). Among the 81 haplotypes, 87 polymorphic sites were identified, 68 of which were parsimony informative. Nucleotide composition was biased toward adenine and thymine, with frequencies of 37.8% thymine (T), 30.0% adenine (A), 16.9% guanine (G), and 15.1% cytosine (C), consistent with typical mitochondrial DNA. No nonsense mutations or internal stop codons were observed in the translated amino acid sequences, supporting the mitochondrial origin of the data. Translation of the COI sequences into 216 amino acids revealed polymorphism at only seven positions. Pairwise nucleotide differences among haplotypes ranged from 1 bp (0.1%) to 37 bp (6.0%). Overall haplotype diversity was high (h = 0.9608 ± 0.0080), and nucleotide diversity was moderate (π = 0.026804 ± 0.013240). Of the 81 haplotypes, six were shared across multiple populations, 22 were private (restricted to a single population), and 53 were singletons (found in only one individual). Notably, all Iranian haplotypes were shared exclusively with the Van population in eastern Türkiye, suggesting a possible historical or contemporary connection between these regions.

The amplified ITS1–5.8S rRNA–ITS2 region from the nuclear genome of T. bifarius totaled 1,339 base pairs (bp). The 5.8S region (123 bp) was invariant across all samples and was therefore excluded from subsequent analyses. Only the ITS1 and ITS2 regions were concatenated and used for allele identification. The resulting combined ITS1–ITS2 sequence was 1,216 bp in length and showed no length variation among the 156 successfully amplified samples. These sequences were collapsed into 59 unique alleles. Among the 1,216 bp sequences, 130 sites were polymorphic, with 127 of these being parsimony-informative. The overall nucleotide composition was 42.14% adenine (A), 35.07% thymine (T), 13.95% guanine (G), and 8.83% cytosine (C). Pairwise differences between alleles ranged from 1 bp to 106 bp. Allele diversity was high (Ad = 0.8883 ± 0.0210), and nucleotide diversity was π = 0.023355 ± 0.011370. Most alleles were unique: 45 were singletons, 6 were private, and only 8 alleles were shared among multiple populations. In the Iranian population, two singleton alleles were identified, along with one private allele and one shared allele (Table SS2).

Neutrality tests were performed to assess deviations from neutrality and potential changes in population size. Fu’s F_S value for all COI haplotype analyses was significantly negative (F_S = −23.4979, P ≤ 0.05), suggesting population expansion. In contrast, Tajima’s D, however, was not significant (Tajima’s D = −0.5032, P ≥ 0.05). Both Harpending’s raggedness index Hri = 0.0036 and SSD = 0.0059 were low, further supporting a recent demographic expansion. When analyzed individually, several Turkish populations such as Çanakkale and Edirne showed signals of past demographic expansion, whereas the remaining populations did not yield statistically significant results. For both the Ardahan and Çanakkale populations, Tajima’s D value was significantly negative (Tajima’s D = −2.0832, P ≤ 0.05; Tajima’s D = −1.6916, P ≤ 0.05, respectively), suggesting deviations from neutrality, possibly due to population expansion (Table 2). In addition, the mismatch distribution analysis of all COI haplotypes revealed a multimodal pattern (Fig. 2).

When analyzed individually, the Çanakkale and Edirne populations exhibited unimodal mismatch distribution curves, indicative of recent population expansion. On the contrary, the Ardahan, Denizli, and Erzincan populations displayed bimodal curves, while the remaining populations showed multimodal patterns. These latter distributions suggest demographic stability or population decline (Fig. S1). In the ITS1–ITS2 region data for all the populations, neither Fu’s F_S nor Tajima’s D value was a significant value (F_S = −1.5243, P ≥ 0.05; Tajima’s D = −0.7292, P ≥ 0.05) while Hri (0.0121) and SSD (0.0152) values calculated for all the alleles were low. When each population was analyzed separately, Fu’s F_S values were non-significant for all populations. Only the Ardahan population generated a significant negative value in Tajima’s D (−1.5900, P ≤ 0.05). The remaining populations produced no significant values (Table 2). Mismatch distribution analysis of all the ITS1–ITS2 datasets produced a multimodal curve, indicating that all populations are likely demographically stable or diminished populations (Fig. 2). When analyzed individually, the Çanakkale, Denizli, Edirne, Muğla and Van populations exhibited unimodal mismatch distribution curves, whereas the remaining populations produced multimodal patterns (Fig. S2). Moreover, no polymorphism was observed in the Erzincan population.

The pairwise F_ST values showed significant genetic differentiation for the COI data set (Table SS3). The highest FST value was observed between Çanakkale and Ardahan (F_ST = 0.9484, P ≤ 0.05), followed by Iran (F_ST = 0.9434, P ≤ 0.05) and Van (F_ST = 0.9304, P ≤ 0.05), respectively. The next highest statistically significant F_ST values involving another population were observed between Edirne and Ardahan (F_ST = 0.9361, P ≤ 0.05), Iran (F_ST = 0.9167, P ≤ 0.05), and Van (F_ST = 0.9101, P ≤ 0.05). The lowest F_ST value involving the Iranian population was with the Van population (F_ST = 0.1318, P ≥ 0.05), even though it was statistically insignificant. Furthermore, for the Iranian population, the highest level of genetic differentiation was observed with the Çanakkale population. Pairwise F_ST analysis using the ITS1–ITS2 dataset supported the presence of considerable genetic differentiation among populations (Table S4). The highest F_ST values were observed in comparisons between Denizli and Erzincan (F_ST = 0.965, P ≤ 0.05), followed by Denizli and Muğla (F_ST = 0.958, P ≤ 0.05), and Denizli and Çanakkale (F_ST = 0.929, P ≤ 0.05). In contrast, the lowest levels of genetic differentiation were found between Antalya and Iran (F_ST = 0.0146, P ≥ 0.05) and Erzincan and Muğla (F_ST = 0.0161, P ≥ 0.05), all of which were statistically insignificant. For the Iranian population, the highest differentiation was detected between Iran and Ardahan (F_ST = 0.808, P ≤ 0.05), whereas the lowest value was found between Iran and Antalya (F_ST = 0.0146, P ≥ 0.05).

An analysis of the hierarchical distribution of genetic diversity based on the COI dataset revealed that the highest level of genetic variation occurred within populations (30.61%, P ≤ 0.01), followed by variation among populations (39.72%, P ≤ 0.01) and among groups (29.68%, P ≤ 0.01), when all populations were divided into two groups (Group 1: IRN, VAN; Group 2: the other Turkish populations. The genetic differentiation between these groups was high (F_ST = 0.69, P ≤ 0.001; Table 3).

On the contrary, analysis based on the ITS1–ITS2 dataset showed that 30.5% of the genetic variation was within populations, 4.93% among populations, and 64.49% among groups, when all populations were divided into three groups (Group 1: ANT, DNZ, MUG; Group 2: ARD, CAN, EDR, ERZ, ESK, HTY, KAY, SNP; Group 3: IRN, VAN). The genetic differentiation between these groups was also high (F_ST = 0.69, P ≤ 0.01) (Table 3).

The ML, MP, and BEAST tree topologies were largely congruent, differing only in bootstrap or posterior probability values at certain branches. Consequently, the BEAST tree for the COI dataset and the ML tree for the ITS dataset are presented here. Molecular analyses based on the COI gene suggest that T. bifarius diverged from the outgroup taxa approximately 5.52 million years ago (MYA) (Fig. 3).

Subsequently, T. bifarius split into two distinct haplogroups, Clades A and B, around 3.63 MYA. In Clade A, a basal haplotype from the Eskişehir population (H45), located in central Anatolia, diverged from its sister group approximately 1.98 MYA, forming a monophyletic group with the remaining haplotypes. This sister group later split into two subclades, both primarily composed of haplotypes from northern Anatolia. Clade B, the second major lineage, is divided into two subclades: B1 and B2. Their divergence is estimated to have occurred around 2.5 MYA. Subclade B1 is characterized by a basal haplotype from the Kayseri population, which diverged approximately 1.25 MYA. The remaining haplotypes in this subclade are mainly from the Hatay and Antalya populations. Subclade B2 shows a more structured organization and is dominated by haplotypes from western and central Anatolia. It is further divided into two geographically distinct subclades: B2-A and B2-B, which split around 1.84 MYA. B2-A includes basal haplotypes (H47, H48, and H49) from the Hatay population, which diverged approximately 1.33 MYA to form a minor distinct clade. The remaining haplotypes within B2-A, originating from Erzincan, Van, and Iran, formed a separate lineage around 1.09 MYA. Despite this diversity, B2-A remains primarily composed of haplotypes from western and central Türkiye, with diversification dating to around 0.7 MYA. B2-B features a basal haplotype (H40) from the Eskişehir population, which diverged approximately 0.92 MYA. This subclade subsequently split around 0.54 MYA into two clades, each consisting of smaller groups of haplotypes predominantly from western Anatolia.

The general pattern observed in the COI gene phylogenetic trees is further supported by the results of the minimum spanning network analysis (Fig. 4), which reveals four main haplogroups.

The first group, labeled as A in the network, corresponds to Clade A in the phylogenetic trees and includes haplotypes from central, northern, and northeastern Anatolia. The remaining haplotypes, forming subclades B1 and B2 (including B2-A and B2-B) in the BEAST tree, are also clearly represented in the network as distinct clusters labeled B1, B2-A, and B2-B. The B1 haplogroup is predominantly composed of haplotypes from central and southern Anatolia. A small star-like phylogeny is observed within this group, with haplotype H62, originating from Kayseri from the central Anatolia. The B2-A haplogroup primarily contains haplotypes from western and central Türkiye, along with another cluster that includes both Iranian and eastern Anatolian (Van, Erzincan) haplotypes. In the B2-B subgroup, haplotypes H25 and H4, shared and among the most frequent, serve as central nodes connected to other haplotypes mostly found in western Anatolia. Several small star-like phylogenies are also evident within this group, which aligns closely with Clade B2-B in the phylogenetic tree.

Both phylogenetic analyses (ML and MP) of the concatenated ITS1–ITS2 region of T. bifarius produced similar topologies, differing only in bootstrap and posterior probability values. Therefore, only the ML tree is presented here, with BI posterior probabilities on the corresponding branches. All branches in the tree are strongly supported (Fig. 5).

Analysis of the ITS1–ITS2 region revealed a predominantly polytomous tree structure, indicating relatively low phylogenetic resolution. The A41 allele from the Kayseri population occupies the most basal position in the tree. The topology is divided into two sister clades. The first clade includes a small group formed by the A16, A15, and A14 alleles from the Denizli population, as well as a larger subclade comprising the A26, A28, A27, A32, A31, and A29 alleles, all sampled from Hatay. The second is characterized by a basal allele from the Eskişehir population and includes additional alleles from Eskişehir, Antalya, and Sinop, suggesting a distinct lineage. Within the broader polytomous structure, the A1 allele (from Antalya) also occupies a basal position. The remainder of the tree is composed of numerous unresolved branches, consisting of alleles or allele groups from populations in eastern, central, northern, and western Türkiye, as well as from Iran.

The minimum spanning network based on the ITS region is presented in Fig. 6.

The most frequent and widely distributed allele, A2 (found in 48 individuals), occupies a central position in a star-like configuration, suggesting it may represent an ancestral or founding allele. A2 is shared among six geographically diverse populations: Antalya, Çanakkale, Edirne, Erzincan, Eskişehir, and Muğla. Another widely shared allele, A4, is present in individuals from the Antalya, Eskişehir, and Sinop populations. Notably, a distinct connection pattern is evident between the Iranian population and the Van population in eastern Türkiye. This pattern may indicate historical gene flow or shared ancestry between these geographically separated populations, supporting the hypothesis of east–west dispersal or connectivity across the region.

This study provides the first comprehensive analysis of mitochondrial and nuclear genetic variation in T. bifarius populations from Türkiye and Iran. The high haplotype diversity (h = 0.9608) and nucleotide diversity (π = 0.0268) observed in the COI gene suggest a long evolutionary history, potentially involving multiple glacial refugia and complex post-glacial recolonization pathways in the region. The A+T content in the COI gene T. bifarius was found to be 67.72%, which falls within the typical range reported for mitochondrial coding genes in many insect groups. Additionally, amino acid polymorphism in the COI gene was relatively low, with only 7 out of 216 residues exhibiting variation. This pattern is commonly observed in mitochondrial genes of insects and other metazoans, largely due to the essential role of mitochondrial-encoded proteins in oxidative phosphorylation. (Clary and Wolstenholme 1985; Crozier and Crozier 1993; Mutun and Atay, 2015). Furthermore, translation of the haplotype sequences into amino acids revealed no stop codons or nonsense mutations. Despite the large number of COI haplotypes identified (N = 81), only six were shared among populations, while a substantial proportion consisted of singletons (53) or private haplotypes (22). Most of these unique haplotypes were detected in Central and Western Anatolian populations. This pattern supports a model of restricted gene flow and pronounced population structure, likely shaped by geographical barriers and historical climatic fluctuations. The high number of private and singleton haplotypes suggests that these haplotypes have been derived relatively recently. Moreover, the presence of such young sequences within expanding populations may reflect the high incidence of derived haplotypes typically associated with recent population growth or range expansion (Crandall and Templeton 1993). In the Iranian populations, only four shared haplotypes were identified, with no singletons or private haplotypes detected. The presence of shared haplotypes between eastern Türkiye (Van population) and Iran suggests recent gene flow between these geographically proximate regions. As observed in many other animal taxa, the high genetic diversity in T. bifarius is consistent with expectations for a species inhabiting Türkiye’s topographically complex and climatically diverse landscape (Dubey et al. 2007; Çınar-Kul 2010; Ansel et al. 2011; Budak et al. 2011; Dinç and Mutun 2011, 2019; Çimen and Mutun 2022; Mutun and Danso 2023). The nuclear ITS region exhibited a similarly high level of genetic variation, with an average diversity (Ad) of 0.8883 and nucleotide diversity (π) of 0.0233. Among 156 individuals analyzed, 59 distinct alleles were identified. Of these, 45 were singletons and 6 were private alleles, while only 8 alleles were shared among individuals, further supporting the patterns observed in the mitochondrial data. Notably, all individuals from the Erzincan population carried unique singleton alleles, suggesting that this population is both young and isolated (Crandall and Templeton, 1993). The prevalence of singleton alleles is indicative of population expansion or recent divergence, a pattern frequently reported in species with ecological plasticity that inhabit fragmented landscapes (Avise, 2000). The concordance between the COI and ITS datasets supports a demographic scenario of recent population expansion.

In the present study, the highest haplotype diversity was observed in the Eskişehir population, whereas the greatest nuclear allelic diversity was detected in Hatay. In contrast, our previous research on T. bromius revealed that both mitochondrial and nuclear diversity peaked in the Sinop population, located along the northern Black Sea coast. These distinct genetic diversity patterns in T. bifarius and T. bromius provide valuable insights into their respective adaptations to the historical and ecological dynamics of Anatolia. Eskişehir, in central-western Türkiye, may act as a long-term refugial area or secondary contact zone for T. bifarius, reinforcing remarkable haplotype diversity. Nuclear allele diversity was highest in Hatay, a biogeographically transitional region connecting Anatolia to the Levant, potentially acting as a corridor for gene flow and historical dispersal (Çıplak 2008; Bilgin 2011). While T. bromius shows a stronger affinity for Palearctic temperate zones and likely persisted in coastal refugia such as Sinop during glacial cycles, T. bifarius, which is more adapted to Mediterranean conditions, appears to have diversified primarily in interior or southern refugial areas. These findings indicate that, despite belonging to the same genus, the two species exhibit markedly different phylogeographic patterns, likely shaped by their distinct ecological requirements and historical dispersal dynamics (Hewitt, 2000; Schmitt, 2007; Çıplak, 2008).

The impact of historical factors can be critically assessed through a range of demographic analyses of population data (Avise, 2004). For T. bifarius, Fu’s F_S value (−23.4979) is significantly negative, indicating an excess of rare haplotypes, which suggests deviation from neutrality and is consistent with population expansion. This inference is further supported by the mismatch distribution analysis of the COI data, where all populations exhibited a multimodal profile (Hri = 0.0036; SSD = 0.0059). Although these values were not statistically significant, the overall pattern aligns with a demographic history of expansion (Xiao et al. 2011). Although Tajima’s D was not significant at the overall level, some individual populations such as those from Ardahan and Çanakkale exhibited significantly negative D values, suggesting localized population expansion or deviations from neutral evolution. Mismatch distribution analyses further supported these demographic trends. Although a multimodal haplotype distribution was observed in several populations, Çanakkale and Edirne exhibited unimodal curves, characteristic of recent population expansion. In contrast, the bimodal profiles observed in populations from Ardahan, Denizli, and Erzincan may reflect historical bottlenecks. Neutrality tests based on the nuclear marker were not significant overall (Fu’s F_S = −1.5243; Tajima’s D = −0.7292); however, the Çanakkale population again showed a significantly negative Tajima’s D (−1.5900, P ≤ 0.05), reinforcing the evidence for localized population growth. The mismatch distribution derived from the ITS dataset was also multimodal overall, suggesting that the nuclear genome may retain more stable genetic signals over time or reflect deeper historical population structure not captured by the more rapidly evolving mitochondrial genome. Notably, several populations (Çanakkale, Edirne, Denizli, Muğla, and Van) exhibited unimodal patterns, indicative of recent demographic expansions. When considered together, demographic analyses of both the mitochondrial COI gene and the nuclear ITS region suggest recent population expansion in certain western populations, while others appear genetically stable, potentially due to historical isolation or environmental barriers restricting gene flow. Furthermore, both pairwise FST values and hierarchical AMOVA results support significant geographic structuring of T. bifarius populations across Türkiye. These findings reveal varying levels of genetic differentiation: some populations exhibit strong genetic structure, whereas others particularly those with low and non-significant F_ST values likely experience historical or ongoing gene flow, possibly due to geographic proximity or shared evolutionary histories (Slatkin 1993; Hutchison and Templeton, 1999).

When the results obtained for T. bifarius in the current study are compared with findings from our previous work, T. bromius displays partially similar yet more structured demographic patterns. This increased structuring may reflect a more complex and regionally heterogeneous demographic history, potentially linked to T. bromius’s broader distribution across the Palearctic region and greater sensitivity to climatic fluctuations during the Pleistocene (Hewitt, 2000; Schmitt, 2007). In contrast, T. bifarius exhibits more recent or incomplete signals of demographic change, likely shaped by its Mediterranean distribution and associated ecological constraints.

Our previous phylogenetic analyses indicated that T. bromius diverged from its outgroup taxa approximately 8.83 MYA, during the Late Miocene (Sanal-Demirci et al. 2021). This period was characterized by dramatic climatic shifts and the widespread expansion of pasture-dominated habitats, highlighting the evolutionary significance of T. bromius in adapting to dynamic environmental conditions. The Miocene epoch witnessed substantial climatic and ecological fluctuations, including the proliferation of open habitats such as grasslands. These environmental changes likely created new ecological opportunities for Tabanid flies, including increased host availability and more favorable breeding conditions. Our current research estimates that T. bifarius diverged from outgroup taxa approximately 5.52 million years ago, during the Messinian period of the Late Miocene. This finding aligns closely with our previous results and strongly supports the hypothesis that the significant diversification within the Tabanus lineage was predominantly driven by transformative environmental events occurring during the Miocene. The consistent results observed for both T. bromius and T. bifarius underscore the pivotal role of the Miocene in the early diversification of the Tabanus genus. Additionally, Morita et al. (2016) have proposed that the divergence of several Tabanidae species occurred during the Miocene. These divergence events correspond with significant ecological shifts throughout the Miocene and early Pliocene, suggesting that environmental instability and the emergence of new geographic barriers were critical drivers of lineage separation and population differentiation in T. bifarius and related taxa (Labandeira 2010). Our analyses suggest that the most ancestral haplotype of T. bifarius (H44 from Eskişehir) originated approximately 1.98 MYA, corresponding to the Early Pleistocene. The estimated divergence time between the two primary clades (Clade A and Clade B) is approximately 3.63 MYA, placing it in the Late Pliocene. These divergence and subsequent diversification events have persisted from the Early Pleistocene through to the present. This timeline underscores the significant climatic fluctuations associated with the glacial and interglacial cycles of the Quaternary period (Bryson 1974; Hewitt 1996; Petit et al. 2005; Çıplak 2008). Climatic oscillations during the Quaternary period drove cycles of species contraction and expansion, profoundly influencing the genetic structure of populations (Hewitt 2001). Isolation in multiple glacial refugia during these periods likely contributed to the preservation of high levels of intrapopulation genetic diversity (Weir and Schluter 2007). Phylogeographic analyses employing hereditary genetic data elucidate these evolutionary dynamics (Avise 2004). Therefore, Quaternary climatic shifts have likely played a central role in shaping the population structure and lineage diversification within T. bifarius.

When molecular data from both markers are considered together, they provide mutually supportive evidence regarding the origin and distribution of the species in Türkiye. Phylogenetic tree analyses indicate that one of the species’ entry routes into Türkiye is from Transcaucasia in the northeast, specifically through the Ardahan area. It is hypothesized that the Ardahan population remained isolated for an extended period due to glacial conditions, subsequently expanding from the northern side of the North Anatolian Mountains towards Sinop. Another probable entry route is from Hatay in the south, with this population dispersing between the Southern Taurus Mountains and the Central Anatolian Diagonal in the direction of Kayseri. The population that reached Central Anatolia likely merged with the Sinop population, facilitating further expansion into the southern and western regions. The detection of the A2 allele (N = 48) in the western region of the Anatolian Diagonal, a prominent mountain range dividing Anatolia into eastern and western halves, suggests that T. bifarius has been significantly influenced by this geographic barrier (Mutun 2010; Mutun 2016). Furthermore, the absence of shared haplotypes and alleles between populations on the eastern and western sides of the Anatolian Diagonal underscores the mountain range’s significant role in shaping the genetic landscape of T. bifarius. It is hypothesized that the populations in Çanakkale and subsequently Edirne represent stages in the westward expansion of the species from Central Anatolia. Notably, the species’ distribution in Europe is restricted to the Mediterranean region and Central Europe (Chvála, 1988). In light of this information, it can be hypothesized that at least some European populations of T. bifarius may have originated from the Thrace region of Türkiye, while lineages dispersing from the northern Black Sea region through Transcaucasia may have also contributed to these populations. However, comprehensive phylogenetic and phylogeographic analyses of European T. bifarius populations are necessary to confirm this hypothesis. Additionally, analysis of Iranian samples indicates that another potential entry point for the species into Türkiye is from the east, specifically through the Van region. This population likely dispersed through the valleys of Eastern Anatolia, advancing from the northern side of the Central Anatolian Diagonal towards Erzincan. Upon reaching Central Anatolia, it is hypothesized to have merged with existing Central Anatolian populations. Phylogenetic and population genetic analyses collectively indicate that T. bifarius entered Anatolia via three distinct pathways: (1) entry from the Levant region into southern Anatolia via Hatay; (2) dispersal from the Caucasus into northeastern Anatolia through Ardahan; and (3) migration from the Iranian Plateau into eastern Anatolia via Van. These findings offer valuable insights into the historical biogeographic processes that have shaped the region’s biodiversity. The ITS1–ITS2 analysis yielded phylogenetic tree topologies similar to those obtained from the COI data. Notably, the allele structure revealed geographically significant groupings within the polytomic clusters. The co-occurrence of individuals from IRN and VAN in certain subgroups suggests a close genetic relationship between the lineage originating from Iran and populations in Eastern Anatolia. This finding is consistent with the results of the COI phylogenetic analysis. Conversely, populations from northwestern (EDR, CAN) and northeastern (ARD) Türkiye clustered separately, further supporting the existence of geographically structured population groups. Distinct clusters formed by populations from Hatay (HTY), Antalya (ANT), and Kayseri (KAY) further underscore the significant role of regional isolation in shaping the genetic structure. The identification of multiple lineages across distinct regions underscores significant historical dispersal events, followed by subsequent lineage diversification. These patterns have likely been shaped by the complex interplay of topographical barriers and the diverse climatic conditions characteristic of Anatolia and its surrounding regions. Overall, the strong bootstrap support observed across both major clades and their respective substructures indicates high reliability of the inferred phylogenetic relationships. Furthermore, the pronounced phylogeographic partitioning evident in the maximum likelihood tree suggests that geographic barriers and historical demographic processes have played a pivotal role in shaping the genetic structure of T. bifarius populations throughout Türkiye. The findings from both our current study and previous phylogenetic analyses of T. bromius (Sanal-Demirci et al., 2021) indicate that horseflies belonging to the same species can exhibit considerable genetic diversity from a geographical point of view. Additionally, the phylogenetic analysis demonstrates that environmental and climatic differences play a significant role in shaping the genetic structure according to the individuals’ native regions. Finally, our findings provide valuable baseline data to inform future biogeographic and ecological studies of Tabanus species in the Middle East. Given the medical and veterinary significance of Tabanus spp. as vectors of disease, elucidating their population structure is essential to enhancing regional control and surveillance strategies.