Skip to main content

The genome of sheep ked (Melophagus ovinus) reveals potential mechanisms underlying reproduction and narrower ecological niches

Abstract

Background

Melophagus ovinus is considered to be of great veterinary health significance. However, little is known about the information on genetic mechanisms of the specific biological characteristics and novel methods for controlling M. ovinus.

Results

In total, the de novo genome assembly of M. ovinus was 188.421 Mb in size (330 scaffolds, N50 Length: 10.666 Mb), with a mean GC content of 27.74%. A total of 13,372 protein-coding genes were functionally annotated. Phylogenetic analysis indicated that the diversification of M. ovinus and Glossina fuscipes took place 72.76 Mya within the Late Cretaceous. Gene family expansion and contraction analysis revealed that M. ovinus has 65 rapidly-evolving families (26 expansion and 39 contractions) mainly involved DNA metabolic activity, transposases activity, odorant receptor 59a/67d-like, IMD domain-containing protein, and cuticle protein, etc. The universal and tightly conserved list of milk protein orthologues has been assembled from the genome of M. ovinus. Contractions and losses of sensory receptors and vision-associated Rhodopsin genes were significant in M. ovinus, which indicate that the M. ovinus has narrower ecological niches.

Conclusions

We sequenced, assembled, and annotated the whole genome sequence of M. ovinus, and launches into the preliminary genetic mechanisms analysis of the adaptive evolution characteristics of M. ovinus. These resources will provide insights to understand the biological underpinnings of this parasite and the disease control strategies.

Peer Review reports

Background

The family Hippoboscidae (Diptera, Hippoboscoidea) is closely associated with mammals and birds in its ecological habits [1] which contain 213 described species until now. Among them, Melophagus ovinus, commonly known as sheep ked, is a serious hematophagous pest with a limited host range, feeding mainly on sheep and goats and occasionally livestock (rabbits and dogs), and wild animals (Tibetan antelope, European bison, and red foxes) [2]. The presence of M. ovinus has been reported in most temperate sheep-rearing areas in Europe, Asia, North America, Africa, and Oceania [3].

M. ovinus is the most studied species of the family Hippoboscidae due to its veterinary importance and significant economic effects. M. ovinus infection results in skin damage, anemia, inflammation, wool loss, and subsequent secondary bacterial infections and cutaneous myiasis [2, 4]. M. ovinus has been reported to serve as potential vectors of zoonotic pathogens including Bluetongue virus, Trypanosoma spp., Bartonella spp., Anaplasma spp., Borrelia burgdorferi, Border disease virus, Rickettsia spp., Theileria spp [5,6,7,8,9,10]. Recent studies on the species have mainly focused on the diversity of pathogenic organisms carried by M. ovinus [3] and gut microbiota and obligate endosymbiont [11, 12]. Host-symbiont interactions between the keds and its obligate endosymbiont seem to influence their vector role [1].

Due to the highly specialized blood-feeding style, M. ovinus has developed a large number of specific adaptions in morphology and physiology. The life cycle of M. ovinus comprises three definite stages: larva, pupa, and wingless adult. Compared with the other hippoboscids such as Glossina morsitans, all life stages of sheep ked occur on the host, being strictly host dependent. Frequent direct contact between host animals makes it relatively easy for the offspring to transfer to a wide variety of hosts. The unique host ecology may explain the complete loss of wings in M. ovinus [13]. The reproductive biology of adenotrophic viviparity is unique to the superfamily Hippoboscoidea. Expansion and adaptation of the female accessory gland and the uterus provide the nutrient synthesis and delivery and the habitat for developing larvae [14, 15]. M. ovinus also adopt such reproductive mode and produce a single 3rd instar larva which immediately molded into pupae. During the whole lifecycle, the female produces only 12–15 progeny. The remarkable reproduction pattern of M. ovinus is considered to be the result of adaptive evolution, making their offspring with higher survival but lower reproductive capacity than other dipteran insects. Similarly, M. ovinus have highly specialized piercing-sucking mouthparts that are used for feeding on host blood and both females and males are strictly blood-feeding [2, 13]. Despite this, it remains unclear which genetic mechanisms form the basis for these long-term adaptive evolution characteristics.

Although the species of Hippoboscoidea are of veterinary and ecological importance, fewer genomic resources have been annotated [15, 16] seriously hindering investigation of the specific biological mechanisms (pathogen vectors and behavior characteristics) from the perspectives at genomic levels. Presently, genome-wide datasets are available for six medical and veterinary important ectoparasites from Glossinidae, which provide insight into the evolutionary biology and vector control strategies [16]. However, little is known about the information on genomes representing most families of Hippoboscidae, including M. ovinus. Here, we assembled a draft genome for M. ovinus and conducted comparative genome analyses with Glossinidae and other dipteran insects. This work gives first insights into the molecular biology of this hematophagous ectoparasite and provides a new resource for future studies of M. ovinus, and comparative genomic and genetic investigations of Hippoboscoidea.

Results and discussion

Genome assembly and annotation

In total, we generated 143.85 Gb of subreads with high sequencing depth (757.09 ×) on the PacBio Sequel System. A total of 21.32 Gb and 8.42 Gb clean data generated from the Illumina HiSeq X Ten platform for genome and transcriptome sequencing were used to correct the primary PacBio assembly. The genome size estimated by genomoscope 2.0 was 183.278 Mb (Figure S2). The de novo genome assembly of M. ovinus was 188.421 Mb in size (330 scaffolds N50 Length: 10.666 Mb), with a mean GC content of 27.74% and a heterozygosity of 0.456% on the basis of k-mer statistics (Table S1). Assembly completeness of M. ovinus (94.5%) was determined by Benchmarking Universal Single-Copy Orthologs (BUSCO) analysis, which demonstrated a highly accurate assembly (Table S2). The M. ovinus genome is slightly larger than that of Drosophila melanogaster (180 Mb) and much smaller than that of Glossina species (315–380 Mb) and Musca domestica (691 Mb) [16,17,18]. Genomic GC content is significantly correlated with biological features of genome organization [19]. The lower GC content in M. ovinus (27.74%) relative to other Glossina species (the GC content ranges from 27 to 35%) [16] could result in lower transcriptional activity of heterochromatic DNA. Furthermore, we assembled a complete mitochondrial (mtDNA) genome sequence of M. ovinus, which showed 99.93% and 99.04% sequences similarity with the previously published sequences (Genebank Accession NO.MH024396 and NC_037368).

A total of 13,372 proteins were identified by the algorithm MAKER with an average number of 3.79 exons and 11.06 introns per gene. Mean length of gene, exon, and intron were 4,477.1, 525.8, and 198.5 bp, respectively (Table 1), smaller than the findings for the genomes of Glossina species. We identified 9,505 genes with homology to proteins in the InterPro database, among which 6,320 were assigned to GO classifications, and 663, 525, and 2,609 genes could be annotated in the KEGG, MetaCyc, and Reactome pathway databases, respectively.

Table 1 Results of de novo genome of M. ovinus

The estimated percentage of masked repeats this draft genome is 27.08%, in contrast to 34.95–39.99% of the Glossina species [16]. The most-abundant repeat types were unclassified repeats (14.28%), DNA/simple repeat (10.53%), DNA/low complexity (1.70%), LINE families (0.44%), LTR retro-elements (0.10%), and DNA/TcMar-Mariner (0.02%) (Table S3). In addition, 330 rRNAs, 57 miRNAs, 28 snRNAs, 3 ribozymes, 153 tRNAs, and 96 other ncRNAs were identified in the assembled genome (Tables 1; S4).

Gene family analysis and evolution

Finally, 17,597 gene families were identified using OrthoFinder, covering 245,501 genes. Among these, 5,185 gene families were shared by all species, 325 were single-copy orthogroups, and 187 gene families containing 571 genes were unique to M. ovinus (Figure S3; Table S5). We then identified 325 shared single-copy orthologues to construct phylogenetic trees. The phylogeny analysis indicated that nine fly species were found to be clustered together into a large branch and strongly supported with high values (bootstrap = 100). M. ovinus was most closely related to Glossina fuscipes (G. fuscipes), and that these two species in turn formed a clade, supporting the sister relationship previously published between Glossinidae and Hippoboscidae [13]. The estimated divergence time between D. melanogaster and other eight fly species was 129.10 Mya, which was consistent with the result of recent studies [20]. To the Hippoboscoidea superfamily, the estimated divergence time was 89.09 Mya. Remarkably, the diversification of M. ovinus and G. fuscipes took place 72.76 Mya within the Late Cretaceous. These results suggested that there were considerable genetic distances between such species.

Gene family expansion and contraction were analyzed using CAFE and gene birth rate was estimated at 0.00120 when accounting for duplications/gene/Mya. This approach revealed a total of 626 gene families for 11 species experienced significant expansion or contraction events (Fig. 1). Among them, M. ovinus has 65 rapidly-evolving families (26 expansion and 39 contractions). The corresponding genes identified from these gene families included 649 expanded genes and a greater degree of 2,983 contracted genes. The expanded families mainly involved DNA metabolic activity, transposases activity, DNA ligase activity, and DNA recombination (Figure S4 and Tables S6-7). In terms of contracted gene families in M. ovinus, these genes encoded odorant receptor 59a/67d-like, IMD domain-containing protein, cuticle protein, troponin, chloride-channel proteins, and fibrinogen C-terminal domain-containing protein, etc. (Table S7).

Fig. 1
figure 1

Comparative genomic analyses among M. ovinus and ten other species

Reproductive genetics of M. ovinus

In tsetse flies, nutrients for intrauterine larval development are synthesized by specialized glands and have been described as a novel milk-specific protein family. To date, thirteen milk proteins in have been characterized in tsetse flies including milk gland proteins 1–10 (mgp1-10), transferrin, acid sphingomyelinase 1, and Peptidoglycan Recognition Protein-LB (PGRP-LB) [16, 21]. A list of milk protein orthologues has been assembled from the genome of M. ovinus, which indicates that milk proteins are universal and tightly conserved in the Hippoboscoidea superfamily.

Search for orthologous sequences to mgp2-10 revealed that M. ovinus lacks an orthologous sequence for the mgp2 gene, while G. brevipalpis lacks the same mgp2 gene. The mgp2-10 gene family of M. ovinus shows a conserved sequential distribution pattern, while localizes to a 20 kb microsyntenic region with gene opposite-direction (Fig. 2A). Phylogenetic analysis showed that these genes can be divided into two distinct groups, one consisting of mgp3, 7, the other of mgp2, 4, 5, 6, 8, 9, 10. Among them, the mgp3, 7, 8, 5, 10 genes of M. ovinus clustered together with the corresponding genes of tsetse flies, but not the mgp6*, 4*, 9* genes, which probably indicated different evolutionary features between M. ovinus and tsetse flies (Fig. 2B). Recent studies showed the mgp2-10 gene family was defined as a highly divergent protein family and lack conserved functional protein domains [16, 21]. Evolutionarily related sequences of M. ovinus reveals a high degree of genetic diversity of milk genes, meanwhile, much more research in other viviparous genera will be helping shed light on their potential origin. Similar to Glossina species, comparative expression analysis of milk genes in male and female specimens showed that female-specific expression profiles (Fig. 2C). Interference with at least two milk proteins has substantial influence over tsetse fecundity and larval development [21,22,23]. In view of important role of these milk proteins in reproduction, functional characterization is warranted, as they might provide novel targets for the control of viviparous insects.

Fig. 2
figure 2

Overview of the conservation of M. ovinus milk protein genes and their expression patterns. A Syntenic analysis of gene structure/conservation in the mgp2-10 genetic locus. B Phylogenetic analysis of orthologs from the mgp2-10 gene family. C Sex-specific expression of M. ovinus milk proteins

In addition, M. ovinus has a single yolk protein gene that is orthologous to Glossina and M. domestica yp2 (Figure S5). In comparison with other fly species such as M. domestica (7 yolk proteins) and D. melanogaster (3 yolk proteins) [24], M. ovinus and Glossina have a marked reduction in yolk genes associated with their reduction in reproductive capacity.

Immune characteristics of M. ovinus

Peptidoglycan recognition proteins (PGRPs) are important pathogen-associated molecular patterns (PAMPs) involved in pathogen molecular recognition in insects and played an important role in natural antibacterial immunity and regulation of intestinal homeostasis [25, 26]. In M. ovinus, only three PGRP genes were identified, two in the long subfamily (PGRP-LB and -LC, and -SA) and one in the short subfamily (PGRP-SA), and Glossina has 6 PGRPs (4 long PGRP-L genes, 2 Short PGRP-S), whereas D. melanogaster has 13 PGRPs (6 long PGRP-L genes, 7 short PGRP-S), and A. gambiae has 7 PGRPs (4 long PGRP-L genes, 3 Short PGRP-S) (Fig. 3A-B). Phylogenetic analysis indicated that the PGRP-LB, PGRP-LC and PGRP-SA of the M. ovinus were clustered with the corresponding genes of the Glossina and D. melanogaster, respectively. Among them, PGRP-LB and PGRP-SA of sheep louse were phylogenically closer to D. melanogaster, while PGRP-LC is closer to Glossina (Fig. 3C). This indicates that PGRPs genes may have different ways in evolution.

Fig. 3
figure 3

Distribution pattern and phylogenetic analysis of PGRPs orthologs. Total number (A), types and distribution (B) of the PGRPs. C Phylogenetic analysis of PGRPs orthologs

The reduced PGRP repertoire of M. ovinus and Glossina might reflect its blood-specific diet, which likely exposes its gut to fewer microbes than other insects. While, the loss of PGRPs such as PGRP-SC1 and PGRP-LE might indicate that the immune system has evolved to protect the symbiosis of M. ovinus. This reduced immune capacity was also observed in tsetse flies and aphids which harbor obligate symbionts [15, 27]. PGRP-LB is a maternally transmitted immune milk protein and plays a pivotal role in tsetse flies by protecting the obligate symbiont Wigglesworthia [28]. In M. ovinus, the PGRP-LB sequence is highly homologous with Glossina PGRP-LB and may have functional relationship. Although phylogenetically distant with the mutualistic Wigglesworthia in Glossina, the obligate symbiont Arsenophonus of M. ovinus displays remarkably similar traits such as transmission via the milk glands and vital function in development and fecundity [29]. Hence, the PGRP-LB-Arsenophonus axis in M. ovinus might an analogy to the PGRP-LB-Wigglesworthia axis in tsetse flies.

Contractions and losses of sensory receptors in M. ovinus

Insects with diversified ecological and biochemical niches impact the size of sensory receptor families with evidence for generalists insects (broad niches) having a larger number of genes compared with specialist (narrower niches) [24, 30]. These gene families encode chemosensory proteins (CSPs), gustatory receptors (GRs), odorant binding proteins (OBPs), odorant receptors (ORs), ionotropic receptors (IRs), and sensory neuron membrane proteins (SNMPs) [31, 32]. Detailed annotation of M. ovinus reveals that they have fewer sensory receptors (22) relative to Glossina (116), D. melanogaster (248), and A. gambiae (305) (Fig. 4). The chemoreceptor of the blood-sucking insect (M. ovinus and Glossina) is lower in number than the polyphagous insect (D. melanogaster and A. gambiae). Through the comparative analysis of homologous genes, it was found that the ORs and IRs of the M. ovinus were significantly lower than other insects. In addition, GRs associated with sweet tastes, present in D. melanogaster and A. gambiae, are missing in M. ovinus and Glossina. These genetic differences are consistent with the combination of a restricted diet of vertebrate blood and their narrow host range [33]. Compared with Glossina, a list of ORs and IRs of the M. ovinus has been lost over evolutionary history. Meanwhile, results from rapidly-evolving families showed that odorant receptor 59a-like and odorant receptor 67d-like were contracted. These results might explain the mechanism of the M. ovinus obligately parasitize sheep, while Glossina have a variety of animal hosts and a diverse living environment.

Fig. 4
figure 4

Comparison of sensory receptors gene homologs between species. Total number (A), types and distribution (B) of the sensory receptors

Vision-associated Rhodopsin genes in M. ovinus

Vision-associated Rhodopsin genes play a critical role in host identification and pursuit and mate-seeking by insects [34, 35]. These color sensing Rhodopsin genes are identified as critical factors in the optimization and development of trap/target technologies [36, 37]. The visual system of Glossina is very similar to that of other Brachycera species, especially detecting colors in the blue wavelength ranges [15]. The search for vision-associated genes across the M. ovinus reveals orthologs to those described in the Glossina species. Unlike Glossina, which has six Rhodopsin genes (Rh1, 2, 3, 5, 6, and 7), M. ovinus has four Rhodopsin genes (Rh1, 3, 6, and 7). Morphology and function of the compound eye retina is relative conserved throughout the Brachycera [15]. Five opsin genes are differentially expressed in the photoreceptors of the D. melanogaster compound eye, four of five in Glossina, M. domestica, and S. calcitrans, three of five in M. ovinus (Fig. 5A). The Rh1 subfamily is deployed in motion vision [16] and shows greater sequence conservation between M. ovinus and Glossina. The loss of one of the two dipteran ultraviolet (UV)-sensitive Rhodopsins (Rh4) is found in M. ovinus, Glossina, M. domestica, and S. calcitrans. The blue sensitive opsin Rh5, which is expressed in color-discriminating inner photoreceptors [16], have been missed in M. ovinus and C. capitata. The common ocellus-specific opsin Rh2 in other insects is not detected because of M. ovinus has no ocellus. In addition, the M. ovinus genome also contains the ortholog of the Rh7, which is recently characterized UV-like and circadian photoreceptor [38, 39]. Phylogenetic analysis suggests the Rhodopsin gene family contracted in M. ovinus after their divergence within the family Hippoboscidae. Besides, the Rhodopsin genes of the M. ovinus were clustered together with the corresponding genes of the Glossina species (Fig. 5B), which were relatively conservative in evolution and more closely related to each other. These genetic differences could account for the attractivity of trap and targets to different insects.

Fig. 5
figure 5

Distribution (A) and phylogenetic analysis (B) of orthologues from the Rhodopsin gene family

Conclusion

In summary, we report a high-quality draft genome of M. ovinus that enabled the molecular landscape of this hematophagous ectoparasite to be characterized for the first time. This 188.421 Mb genome, with a repeat content of 27.08%, encodes 13,372 protein-encoding genes. Conspicuous in the predicted proteome are milk-specific protein family, PGRPs, sensory receptors, and Rhodopsin gene family, some of which appear worthy of investigation as novel targets for vector control strategies. In addition, comparative genomics allowed a re-evaluation of the genetic relationships of M. ovinus and Glossina. In conclusion, this study provides some insight into the genetics underlying the evolution of unique adaptive traits of M. ovinus and a new resource for comparative genomic and genetic explorations of the parasites in the family Hippoboscidae.

Materials and methods

Sample collection, sequencing, and quality control

A total of 15 adult female specimens of M. ovinus were collected from a sheep farm located at Zhidoi County, Yushu Tibetan Autonomous Prefecture, Qinghai province, China (33°37’N, 95°58’E). The species was identified by morphological approaches (Figure S1) and 18S rRNA gene sequence analysis according to previous studies [3]. The specimens were washed extensively in physiological saline and every 5 specimens were pooled with their gut contents removed. The remaining specimens were deposited at the Institute of Zoology, Chinese Academy of Sciences (IZCAS).

Genomic DNA/RNA extraction, library preparation and sequencing were conducted by conducted by the company Berry Genomics (Beijing, China). Total genomic DNA/RNA extraction from pooled keds was conducted according to the instruction of Blood & Cell Culture DNA Kit (Qiagen, USA) and Trizol reagent (Invitrogen, USA), respectively. For long-read sequencing, high quality gDNA (≥ 10 μg, ≥ 100 ng/μL) was prepared for a 20-kb insert SMRTbell library construction using SMRTbell® Express Template Prep Kit 2.0 (Pacific Biosciences, USA) from the pooled DNA isolates and sequenced on the PacBio Sequel platform employing P6-C4 chemistry. For short-read and RNA sequencing, paired-end libraries were constructed with an insert size of 400 bp and sequenced on the Illumina HiSeq X Ten platform (150-nt paired-end reads).

We use clumpify.sh (BBTools suite v37.93, Bushnell) to compressed raw Illumina short reads into clumps and removed duplicates. Quality control was performed with bbduk.sh (BBTools): both sides were trimmed to Q20 based on Phred scores, reads shorter than 15 bp or with more than 5 Ns were discarded, poly-A or poly-T tails of at least 10 bp were trimmed, and overlapping paired reads were corrected.

Genome size estimation

We employed the strategy of short-read k-mer distributions to estimate the genome size using short sequencing data. The histogram of k-mer frequencies was computed with 17-mers using jellyfish [40]. The genome size was then estimated using GenomeScope v2.0 [41].

Genome, mitochondrion and transcriptome assembly

We performed de novo genome assembly with long reads using Canu v2.1.1 [42] and wtdbg2 v2.5 [43]. Both assemblies were first polished by Flye v2.4.2 (–polish_target) [44] on raw Pacbio sequences. To improve assembly contiguity, the two assemblies were merged into one assembly after two rounds of quickmerge v0.3 of USAGE 2 [45], and further polished with Illumina short reads using Pilon v1.22 [46]. Subsequently, we filtered possible contaminants by HS-BLASTN [47] and VecScreen against the NCBI and UniVec database, respectively.

We assembled and annotated the mitochondrial genome of M. ovinus based on Illumina short reads using MitoZ v2.3 [48] and MITOS webserver [49]. Transcriptome assembly was performed under a genome-guied method, RNA-seq reads to our assembled genome using HISAT2 v2.1.0 [50] and then assembled with StringTie v1.3.4 [51]. Redundant isoforms were removed with Redundans v0.13c [52] under default parameters. Finally, Benchmarking Universal Single-Copy Orthologs (BUSCO) [53] was used to evaluate the completeness of all assemblies.

Gene prediction and functional annotation

A custom library combining a de novo species-specific repeat library with RepBase-20181026 and Dfam 3.0 databases was generated by RepeatModeler v1.0.11 [54, 55]. Repetitive sequences were predicted by RepeatMasker v4.0.9 together with the custom library [56]. Gene prediction was conducted with the MAKER v2.31.10 pipeline [57] by integrating ab initio predictions, RNA-seq evidence, and protein homology-based searches. For ab initio annotation, Augustus v3.3 [58] and GeneMark-ET v4.33 [59] were simultaneously performed to predict gene coding regions. Two predictors were trained using BRAKER v2.1.0 [60] with RNA-seq data, and previously genome-guided assembled transcripts were used as transcriptome-based evidence. Protein sequences of Drosophila melanogaster, Glossina fuscipes and Glossina morsitans were downloaded as protein homology-based evidence.

Homology-based gene functions were assigned using Diamond v0.9.18 [61] and the UniProtKB (SwissProt + TrEMBL). Protein domains, Gene Ontology (GO) and pathway annotations, were searched with InterProScan 5.34–73.0 [62] against the Pfam [63], PANTHER [64], Gene3D [65], Superfamily [66], and CDD [67] databases. For non-coding genes, ribosomal RNA (rRNA), microRNA (miRNA), and small nuclear RNA (snRNA) genes were predicted by Infernal v1.1.2 [68] against the Rfam v14.0 [69] database. Transfer RNAs were further refined using tRNAscan-SE v2.0 [70].

Gene family analysis and evolution

We identified gene families using 10 public genome protein sequences of insect species, including Aedes aegypti, Anopheles gambiae, Bactrocera oleae, Ceratitis capitate, Drosophila melanogaster,Glossina fuscipes,Lucilia cuprina, Musca domestica, Sarcophaga bullata, and Stomoxys calcitrans. Aedes aegypti and Anopheles gambiae were selected as the out-group. OrthoFinder v2.2.7 [71] was used to infer orthogroups with Diamond [61] as the sequence aligner. Gene family evolution (gain and loss) was analyzed using CAFE v4.2 [72] with the lambda parameter used to calculate birth and death rates. The ultrametric tree generated from OrthoFinder was transformed using r8s [73] and time-calibrated by the divergence time (35 mya) of the most recent common ancestor between Musca domestica and Stomoxys calcitrans from the TimeTree database [74].

Phylogenetic analysis

The identified sequences of mgp genes, PGRP genes, sensory receptors gene, and vision-associated rhodopsin genes were combined with the corresponding homologous genes of Drosophila melanogaster, Musca domestica, Glossina fuscipes and Glossina morsitans and used as template sequences to homologous blast with the M. ovinus genome database (with e-value as 1e-5) to determine the annotation results of gene families [75]. The protein sequences of the target genes were aligned using MAFFT software [76]. The phylogenetic tree reconstructions were made using maximum likelihood (ML) based on IQ-TREE v2 [77] with 1,000 ultrafast bootstraps. Dendrograms were edited using Interactive Tree of Life (iTOL) v5 [78].

Reverse transcription quantitative PCR (RT-qPCR)

Total RNAs were isolated from female and male keds with Trizol reagent (Invitrogen, USA). cDNA templates were synthesized from 1 μg RNA using the GoScript Reverse Transcription System (Promega, USA). Specific primers for mgp genes and Tubulin were designed with the identified gene sequences by genome analysis (Table S8). The relative expression levels were quantified using the SYBR Green qPCR kit (Cwbio, CHINA) and Applied Biosystems 7500 Fast Real-Time PCR system (Thermo, USA). Expression of target genes were normalized to the reference endogenous Tubulin and calculated by using the 2−ΔΔCt method. Statistical analysis was analyzed by Student’s t-test using GraphPad Prism 7 (GraphPad).

Availability of data and materials

Raw sequencing data are deposited in the NCBI Sequence Read Archive under the accessions SRR17267914-SRR17267916. The genome and transcriptome assemblies are available at GenBank with the accession number SRX13445939, SRX13445940, respectively, corresponding to the BioProject PRJNA789604 and the BioSample SAMN24146076. All supplementary figures and tables are provided as Additional File.

Abbreviations

CSPs:

Chemosensory proteins

GRs:

Gustatory receptors

OBPs:

Odorant Binding Proteins

ORs:

Odorant Receptors

IRs:

Ionotropic Receptors,

SNMPs:

Sensory Neuron Membrane Proteins

References

  1. Marcos ABS, Domenico O. Keds, the enigmatic flies and their role as vectors of pathogens. Acta Trop. 2020;209:105521.

    Article  Google Scholar 

  2. Small RW. A review of Melophagus ovinus, (L.), the sheep ked. Vet Parasitol. 2005;130(1–2):141–55.

    Article  Google Scholar 

  3. Zhang QX, Wang Y, Li Y, Han SY, Wang B, Yuan GH, et al. Vector-Borne pathogens with veterinary and public health significance in Melophagus ovinus (sheep ked) from the Qinghai-Tibet Plateau. Pathogens. 2021;10(2):249.

    Article  CAS  Google Scholar 

  4. Rudolf I, Betasova L, Bischof V, Venclikova K, Blazejova H, Mendel J, et al. Molecular survey of arthropod-borne pathogens in sheep keds (Melophagus ovinus), Central Europe. Parasitol Res. 2016;115(10):3679–82.

    Article  Google Scholar 

  5. Chu CY, Jiang BG, Qiu EC, Zhang F, Zuo SQ, Yang H, et al. Borrelia burgdorferi sensu lato in sheep keds (Melophagus ovinus), Tibet. China Vet Microbiol. 2011;149(3–4):526–9.

    Article  Google Scholar 

  6. Gibson W, Pilkington JG, Pemberton JM. Trypanosoma melophagium from the sheep ked Melophagus ovinus on the island of St Kilda. Parasitology. 2010;137(12):1799.

    Article  CAS  Google Scholar 

  7. Hao LL, Yuan DB, Li SH, Jia T, Guo L, Hou W, et al. Detection of Theileria spp. in ticks, sheep keds (Melophagus ovinus), and livestock in the eastern Tibetan Plateau. China Parasitol Res. 2020;119(8):2641–8.

    Article  Google Scholar 

  8. Liu YH, He B, Li KR, Li F, Zhang LY, Li XQ, et al. First report of border disease virus in Melophagus ovinus (sheep ked) collected in Xinjiang, China. PLoS ONE. 2019;14(8):e0221435.

    Article  CAS  Google Scholar 

  9. Zhao L, Wang JL, Ding YL, Li KR, He B, Li F. Theileria ovis (Piroplasmida: Theileriidae) detected in Melophagus ovinus (Diptera: Hippoboscoidea) and Ornithodoros lahorensis (Ixodida: Argasidae) removed from sheep in Xinjiang. China J Med Entomol. 2020;57(2):631–5.

    Article  CAS  Google Scholar 

  10. Lu M, Tian J, Zhao H, Jiang H, Qin X, Wang W, et al. Molecular survey of vector-borne pathogens in ticks, sheep keds, and domestic animals from Ngawa, Southwest China. Pathogens. 2022;11(5):606.

    Article  CAS  Google Scholar 

  11. Duan DY, Zhou HM, Cheng TY. Comparative analysis of microbial community in the whole body and midgut from fully engorged and unfed female adult Melophagus ovinus. Med Vet Entomol. 2020;34(2):215–24.

    Article  Google Scholar 

  12. Husnik F, Hypsa V, Darby A. Insect-symbiont gene expression in the midgut bacteriocytes of a blood-sucking parasite. Genome Biol Evol. 2020;12(4):429–42.

    Article  CAS  Google Scholar 

  13. Petersen FT, Meier R, Kutty SN, Wiegmann BM. The phylogeny and evolution of host choice in the Hippoboscoidea (Diptera) as reconstructed using four molecular markers. Mol Phylogenet Evol. 2007;45(1):111–22.

    Article  CAS  Google Scholar 

  14. Benoit JB, Attardo GM, Baumann AA, Michalkova V, Aksoy S. Adenotrophic viviparity in tsetse flies: potential for population control and as an insect model for lactation. Annu Rev Entomol. 2015;60:351–71.

    Article  CAS  Google Scholar 

  15. IGGI. Genome sequence of the tsetse fly (Glossina morsitans): vector of African trypanosomiasis. Science. 2014;344(6182):380–6.

    Article  Google Scholar 

  16. Attardo GM, Abd-Alla AMM, Acosta-Serrano A, Allen JE, Bateta R, Benoit JB, et al. Comparative genomic analysis of six Glossina genomes, vectors of African trypanosomes. Genome Biol. 2019;20(1):1–31.

    Article  Google Scholar 

  17. Adams MD, Celniker SE, Holt RA, Evans CA, Gocayne JD, Amanatides PG, et al. The genome sequence of Drosophila melanogaster. Science. 2000;287(5461):2185–95.

    Article  Google Scholar 

  18. Scott JG, Warren WC, Beukeboom LW, Bopp D, Clark AG, Giers SD, et al. Genome of the house fly, Musca domestica L., a global vector of diseases with adaptations to a septic environment. Genome Biol. 2014;15(10):466.

    Article  Google Scholar 

  19. Bernardi G. The neoselectionist theory of genome evolution. Proc Natl Acad Sci USA. 2007;104(20):8385–90.

    Article  CAS  Google Scholar 

  20. Ren LP, Shang YJ, Yang L, Wang SW, Wang X, Chen S, et al. Chromosome-level de novo genome assembly of Sarcophaga peregrina provides insights into the evolutionary adaptation of flesh flies. Mol Ecol Resour. 2021;21(1):251–62.

    Article  CAS  Google Scholar 

  21. Benoit JB, Attardo GM, Michalkova V, Krause TB, Bohova J, Zhang Q, et al. A novel highly divergent protein family identified from a viviparous insect by RNA-seq analysis: a potential target for tsetse fly-specific abortifacients. PLoS Genet. 2014;10(4):e1003874.

    Article  Google Scholar 

  22. Benoit JB, Attardo GM, Michalkova V, Takac P, Bohova J, Aksoy S. Sphingomyelinase activity in mother’s milk is essential for juvenile development: a case from lactating tsetse flies. Biol Reprod. 2012;87(1):1–10.

    Article  Google Scholar 

  23. Yang G, Attardo GM, Lohs C, Aksoy S. Molecular characterization of two novel milk proteins in the tsetse fly (Glossina morsitans morsitans). Insect Mol Biol. 2010;19(2):253–62.

    Article  CAS  Google Scholar 

  24. Olafson PU, Aksoy S, Attardo GM, Buckmeier G, Chen X, Coates CJ, et al. The genome of the stable fly, Stomoxys calcitrans, reveals potential mechanisms underlying reproduction, host interactions, and novel targets for pest control. BMC Biol. 2021;19(1):41.

    Article  CAS  Google Scholar 

  25. Dyer NA, Rose C, Ejeh NO, Acosta-Serrano A. Flying tryps: survival and maturation of trypanosomes in tsetse flies. Trends Parasitol. 2013;29(4):188–96.

    Article  Google Scholar 

  26. Lu Y, Su FH, Li QL, Zhang J, Li YJ, Tang T. Pattern recognition receptors in Drosophila immune responses. Dev Comp Immunol. 2020;102:103468.

    Article  CAS  Google Scholar 

  27. Elsik CG. The pea aphid genome sequence brings theories of insect defense into question. Genome Biol. 2010;11(2):106.

    Article  Google Scholar 

  28. Wang J, Aksoy S. PGRP-LB is a maternally transmitted immune milk protein that influences symbiosis and parasitism in tsetse’s offspring. Proc Natl Acad Sci USA. 2012;109(26):10552–7.

    Article  CAS  Google Scholar 

  29. Novakova E, Husnik F, Sochova E, Hypsa V. Arsenophonus and Sodalis symbionts in louse flies: an analogy to the Wigglesworthia and Sodalis system in tsetse flies. Appl Environ Microbiol. 2015;81(18):6189–99.

    Article  CAS  Google Scholar 

  30. Robertson HM. Molecular evolution of the major arthropod chemoreceptor gene families. Annu Rev Entomol. 2019;64:227–42.

    Article  CAS  Google Scholar 

  31. Liu R, He X, Lehane S, Lehane M, Hertz-Fowler C, Berriman M. Expression of chemosensory proteins in the tsetse fly Glossina morsitans morsitans is related to female host-seeking behaviour. Insect Mol Biol. 2012;21(1):41–8.

    Article  Google Scholar 

  32. Obiero GF, Mireji PO, Nyanjom SR, Christoffels A, Robertson HM, Masiga DK. Odorant and gustatory receptors in the tsetse fly Glossina morsitans morsitans. PLoS Negl Trop Dis. 2014;8(4):e2663.

    Article  Google Scholar 

  33. Kabaka JM, Wachira BM, Mangera CM, Rono MK, Hassanali A, Okoth SO. Expansions of chemosensory gene orthologs among selected tsetse fly species and their expressions in Glossina morsitans morsitans tsetse fly. PLoS Negl Trop Dis. 2020;14(6):e0008341.

    Article  CAS  Google Scholar 

  34. Gibson G, Torr SJ. Visual and olfactory responses of haematophagous Diptera to host stimuli. Med Vet Entomol. 1999;13(1):2–23.

    Article  CAS  Google Scholar 

  35. Yuval B. Mating systems of blood-feeding flies. Annu Rev Entomol. 2006;51:413–40.

    Article  CAS  Google Scholar 

  36. Allan SA, Day JF, Edman JD. Visual ecology of biting flies. Annu Rev Entomol. 1987;32:297–316.

    Article  CAS  Google Scholar 

  37. Lindh JM, Goswami P, Blackburn RS, Arnold SE, Vale GA, Lehane MJ. Optimizing the colour and fabric of targets for the control of the tsetse fly Glossina fuscipes fuscipes. PLoS Negl Trop Dis. 2012;6(5):e1661.

    Article  Google Scholar 

  38. Leung NY, Thakur DP, Gurav AS, Kim SH, Di Pizio A, Niv MY. Functions of opsins in Drosophila taste. Curr Biol. 2020;30(8):1367–79.

    Article  CAS  Google Scholar 

  39. Zhang SF, Kong XB, Liu F, Zhang Z. Identification and expression patterns of opsin genes in a forest insect, Dendrolimus punctatus. Insects. 2020;11(2):116.

    Article  Google Scholar 

  40. Marcais G, Kingsford C. A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. Bioinformatics. 2011;27(6):764–70.

    Article  CAS  Google Scholar 

  41. Ranallo-Benavidez TR, Jaron KS, Schatz MC. GenomeScope 2.0 and Smudgeplot for reference-free profiling of polyploid genomes. Nat Commun. 2020;11(1):1432.

    Article  CAS  Google Scholar 

  42. Koren S, Walenz BP, Berlin K, Miller JR, Bergman NH, Phillippy AM. Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. Genome Res. 2017;27(5):722–36.

    Article  CAS  Google Scholar 

  43. Ruan J, Li H. Fast and accurate long-read assembly with wtdbg2. Nat Methods. 2020;17(2):155–8.

    Article  CAS  Google Scholar 

  44. Kolmogorov M, Yuan J, Lin Y, Pevzner PA. Assembly of long, error-prone reads using repeat graphs. Nat Biotechnol. 2019;37(5):540–6.

    Article  CAS  Google Scholar 

  45. Chakraborty M, Baldwin-Brown JG, Long AD, Emerson JJ. Contiguous and accurate de novo assembly of metazoan genomes with modest long read coverage. Nucleic Acids Res. 2016;44(19):e147.

    Google Scholar 

  46. Walker BJ, Abeel T, Shea T, Priest M, Abouelliel A, Sakthikumar S, et al. Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. PLoS ONE. 2014;9(11):e112963.

    Article  Google Scholar 

  47. Chen Y, Ye W, Zhang Y, Xu Y. High speed BLASTN: an accelerated MegaBLAST search tool. Nucleic Acids Res. 2015;43(16):7762–8.

    Article  CAS  Google Scholar 

  48. Meng GL, Li YY, Yang CT, Liu SL. MitoZ: a toolkit for animal mitochondrial genome assembly, annotation and visualization. Nucleic Acids Res. 2019;47(11):e63.

    Article  CAS  Google Scholar 

  49. Bernt M, Donath A, Juhling F, Externbrink F, Florentz C, Fritzsch G, et al. MITOS: improved de novo metazoan mitochondrial genome annotation. Mol Phylogenet Evol. 2013;69(2):313–9.

    Article  Google Scholar 

  50. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    Article  CAS  Google Scholar 

  51. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.

    Article  CAS  Google Scholar 

  52. Pryszcz LP, Gabaldon T. Redundans: an assembly pipeline for highly heterozygous genomes. Nucleic Acids Res. 2016;44(12):e113.

    Article  Google Scholar 

  53. Waterhouse RM, Seppey M, Simao FA, Manni M, Ioannidis P, Klioutchnikov G, et al. BUSCO applications from quality assessments to gene prediction and phylogenomics. Mol Biol Evol. 2018;35(3):543–8.

    Article  CAS  Google Scholar 

  54. Bao W, Kapitonov VV, Pavlicek A, Klonowski P, Kohany O, Walichiewicz J. Repbase update, a database of repetitive elements in eukaryotic genomes. Mob DNA. 2015;6(1–4):1–6.

    Google Scholar 

  55. Hubley R, Finn RD, Clements J, Eddy SR, Jones TA, Bao W, et al. The Dfam database of repetitive DNA families. Nucleic Acids Res. 2016;44(D1):D81–9.

    Article  CAS  Google Scholar 

  56. Smit AFA, Hubley R, Green P. RepeatMasker Open-4.0. http://www.repeatmasker.org. (2013–2015).

  57. Holt C, Yandell M. MAKER2: an annotation pipeline and genome-database management tool for second-generation genome projects. BMC Bioinformatics. 2011;12:491.

    Article  Google Scholar 

  58. Stanke M, Steinkamp R, Waack S, Morgenstern B. AUGUSTUS: a web server for gene finding in eukaryotes. Nucleic Acids Res. 2004;32:W309–12.

    Article  CAS  Google Scholar 

  59. Lomsadze A, Ter-Hovhannisyan V, Chernoff YO, Borodovsky M. Gene identification in novel eukaryotic genomes by self-training algorithm. Nucleic Acids Res. 2005;33(20):6494–506.

    Article  CAS  Google Scholar 

  60. Hoff KJ, Lange S, Lomsadze A, Borodovsky M, Stanke M. BRAKER1: unsupervised RNA-Seq-based genome annotation with GeneMark-ET and AUGUSTUS. Bioinformatics. 2016;32(5):767–9.

    Article  CAS  Google Scholar 

  61. Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2015;12(1):59–60.

    Article  CAS  Google Scholar 

  62. Finn RD, Attwood TK, Babbitt PC, Bateman A, Bork P, Bridge AJ, et al. InterPro in 2017-beyond protein family and domain annotations. Nucleic Acids Res. 2017;45(D1):D190–9.

    Article  CAS  Google Scholar 

  63. Finn RD, Miller BL, Clements J, Bateman A. iPfam: a database of protein family and domain interactions found in the Protein Data Bank. Nucleic Acids Res. 2014;42:D364–73.

    Article  CAS  Google Scholar 

  64. Mi H, Huang X, Muruganujan A, Tang H, Mills C, Kang D, et al. PANTHER version 11: expanded annotation data from Gene Ontology and Reactome pathways, and data analysis tool enhancements. Nucleic Acids Res. 2017;45(D1):D183–9.

    Article  CAS  Google Scholar 

  65. Lewis TE, Sillitoe I, Dawson N, Lam SD, Clarke T, Lee D, et al. Gene3D: extensive prediction of globular domains in proteins. Nucleic Acids Res. 2018;46(D1):D435–9.

    Article  CAS  Google Scholar 

  66. Wilson D, Pethica R, Zhou Y, Talbot C, Vogel C, Madera M, et al. SUPERFAMILY-sophisticated comparative genomics, data mining, visualization and phylogeny. Nucleic Acids Res. 2009;37:D380–6.

    Article  CAS  Google Scholar 

  67. Marchler-Bauer A, Bo Y, Han L, He J, Lanczycki CJ, Lu S, et al. CDD/SPARCLE: functional classification of proteins via subfamily domain architectures. Nucleic Acids Res. 2017;45(D1):D200–3.

    Article  CAS  Google Scholar 

  68. Nawrocki EP, Eddy SR. Infernal 1.1: 100-fold faster RNA homology searches. Bioinformatics. 2013;29(22):2933–5.

    Article  CAS  Google Scholar 

  69. Kalvari I, Argasinska J, Quinones-Olvera N, Nawrocki EP, Rivas E, et al. Rfam 13.0: shifting to a genome-centric resource for non-coding RNA families. Nucleic Acids Res. 2018;46(D1):D335–42.

    Article  CAS  Google Scholar 

  70. Lowe TM, Eddy SR. tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res. 1997;25(5):955–64.

    Article  CAS  Google Scholar 

  71. Emms DM, Kelly S. OrthoFinder: solving fundamental biases in whole genome comparisons dramatically improves orthogroup inference accuracy. Genome Biol. 2015;16(1):1–14.

    Article  CAS  Google Scholar 

  72. Han MV, Thomas GW, Lugo-Martinez J, Hahn MW. Estimating gene gain and loss rates in the presence of error in genome assembly and annotation using CAFE 3. Mol Biol Evol. 2013;30(8):1987–97.

    Article  CAS  Google Scholar 

  73. Sanderson MJ. r8s: inferring absolute rates of molecular evolution and divergence times in the absence of a molecular clock. Bioinformatics. 2003;19(2):301–2.

    Article  CAS  Google Scholar 

  74. Kumar S, Stecher G, Suleski M, Hedges SB. TimeTree: a resource for timelines, timetrees, and divergence times. Mol Biol Evol. 2017;34(7):1812–9.

    Article  CAS  Google Scholar 

  75. He P, Wang MM, Wang H, Ma YF, Yang S, et al. Genome-wide identification of chemosensory receptor genes in the small brown planthopper, Laodelphax striatellus. Genomics. 2020;112(2):2034–40.

    Article  CAS  Google Scholar 

  76. Rozewicki J, Li S, Amada KM, Standley DM, Katoh K. MAFFT-DASH: integrated protein sequence and structural alignment. Nucleic Acids Res. 2019;47(W1):W5–10.

    CAS  Google Scholar 

  77. Minh BQ, Schmidt HA, Chernomor O, Schrempf D, Woodhams MD, von Haeseler A, et al. IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Mol Biol Evol. 2020;37(5):1530–4.

    Article  CAS  Google Scholar 

  78. Letunic I, Bork P. Interactive Tree of Life (iTOL) v5: an online tool for phylogenetic tree display and annotation. Nucleic Acids Res. 2021;49(W1):W293–6.

    Article  CAS  Google Scholar 

Download references

Funding

This study was financially supported by the National Key R&D Program of China (2022YFC2601602), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA19050204), the Regular Assistance Project of International Department of the Ministry of Science and Technology of China (KY201904013), the Major Program of National Natural Science Foundation of China (No. 32090023), National Forestry and Grassland Administration, China, and BJAST Budding Talent Program.

Author information

Authors and Affiliations

Authors

Contributions

HX. H. designed the study and reversed the manuscript. QX. Z. and QS. Z. analyzed data and wrote the manuscript. SY. H, Y. L, and Y. W. performed sample collections and methodology experiments. The author(s) read and approved the final manuscript.

Corresponding author

Correspondence to Hongxuan He.

Ethics declarations

Ethics approval and consent to participate

Not applicable. All samples had been collected in the context of previous studies. The authors declare that all methods were carried out in accordance with relevant guidelines and regulations.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1:

Figure S1. Morphological characteristics of M. ovinus. A: Dorsal view of the female; B: Ventral view of the female; C: Dorsal view of the male; D: Ventral view of the male. Figure S2. GenomeScope results for M. ovinus. A: untransformed linear plot; B: untransformed log plot; C: transformed linear plot; D: transformed log plot. Figure S3. A comparative representation of orthologous and paralogous genes aligned with other ten insect genomes. Figure S4. Rapidly evolving gene families (gene family expansion). Figure S5. Phylogenetic analysis of yolk protein genes from M. ovinus and ten other species.

Additional file 2:

Table S1. Summary of each assembly at each step for M. ovinus. Table S2. GenomeScope genome size estimates for M. ovinus with different combinations of k-mer length and maximum k-mer coverage. Table S3. Repeat annotation in M. ovinus. Table S4. Summary statistics for non-coding RNAs in M. ovinus. Table S5. Summary statistics for gene family inference and gene counts for each family and each species. Table S6. Gene family evolution reported from CAFE. Table S7. Rapidly evolving gene families in M. ovinus. Table S8. Primers used in the present study.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhang, Q., Zhou, Q., Han, S. et al. The genome of sheep ked (Melophagus ovinus) reveals potential mechanisms underlying reproduction and narrower ecological niches. BMC Genomics 24, 54 (2023). https://doi.org/10.1186/s12864-023-09155-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-023-09155-1

Keywords