Insights into the phylogenetic relationship and drug targets of Babesia infective to small ruminants from its mitochondrial genomes

Background: Babesiosis, a tick-borne disease caused by protozoa of the genus Babesia , is widespread in subtropical and tropical countries. Mitochondrion is an essential organelle that is responsible for energy transduction and metabolism, calcium homeostasis and cell signaling. Mitochondrial genomes could provide a new insight to comprehend and investigate the biological features, genetic evolution and classification of the causative agent. Nevertheless, there are limited data on the mitochondrial genomes of ovine Babesia spp. in China. Methods: Herein, we sequenced, assembled and annotated the mitochondrial genomes of six ovine Babesia isolates, and analyzed genome size, gene content, genome structure and cytochrome b ( cob ) amino acid sequences, and performed comparative mitochondrial genomics and phylogenomic analyses among apicomplexan parasites. Results: The mitochondrial genomes range from 5767 to 5946 bp in length with a linear form and contain three protein-encoding genes, cytochrome c oxidase I ( cox1 ), cytochrome c oxidase III ( cox3 ) and cob , six large subunit rRNA gene (LSU), and two terminal inverted repeats (TIR) on both ends. The cob gene sequence analysis indicated that the binding site of anti- Babesia drugs targeted on the cytochrome bc 1 complex. Babesia microti and Babesia rodhaini have a dual flip-flop inversion of 184-1082 bp, whereas other Babesia spp. and Theileria spp. have one pair of TIR, 25-1563 bp. Phylogenetic analysis indicated that six ovine Babesia isolates were divided into two clades, Babesia sp. and Babesia motasi . B. motasi isolates were further separated into two small clades ( B. motasi Hebei/Ningxian and B. motasi Tianzhu/Lintan). Conclusions: The data provided new insights into the taxonomic relationships and drug targets of apicomplexan parasites.

Mitochondrion is a crucial organelle that is responsible for energy transduction and metabolism, calcium homeostasis and cell signaling [16]. The mitochondrial cytochrome bc1 complex, known as complex III, is a multimeric enzyme that is an indispensable element of respiratory chain and energy conversion [17]. The complex III comprises three redox active subunits, cytochrome b (encoded by mitochondrial genome), cytochrome c1 and the Rieske iron-sulfur protein. The catalytic cycle of cytochrome bc1 (Q cycle) is in quinone reduction (Qi) and quinol oxidation (Qo) sites that are mainly constituted by cytochrome b. Inhibitors of bc1 complex may be divided into two types that act on Qi and Qo sites according to the binding site. Therefore, the bc1 complex has been considered as a promising target to develop the anti-microbe compounds. Some studies suggest that yeast bc1 complex structure could be used as a model for discovering new antimalarial drugs [18,19]. However, the mutations in the cytochrome b (cob) gene have been detected to be the molecular basis of drug resistance of some anti-microorganisms [19][20][21][22]. It has demonstrated that several anti-Babesia drugs are ineffective owing to drug resistance, for instance diminazene aceturate and atovaquone that possibly inhibit mitochondrial respiratory activity and electron transport [23][24][25]. Therefore, it needs to well understand the mechanism of resistance of these anti-Babesia drug and develop new drugs.  [26][27][28][29][30][31]. However, there are limited data on the mitochondrial genomes of ovine Babesia spp. in China. Though the mitochondrial genomes of BmLT and BspXJ were sequenced by using Illumina sequencing technology, they have not been sequenced using Sanger dideoxy chain-termination method to verify [32]. In this study, we sequenced six mitochondrial genomes of Babesia isolates infective to small ruminants, and the assembled and annotated sequences were submitted to GenBank after comparing with those of other piroplasms. And they were used to clarify phylogenetic relationships and classification of the babesiae, and to determine novel molecular markers for identification studies of Babesia species. The current study will provide valuable information for understanding of mitogenome evolution among apicomplexan parasites, identifying diagnostic markers and screening drug target.

Parasites and isolation of genomic DNA
The purified merozoites of six ovine Babesia were provided by our laboratory [33]. Genomic DNA extraction, concentration measure and quality evaluation were according to the previously described by Wang et al. [34]. All the genomic DNAs were kept at -20 until using.

Amplification and sequencing
The PCR primers were designed by aligning with the reported genomic sequences of B. bovis and B. bigemina (AB499088 and AB499085) (Table S1). Mitochondrial genome fragments were amplified and cloned into pGEM ® -T Easy Vector (Promega, Beijing, China) for subsequent sequencing by Sangon Biotech company.

Genome assembly and annotation
Mitochondrial genomic fragments were assembled with the CLC Genomics Workbench v.7.5.1. The mitogenome annotation was performed using the DOGMA [35] and Artemis [36], followed by application of BLAST (http://blast.ncbi.nlm.nih.gov/Blast.cgi) to identify homologous proteins in apicomplexan parasites in the GenBank database. The tRNA genes were searched using tRNAscan-SE v.2.0 with the default search mode and other mitochondrial sequence sources [37]. The rRNA genes were annotated by searching previously reported rRNA sequences of B. bovis, B. microti, B. orientalis, and T. parva. The genome comparisons were aligned using Mauve [38]. The sequences were deposited in GenBank (accession numbers: MK962313, MK962314 and MN605889-MN605892).

Sequence alignment and amino acid conservation of cob gene
We referred to the previous researches work [18][19][20][21][22][39][40][41][42][43]] that use mitochondrial cytochrome bc1 complex from S. cerevisiae to explicate the binding residues of the inhibitor (acting on bc1 complex) by X-ray crystallography, spectroscopy and cob sequence alignment and analysis. In the previous reported, sequence alignments of the cob amino acid residue from Saccharomyces cerevisiae, Plasmodium falciparum, Toxoplasma gondii, B. microti and Bos taurus indicated that the drugs binding residues of COB is conservation. Therefore, we used MegAlign software to compare the COB sequences of S. cerevisiae, Bos taurus, T. gondii, Theileria parva, B. microti, B. duncani, P. falciparum, BspXJ/DH, BmLT/TZ and BmNX/HB to found main resistance mutations and drug binding residues in the genus Babesia.

Phylogenetic analysis
The concatenated sequences of COB and cytochrome c oxidase I (cox1) amino acid residues of 26 apicomplexan parasites (Table 1) were aligned using Clustal W of MEGA v.6.06 (http://www.megasoftwa re.net/) software. Subsequently, phylogenetic tree was constructed using MEGA v.6.06 with the maximum likelihood (ML) or neighbour-joining (NJ) analysis based on the JTT with Freqs model. T. gondii and Plasmodium spp. were used as the outgroup. Furthermore, the phylogenetic analysis of the whole mitochondrial nucleotide sequence was conducted by the ML method using the Kimura 2-parameter model. The consensus trees were created after bootstrap analyses with 1000 replications.

Sequence analysis
Sequence analysis revealed that the ovine Babesia mitochondrial genomes were linear DNA of 5767 to 5946 bp, with 70.05-70.87% of A + T content (Table 1). Mitochondrial genomes of six ovine Babesia contained three protein-encoding genes, cox1, cytochrome c oxidase III (cox3), cob, six large subunit rRNA gene (LSU), and two terminal inverted repeats (TIR). The transcriptional direction of cox3, LSU3, LSU6, LSU2, cob and LSU5 was from 3′ to 5′, whereas the direction of cox1, LSU1 and LSU4 was from 5′ to 3′ ( Figure 1; Figure S1). The start codons of cox3 and cob genes of ovine Babesia were ATA and ATG, respectively. The initiation codon of cox1 gene of BspXJ/DH was ATA, whereas that of B. motasi was ATG. Most of protein-coding genes had TAA as a termination codon, followed by TGA ( Table 2).

Comparison of mitochondrial genomes sequenced by the Sanger and Illumina technology
The aligment of BspXJ and BmLT mitochondrial genomes sequenced using Illumina method in previous study [32]and Sanger technique in this study were performed. The results showed that two sets of data from Sanger and Illumina have some differences in the size of the mitochondrial genome, the A + T contents, the number of LSU, and the start and stop codons ( Table 1, Table 2) [32]. Furthermore, there were base differences at several positions between the two sets of sequences (Table 3).

Cob gene sequence analysis
The results of amino acid sequence alignment indicated that all COB contain the highly conserved PEWY motif. The residue at position Leu275 of S. cerevisiae is a key determinant of efficacy of atovaquone and myxothiazol binding to the bc1 complex. This position is occupied by a Leu in the T. parva, B. microti, B. duncani, BspXJ/DH, BmLT/TZ and BmHB/NX sequence, whereas a Phe is present in T. gondii, bovine and P. falciparum sequence ( Figure 2). The inhibitors of Qi and Qo sites include atovaquone, stigmatellin, myxothiazol, endochin-like quinolone (ELQ), antimycin A and NQNO. The amino acid changes conferring atovaquone resistance in the yeast numbering system included five mutations (I269M, F278I/A, Y279C/S, L275F and L282V). The drugs (target protein is Qo site of bc1 complex) binding residues of COB of ovine Babesia included Met128, Gly132, Glu259, Leu262, Phe265 and Tyr266. And drugs (target protein is Qi site of bc1 complex) binding residues of COB of ovine Babesia involved His187, Ser191 and Asp214 (Table 4).

Phylogenetic analysis
Phylogenetic trees were constructed with the concatenated amino acid sequences of cob and cox1 using ML and NJ methods. The two approaches showed no significant changes in the topology. The piroplasms were divided into seven groups  Figure 3). The phylogenetic tree using the whole mitochondrial nucleotide sequence was constructed by the ML method based on the Kimura 2-parameter model. The result was similar to that using concatenated amino acid sequences of cob and cox1 with the exception of T. orientalis that was located in the B. conradae clade ( Figure  S2).

Discussion
In the present study, we assembled and annotated the mitochondrial genomes of six ovine Babesia isolates and performed a mitochondrial genomic analysis with published mitochondrial genomes of apicomplexan parasites. The mitochondrial genomes of six Babesia infective to small ruminants are highly similar to most of Babesia spp. in respect of genome size, high A + T content, genome form, and gene content. However, they are smaller than those of B. microti, B. rodhaini and Theileria equi, and larger than Toxoplasma gondii [28,29,44]. It's consistent with other apicomplexan parasites that the tRNA genes were not found in the mitochondrial genomes. We speculate that they may be directly encoded by the nuclear genome. The order and transcriptional direction of three protein-encoding genes are the same in B. bovis, B. bigemina, B. gibsoni [28], but different from B. microti, T. equi, T. orientalis and P. falciparum [28,29].
In the reported mitochondrial genomes of piroplasma, B. microti and B. rodhaini have a dual flip-flop inversion system, range from 184-1082 bp in length [29]. However, one pair of TIR was found at the mitochondrial genomes of other Babesia spp. and Theileria spp., range from 25-1563 bp in length [26-28, 30, 32]. These findings indicated that the number and size of TIRs are one of the main causes of different mitochondrial genome sizes. We also found that different numbers of LSUs are responsible for the size of the mitochondrial genomes. The TIRs are considered to play a crucial role in the replication and stabilization of the linear mitochondrial genomes [45]. In the published mitochondrial genome of apicomplexan parasites, the sequences of the coding genes and LSU are basically the same. Therefore, we speculated that differences in lengths and sequences of TIRs may be responsible for divergence in host-specific and in vitro culture characteristics of protozoa.
The difference between the Illumina and Sanger sequencing data is mainly caused by the nucleotide substitutions, deletions and insertions, which results from the using of different techniques of sequencing and software of assembly and annotation. And settings of parameters are different, which is one of the reasons for the difference between the two sets of data. The results of Illumina sequencing method are more prone to errors. The Sanger approach is more accurate and thus more appropriate for further studies of small genome sequencing.
Mitochondria are essential organelles and play an important role of energy metabolism, growth and development of apicomplexan protozoa [46]. Cytochrome bc1 is an integral membrane protein complex, which is vital to cellular respiration. The highly conserved binding site of inhibitors in cytochrome bc1 is the molecular basis of the drug effect on yeast, fungi and parasites [43]. The drug binding residues in COB of six ovine Babesia are completely consistent. In the conserved PEWY region, Phe278 (yeast number of cob) is present in most organism, whereas T. gondii and B. taurus at this position are Tyr and Ala, respectively. Studies have reported that the L275F mutation in yeast has no effect on enzyme activity, but IC50 has increased ten-fold [22]. Compared with the wild type of yeast, the mutation Y279S indicated a B40-fold increased IC50 (B1.7 mM) for atovaquone [22], and Y268S in the P. falciparum numbering system results in a 3,000-fold loss of atovaquone sensitivity. In addition, the mutations M133I, L271V and Y268N of P. berghei confer resistance to atovaquone [43]. Therefore, we conclude that the mutations of COB are largely responsible for the efficacy of drugs (target protein is bc1 complex) in apicomplexan parasites.
Currently, atovaquone and ELQs have been reported for the treatment of human babesiosis and malarial by modifying of the drug target by disruption of cytochrome bc1 complex [19,21,22,41,43]. In 2019, a B. motasilike parasite has been detected in human blood in Korea [47], which suggests that B. motasi may be potentially zoonotic. Therefore, we should investigate the infection of B. motasi in human in China, evaluate the zoonotic potential of B. motasi, and the effect of inhibitors binding to cytochrome bc1 complex. Our data could show that atovaquone, stigmatellin, myxothiazol, endochin-like quinolone (ELQ), antimycin A and NQNO drugs can be used in the treatment of babesiosis in future. The molecular mechanism of the resistance of these drugs is the mutation of the COB, which it suggests that a combined drug strategy is possible for avoiding drug resistance process during babesiosis therapy.
In this study, the taxonomical relationships of Babesia infective to sheep and goats are consistent with the reported phylogenetic analyses based on cob, cox1, cox3, nuclear small subunit (SSU) and internal transcribed spacer (ITS) [5,6,8,32,48]. The ovine Babesia are divided into two species: Babesia sp. and B. motasi. B. motasi further fell into two small clades, named BmLT/TZ and BmNX/HB. With the exception of B. conradae, piroplasm infective to the same host fell into one clade. These findings are consistent with the phylogenetic tree of B. gibsoni, B. duncani and B. orientalis based on the amino acid sequences of cox1 and cob [26,27,30].

Conclusions
In conclusion, we reported the mitochondrial genome of six ovine Babesia infective to small ruminants in China.

Ethics approval and consent to participate
Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
The datasets are included in the article and its additional files. The sequences were submitted to the GenBank database, accession numbers MK962313, MK962314 and MN605889-MN605892.

Competing interests
The authors declare that they have no competing interests. Tables Table 1 Mitochondrial genome sequences of apicomplexan parasites used in the present study b b1 and b2 is the same sample, which sequenced using Sanger method and Illumina method, respectively. Table 2 Gene contents and initiation and termination codons of encoding genes of the piroplasma mitochondrial genomes Table 3 Comparison of Babesia sp. Xinjiang and Babesia motasi Lintan mitochondrial genome were sequenced using Illumina (BspXJ-Illumina and BmLT-Illumina) and the Sanger method (BspXJ-Sanger and BmLT-Sanger)  Table 4 Main resistance mutations and drug binding residues of cytochrome b in Saccharomyces cerevisiae, Plasmodium falciparum, B. microti, B. duncani and six ovine Babesia Figure 1 Comparison of the mitochondrial genomes of six ovine Babesia isolates and other apicomplexan parasites. It was performed using SnapGene software and Adobe Photoshop. The different shades of gray represent different gene types. In detail, white represents TIR, 40% grey represents protein-encoding genes, and 80% grey represents LSU  Phylogenetic relationships of Babesia infective to small ruminants in China and other apicomplexan parasites. Phylogeny was inferred with a maximum likelihood analysis of amino acid sequence of cox1 and cob genes