Babesia ovis–infected blood sample was obtained from the natural case (50% parasitemia levels) located in an area in the Central Anatolian Region of Turkey in which ovine babesiosis is enzootic in 2015. The parasite was monitored by a microscopic examination of Giemsa-stained thin blood smears and verified by the PCR analysis. B. ovis–infected blood was washed five times with phosphate-buffered saline (PBS) by centrifugation at 2,000 rpm for 15 min to separate RBCs from other cells, and, then, infected RBCs were diluted 10 times with 0.1% saponin (in PBS) to lyse. It was kept at room temperature for 5 min and centrifuged at 4,000 rpm for 20 min to collect merozoites. The merozoite pellet was diluted with PBS, and DNA extraction was performed by using Phenol Chloroform Isoamyl Alcohol (Sigma P2069) (Green and Sambrook, 2012). MinION libraries were constructed with the Ligation Sequencing Kit 1D, LSK108 (Oxford Nanopore Technologies), and then sequenced with three FLO-MIN106 flow cells. A PacBio library was sequenced using the PacBio RS II sequencer with two SMRT cells. An Illumina library was constructed with a TruSeq DNA PCR-Free kit (Illumina), and 300-bp paired-end sequencing was performed using a MiSeq sequencer (Illumina). Reads with low quality were removed using Trimmomatic v0.36 based on the following parameters: LEADING:20 TRAILING:20 SLIDINGWINDOW:4:15 MINLEN:36.

Infected blood from the same sheep was dissolved into a PAXgene blood RNA tube, and RNA was purified using the PAXgene Blood RNA Kit (QIAGEN). The quality and quantity of the purified RNA were validated with Bioanalyzer (Agilent). Library construction for RNAseq was performed as per the instruction manual with the TruSeq Stranded mRNA LT Sample Prep Kit (Illumina), and 300-bp paired-end sequencing was performed using a MiSeq sequencer (Illumina).

Obtained fast5 data were base-called to fastq using poretools 0.6.0 (Loman and Quinlan, 2014) or albacore.1.2.6. The fastq derived from MinION and Pacbio preassembled reads and Illumina paired-end reads were assembled using SPAdes-3.11.1 with -k 127 option (Bankevich et al., 2012). The resulting scaffolds were examined to subtract possible host contamination and artifacts based on the following criteria. First, scaffolds with equal or more than 1,000 bp were selected. Second, redundancy among the scaffolds was excluded if more than 90% of the length of each scaffold had more than 80% identity against the other scaffolds. The identity was evaluated by BLASTN (Altschul et al., 1990). Third, scaffolds were excluded if the region showing the highest identity against B. bovis or B. bigemina was less than 10% of each scaffold. The identity was evaluated by BLASTX (Altschul et al., 1990) using B. bovis (Brayton et al., 2007), B. bigemina (Jackson et al., 2014), and Ovis aries (Archibald et al., 2010) coding sequences and non-redundant protein database provided by National Center for Biotechnology Information (NCBI). Finally, sequencing depth was calculated using BamDeal (https://github.com/BGI-shenzhen/BamDeal) with aligned reads made by bowtie2 version 2.3.4.1 (Langmead and Salzberg, 2012) for the Illumina reads and BWA version 0.7.17 (Li and Durbin, 2010) for the nanopore and PacBio reads, and, then, scaffolds with too low (less than 10) or too high (more than 100) coverage were excluded using the Illumina reads. Scaffolds corresponding to both apicoplast and mitochondria genomes were additionally specified on the basis of sequence similarity. The qualified scaffolds were polished using the Illumina reads with Pilon with 15 iterations until a plateau was reached (Walker et al., 2014). To eliminate redundancy, further quality filtering was performed using Redundans (Pryszcz and Gabaldón, 2016). The qualified contigs were validated using BlobTools v1.1.1 (Laetsch and Blaxter, 2017). The translated amino acids sequences of each contigs were examined using Diamond version 2.1.6 (Buchfink et al., 2021) against NCBI non-redundant protein (nr) database and then subjected to BlobTools as inputs. Integrity of the assembly was validated using BUSCO using the apicomplexa_odb10 dataset (Manni et al., 2021).

The gene model was estimated by AUGUSTUS version 3.1.0 (Stanke et al., 2008). Trained parameters for AUGUSTUS were acquired by webAugustus (Hoff and Stanke, 2013) using the B. bovis gene model obtained from PiroplasmaDB version 5.1 (Ma et al., 2013). The RNAseq data were mapped onto the B. ovis–assembled draft genome with Tophat2 (Trapnell et al., 2009). The first round of AUGUSTUS was performed to obtain an exon–exon junction database. The second round of AUGUSTUS was performed using the database to obtain the final gene model. Our previous report describes the procedure in detail (Yamagishi et al., 2017). Functional annotation of B. ovata genes, including gene ontology (GO) terms, was conducted by Blast2GO (Conesa et al., 2005). tRNA and rRNA genes in the genome were predicted by tRNAscan (Lowe and Eddy, 1997) and RNAmmer 1.2 web servers (Lagesen et al., 2007) with default parameters, respectively. Gene model in the apicoplast and mitochondria genomes was estimated by MiGAP considering their bacterial origin.

For the comparative genomic analysis among apicomplexan parasites, we used genome assembly and annotation released in PiroplasmaDB-37 for B. bovis T2Bo strain, PiroplasmaDB-36 for B. bigemina BOND strain, PiroplasmaDB-40 for B. ovata Miyake strain, PiroplasmaDB-61 for B. caballi USDA-D6B2 strain, PiroplasmaDB-40 for B. microti RI strain, PiroplasmaDB-46 for B. divergens 1802A strain, PiroplasmaDB-55 for Babesia sp. Xinjiang, ToxoDB-40 for the T. gondii ME49 strain, and PlasmoDB-56 for the P. falciparum 3D7 strain.

Orthologous genes conserved among these parasites were specified using OMA version 2.5.0 with the default parameters (Altenhoff et al., 2015). The orthologs coded in the B. ovis genome and all the nine parasites above were selected. Amino acid sequences of the orthologs were aligned using MAFFT with the default parameters (Katoh et al., 2002). The resulting gaps were trimmed, and the remaining sequences were concatenated. A phylogenetic tree by the maximum likelihood method with 100 times of bootstrap replications was constructed using MEGA version 10.0.5 (Kumar et al., 2018)

Synteny between B. ovis and B. bovis genomes were examined using MCScanX (Wang et al., 2012). The identity between the amino acid sequences of the genes was calculated using BLASTP, and gene pairs with e-values less than 1e-100 were selected and then the synteny was visualized by AccuSyn (https://accusyn.usask.ca/). Structural diversity among apicoplast genomes of B. ovis (BaOVIS_AP_20230606), B. bovis (NC_011395.1), and Babesia sp. Xinjiang (KX881914.1) was examined by dotplot using YASS (Noe and Kucherov, 2005). Structural diversity among mitochondrial genomes of B. ovis (BaOVIS_MT_20230606), B. bovis (EU075182.1), and Babesia sp. Xinjiang (KX698108.1) was examined in the same way.

The uniquely gained genes were defined as the B. ovis genes that have no orthologs in either species, as shown below. The first set was closely related to B. bovis and Babesia sp. Xinjiang. The other set was wider range of Babesia spp., B. bovis, B. bigemina, B. ovata, B. caballi, B. divergens, and Babesia sp. Xinjiang. The uniquely lost genes were defined as the genes conserved in Babesia spp. mentioned above but lack in B. ovis. Those in B. bovis were retrieved as representatives. GO enrichment analysis for both gained and lost genes was performed using agriGO (Tian et al., 2017).

The ves and smorf multigene families in the B. ovis genome were searched using hidden Markov model. The training dataset was constructed on the basis of the public genomes and annotations mentioned above. In detail, ves1α, ves1β, and smorf in B. bovis were retrieved from the gff by searching “variant erythrocyte surface antigen-1, alpha,” “variant erythrocyte surface antigen-1, beta,” and “smorf,” respectively. Babesia divergens ves1 genes were retrieved from the gff by searching “variant erythrocyte surface antigen.” Babesia bigemina ves1a, ves1b, ves1ba, and ves2 were retrieved from the website by Sanger institute (https://www.sanger.ac.uk/resources/downloads/protozoa/babesia-bigemina.html). Babesia sp. Xinjiang ves was obtained from GenBank sequence MBFZ00000000.1 (Guan et al., 2016). Hidden Markov models for each gene set were constructed, and, then, ves-related genes in B. ovis were identified using HMMER version 3.3.2 (Mistry et al., 2013). The genes in B. ovis and the abovementioned ves and smorf genes in B. bovis, B. divergens, B. bigemina, and Babesia sp. Xinjiang were aligned using BLASTP (Altschul et al., 1990). Those with more than 200-bit score were selected, and mutual similarity was visualized with Gephi using a Fruchterman–Reingold layout (Bastian et al., 2009).

The mtm multigene family in B. ovis was examined using an algorithm predicting transmembrane domain, TMHMM version 2.0c (Krogh et al., 2001), if there were eight or more predicted transmembrane domains and at least one TM domain average per 100 amino acids. Those in B. caballi, B. divergens, and Babesia sp. Xinjiang identified with the same method. Those in B. bovis and B. bigemina were retrieved from the previous study (Hakimi et al., 2020). Mutual similarity among them was examined and visualized as mentioned above.

To find other multigene families, the genes in B. ovis were aligned to themselves using BLASTP. Those with more than 200-bit score were selected, and mutual similarity was visualized with Gephi using a Fruchterman–Reingold layout (Bastian et al., 2009). Clusters were assigned by Gephi function, and their functions were estimated using BLASTP with nr database.

Using MiSeq, we obtained 1.9 million paired-end reads consisting of 1.1 Gbp. In parallel, a total of 196,077 reads and 295,660 subreads, consisting of 0.6 and 1.8 Gbp, were obtained using MinION FLO-MIN106 and PacBio RS II, respectively.

Using SPAdes assembler, these reads were assembled into 2,920 scaffolds comprising approximately 9.2 Mbp. Following subtraction based on redundancy, possible host contamination, and coverage, 80 scaffolds were qualified. Sequentially, we conducted error correction using the Illumina reads with Pilon. In 15 iterations, the modification reached a plateau resulting in the correction of substitution of 353 nucleotides, 11 insertions, and six deletions, respectively. To eliminate redundancy, further quality filtering was performed using Redundans (Pryszcz and Gabaldón, 2016), and 39 redundant contigs were removed. Coverage, guanine-cytosine (GC) contents, and coding sequences of the qualified contigs were validated using BlobTools (Laetsch and Blaxter, 2017) and confirmed that the coverage and GC contents were consistent except for apicoplast and mitochondria contigs. In addition, the majority of the amino acid sequences encoded by each contig were related to Apicomplexan parasite genes ( Figure S1 ), suggesting that there was no migration of sequences other than B. ovis in the qualified contigs. The sequencing depth was 27.1, 9.3, and 1.9 for the Illumina, Nanopore and PacBio reads, respectively. Finally, the draft genome consisting of 7,811,629 bp was successfully obtained and designated as Selcuk assembly ( Table 1 ). The longest scaffolds, N₅₀, and L₅₀ were 1,533,616 bp, 525,571 bp, and the fifth longest contig out of the 41 scaffolds, respectively. The integrity of the assembly was validated by BUSCO, and percentages of completeness, fragmented, and missing core genes were 96.6%, 0.9%, and 2.5%, respectively.

In this study, we applied AUGUSTUS, which provides a more reliable gene model estimation supported by experimentally obtained transcriptome (Stanke et al., 2008). A total of 3.8-Mbp paired-end reads were obtained by RNAseq. Of these, 10.7% were successfully mapped to the assembled draft genome using Tophat2. Considering that the blood specimen included host sheep cells together with B. ovis–infected RBC, the mapped ratio was reasonable. Two scaffolds had high sequence similarity with significant e-value (“0.0” for both) against the apicoplast and mitochondria genomes of B. bovis, respectively, and, then, MiGAP was utilized for gene modeling, considering that they evolved from bacterial symbionts. Collectively, 3,419 coding regions (CDS) were predicted ( Table 1 ). It corresponded to the number of CDS in B. bovis. tRNAs, 5.8s rRNA, and 18S rRNA were estimated in 64, five, and one locus, respectively. Functional annotation of the CDS was performed with Blast2GO. On the basis of the estimation, 2,772 CDS were functionally annotated, whereas the remaining 647 CDS were then simply assigned as hypothetical proteins.

To create an orthology-based phylogenetic tree, we identified 259 single-copy orthologous genes among B. bovis, B. bigemina, B. caballi, B. microti, B. divergens, B. ovata, Babesia sp. Xinjiang, P. falciparum, T. gondii, and B. ovis using OMA. Their amino acid sequences were aligned; then, gaps were removed and concatenated; and a phylogenetic tree was constructed ( Figure 1 ). This genome-wide analysis revealed that B. ovis was phylogenetically most close to B. bovis.

We then examined synteny between B. ovis and the closest Babesia species, B. bovis ( Figure 2 ). It implied that synteny blocks among them were generally conserved although better genome assembly is needed to conclude this. In addition, the contig 17 of the B. ovis genome diverges from chromosomes 2 and 3 of the B. bovis genome, indicating the presence of structural recombination. The synteny of apicoplast and mitochondrial genome among B. ovis, B. bovis, and Babesia sp. Xinjiang was examined as well. Dotplot analyses were applicable due to high sequence identity among them, and their synteny was well conserved ( Figures S2 , S3 ).

We have identified 281 uniquely gained genes compared to closely related Babesia spp., B. bovis, and Babesia sp. Xinjiang. To evaluate systematic gain in function, we analyzed the 281 genes using gene ontology; however, there was no enrichment compared to the whole ontology in B. ovis. We also identified 19 uniquely gained genes in comparison to wider Babesia spp. However, there was no specific enrichment as well. We have identified 229 and 57 uniquely lost genes against a close and wide range of Babesia spp., but there was not any enrichment as well.

Using hmm model, 43 ves1-like genes were identified in B. ovis ( Table S1 ). Mutual homologies among the B. ovis ves1-like genes, smorf, and ves of B. bovis, B. bigemina, B. caballi, B. divergens, and Babesia sp. Xinjiang were visualized. It demonstrated that the B. ovis ves1-like genes were differentiated from ves1α and ves1β of B. bovis ( Figure 3 ). A part of Babesia sp. Xinjiang ves and B. ovis ves1-like genes had similarity; however, they clustered separately from each other.

Genes that encode proteins with multi-transmembrane domains were examined, and 60 genes with eight or more predicted transmembrane domains were identified ( Figure S4 ). It is known that B. bovis and B. bigemina have unique duplicated genes with the multi-transmembrane domains, mtm and msf, respectively. However, B. ovis did not have such duplicated genes.

Other multigene families in B. ovis were examined using a similarity search among B. ovis genes ( Figure S5 ). The biggest cluster was cluster 1, consisting of 16 genes annotated as DEAD box ATP-dependent RNA helicase, followed by clusters 2 and 3 consisting of 10 and 8 genes, respectively. They were annotated as 26S protease regulatory subunit 8 and T-complex protein 1, respectively.

Regarding the expanded B. ovis ves1-like genes, those number could be under estimated because the redundancy of the assembled contigs was enough removed. On the other hand, it is also less likely that the expansion of the mtm family genes had been overlooked at all. Because the expanded mtm genes of B. bovis and B. bigemina are scattered over the contigs regardless of the size of contig, however, there is no expansion of mtm in the B. ovis genome. It is true that gene expansion and contraction can only be fully concluded by the chromosome level of the genome; however, it is possible to qualitatively assert the presence or absence of expansions.