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

Background Babesiosis, a tick-borne disease caused by protozoans of the genus Babesia, is widespread in subtropical and tropical countries. Mitochondria are essential organelles that are responsible for energy transduction and metabolism, calcium homeostasis and cell signaling. Mitochondrial genomes could provide new insights to help elucidate and investigate the biological features, genetic evolution and classification of the protozoans. 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; analyzed the genome size, gene content, genome structure and cytochrome b (cytb) 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 subunit 1 (cox1), cytochrome c oxidase subunit 3 (cox3) and cytb, six large subunit rRNA genes (LSU) and two terminal inverted repeats (TIR) on both ends. The cytb gene sequence analysis indicated the binding site of anti-Babesia drugs that targeted the cytochrome bc1 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 TIRs, 25–1563 bp. Phylogenetic analysis indicated that the six ovine Babesia isolates were divided into two clades, Babesia sp. and Babesia motasi. Babesia 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.

Mitochondria are crucial organelles that are 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 the respiratory chain and energy conversion [17]. Complex III comprises three redox active subunits, cytochrome b (encoded by the mitochondrial genome), cytochrome c oxidase subunit 1 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 form by cytochrome b. Inhibitors of the 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 a promising target for the development of antimicrobial compounds. Some studies suggest that the yeast bc1 complex structure could be used as a model for discovering new antimalarial drugs [18,19]. However, mutations in the cytochrome b (cytb) gene have been shown to be the molecular basis of drug resistance of some microorganisms [19][20][21][22]. It has been demonstrated that several anti-Babesia drugs, such as diminazene aceturate and atovaquone, which possibly inhibit mitochondrial respiratory activity and electron transport, are ineffective owing to drug resistance [23][24][25]. Therefore, it is necessary to understand the mechanism of resistance of these anti-Babesia drugs 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. Although the mitochondrial genomes of BmLT and BspXJ were sequenced by using Illumina sequencing technology, they have not been sequenced using the Sanger dideoxy chain-termination method for verification [32]. In this study, we sequenced six mitochondrial genomes of Babesia isolates that can infect small ruminants, and the assembled and annotated sequences were submitted to GenBank after comparison with those of other piroplasms. These data were used to clarify the phylogenetic relationships and classification of the babesiae and to determine novel molecular markers for identification of Babesia species. The present study provides valuable information for understanding mitogenome evolution among apicomplexan parasites, identifying diagnostic markers and screening drug targets.

Parasites and isolation of genomic DNA
The purified merozoites of six ovine Babesia isolates were provided by our laboratory [33]. Genomic DNA extraction, concentration measurement and quality evaluation were performed as previously described by Wang et al. [34]. All genomic DNA samples were kept at − 20 °C until use.

Amplification and sequencing
The PCR primers were designed based on the reported genomic sequences of B. bovis and B. bigemina (Gen-Bank: AB499088 and AB499085) (Additional file 1: Table S1). Mitochondrial genome fragments were amplified and cloned into the pGEM ® -T Easy Vector (Promega, Beijing, China) for subsequent sequencing by Sangon Biotech Company.

Genome assembly and annotation
Mitochondrial genomic fragments were assembled with CLC Genomics Workbench v.7.5.1. The mitogenome annotation was performed using 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 Gen-Bank 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 the GenBank database under the accession numbers MK962313, MK962314 and MN605889-MN605892.

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

Phylogenetic analysis
The concatenated sequences of cytb and cytochrome c oxidase subunit 1 (cox1) amino acid residues of 26 apicomplexan parasites (Table 1) were aligned using Clustal W implemented in MEGA v.6.06 (http://www.megas oftwa re.net/) software. Subsequently, a phylogenetic tree was constructed using MEGA v.6.06 with maximum likelihood (ML) or neighbor-joining (NJ) analysis based on the JTT with the Freqs model. T. gondii and Plasmodium spp. were used as the outgroup. Furthermore, phylogenetic analysis of the whole mitochondrial nucleotide sequence was conducted by the ML method using the Kimura 2-parameter nucleotide substitution model. 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% A + T content (Table 1). Mitochondrial genomes of six ovine Babesia isolates contained three protein-encoding genes, cox1, cytochrome c oxidase subunit 3 (cox3), cytb, six large subunit rRNA genes (LSU) and two terminal inverted repeats (TIRs). The transcriptional direction of cox3, LSU3, LSU6, LSU2, cytb and LSU5 was from 3′ to 5′, whereas the direction of cox1, LSU1 and LSU4 was from 5′ to 3′ (Fig. 1, Additional file 2: Figure S1). The start codons of the cox3 and cytb genes of ovine Babesia were ATA and ATG, respectively. The initiation codon of the cox1 gene of BspXJ/DH was ATA, whereas that of B. motasi was ATG. Most of the protein-coding genes had TAA as a termination codon, followed by TGA ( Table 2).

Comparison of mitochondrial genomes sequenced by Sanger and Illumina technology
The alignment of BspXJ and BmLT mitochondrial genomes sequenced using the Illumina method in a previous study [32] and the Sanger technique in this study was 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 LSUs, and the start and stop codons (Tables 1, 2, Additional file 3: Table S2) [32]. Furthermore, there were base differences at several positions between the two sets of sequences (Table 3).

cytb gene sequence analysis
The results of amino acid sequence alignment indicated that all cytb sequences contain the highly conserved PEWY motif. Leu 275 of S. cerevisiae is a key determinant of the 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 sequences, whereas a Phe is present in the T. gondii, bovine and P. falciparum sequences (Fig. 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 drug (target protein is Qo site of bc1 complex) binding residues of cytb of ovine Babesia included Met 128 , Gly 132 , Glu 259 , Leu 262 , Phe 265 and Tyr 266 . Drug (target protein is Qi site of bc1 complex) binding residues of cytb of ovine Babesia included His 187 , Ser 191 and Asp 214 (Table 4).

Phylogenetic analysis
Phylogenetic trees were constructed with the concatenated amino acid sequences of cytb and cox1 using the ML and NJ methods. The two approaches showed no significant changes in the topology. The piroplasms were divided into seven groups: (i) classical Babesia species that could infect ruminants, canines and equines; (ii) classical Theileria species that could infect bovines; (iii) Theileria equi; (iv) Cytauxzoon felis;   (Fig. 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 cytb and cox1, with the exception of T. orientalis, which was located in the B. conradae clade (Additional file 4: 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 the six Babesia isolates infective to small ruminants are highly similar to those of most Babesia spp. with respect to genome size, high A + T content, genome form, and gene content. However, they are smaller than those of B. microti, B. rodhaini and T. equi and larger than that of T. gondii [28,29,44]. Consistent with the results from other apicomplexan parasites, 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 and B. gibsoni [28] but different from those in B. microti, T. equi, T. orientalis and P. falciparum [28,29].
In the reported mitochondrial genomes of piroplasms, B. microti and B. rodhaini have a dual flip-flop inversion system, ranging from 184 to 1082 bp in length [29]. However, one pair of TIRs was found in the mitochondrial genomes of other Babesia spp. and Theileria spp., ranging from 25 to 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. TIRs are considered to play a crucial role in the replication and stabilization of linear mitochondrial genomes [45]. In the published mitochondrial genomes of apicomplexan parasites, the sequences of the coding genes and LSU are basically the same. Therefore, we speculated that differences in the lengths and sequences of TIRs may be responsible for divergence in the host-specific and in vitro culture characteristics of protozoans.
The difference between the Illumina and Sanger sequencing data is mainly caused by nucleotide substitutions, deletions and insertions, which result from the use of different sequencing techniques assembly and annotation software. The settings of the parameters are different, which is one of the reasons for the difference Fig. 1 Comparison of the mitochondrial genomes of six ovine Babesia isolates and other apicomplexan parasites. The analysis was performed using SnapGene software and Adobe Photoshop. The different shades of gray represent different gene types. In detail, white represents TIR, 40% gray represents protein-encoding genes and 80% gray represents LSU Table 2 Gene contents and initiation and termination codons of genes encoding the piroplasm mitochondrial genomes between the two datasets. The results of the Illumina sequencing method are more prone to errors than those of the Sanger method. 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 in energy metabolism, growth and development of apicomplexan protozoa [46]. Cytochrome bc1 is an integral membrane protein complex that 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 the cytb sequences of six ovine Babesia isolates are completely consistent. In the conserved PEWY region, Phe 278 (yeast number of cytb) is present in most organisms, whereas T. gondii and B. taurus have Tyr and Ala, respectively, at this position. Studies have reported that the L275F mutation in yeast has no effect on enzyme activity, but the IC50 increased ten-fold [22]. Compared with the wild-type yeast, the Y279S mutant had a 40-fold increased IC50 (B1.7 mM) for atovaquone [22], and Y268S in the P. falciparum numbering system resulted in a 3000-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 mutations in cytb are largely responsible for the efficacy of drugs (the target protein is the bc1 complex) in apicomplexan parasites.
Currently, atovaquone and ELQs have been reported for the treatment of human babesiosis and malaria by modifying the drug target through disruption of the cytochrome bc1 complex [19,21,22,41,43]. In 2019, a B. motasi-like parasite was 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 humans in China and evaluate the zoonotic potential of B. motasi and the effect of inhibitors binding to the cytochrome bc1 complex. Our data showed that atovaquone, stigmatellin, myxothiazol, endochin-like quinolone (ELQ), antimycin A and NQNO drugs can be used in the treatment of babesiosis in the future. The molecular mechanism of the resistance of these drugs is the mutation of cytb, which suggests that a combined drug strategy is possible to avoid drug resistance during treatment of babesiosis.
In this study, the taxonomical relationships of Babesia isolates infective to sheep and goats are consistent with the reported phylogenetic analyses based on cytb, cox1, cox3, nuclear small subunit (SSU) and internal transcribed spacer (ITS) [5,6,8,32,48]. The ovine Babesia isolates are divided into two species: Babesia sp. and B. motasi. Babesia 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 position of B. gibsoni, B. duncani and B. orientalis   based on the amino acid sequences of cox1 and cytb [26,27,30].

Conclusions
In conclusion, we reported the mitochondrial genome of six ovine Babesia isolates that infect small ruminants in China. The phylogeny based on the concatenated amino acid sequences indicated that there are two Babesia species (Babesia sp. and B. motasi) infective to sheep and goats in China and that the B. motasi isolates possibly belong to two subspecies (BmHB/NX and BmTZ/LT). The possible efficacy of the inhibitor (target protein is bc1 complex) should be evaluated on these six ovine Babesia isolates in future work. Further studies are needed to analyze the TIR and protein functions, which could provide new insights into the phylogenetic relationships, biology and therapy of Babesia.