Phylogenetic Analysis of Mitochondrial Genome of Tabanidae (Diptera: Tabanidae) Reveals the Present Status of Tabanidae Classification

Simple Summary Tabanidae suck the blood of humans and animals, are important biological vectors for the transmission of diseases, and are of considerable economic and medical significance. However, current knowledge about the mitochondrial genome of this family is limited. Therefore, six newly completed mitochondrial genomes of four genera of Tabanidae (Haematopota turkestanica, Chrysops vanderwulpi, Chrysops dissectus, Tabanus chrysurus, Tabanus pleskei, and Hybomitra sp. species) were sequenced and analyzed. The results show that the six newly mitochondrial genomes have quite similar structures and features. Phylogeny was inferred by analyzing the 13 amino acid sequences coded by mitochondrial genes of 22 mitogenomes (all available complete mitochondrial genomes of tabanidae). Bayesian inference, maximum likelihood trees, and maximum parsimony inference analyses all showed consistent results. This study supports the concept of monophyly of all groups, ratifies the current taxonomic classification, and provides useful genetic markers for studying the molecular ecology, systematics, and population genetics of Tabanidae. Abstract Tabanidae suck the blood of humans and animals, are important biological vectors for the transmission of diseases, and are of considerable economic and medical significance. However, current knowledge about the mitochondrial genome of this family is limited. More complete mitochondrial genomes of Tabanidae are essential for the identification and phylogeny. Therefore, this study sequenced and analyzed six complete mitochondrial (mt) genome sequences of four genera of Tabanidae for the first time. The complete mt genomes of the six new sequences are circular molecules ranging from 15,851 to 16,107 base pairs (bp) in size, with AT content ranging from 75.64 to 77.91%. The six complete mitochondrial genomes all consist of 13 protein-coding genes (PCGs), 2 ribosomal RNA genes (RRNA), 22 transfer RNA genes (tRNAs), and a control region, making a total of 37 functional subunits. ATT/ATG was the most common start codon, and the stop codon was TAA of all PCGS. All tRNA except tRNA Ser1 had a typical clover structure. Phylogeny was inferred by analyzing the 13 concatenated amino acid sequences of the 22 mt genomes. Bayesian inference, maximum-likelihood trees, and maximum-parsimony inference analyses all showed consistent results. This study supports the concept of monophyly of all genus, ratifies the current taxonomic classification, and provides effective genetic markers for molecular classification, systematics, and genetic studies of Tabanidae.


Introduction
Tabanidae is commonly known as the horsefly. Adult Tabanidae, whole metamorphosis insects, are stout, good flyers, and have scraping and licking mouthparts [1,2]. Brachycera are divided into over 20 families, with Tabanidae being one of the largest, The samples were thoroughly washed three times in physiological saline [11,23], then fixed in 75% (v/v) ethanol and stored at −80 °C until further use. TIANamp Genomic DNA Kit (Tiangen) was used to extract insect cox 1 sequence and amplify the cox 1 sequence by PCR for molecular identification, using forward and reverse primers (5′-ATT CAA CCA ATC ATA AAG ATA TTG G-3′) and (5′-TAA ACT TCT GGA TGT CCA AAA AAT CA-3′), respectively. The final PCR reaction mix (25 μL) consisted of 18.3 μL of distilled water, 2.5 μL of 10× Ex Taq buffffer, 2 μL of dNTP mixture, 0.5 μL of each primer, 1 μL of DNA template, and 0.2 μL of Ex Taq DNA polymerase. The PCR reaction mix, in 0.2 Ml PCR tubes, was placed in a Takara TP600 thermocycler (TAKARA, Kusatsu, Japan) and the following thermocycling conditions were used: 94 °C for 5 min (initial denaturation), then 94 °C for 30 s (denaturation), 40 °C for 1 min (annealing), and 72 °C for 1 min/kb (extension) for 40 cycles, and a final extension at 72 °C for 10 min. The PCR products were sent to Sangon Biotech Company (Shanghai, China) for sequencing in both directions using the same primers.

Construction of the Genomic Library and Sequencing
The standard Illumina TruSeq Nano DNA LT library Preparation Guide was used to construct the required on-machine libraries. The main process was as follows (Kit: TruSe-qTM DNA Sample Prep Kit, San Diego, CA, USA): (1) DNA fragmentation: Covaris is used to interrupt DNA and make it fragmented, and the target fragment is purified and sorted by magnetic beads; (2) DNA double-end Repair: the DNA fragment with protruding End was repaired by the joint action of 3'-5' exonuclease and polymerase in End Repair Mix; (3) Introduction of "A" base at the 3'end: introduction of A single base "A" at the 3' end of the repaired and flat DNA fragment; The 3'end of the joint contains A single base "T", so as to ensure that the DNA fragment and the joint can be connected by "A" and "T" complementary pairing, and to prevent the DNA insertion fragments from connecting to each other in the process of connecting the DNA fragment; (4) Connector con-

Construction of the Genomic Library and Sequencing
The standard Illumina TruSeq Nano DNA LT library Preparation Guide was used to construct the required on-machine libraries. The main process was as follows (Kit: TruSeqTM DNA Sample Prep Kit, San Diego, CA, USA): (1) DNA fragmentation: Covaris is used to interrupt DNA and make it fragmented, and the target fragment is purified and sorted by magnetic beads; (2) DNA double-end Repair: the DNA fragment with protruding End was repaired by the joint action of 3'-5' exonuclease and polymerase in End Repair Mix; (3) Introduction of "A" base at the 3'end: introduction of A single base "A" at the 3' end of the repaired and flat DNA fragment; The 3'end of the joint contains A single base "T", so as to ensure that the DNA fragment and the joint can be connected by "A" and "T" complementary pairing, and to prevent the DNA insertion fragments from connecting to each other in the process of connecting the DNA fragment; (4) Connector connection: under the action of ligase, the connector containing the tag is incubated with the DNA fragment to make it connected; (5) Purification of ligation products: magnetic beads purify ligation products to remove free and self-ligation sequence; (6) Verification library: quantitative library using Pico Green; Agilent Bioanalyzer 2100 (State of California) was used for quality control of PCR-enriched fragments to verify the size and distribution of DNA library fragments; (7) Homogenized and mixed library: multiplexed DNA libraries (multiplexed DNA libraries) homogenized to 10 nM after equal volume mixing; and (8) Sequencing: the mixed library (10 nM) was gradually diluted and quantified to 4-5 pM before computer sequencing. The raw data generated were transferred to a computer workstation for analysis and genome characterization. In order to obtain high-quality sequences, FastQC V.0.11.9 [24] software was used to verify the quality index of the data obtained, and Trim Galore V.0.6.5 [25] software was used to remove the adaptation sequences. Finally, Fast QC software was used for further quality inspection and data verification.

Annotation and Bioinformatics Analysis
The complete mitochondrial genome sequences obtained by splicing were uploaded to MITOS Web for function annotation (http://mitos2.Bioinf.unileipzig.de/index.py accessed on 1 July 2022) [30]. TRNAscan-SE was used to identify the secondary structure of tRNA (http://lowelab.ucsc.edu/tRNAscan--SE accessed on 1 July 2022) [31]. DNAStar (V.5.0) was used to calculate the boundaries between genes of A + T, G + C and content, AT-skew = (A − T)/(A + T) and GC-skew = (G − C)/(G + C) [32]. An online open reading frame finder was used to analyze and translate PCGs (https://www.ncbi.nlm.nih. gov/orffinder/ accessed on 1 July 2022). The relative synonymous codon usage (RSCU) of PCGs was determined using the invertebrate mitochondrial genetic code of MEGA X [33]. Comparisons in gene lengths and nucleotides were made between Ha. turkestanica, C. vanderwulpi, C. dissectus, T. chrysurus, T. pleskei, and Hy. sp., all of which belong to the Tabanidae family. Differences in the nucleotide and amino acid sequences were calculated using the MEGAX software. Sliding window analysis was also used to estimate nucleotide diversity (π) between the mitogenomes of Tabanidae species by overlapping at 25 bp intervals for every 200 bp and π was drew at the midpoint position. DnaSP software (V.6) was used to compare the proportion of nonsynonymous (dN) and synonymous (dS) substitutions (dN/dS) in the obtained sequences [34].

Phylogenetics Analysis
Concatenated amino acid sequences of the complete Ha. turkestanica, C. vanderwulpi, C. dissectus, T. chrysurus, T. pleskei, and Hy. sp. mt genome were aligned with the corresponding amino acid sequences of 15 horseflies from a suborder all available in GenBank, with Anopheles sacharovi (MZ382545) as the outgroup. The MAFFT algorithm was used to align all 13 PCGs of each sequence obtained in this study and GenBank database [35].
Phylogenetic relationships between the analyzed species were reconstructed using three methods: BI, MP, andML. MrBayes 3.1 was used to reconstruct the BI tree, and four independent Markov chain runs were performed for 1,000,000 metropolis-coupled (MCMC) generations, sampling a tree every 100 generations. The first 25% (2500) of the generated trees were omitted as burn-in, with the remaining trees being used to calculate Bayesian posterior probabilities. Phylogenetic trees were visualized using the FigTree software. MP methods were performed using the Fitch criterion (1000 bootstrap replicates) within PAUP 4.0 Beta 10. The ML methods were performed using MEGA X software (https: //www.megasoftware.net accessed on 1 July 2022), and bootstrapping was performed with 1000 replicates. Phylograms were drawn using FigTree (v. 1.42) [36].
The cox 1 sequences of the six horseflies in this study were compared with 31 horseflies of five genera of Tabanidae in GenBank, using Rhamphomyia insignis as the outgroup (KT225299.1). A phylogenetic evolutionary tree was constructed using the ML method, as described above.

Acquisition of Mitochondrion cox1 Genes
The specimens were identified as Ha. turkestanica, C. vanderwulpi, C. dissectus, T. chrysurus, T. pleskei, and Hy. sp., which was verified by morphological characteristics. Hy. sp. identification was not species-specific as only the genera was identified. The amplified cox1 sequences were observed to be 93.5%, 93.8%, 96.2%, 95.62%, 96.48%, and 96.15% identical to the sequences of species available in GenBank (MT231188, OM991886, OM991887, NC062705, NC062705, and MT410834), respectively. The cox 1 sequences of the six horseflies in this study were compared with 31 horseflies of five genera of Tabanidae in GenBank, using Rhamphomyia insignis as the outgroup (KT225299.1). A phylogenetic tree was constructed using an ML analytical approach. The results of the phylogenetic trees were consistent with those of the six branches. In this study, Hy. sp. was noted to be in the same branch as Hybomitra astur; Hybomitra lurida, and Hybomitra bimaculata, all belonging to the genus Hybomitra. Furthermore, T. pleskei and T. chrysurus were observed to belong to the genus Tabanus and Ha. turkestanica was observed in the genus Haematopota. The six horseflies of the genus Atylotus were noted to belong to the same branch and both C. dissectus and C. vanderwulpi were located in the same branch as the horseflies of the genus Chrysops ( Figure 2). The results of cox 1 sequence-based evolution analysis were consistent with the results of homology and morphology analyses, indicating that our identification methods are accurate.
flies of five genera of Tabanidae in GenBank, using Rhamphomyia insignis as the outgroup (KT225299.1). A phylogenetic evolutionary tree was constructed using the ML method, as described above.

Acquisition of Mitochondrion cox1 Genes
The specimens were identified as Ha. turkestanica, C. vanderwulpi, C. dissectus, T. chrysurus, T. pleskei, and Hy. sp., which was verified by morphological characteristics. Hy. sp. identification was not species-specific as only the genera was identified. The amplified cox1 sequences were observed to be 93.5%, 93.8%, 96.2%, 95.62%, 96.48%, and 96.15% identical to the sequences of species available in GenBank (MT231188, OM991886, OM991887, NC062705, NC062705, and MT410834), respectively. The cox 1 sequences of the six horseflies in this study were compared with 31 horseflies of five genera of Tabanidae in GenBank, using Rhamphomyia insignis as the outgroup (KT225299.1). A phylogenetic tree was constructed using an ML analytical approach. The results of the phylogenetic trees were consistent with those of the six branches. In this study, Hy. sp. was noted to be in the same branch as Hybomitra astur; Hybomitra lurida, and Hybomitra bimaculata, all belonging to the genus Hybomitra. Furthermore, T. pleskei and T. chrysurus were observed to belong to the genus Tabanus and Ha. turkestanica was observed in the genus Haematopota. The six horseflies of the genus Atylotus were noted to belong to the same branch and both C. dissectus and C. vanderwulpi were located in the same branch as the horseflies of the genus Chrysops ( Figure 2). The results of cox 1 sequence-based evolution analysis were consistent with the results of homology and morphology analyses, indicating that our identification methods are accurate.

mtDNA Features of Six Horseflies in This Study
Complete mitogenomes of Ha. turkestanica, C. vanderwulpi, C. dissectus, T. chrysurus, T. pleskei, and Hy. sp. ranged from 15,851-16,107 bp in length, and were arranged in a classic double-stranded circular DNA molecules. The complete mt genome of Hy. sp. obtained in this study is the first of the genus Hybomitra. Detailed annotations of the mt genomes of the six newly sequenced, including the position and length of each gene and direction of the sequence, are reported in Figure 3.
Complete mitogenomes of Ha. turkestanica, C. vanderwulpi, C. dissectus, T. chrysurus, T. pleskei, and Hy. sp. ranged from 15,851-16,107 bp in length, and were arranged in a classic double-stranded circular DNA molecules. The complete mt genome of Hy. sp. obtained in this study is the first of the genus Hybomitra. Detailed annotations of the mt genomes of the six newly sequenced, including the position and length of each gene and direction of the sequence, are reported in Figure 3. All mitogenomes comprised 13 PCGs (cox 1-3, nad 1-6, nad 4L, atp 6, atp 8, and cyt b), 2 rRNA (rrn L and rrn S), 22 tRNAs and a control region, a total of 37 functional subunits, organized along the N (forward) and J (reverse) strands (Table 1). It has the same structure as other horseflies. Each of the six newly sequenced mitochondrial genomes had one or only one control region, ranging in length from 946 to 1122 bp. The control region contains a large number of AT-rich regions and repeats, which have been shown to be the main reason for the failure of PCR amplification and Sanger sequencing of complete split genomes in invertebrates [37,38]. All mitogenomes comprised 13 PCGs (cox 1-3, nad 1-6, nad 4L, atp 6, atp 8, and cyt b), 2 rRNA (rrn L and rrn S), 22 tRNAs and a control region, a total of 37 functional subunits, organized along the N (forward) and J (reverse) strands (Table 1). It has the same structure as other horseflies. Each of the six newly sequenced mitochondrial genomes had one or only one control region, ranging in length from 946 to 1122 bp. The control region contains a large number of AT-rich regions and repeats, which have been shown to be the main reason for the failure of PCR amplification and Sanger sequencing of complete split genomes in invertebrates [37,38].
The sequences of this study were compared with those of Diptera families, and the results showed that the six newly sequenced horseflies had the same gene order as Tanbanidae, Athericidae, and Rhagionidae species are all from Brachycera, but were different from the Culicidae, Sciaridae, and Trichoceridae species all from Nematocera. Tabanidae mitochondrial genomes are conserved and characterized by the sequence in Diptera. In their mt genome, 23 genes were located on the plus strand, and the remaining 14 genes were located on the minus strand ( Figure S1) [18][19][20][21].
All 13 PCGs (except Hy. sp. and C. vanderwulpi for nad 5 with GTG) used ATN as the start codon, with ATG being the most frequently used (cox 2, atp 6, cox 3, nad 4, nad 4L, and cyt b), followed by ATT (nad 2, cox1, atp 8, and nad 3). All genes used TAA as the standard stop codon (except T. chrysurus and T. pleskei for nad 5, T. chrysurus for cytb, Ha. turkestanica for cox3 and nad 6 using TAG as stop codon) ( Table 1). The canonical start codons most commonly used for invertebrate mitogenomes are ATN, TTG, GTT, and GTG. Generally, horsefly PCGs use ATN as the start codon, whereas the nad 5 gene in
All 13 PCGs (except Hy. sp. and C. vanderwulpi for nad 5 with GTG) used ATN as the start codon, with ATG being the most frequently used (cox 2, atp 6, cox 3, nad 4, nad 4L, and cyt b), followed by ATT (nad 2, cox1, atp 8, and nad 3). All genes used TAA as the standard stop codon (except T. chrysurus and T. pleskei for nad 5, T. chrysurus for cytb, Ha. turkestanica for cox3 and nad 6 using TAG as stop codon) ( Table 1). The canonical start codons most commonly used for invertebrate mitogenomes are ATN, TTG, GTT, and GTG. Generally, horsefly PCGs use ATN as the start codon, whereas the nad 5 gene in some species use GTG as the start codon, which is considered common across various organisms [17,42,43].

Analysis of the RNA (2 rRNAs and 22 tRNAs)
The location of the two ribosomal RNA genes (rRNAs, rrn L and rrn S) are conse in all mitogenomes, with rrn L being flanked by trn L1 and trn V, and rrn S being flan by the trn V-and AT-rich regions. The AT content was similar, despite almost all sp belonging to different genera. Individually, the variation in length of the six rRNAs minimal. The large subunit of mt RNA of the six horseflies varied from 1293 bp in vexativa to 1332 bp in T. pleskei, and the small subunit was 789 bp long in Ha. vexativa

Analysis of the RNA (2 rRNAs and 22 tRNAs)
The location of the two ribosomal RNA genes (rRNAs, rrn L and rrn S) are conserved in all mitogenomes, with rrn L being flanked by trn L1 and trn V, and rrn S being flanked by the trn V-and AT-rich regions. The AT content was similar, despite almost all species belonging to different genera. Individually, the variation in length of the six rRNAs was minimal. The large subunit of mt RNA of the six horseflies varied from 1293 bp in Ha. vexativa to 1332 bp in T. pleskei, and the small subunit was 789 bp long in Ha. vexativa and 796 bp in Tabanus spp. mitogenome. These values were similar to those of previous horsefly mitogenome studies.
The total length of the 22 tRNAs ranged from 2076 to 2444 bp in the six newly sequenced mt genomes, and the individual gene lengths varied from 64 to 72 bp (Table 1). All six horseflies characteristically lacked a DHU arm. The lack of a DHU arm was observed in trnS1 of several other metazoans [44], including tabanid flies [23]. The most common nucleotide mismatch was G-U, followed by U-U, which plays an important role in maintaining the stability of tRNA secondary structure ( Figure S2).

Evolutive Analysis
In genetics, dN/dS represents the ratio between non-synonymous (dN) and synonymous replacement rates (dS). This ratio determines whether there is selection pressure on the protein-coding gene. The corresponding sequences of each PCG in the studied species were paired based on observations of non-synonymous and synonymous replacement rates (dN/dS). The results are shown in two images. The results obtained indicated that the different regions evaluated are evolving globally under the effect of negative pressure, with dN/dS values < 1 and with ratios ranging from cox 1 (0.02 to 0.06) to atp 8 (0.12 to 0.36). Thus, the order of influence of evolutionary pressure was depicted according to the averages obtained: cox 1 < atp 6 < cyt b < cox 3 < cox 2 < nad 3 < nad 1 < nad 4 < nad 4L < nad 2 < nad 5 < nad 6 < atp 8 ( Figure 6). In addition, as reported in other studies, PCGs belonging to mitochondrial complexes III and IV had strong purification pressures and PCGs belonging to the complex I region were weaker. Although the dN/dS rates of atp 8 and nad 6 were < 1, there were signs of weak purification pressure and transient positive evolutionary pressure [16,45,46]. quenced mt genomes, and the individual gene lengths varied from 64 to 72 bp (Table All six horseflies characteristically lacked a DHU arm. The lack of a DHU arm was served in trnS1 of several other metazoans [44], including tabanid flies [23]. The m common nucleotide mismatch was G-U, followed by U-U, which plays an important in maintaining the stability of tRNA secondary structure ( Figure S2).

Evolutive Analysis
In genetics, dN/dS represents the ratio between non-synonymous (dN) and syno mous replacement rates (dS). This ratio determines whether there is selection pressure the protein-coding gene. The corresponding sequences of each PCG in the studied spe were paired based on observations of non-synonymous and synonymous replacem rates (dN/dS). The results are shown in two images. The results obtained indicated the different regions evaluated are evolving globally under the effect of negative press with dN/dS values < 1 and with ratios ranging from cox 1 (0.02 to 0.06) to atp 8 (0.1 0.36). Thus, the order of influence of evolutionary pressure was depicted according to averages obtained: cox 1 < atp 6 < cyt b < cox 3 < cox 2 < nad 3 < nad 1 < nad 4 < nad 4 nad 2 < nad 5 < nad 6 < atp 8 ( Figure 6). In addition, as reported in other studies, PC belonging to mitochondrial complexes III and IV had strong purification pressures PCGs belonging to the complex I region were weaker. Although the dN/dS rates of a and nad 6 were < 1, there were signs of weak purification pressure and transient posi evolutionary pressure [16,45,46]. Three additional analyses were performed to assess π among the mitochondria quences obtained in this study and other horsefly sequences of Tabanidae. Only one tochondrial whole genome sequence was found in Hybomitra, and thus the nucleotide versity analysis of Hybomitra could not be performed. The first analysis evaluated th among mitochondrial sequences of horseflies belonging to the genera of Tabanus, the ond analysis evaluated the π belonging to the Chrysops, and the third analysis evalua the π belonging to the Haematopota. Figure 7 showed the values of π, considering groups of evaluated sequences. The π values ranged from 0.01 to 0.12 in Tabanus (red li Three additional analyses were performed to assess π among the mitochondrial sequences obtained in this study and other horsefly sequences of Tabanidae. Only one mitochondrial whole genome sequence was found in Hybomitra, and thus the nucleotide diversity analysis of Hybomitra could not be performed. The first analysis evaluated the π among mitochondrial sequences of horseflies belonging to the genera of Tabanus, the second analysis evaluated the π belonging to the Chrysops, and the third analysis evaluated the π belonging to the Haematopota. Figure 7 showed the values of π, considering the groups of evaluated sequences. The π values ranged from 0.01 to 0.12 in Tabanus (red line), 0.01 to 0.12 in Chrysops (green line); 0.00 to 0.11 in Haematopota (blue line) (Figure 7). The concatenated sequences revealed that the lowest nucleotide diversity was observed in nad 1 in Tabanus; Chrysops and Haematopota. The nad 2, nad 6, and cytb showed higher nucleotide diversity in all three data sets [41]. 0.01 to 0.12 in Chrysops (green line); 0.00 to 0.11 in Haematopota (blue line) (Figure 7) concatenated sequences revealed that the lowest nucleotide diversity was observed in 1 in Tabanus; Chrysops and Haematopota. The nad 2, nad 6, and cytb showed higher n otide diversity in all three data sets [41].

Phylogenetic Analyses
Phylogenetic analyses of the amino acid concatenated coding regions were ali using the 13 PCGs from 22 taxa (six from this study, the remaining from the GenB database) using three analytical approaches (BI, MP, and ML). All the methods prod nearly identical tree topologies (Figure 8). A monophyletic large group included 21 corresponding to the suborder Brachycera, with Anopheles sacharovi (Diptera: Culic serving as the outgroup. The suborder Brachycera presented a topology with two w supported clades, with the top clade corresponding to the subfamilies Tabanomorpha the lower clades corresponding to Muscomorpha and Stratiomyomorpha. The Tab morpha clade, contained 13 species (including the six gadflies sequenced in this stu and five subclades, representing the genera Chrysops, Haematopota, Tabanus, Hybom and Atylotus, were recovered as a sister group to the Tabanidae tribe, as previously cated. The phylogenetic trees revealed that C. dissectus and C. vanderwulpi belonged t genus Chrysops in the same clades as Chrysops silvifacies; Ha. turkestanica and Ha. vex belonged to the genus Haematopota in the same clade; T. chrysurus, in the same clad T. pleskei, belonged to the genus Tabanus. Hy. sp., and T. haysi were closely related in same clade, followed by A. miser.

Phylogenetic Analyses
Phylogenetic analyses of the amino acid concatenated coding regions were aligned using the 13 PCGs from 22 taxa (six from this study, the remaining from the GenBank database) using three analytical approaches (BI, MP, and ML). All the methods produced nearly identical tree topologies (Figure 8). A monophyletic large group included 21 taxa corresponding to the suborder Brachycera, with Anopheles sacharovi (Diptera: Culicidae) serving as the outgroup. The suborder Brachycera presented a topology with two wellsupported clades, with the top clade corresponding to the subfamilies Tabanomorpha and the lower clades corresponding to Muscomorpha and Stratiomyomorpha. The Tabanomorpha clade, contained 13 species (including the six gadflies sequenced in this study), and five subclades, representing the genera Chrysops, Haematopota, Tabanus, Hybomitra, and Atylotus, were recovered as a sister group to the Tabanidae tribe, as previously indicated. The phylogenetic trees revealed that C. dissectus and C. vanderwulpi belonged to the genus Chrysops in the same clades as Chrysops silvifacies; Ha. turkestanica and Ha. vexativa belonged to the genus Haematopota in the same clade; T. chrysurus, in the same clade as T. pleskei, belonged to the genus Tabanus. Hy. sp., and T. haysi were closely related in the same clade, followed by A. miser.   2021) reported that ML and BI phylogenies showed nonmonophyletic relationships among species of Tabanus and/or Atylotus, where A. miser clustered with T. formosiensis, while Tabanus did not form a single clade. Previous phylogenetic analyses of the nuclear and mitochondrial genomes based on conventional barcoding have also suggested paraphyletism in the genera Tabanus and/or Atylotus [1,5,47,48].
The evolutionary tree constructed in this study showed that Hy. sp. and T. haysi were in the same branch. The evolutionary tree constructed with the cox 1 sequence showed that Hy. sp. was in the same branch as Hybomitra astur; Hybomitra lurida, and Hybomitra bimaculata, and they all belong to the genus Hybomitra. However, T. haysi was in the genus Tabanus, with each genus having its own separate branch, which was consistent with the method of using cox 1 as molecular marker.
Therefore, it is considered that such classification results can be formed because there was only one sequence of Hybomitra in this study, which is the first reported for the genera. Therefore, Hy. sp. was in the same clade as T. haysi, suggesting that phylogenetic analysis using only 21 protein-coding sequences of these species is not representative of the 4455 clades of these species. The evolutionary relationship between the genera Atylotus and Hybomitra within Tabanidae remains unclear, owing to the scarcity of complete mitogenome sequences. Due to its higher mutation rate, mt DNA is generally more prone to saturation than the nuclear genome. Thus, in order for a more accurate and comprehensive analysis of the classification and evolution of horseflies, it is necessary to sequence more mitogenomes of horseflies and include them in future analyses to solve the problem of horsefly classification, as complete mitochondrial genome sequences have been shown to resolve the phylogenetic relationships of many other surface parasites [49][50][51][52].

Conclusions
In this study, we first determined the complete mitochondrial genomes of six horseflies, and secondly, the mitochondrial genomes characteristics and phylogenetic analysis of six newly sequenced were carried out. All mitogenomes evaluated were similar to the   2021) reported that ML and BI phylogenies showed nonmonophyletic relationships among species of Tabanus and/or Atylotus, where A. miser clustered with T. formosiensis, while Tabanus did not form a single clade. Previous phylogenetic analyses of the nuclear and mitochondrial genomes based on conventional barcoding have also suggested paraphyletism in the genera Tabanus and/or Atylotus [1,5,47,48].
The evolutionary tree constructed in this study showed that Hy. sp. and T. haysi were in the same branch. The evolutionary tree constructed with the cox 1 sequence showed that Hy. sp. was in the same branch as Hybomitra astur; Hybomitra lurida, and Hybomitra bimaculata, and they all belong to the genus Hybomitra. However, T. haysi was in the genus Tabanus, with each genus having its own separate branch, which was consistent with the method of using cox 1 as molecular marker.
Therefore, it is considered that such classification results can be formed because there was only one sequence of Hybomitra in this study, which is the first reported for the genera. Therefore, Hy. sp. was in the same clade as T. haysi, suggesting that phylogenetic analysis using only 21 protein-coding sequences of these species is not representative of the 4455 clades of these species. The evolutionary relationship between the genera Atylotus and Hybomitra within Tabanidae remains unclear, owing to the scarcity of complete mitogenome sequences. Due to its higher mutation rate, mt DNA is generally more prone to saturation than the nuclear genome. Thus, in order for a more accurate and comprehensive analysis of the classification and evolution of horseflies, it is necessary to sequence more mitogenomes of horseflies and include them in future analyses to solve the problem of horsefly classification, as complete mitochondrial genome sequences have been shown to resolve the phylogenetic relationships of many other surface parasites [49][50][51][52].

Conclusions
In this study, we first determined the complete mitochondrial genomes of six horseflies, and secondly, the mitochondrial genomes characteristics and phylogenetic analysis of six newly sequenced were carried out. All mitogenomes evaluated were similar to the mtDNA molecular pattern for the Tabanidae family: 37 subunits were subdivided into 13 PCGs, 22 tRNAs, 2 rRNAs, and a control region. All tRNAs had a typical leaf clover structure, except tRNA Ser 1. The three analysis methods, BI, MP, and ML, yielded similar topologies. Our mitogenomes phylogenetic analysis supports the paraphyly of genus Tabanus; Chrysops and Haematopota. Only 21 sequences of these species is not representative of the whole family Tabanidae. Due to a lack of complete genome sequence, Hybomitra and Atylotus genera evolutionary relationship are unclear. Thus, in order for a more accurate and comprehensive analysis of the classification and evolution of horseflies, it is necessary to sequence more mitogenomes of horseflies to solve the problem of horsefly classification to enrich the knowledge of their molecular aspects. These findings provide new genetic markers for further studies of taxonomic, systematics, and genetics of Tabanidae.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/insects13080695/s1. Figure S1: Comparison of linearized mitogenomes of Diptera species. The tRNA genes are shown by a single-letter code of corresponding amino acids. Newly sequenced mitogenomes are indicated in bold. Asterisks indicate the incomplete mitogenomes. Figure