Tabanids were captured near horse and cattle farms in Daqing City, Heilongjiang Province, Northeastern China, using the net method from July to September 2021. Collectors held the end of the net handle and waved it in an “X” pattern, allowing for easy maneuvering while walking. Tabanids of the same genus were placed in specialized cages to shield them from direct sunlight and ensure their activity. Species identification was based on morphological characteristics, including the shape of the basal scapula of the head, as well as the length, width, and degree of dorsal protrusion of the antennae. We also examined the colour of the back plate on the abdomen and wing pattern (China fauna editorial committee of the Chinese Academy of Sciences, 1998). The cleaning procedures were as follows: (1) The tabanids were placed into a sterile centrifuge tube, followed by the addition of 75% ethanol solution. The tube was then gently shaken for 1 min; (2) the tabanid samples were transferred to a new centrifuge tube, and 0.1% sodium hypochlorite solution was added. The tube was again shaken gently for 1 min; (3) the samples were transferred once more to a new centrifuge tube and cleaned twice with distilled water; (4) distilled water was added again, and after gentle agitation, the liquid was poured out for preservation; and (5) and the cleaned tabanids were fixed onto wax plates and dissected under sterile conditions, and tissue samples were collected from the Malpighian tubules, midgut, and ovary. The tissue samples were stored at −80°C for subsequent use. A total of 432 tabanids were collected. Among them, 135 tabanids were used in the study, including 40 T. griseinus, 45 Atylotus sp., 25.

H. turkestanica and 25 Hybomitra sp. After dissecting the tabanids five samples was combined into one pool and submitted to sequencing. All tabanids were adults, and the detailed sample information is shown in Table 1.

DNA was extracted from the Malpighian tubules, midgut, and ovary tissues using the DNA Stool Mini Kit, following the manufacturer’s protocol. Subsequently, the extracted DNA was utilized for Illumina sequencing using primers 338F: 5′- ACT CCT ACG GGA GGC AGC AG-3′ and 806R: 5′- GGA CTA CCA GGG TATC-TAA TCC-3′, targeting the 16S (V3–V4) region. PCR products were subjected to 250 bp paired-end sequencing on an Illumina NovaSeq 6,000 platform.

Raw reads were initially filtered using Trimmomatic (v0.33) (Bolger et al., 2014). Subsequently, the primer sequences were identified and eliminated using Cutadapt (v1.9.1) (Martin, 2011), resulting in the generation of high-quality reads without primer sequences. The clean reads were then spliced and clustered on the QIIME2 platform (Bolyen et al., 2019) using the DADA2 algorithm (Callahan et al., 2016) to obtain amplicon sequence variants (ASVs). Taxonomic annotation of the feature sequences was processed using a Bayes classifier with SILVA (release 138.1) as the reference database (Quast et al., 2012). The Shannon diversity index and Bray-Curtis distance of the samples were calculated using the “diversity” and “vegdist” functions in the R package: vegan. Permutational multivariate analysis of variance (PERMANOVA) was performed using the adonis function in vegan. The ggplot2 package was used to visualize data such as species composition and principal component analysis. Statistics on the composition of each sample were calculated at various taxonomic levels, including phylum, class, order, family, genus, and species. The abundance of each species in the samples and distribution histograms at each taxonomic level were obtained using QIIME2.

All analyses were conducted using R version 4.2.1. Rarefaction curves based on the number of feature sequences were generated using “estimate R” function in the package vegan v2.5–7 (Dixon, 2003). The profile was transformed into a table of relative abundances prior to statistical analyses. Wilcoxon rank sum test was employed to assess differences in species between groups, and a p-value of <0.05 was considered statistically significant. The Shannon index was used to evaluate microbiota diversity. Principal coordinate analysis (PCoA) was conducted based on the Bray-Curtis distance and permutational multivariate analysis of variance (PERMANOVA) with a permutation of 999. Linear discriminant analysis of effect size (LEfSe) was performed to identify the key species that most likely explained the differences between the groups. Multiple bacterial pathogen detection (MBPD), which can detect a broad range of animal, plant, and zoonotic pathogens based on 16S rRNA gene sequencing (Yang et al., 2023), was employed to explore putative pathogens in horsefly microbial communities.

Phenotypic sex determination relies on secondary sex characteristics, particularly eye shape. Morphological identification was conducted individually for each collected specimen, using morphological descriptions to ascertain the genus or species level. The captured tabanids were classified into four genera: Atylotus: Usually small sized fly, frons with spotted calli or without it. Colour of eyes in living specimen green or yellow. Basal plate of flagellum broad and obtuse dorsal angle. Basicosta pale to brown setulose. Haematopota: Generally small and slender flies of brownish to blackish grey in colour; eyes with several wavy bands in live condition; frons with velvety black spot on each side above the frontal callus and often a mid-frontal spot above these; picture wing pattern, i.e., dark wing with pattern of pale spots; mid tibiae and hind tibiae often with pale rings. Tabanus: Frons with prominent callus. Hybomitra: Medium sized fly; vertex with ocellar tubercle; eyes with dense pubescence and 3 green or purple band in live condition; basal and median callus usually broad; body blackish to dark greyish often with orangish side markings in at least 2 to 3 anterior abdominal segments. Two of these were identified at the genus level (Atylotus sp. and Hybomitra sp.), while the other two were identified at the species level (Haematopota turkestanica and Tabanus griseinus). The morphological characteristics of the cells are depicted in Figure 1A.

After dissecting the tabanids We characterized different organs bacterial community composition by sequencing the V3–V4 region of the 16S rRNA (16S) gene. In total, 7,516,786 pairs of raw reads were obtained. After quality control, 7,486,165 clean reads were obtained, with an average of 77,981 clean reads per sample. A total of 1,338 ASVs were obtained following the splicing and clustering of clean reads. Notably, more than 1,000 ASVs were shared among the different organs of the four tabanid species, indicating the presence of a relatively stable core community within each organ. Taxonomic classification revealed that nearly all ASVs could be classified above the genus level, indicating that the majority of the tabanid microbiota comprised known bacteria (Figures 1B–D).

To investigate the diversity of the microbiota in Tabanidae, the microbiota present in the midgut, Malpighian tubules, and ovary of four genera of Tabanidae were analyzed. The rarefaction analysis indicated that the sampling number curve for each organ tended to saturate at 20 (Figure 2A). Moreover, the number of species observed in the ovaries exceeded that in the other two organs with the same sample size, indicating higher species richness in the ovaries.

Subsequently, the Shannon index was used to evaluate the microbial diversity in each organ. The midgut exhibited significantly lower microbial diversity than the other organs, whereas the ovaries displayed the highest microbial diversity (p < 0.05). The rank-sum test revealed that Tabanus griseinus, Hybomitra sp., and Haematopota turkestanica exhibited significantly higher microbiota diversity in their ovaries than in the other organs (p < 0.05; Figure 2B). Furthermore, microbiota diversity within the same organ among different tabanids was compared. The microbiota diversity of Atylotus sp. ovaries was low (p < 0.05, compared with other groups), whereas the differences among the other groups were not statistically significant (Figure 2C).

PCoA was conducted using the Bray-Curtis distance metric revealed significant effects of different organs and tabanids on the microbiota composition (Figure 2D). Notably, the organs exhibited stronger explanatory power for microbiota variation (R2 = 0.2264, p < 0.001) than the tabanids. Further analysis of the characteristics of microbial diversity among the different tabanids involved the calculation of the distance between samples at the species and organ levels. We observed significant differences in the β-diversity of the microbiota across different organs within the same species (Figure 2E). Among different species, the β-diversity of the Malpighian tubules microbiota is relatively similar, whereas the midgut and ovary microbiota exhibit substantial differences (Figure 2F).

A total of 1,338 ASVs were classified into 20 phyla, 225 families, and 506 genera. Firmicutes (37.6% ± 13.7%), Proteobacteria (32.7% ± 20.4%) and Bacteroidetes (11.6% ± 6.3%) were the dominant phyla. Other phyla, including Tenericutes (7.8% ± 17.4%), Epsilonbacteraeota (2.3% ± 3.5%), Cyanobacteria (2.3 ± 2.0%), Actinobacteria (1.9% ± 1.0%), and Acidobacteria (1.4% ± 1.0%) were also observed (Figure 3A). At the family level, Enterobacteriaceae (22.8% ± 20.0%), Lachnospiraceae (9.5% ± 8.8%), Ruminococcaceae (8.7% ± 4.2%), and Spiroplasmataceae (7.4% ± 17.5%) were predominant (Figure 3B).

The characteristics of the microbial composition in different organs of the four Tabanidae species were examined to analyze the variation in microbiota in tabanids. Firmicutes, Proteobacteria, and Bacteroidetes were the most dominant bacterial phyla in the various organs (Figure 3C). Notably, a high abundance of Tenericutes (20.9% ± 24.9%) was observed in the midguts of most tabanids. However, Haematopota turkestanica exhibited a lower abundance of Tenericutes (1.2% ± 0.7%) in its midgut than other tabanids. In Malpighian tubules, Tabanus griseinus and Haematopota turkestanica displayed higher abundances of Epsilonbacteraeota(5.0 ± 5.6 and 5.6% ± 7.5%). Furthermore, Tabanus griseinus exhibited a higher abundance of Tenericutes (11.3% ± 14.8%) and a lower abundance of Proteobacteria (16.8% ± 10.1%) in the Malpighian tubule. However, the overall structure of ovarian microbiota across the four tabanid species remained relatively consistent.

At the family level, Enterobacteriaceae (19.8% ± 22.9%), Lachnospiraceae (10.9% ± 11.1%), and Ruminococcaceae (9.7% ± 4.3%) emerged as the dominant families across the organs (Figure 3D). Notably, the abundance of Spiroplasmataceae (11.1% ± 14.8%) in the Malpighian tubules of Tabanus griseinus surpassed that in other tabanids. Additionally, higher levels of Enterobacteriaceae (23.7% ± 14.3%) and lower levels of Lachnospiraceae (7.1% ± 1.9) were observed in the Malpighian tubules of Atylotus sp. Spiroplasmataceae (29.1 ± 21.1 and 41.3% ± 38.3%) was notably abundant in the midgut of Tabanus griseinus and Hybomitra sp. Interestingly, the midgut of Haematopota turkestanica contains more Lachnospiraceae (6.7% ± 4.1%) but almost no Spiroplasmataceae (0.6% ± 0.6). It is noteworthy that compared to other tabanids, Atylotus sp. exhibited a higher enrichment of Enterobacteriaceae (33.2 ± 26.8 and 39.0% ± 15.9%) in both the midgut and ovary. The structures of the ovary microbiota of Tabanus griseinus, Haematopota turkestanica, and Hybomitra sp. were consistent, suggesting that the composition of the ovary microbiota may not be significantly influenced by genetic background.

The main changes observed at the family level extended downstream to the hierarchical groups; therefore, we further examined the characteristics at the genus level. The relative abundances of dominant bacterial genera were not evenly distributed among each group (Supplementary Table S1), therefore, LEfSe analysis was used to detect characteristic genera in each group. This study confirmed consistent enrichment of Lachnospiraceae NK4A136 group and Helicobacter (LDA > 3, p < 0.05) in the Malpighian tubules of Tabanus griseinus, Atylotus sp., and Hybomitra sp. (Figures 4A,B,D). Notably, there was minimal difference in the microbiota of each organ in Haematopota turkestanica, with only three enriched genera (Figure 4C) found in the Malpighian tubules after lowering the threshold of the LDA score (LDA > 2, p < 0.05). A few bacteria were enriched in the midgut; however, Photobacterium was uniformly enriched in Atylotus spp., Haematopota turkestanica, and Hybomitra spp. In addition, only Photobacterium was significantly enriched in the midguts of Hybomitra spp. (Figure 4D).

The ovaries were enriched with the most characteristic bacteria, and similar to the Malpighian tubules, Tabanus griseinus, Atylotus spp., and Hybomitra spp. showed significant consistency (Figures 4A,B,D). Enterobacter, Lactococcus, Caproiciproducens, Faecalibacterium, Acinetobacter, Pseudomonas, Klebsiella, Subdoligranulum, and Stenotrophomonas were enriched in the ovaries (LDA > 3, p < 0.05). Interestingly, the ovarian microbiota of Haematopota turkestanica differed from that of the other tabanids, with bacteria such as Christensenellaceae R 7 Group, Blautia, and Ruminococcaceae UCG 013 being significantly enriched in the ovaries (Figure 4C).

LEfSe was employed to analyze how the microbiota varied across the organs of different species, revealing bacterial composition differences in the Malpighian tubules, midgut, and ovaries of each tabanid species. In the Malpighian tubules, Atylotus sp. was the most prominent genus (Figure 4E), with Streptococcus displaying the highest LDA value and was identified as the most significantly enriched genus (LDA > 3, p < 0.05). Furthermore, the genera most significantly enriched in Tabanus griseinus, Haematopota turkestanica, and Hybomitra sp. were Millionella, Loktanella, and RB41 (LDA >3, p < 0.05). In the midgut samples collected from the Hybomitra sp., yet in the context of an expected microbiome diversity, we observed Spiroplasma being the most significantly enriched (Figure 4G). The uncultured F_Xanthomonadaceae, Saccharibacter, and Lachnospiraceae AC2044 group were the most significantly enriched genera in Tabanus griseinus, Atylotus sp., and Haematopota turkestanica, respectively.

Consistent with the organ-level analysis, Haematopota turkestanica exhibited the most specialized ovarian microbiota, thus having the highest number of marked genera (Figure 4F). F_Lachnospiraceae uncultured and Escherichia/Shigella were the most abundant bacteria in the ovaries. The most significantly enriched genera in the ovaries of Tabanus griseinus, Atylotus spp. and Hybomitra spp. were Ruminococcaceae UCG 010, Ensifer, and Erysipelothrix, respectively.

Putative pathogens within the microbiota of tabanids were identified using a pathogen detection pipeline to assess the potential public health risks associated with tabanids. The MBPD is known for its ability to detect a broad spectrum of animal, plant, and zoonotic pathogens. The results revealed over 100 genera of putative animal pathogens were the most abundant (Figure 5). Among these, 20 genera were annotated as zoonotic agents, with Citrobacter exhibiting the highest abundance among the zoonotic pathogens. Interestingly, Brucella species had varying abundance distributions among the four fly species, highlighting the presence of this important group of zoonotic pathogens (Figure 5).

Research on vector-microbiota interactions has revealed the significant impact that microbiota can have on various biological traits in disease-transmitting vectors. This change in the composition of the vector microbiota can potentially influence pathogenic infections through diverse mechanisms. Microbiota may directly interact with vector pathogens or modulate pathogen infection by regulating host immune defenses and nutritional status. The influence of the microbiota on other aspects of vector physiology and pathogen transmission, as observed in the microbiome of tabanids, is summarized in Table 2.