Introduction
The Streblidae is a group of dipteran insects belonging to the superfamily Hippoboscoidea, comprising five subfamilies, 33 genera, and 229 species as of current knowledge (Dick et al., Reference Dick, Graciolli and Guerrero2016). These dipteran flies are important specialized ectoparasites of bats and feed on bat blood. They are commonly referred to as bat flies, along with Nycteribiidae (da Silva et al., Reference da Silva, Ferreira, Hrycyna, Eriksson, Graciolli and Canale2023). Bat flies exhibit a range of unique and significant morphological and physiological adaptations as a result of their long-term parasitic lifestyle. One of the most notable is their adenotrophic viviparity reproductive mode, wherein larvae develop in the female oviduct and are nourished by secretions from accessory glands until maturation into third instar larva (Yang et al., Reference Yang, Huang, Wang, Yang, Zhang and Zheng2023b). Nycteribiidae have wholly vestigial wings that are flattened dorsally and ventrally. They are also inserted dorsally on legs, with the head folded backwards and placed on the thorax, giving them a resemblance to spiders. In contrast, Streblidae have distinctively shaped wings and extensive membranous abdomens (Han et al., Reference Han, Li, Li, Liu, Peng, Wang, Gu, Jiang, Zhou and Li2022). Recent global outbreaks of bat-associated pathogens, such as Ebola viruses, severe acute respiratory syndrome coronavirus (SARS-CoV), SARS-CoV-2, and Middle East respiratory coronavirus (MERS-CoV), have brought bats into the spotlight (Letko et al., Reference Letko, Seifert, Olival, Plowright and Munster2020). A substantial body of research has detected human and animal pathogens in bat flies, including the méjal virus, Amate virus (Ramîrez-Martínez et al., Reference Ramírez-Martínez, Bennett, Dunn, Yuill and Goldberg2021), Hemoplasmas (Wang et al., Reference Wang, Li, Peng, Gu, Zhou, Xiao, Han and Yu2023), Bartonella (Yang et al., Reference Yang, Wang, Yang, Zhang, Zheng and Huang2024) and others. Additionally, Kuang et al. (Reference Kuang, Xu, Wang, Gao, Yang, Wu, Liang, Shi and Feng2023) successfully isolated Nelson Bay reovirus (NBV) from bat flies collected in Yunnan Province, with the S1, S2, and M1 segments showing high similarity to human pathogens. Streblidae have been known to occasionally bite humans (Szentiványi et al., Reference Szentiványi, Heintz, Markotter, Wassef, Christe and Glaizot2023), and their flight ability increases the possibility of pathogen spillover, which has generated widespread interest among researchers. However, conducting research in this area remains challenging due to the small size of Streblidae, the absence of distinct morphological characteristics for identification, as well as morphological identification reference material.
Mitochondrial genomes are circular double-stranded DNA molecules that typically consist of 13 protein-coding genes (PCGs), 22 transfer RNAs (tRNAs), two ribosomal RNAs (rRNAs), and a control region (CR) (Clary and Wolstenholme, Reference Clary and Wolstenholme1985). They are the only extranuclear genetic information carrier in animals and are characterized by maternal inheritance, conserved gene content, small molecular weight (approximately 14–21 kb), high mutation rates, and rapid evolution. Therefore, mitochondrial DNA (mtDNA) is often employed in comparative evolutionary genomics and phylogenetic analysis research (Crozier and Crozier, Reference Crozier and Crozier1993). However, the GenBank database only currently contains three complete mitochondrial genomes of Streblidae. This is significantly disproportionate to the diversity of Streblidae.
In this study, we report on the first complete mitochondrial genome of Brachytarsina amboinensis. We provide a detailed analysis of the features of the new mitochondrial genome, including nucleotide composition and each tRNA's secondary structure. We also compare the differences in nucleotide composition, codon usage, relative synonymous codon usage (RSCU), and AT-skew and GC-skew among different Streblidae mitochondrial genomes. Finally, we investigate the phylogenetic relationships between 17 species in the Hippoboscoidea, 2 species in the Muscoidea, and 2 species in the Oestroidea, based on PCGs. The findings not only enrich the molecular database for bat flies but also provide a reference for the identification and phylogenetic analysis of Streblidae.
Materials and methods
Sample collection and processing
In July 2022, a total of 22 Streblidae specimens were collected from the surface of the eastern bent-wing bat (Miniopterus fuliginosus) in Binchuan County (100.58′E, 25.83′N), Dali Bai Autonomous Prefecture, Yunnan Province, China, and were subsequently preserved in 95% anhydrous ethanol in the laboratory. The samples were identified as B. amboinensis using a SZ2-ILST dissection microscope (Olympus, Tokyo, Japan) and following the identification characteristics described in Wenzel (Reference Wenzel1976). One of the samples was sent to Novogene Co., Ltd. (Beijing, China) for DNA extraction and sequencing, while the remaining specimens were stored in a −20°C freezer at the Institute of Pathogens and Vectors, Dali University. Upon completion of DNA extraction, a portion of the DNA is reserved for the purpose of validating CR. An Illumina paired-end (PE) library was prepared and DNA sequencing was performed on the Illumina Novoseq 6000 sequencing platform, generating a total data volume of 3 G (150 bp reads).
Sequence assembly, annotations and analysis
The raw sequence data of B. amboinensis were assembled using NOVOPLASTY 4.2.1 (Dierckxsens et al., Reference Dierckxsens, Mardulyn and Smits2017), and the CR was verified by Sanger sequencing. The specific primers were designed to amplify the CR sequence based on the 16s rRNA and nad2 genes in the genome. The forward primer sequence is AACAGGCGAACATTATATTTGCCG, and the reverse primer sequence is AATCAGTAATGAAACGGAGCAG. The PCR amplification system outlined in the instructions of TaKaRa LA Taq (TaKaRa, Dalian, China) was employed, with the PCR reaction conditions involving an initial denaturation at 94°C for 4 min, followed by 37 cycles of denaturation at 94°C for 30 s, annealing at 50°C for 30 s, extension at 72°C for 1 min, and a final extension at 72°C for 10 min. The resulting amplified products were analyzed by electrophoresis on a 1% agarose gel, and then sent to Sangon Biotech (Shanghai) Co., Ltd. for bidirectional sequencing. The Sanger trace files were viewed using Geneious Primer 2.2.0, followed by removal of the poorly sequenced parts and manual editing of the sequence. The edited sequence was then aligned with the assembled complete mitochondrial genome. The assembled CR was found to be a perfect match with the CR obtained by PCR. Annotation of the PCGs, tRNAs, and rRNAs of the assembled circular sequence was performed using the MITOS Web Server (Donath et al., Reference Donath, Jühling, Al-Arab, Bernhart, Reinhardt, Stadler, Middendorf and Bernt2019), and the annotated tRNAs were validated using the tRNA-scan Web Server (Chan and Lowe, Reference Chan and Lowe2019). The 13 PCGs were compared with known sequences in the database and manually corrected. The starting positions of the genes cox1, nad1, nad5, and nad6 were modified, while the ending positions of cox2, nad4, and nad5 were altered. The secondary structure of tRNAs was predicted using MITOS, and Adobe Illustrator CC 2019 was used for image processing. The circular mitochondrial genome was visualized using the online tool GENOMEVX (Conant and Wolfe, Reference Conant and Wolfe2008). The mitochondrial genome data has been deposited in GenBank with the accession number OQ247894.
Comparative analysis
MEGA 11 (Tamura et al., Reference Tamura, Stecher and Kumar2021) was utilized for calculating nucleotide composition statistics, pairwise genetic distances, and Relative Synonymous Codon Usage (RSCU) for four Streblidae sequences (GenBank accession numbers: MK896865, MK896866, OQ301747, and OQ247894). RSCU was calculated by excluding incomplete codons of PCGs. Moreover, AT-skew and GC-skew were computed using the following formulae: AT-skew = (A − T)/(A + T) and GC-skew = (G − C)/(G + C). The Ka/Ks ratio of non-synonymous (Ka) and synonymous (Ks) substitutions between B. amboinensis and the other three Streblidae (Paradyschiria parvula, Paratrichobius longicrus, and Raymondia sp.) was analyzed using the mlwl-based Ka/Ks calculator 2.0 (Wang et al., Reference Wang, Zhang, Zhang, Zhu and Yu2010). Additionally, the nucleotide diversity (Pi) of PCGs in the four Streblidae was examined using DnaSP 6.0 (Rozas et al., Reference Rozas, Ferrer-Mata, Sánchez-DelBarrio, Guirao-Rico, Librado, Ramos-Onsins and Sánchez-Gracia2017) with a sliding window of 100 bp and a step size of 20 bp.
Phylogenetic analyses
We selected a total of 17 species, representing all four families of the superfamily Hippoboscoidea, two species from the superfamily Muscoidea, and two species from the superfamily Oestroidea, with Chironomus tepperi and Dixella aestivalis as outgroups (table 1). The 13 PCGs of the mitochondrial genome were extracted for building the phylogenetic trees using the PhyloSuite platform (Zhang et al., Reference Zhang, Gao, Jakovlić, Zou, Zhang, Li and Wang2020). The PCGs were aligned using MAFFT 7 (Rozewicki et al., Reference Rozewicki, Li, Amada, Standley and Katoh2019), and the fuzzy alignments were removed using Gblocks 0.91b (Castresana, Reference Castresana2000). The aligned sequences were manually concatenated using PhyloSuite. Maximum Likelihood (ML) and Bayesian inference (BI) analyses were performed using IQ-TREE 2.2.0 (Minh et al., Reference Minh, Schmidt, Chernomor, Schrempf, Woodhams, Von Haeseler and Lanfear2020) and MrBayes 3.2.7a (Ronquist et al., Reference Ronquist, Teslenko, Van Der Mark, Ayres, Darling, Höhna, Larget, Liu, Suchard and Huelsenbeck2012), respectively. The best partitioning schemes and corresponding nucleotide substitution models were inferred using ModelFinder (Kalyaanamoorthy et al., Reference Kalyaanamoorthy, Minh, Wong, Von Haeseler and Jermiin2017) with the corrected Akaike information criterion (AICC) and greedy search algorithm. The best-fit partitioning scheme and corresponding models used on the BI phylogenetic tree are shown in table S1. The ML analysis was evaluated using the ultrafast bootstrap approximation approach for 1000 replicates. BI analyses were performed with two Markov chain Monte Carlo (MCMC) chains runs of 10,000,000 generations, with sampling every 1000 generations. Convergence was considered to be reached when the average standard deviation of the split frequencies was lower than 0.01. The first 25% of generations from each MCMC chain run to account for potential lack of convergence were discarded as burn-in, and the remaining samples were used to generate the majority consensus trees and to estimate the posterior probabilities. The generated phylogenetic trees were viewed and edited using FigTree 1.44 (http://tree.bio.ed.ac.uk/software/figtree/).
Mitochondrial genome information used in this study
Results
Genomic organization and base compositions
The mitochondrial genome of B. amboinensis is circular and double-stranded, with a length of 15,693 bp. It contains 13 PCGs, 22 tRNAs, two rRNAs, and one CR (fig. 1 and table 2), exhibiting a relatively high degree of conservation, with an identical gene count, arrangement, and orientation as the ancestral insect's genome (Boore, Reference Boore1999). The nucleotide composition of B. amboinensis mitochondrial DNA is 42.46% A, 38.67% T, 7.25% G, and 11.62% C, with an overall AT content of 81.13%.
The circular mitochondrial genome map of Brachytarsina amboinensis.
Mitochondrial genome organization of Brachytarsina amboinensis
The protein coding genes
The 13 PCGs of B. amboinensis range in length from 156 to 1696 bp, with nad5 being the longest and atp8 being the shortest. The total length of the 13 PCGs is 11,083 bp, and the AT content is 80.46%. All PCGs start with ATN (5 ATG, 6 ATT) start codons, except for cox1 (starting with TCG) and nad1 (starting with TTG), which use non-standard start codons. All PCGs, except cox1, nad4, and nad5, end with TAN stop codons (7 TAA, 3 TAG), while cox1, nad4, and nad5 have single T stop codons. The use of TCG as a start codon for cox1 is common in Diptera (Ma et al., Reference Ma, Wang, Wu, Li, Sun, Wang and Chang2022), as non-standard start codons for cox1 are frequently observed in holometabolous insects (Li et al., Reference Li, Yan, Pape, Gao and Zhang2020). Incomplete stop codons are fairly common in the Diptera mitogenomes and can be converted into a potential stop codon via polyadenylation to TAA (Nelson et al., Reference Nelson, Lambkin, Batterham, Wallman, Dowton, Whiting, Yeates and Cameron2012).
The tRNAs, rRNAs and control region
The complete mitochondrial genome of B. amboinensis contains 22 tRNAs, ranging in length from 62–72 bp (table 2). Nineteen of them can be folded into typical cloverleaf structures (fig. 2), while trnS1 lacks the dihydrouridine arm (DHU), and trnG and trnT lack the thymidine–pseudouridine–cytidine (TΨC) loop. There are 17 non-matching base pairs in the tRNAs, mostly consisting of GU mismatches (15 in total). Additionally, there are 2 UU mismatches in the acceptor stem of trnA and trnL2 (fig. 2).
Secondary structure of the tRNA genes from the mitochondrial genome of Brachytarsina amboinensis.
The combined length of the two rRNAs is 2150 bp, with an overall A + T content of 82.22%. The 16S rRNA is 1367 bp in length and located between trnL and trnV, while the 12S rRNA is 783 bp in length and located between trnV and the CR.
The length of the CR is 968 bp, and is located between the 12S rRNA and the trnI. It possesses a very high AT content of 85.43%. This region is significant in the mitochondrial genome, as it plays a critical role in DNA replication and the initiation of transcription (Zhang and Hewitt, Reference Zhang and Hewitt1997).
Comparative analysis
It is commonly recognized that the analysis of base composition and strand asymmetry, including AT- and GC-skew, provides insights into base composition bias in mitochondria (fig. 3). Based on the analysis of the four Streblidae mitochondrial genomes, an AT content ranging from 78.63 to 81.13% was observed, indicative of a preference for A and T nucleotides, which is typical of insect mitochondrial genomes (Clary and Wolstenholme, Reference Clary and Wolstenholme1985). Furthermore, the AT content in CR of the four Streblidae is notably higher, ranging from 83.09 to 92.3%, compared to the entire genome. The observed strand asymmetry pattern in the Streblidae mitochondrial genomes aligns with the characteristic strand bias of dipteran insects, with an AT-skew range of 0.02 to 0.08 and a GC-skew range of −0.28 to −0.23 (Ren et al., Reference Ren, Shang, Yang, Shen, Chen, Wang, Cai and Guo2019). These skewed strand compositions are likely attributed to mutational and selective pressures (Kono et al., Reference Kono, Tomita and Arakawa2018), with the GC-skew value associated with the direction of replication within the genome (Chen et al., Reference Chen, Lu, Bork, Hu and Lercher2016).
AT content (A), AT-skew (B) and GC-skew (C) of the 13 PCGs of four Streblidae.
For evaluating synonymous codon bias, the Relative Synonymous Codon Usage (RSCU) serves as a key parameter, and the obtained values directly reflect the codon usage preference (fig. 4). In B. amboinensis' complete mitochondrial genome, there are a total of 3684 codons encoding proteins, excluding the stop codons. The most commonly used codons are UUA, AUU, UUU, AUA, and AAU, whereas the least commonly used codons are CUC, CGC, AGG, and UGG. Among all 62 codons, 26 of them are considered to have a usage preference, as these synonymous codons have a positive codon bias (RSCU value > 1). The codon usage preference among Streblidae is similar, with the top five most commonly used codons across all four species being UUA, AUU, UUU, AUA, and AAU (fig. 4). Additionally, the RSCU values of the third position nucleotides being A or T are significantly higher than those of G or C, indicating that the high A + T content in the mitochondrial genome leads to codon usage bias.
Heat map of RSCU value of each codon among four species of Streblidae.
The evolutionary patterns for each PCG were assessed using pairwise Ka/Ks and genetic distances (fig. 5). The highest Ka/Ks values were observed for nad5 (0.59) and nad6 (0.59), suggesting that they may have undergone faster evolutionary rates compared to other PCGs. This implies that nad5 and nad6 may have experienced more relaxed selection constraints and accumulated a greater number of mutations over time. In contrast, cox1 (0.12) displayed the smallest Ka/Ks value, indicating a slower evolutionary rate. This suggests that cox1 faces greater evolutionary pressure and evolves at a slower rate compared to the other genes. The Ka/Ks values of all 13 PCGs were less than 1, indicating that purifying selection may have dominated the evolution of the mitochondrial genome (Hurst, Reference Hurst2002; Chang et al., Reference Chang, Qiu, Yuan, Wang, Li, Sun, Guo, Lu, Feng and Majid2020). Genetic distances showed similar results: cox1 (0.21) was evolving comparatively slowly, while nad2 (0.43), nad3 (0.42), and nad6 (0.43) were evolving comparatively quickly. The Pi values of the 13 PCGs in the four Streblidae ranged from 0.178 to 0.328 (fig. 6). Among these genes, nad2 showed the highest variability (0.328), followed by nad6 (0.313), while cox1 (0.178) showed the lowest variability.
Substitution ratios in the individual PCG of mitochondrial genomes of four Streblidae (the mean value of Ka/Ks and genetic distance under figure).
The nucleotide diversity (Pi) of 13 PCGs of five Streblidae.
Phylogenetic analysis
Using BI and ML methods, the phylogenetic relationships of the Hippoboscoidea were analyzed based on 13 PCGs (fig. 7). The constructed BI and ML trees both supported the monophyly of Hippoboscoidea, Glossinidae, Hippoboscidae, and Nycteribiidae, and had high node support values. Streblidae were found to be paraphyletic, with differing relationships to the Nycteribiidae in the two trees. In the BI tree, New World Streblidae were grouped with Nycteribiidae with a high posterior probability (BPP = 1), providing strong support for this relationship. In the ML tree, Old World Streblidae were grouped with Nycteribiidae with moderate bootstrap support (BS = 70), indicating that this relationship is less certain. The Hippoboscidae was found to be the sister group to the bat flies (Streblidae and Nycteribiidae), while the Glossinidae was sister to a clade comprising all three of these families.
Phylogenetic analysis based on nucleotide sequences of 13 PCGs from the mitogenome. Numbers at nodes indicate bootstrap (ML) and posterior probability (BI) support values, respectively. ‘■’ indicates the specie in this study.
Discussion
Designing universal markers is indeed crucial for resolving species identification issues. Universal markers, such as specific DNA sequences or genetic loci that are conserved across related species, are valuable tools for distinguishing between different organisms. The cox1 gene has historically been widely used as a universal DNA marker for insect species identification, but its high conservation within the insect mitogenome can limit its efficiency (Ožana et al., Reference Ožana, Dolný and Pánek2022; Yang et al., Reference Yang, Chen and Dong2023a). Due to the conservativeness of the cox1 gene, it may not be possible to reach the critical threshold for species identification of closely related species. Before determining the species identity, it may only be possible to assign it to a higher taxonomic group, such as a phylum or order. However, the nad2 and nad6 genes display significant nucleotide variations and a high Ka/Ks ratio, which can be utilized for distinguishing species with close phylogenetic relationships, thus establishing them as potential markers for taxonomic categorization within the family Streblidae. Furthermore, combining cox1 with other genes has been demonstrated to enhance the accuracy and reliability of species identification in some studies (Alberfkani et al., Reference Alberfkani, Albarwary, Jaafar, Zubair and Abdullah2022). However, it is important to emphasize that further experimental validation and cross-validation are necessary to ensure the effectiveness and stability of newly designed markers. This validation process is critical for establishing the reliability and utility of these markers for species identification and delimitation.
The formation of the concept of the superfamily Hippoboscoidea was a lengthy process. Initially, it only included three core families, which are the ‘Pupipara’ (Nycteribiidae, Streblidae, and Hippoboscidae). The current understanding of Hippoboscoidea includes the Glossinidae and the Pupipara, and this was proposed by Hennig (Reference Hennig1971) in 1971, and authors such as McAlpine (Reference McAlpine1989) and Ziegler (Reference Ziegler2003) have since defended this view. But Nycteribiidae and Streblidae have even been classified under the family Hippoboscidae at times (McAlpine, Reference McAlpine1989). At the same time, different opinions and conclusions existed regarding the evolutionary relationships within the Hippoboscoidea. Hennig (Reference Hennig1971) suggested that Glossinidae forms the sister group of the other three families within Hippoboscoidea, and that the families Nycteribiidae and Streblidae constitute a monophyletic group. Conversely, McAlpine (Reference McAlpine1989) supported the concept that Glossinidae and Hippoboscidae together form the sister group of Nycteribiidae and Streblidae. Molecular data and ML analyses by Nirmala et al. (Reference Nirmala, Hypša and Žurovec2001) and Dittmar et al. (Reference Dittmar, Porter, Murray and Whiting2006) supported the latter view, but the Streblidae were found to be a paraphyletic group. Subsequent research by Petersen et al. (Reference Petersen, Meier, Kutty and Wiegmann2007) conducted a phylogenetic analysis using combined sequence data from CAD, cox1, 16S, and 28S, and supported the view that Glossinidae, Hippoboscidae, and Nycteribiidae were monophyletic. However, upon analyzing Old and New World Streblidae together, the Streblidae displayed paraphyly, with those from the Old World clustering with the Nycteribiidae. Poon et al. (Reference Poon, Chen, Tsang, Shek, Tsui, Zhao, Guénard and Sin2023) constructed a phylogenetic tree based on the cox1 gene and demonstrated that the Streblidae from the Old World were closely related to the Nycteribiidae and formed a clade which was sister to the New World Streblidae. Kutty et al. (Reference Kutty, Pape, Wiegmann and Meier2010) conducted a phylogenetic analysis based on PCGs, including cox1, cytb, Ef1α, and CAD, and showed that Glossinidae, along with three other families within the superfamily Hippoboscoidea, formed a sister group. Nycteribiidae and Streblidae formed a clade, with Nycteribiidae being a monophyletic group, and Streblidae being paraphyletic. Moreover, the previous studies on the Hippoboscoidea predominantly used partial PCGs, which reflects the relatively limited phylogenetic relationships (Cameron et al., Reference Cameron, Lambkin, Barker and Whiting2007).
The present study employed all 13 PCGs of the mitogenomes and encompassing all major lineages to investigate the systematic relationships within the Hippoboscoidea. Both the ML and BI trees indicated that Glossinidae, Hippoboscidae, and bat flies (Nycteribiidae and Streblidae) were monophyletic groups. The Glossinidae were observed to form a sister group with the other three families, while Hippoboscidae formed a sister group with the bat flies. Interestingly, the New World Streblidae and Nycteribiidae formed a clade, representing a divergence from previous studies and providing new insights into the evolutionary relationships of the Hippoboscoidea. However, the limited number of complete mitochondrial genomes for this group still poses a challenge in understanding their systematic evolution, highlighting the need for additional genomic data to facilitate a comprehensive understanding of their relationships.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0007485324000762.
Author contributions
J. T. Y., X. B. H., Y. J. W. and H. J. Y. designed the study, J. T. Y., X. B. H., Y. J. W. and H. J. Y. collated the data, J. T. Y. analysed the data and wrote the first draft of manuscript, X. B. H. completed the final draft of the manuscript and all authors contributed to the final version.
Financial support
This work was supported by the National Natural Science Foundation of China [No. 32001096 to Xiaobin Huang] and the Science and Technology Planning Project of the Yunnan Provincial Department of Science and Technology (No. 202001AT070025 to Xiaobin Huang).
Competing interests
None.
Ethical standards
No specific permits were required for the insects collected for this study. The capture of bats according to the standards and procedures set by the Animal Ethics Committee of Dali University. (Name: Dali University Ethics Committee; approval ID: MECDU-202104-27).
