The genomic complement of T. equi contains six molecules, including four chromosomes (Figure 1) of ~11.6 Mbp, an apicoplast genome of 47.8 kbp and a mitochondrial genome of 9 kbp. Chromosomes 1 (3,677,484 bp) and 3 (2,338,319 bp) were completely assembled, while chromosome 2 (2,060,349 bp) contains one assembly gap. Chromosome 4 (3,480,987 bp) is fragmented into six pieces with four physical gaps and one assembly gap. The chromosome assemblies agree with the sizes of the chromosomes seen on pulsed field gel electrophoresis, suggesting that there are no large gaps in the genome sequence.

The larger genome size of T. equi as compared to other piroplasms (Theileria and Babesia) (Table 1) is also reflected in the number of predicted protein coding genes (5330), ~ 25% greater than found for T. parva, T. annulata and B. bovis. While T. equi contains homologs of genes found only in the two Theileria spp. (366), or B. bovis (137), it has far more unique genes (1985), which account for the increased size of the genome. The four species have 2,482 genes in common. T. equi has two rRNA operons, and 46 tRNA genes.

^a Pf = P. falciparum, Tp = T. parva, Ta = T. annulata, Bb = B. bovis, Te = T. equi.

^b L indicates the presence of a linear genome.

^c not including introns.

^d genome size/number of protein coding genes.

^e not reported.

The apicoplast is a plastid-like organelle thought to be derived from a secondary endosymbiotic event with green algae 20 21 . Like plastids, many of the genes for metabolic processes in the apicoplast have migrated to the nuclear genome leaving a remnant genome 22 . In T. equi, the A+T rich (71%) 47.8 kb apicoplast genome is larger than those of other piroplasms, due primarily to expansion of the repertoires of three hypothetical genes (Additional file 1). There are 43 unidirectionally encoded CDSs in the T. equi molecule, which includes 11 ribosomal protein coding sequences. Additionally, each of the 20 tRNA and two rRNA genes are present.

T. equi contains 509 nuclear-encoded proteins potentially targeted to the apicoplast as predicted by PlasmoAP 23 , ApicoAP 24 and/or by homology to genes for pathways predicted to occur in the apicoplast (Additional file 2). Similar to other apicomplexans, T. equi has a complete set of nuclear-encoded enzymes for isoprenoid precursor biosynthesis via the methylerythritol phosphate pathway and these activities are predicted to occur in the apicoplast.

The mitochondrial (mt) genome is 9001 bp in length, longer than in other piroplasms and shows evidence of gene duplication and rearrangement. The T. equi mt sequence was recently reported and the authors suggest that a duplicated CDS (BEWA_044660 and BEWA_044650) is actually cytochrome c oxidase subunit III (cox3) 25 . Although there is no sequence similarity to known cox3 sequences, T. equi could not conduct respiration without cox3, and these genes are the only candidate for this function. Our findings corroborate the linear structure of the mt genome and the long inverted terminal repeat structure 25 .

Figure 1 shows blocks of synteny shared between T. equi and B. bovis, T. parva, or T. annulata. There are relatively few large regions that do not have syntenic matches in the other piroplasm genomes. Analysis of chromosome 1 demonstrated that these regions contain unique proteins (184), without Pfam hits or other functional assignment. The T. parva and T. annulata chromosomes are highly syntenic, and a much more fragmented pattern of synteny is seen between Theileria spp. and B. bovis 29 37 . Large regions of synteny are observed between the Theileria spp. and T. equi, suggesting more recent shared ancestry than for T. equi and B. bovis, where the blocks of synteny are more fragmented. T. equi chromosome 1 appears to be evolutionarily related to Thelieria spp. chromosomes 1 and 2. The ancestral Theileria chromosome 2 has split to provide elements of T. equi chromosomes 1 and 2. Approximately 2/3 of each of T. equi chromosomes 2 and 3 share synteny with Theileria spp. chromosomes 1 and 4, respectively. The remaining 1/3 of each of these chromosomes contains more of the B. bovis specific lineage sequences or unique sequences. With a few exceptions, T. equi chromosome 4 and is syntenic to both T. parva and T. annulata chromosome 3. Notably, the orthologous gene matches extend to the ends of each of the T. equi chromosomes demonstrating the lack of telomerically located species-specific gene families and repeats seen in other hemoparasites 26 28 .

TribeMCL placed 2614 proteins in 334 families ranging in size from 2 to 356 members (Additional file 5). Many of the families, including the two largest, have little significant sequence identity, for any two proteins, in pairwise comparisons; however, this result did not change when we used more stringent parameters for clustering. Therefore, we analyzed the families for common functional annotation and show that most members of a given family share a common functional domain or feature (Additional file 5). Exceptions are family 1 and family 17. A few families have high levels of sequence conservation and readily identifiable functional attributes found through Pfam hits. For example, family 11 contains 42 members with hits to Pfam PF00005, a family of ABC transporters. A survey of T. equi protein families that have at least 20 members revealed that most of these families are comprised of functional attributes similar to those represented in similar scans of B. bovis and Theileria spp. None of the T. equi families appear to contain a family of known immunodominant antigenically variable genes, nor are any of the families telomerically associated. Notable families are presented below.

Although host cell transformation has not been reported for T. equi, it has been reported to invade the lymphocyte and develop to a macroshizont stage 13 . In T. parva and T. annulata, the macroschizont life cycle can induce transformation of the infected cell 46 . To assure that the T. equi isolate used in this study infected equine lymphocytes/PBMCs, equine PBMCs were infected with sporozoites obtained from infected adult male Rhipicephalus microplus ticks. To verify infection of equine PBMC cultures, EMA-1 and EMA-2 monoclonal antibodies were tested by immunofluorescence and bound corresponding T. equi antigens within equine PBMCs (Figure 2). The lack of observed transformation in this close relative of pathogens (T. parva and T. annulata) that do transform their host led us to search the T. equi genome for putative host cell transformation genes to see if these genes were retained in the genome. A set of genes 46 identified as potential manipulators of host cell phenotype in T. annulata and T. parva was used to explore the T. equi genome for orthologous genes, along with the EST expression data. T. parva and T. annulata express prohibitin in macro- and microschizonts, while T. equi expresses the ortholog (BEWA_014000) in erythrocyte-stage parasite ESTs. Prohibitin is reported to be a tumor suppressor that is involved in pathways that lead to immortalization, however, this gene was also identified in B. bovis, which does not undergo transformation 46 . Similarly, T. equi has several homologs of cyclophilins, which have been linked to transformation. This gene is found in B. bovis, T. parva and T. annulata, again, limiting the implication of this gene product as a host cell transforming factor. T. equi does not have a homolog of the Theileria schizont AT hook (TashAT) family of proteins. The AT hooks function as nuclear localization signals and DNA binding domains, and may function in cell transformation. Likewise, TashHN (host nucleus) and SuAT1 do not have T. equi orthologs. TashHN is expressed by Theileria macroschizonts that have lost the ability to form merozoites, while SuAT1 is a schizont protein with characteristics of the TashAT family, that can modulate the phenotype of cultured bovine macrophages by altering the expression profile of cytoskeletal proteins in a manner similar to that seen during infection.

Mehlhorn and Schein renamed Babesia equi as T. equi in 1998 due to the description of a lymphocytic stage for this parasite 8 ; however, while this classification has been well adopted, it is not borne out by phylogenetic studies. The phylogenetic position of T. equi was explored using 150 deduced amino acid sequences (Additional file 7) from eight fully sequenced apicomplexan genomes. Results of the Bayesian analysis (Figure 3) support two possible scenarios: the four piroplasm species, T. parva, T. annulata, B. bovis and T. equi, represent a single genus or alternatively, these species represent three separate genera. Maximum Parsimony (MP) analyses of individual sequences resulted in many different tree topologies; however, the four most commonly recovered topologies (Additional file 8) represent 78% of all trees. Consistent with Bayesian inference, the three most common topologies (Additional file 8, topologies a-c) suggest that T. equi shares a more recent common ancestor with the two Theileria species than with B. bovis. MP analysis of the concatenated datasets unambiguously recovered this same relationship of the four piroplasm species, and completely corroborated the tree estimated using Bayesian inference. Finally, maximum likelihood analyses of the concatenated datasets using a variety of substitution models (JTT, MTREV, WAG, RTREV, CPREV, VT) consistently yielded the same estimates, which completely corroborate the trees generated by Bayesian inference and MP analyses of the concatenated datasets (Additional file 8). Collectively, robust phylogeny estimation, which utilized three separate optimality criteria and various models of protein evolution, consistently placed T. equi as a sister taxon to the Theileria clade, with all three of these piroplasm species subtended by B. bovis.

Previous studies employing 18S rRNA sequences across many piroplasmid taxa also revealed uncertainty in the phylogenetic position of T. equi and several other piroplasms 61 62 . Our robust phylogenetic analyses, which strongly imply the placement of T. equi as sister to the Theileria clade provide stability to piroplasm systematics. Our goal was to determine if a selection of orthologous polypeptides from completed Babesia and Theileria genomes could be used to discriminate the appropriate phylogenetic position of T. equi. However, species sampling was restricted by available complete genome sequences and is too limited to completely resolve this issue. A more appropriate solution may be to place T. equi in a genus distinct from both Theileria and Babesia, as has been previously suggested 9 . This option would likely result in the renaming of some related taxa including Cytauxzoon felis, which groups as the closest relative of T. equi in the 18S rRNA study 61 . Importantly, our analysis revealed a set of 50 informative genes that could be analyzed in a broader sampling of Piroplasmida taxa to gain a greater understanding of the evolutionary relationships of the piroplasms (Boxed sequences in Additional file 7).

Animal experiments were approved by the Institutional Animal Care and Use Committee at University of Idaho, USA, in accordance with institutional guidelines based on the U.S. National Institutes of Health (NIH) Guide for the Care and Use of Laboratory Animals.

DNA from T. equi USDA (Florida) strain was used to construct a small and large insert plasmid library as well as a bacterial artificial chromosome library. A total of 100,835 high quality sequence reads (822 bp average read length) were generated and assembled using the Celera Assembler ( http://sourceforge.net/projects/wgs-assembler/), resulting in 119 scaffolds consisting of 495 contigs. A BAC library was end sequenced to generate an additional 3023 reads which were used to confirm the assembly and for targeted sequencing in the closure phase. Gaps were closed by a combination of primer walking and transposon based or shotgun sequencing of medium insert clones, BAC clones or PCR products. The genome sequence has been deposited in GenBank under accession number ACOU00000000. A cDNA library was constructed using cultured infected erythrocytes (Florida strain) using the Creator SMART cDNA library construction kit (Clontech). cDNA in the size range of 0.3 -2kb was normalized using the Trimmer normalization kit (Evrogen). 7629 sequences (GenBank Accession #s HS032712 - HS040340) were generated and assembled into 2355 contigs using PASA 63 .

Genes encoding tRNA’s were identified using tRNA scan-SE 64 . Gene models were predicted using the ab initio gene finding programs GeneZilla, GlimmerHMM 65 , Phat 66 and Snap 67 that used 395 partial and full length high confidence genes in the training set. The training data was manually constructed and inspected for its alignment against highly conserved protein sequences using the AAT package 68 and PASA 63 to align the ESTs to the genomic sequence with a stringent criteria of 95% identity over 90% length using gmap 69 . Data were combined and consensus gene models were derived using Evidence Modeler EVM 70 . The consensus gene models were manually checked for obvious errors. Such models were corrected using a Java based tool called Neomorphic Annotation Station 71 .

Functional annotation was as for B. bovis 29 . In addition, EC numbers were assigned in an automated fashion using PRIAM 72 and metabolic pathways constructed using SRI’s pathway tools 73 . TransportDB is a relational database that was used to determine transporter complement 74 . TribeMCL 75 was used with default parameters to construct a database of genes that are part of paralogous families of proteins represented in the genome. The method is a sequence similarity matrix-based Markov clustering method.

Sybil ( http://sybil.sourceforge.net) was used to create an all-versus-all BLASTP search using the proteomes of T. equi, B. bovis, T. parva, T.annulata, P. falciparum, P. vivax, C. parvum and T. gondii. These outputs were subjected to Jaccard clustering, placing proteins into distinct clusters for each proteome. Clusters from different proteomes were linked based on best bidirectional BLASTP hits between them to provide Jf-COGs. A minimum block size of five with one gap was allowed in the analyses.

To verify PBMC cultures were infected with T. equi, cytospin preparations were made for immunofluorescence antibody microscopy on day nine after cell culture inoculation. After drying, cytospin preparations were fixed in 1% formalin for 2 min and primarily labeled with T. equi EMA-1 and EMA-2-specific monoclonal antibody 36/133.97. Bound antibody was detected using fluorescein-conjugated goat anti-mouse IgG1 and subsequently visualized by fluorescence microscopy.

The phylogenetic position of T. equi relative to other piroplasm species with published genomes, namely T. parva, T. annulata, and B. bovis was determined based on 150 polypeptide sequences, and rooted with four other apicomplexan taxa: P. vivax, P. falciparum, T. gondii, and C. parvum. Polypeptides were selected from single copy genes identified in all eight species. COGs (clusters of orthologous genes) were defined using an in-house comparative pipeline, which starts with BLASTP among protein sequences within and across species. We used the BLOSUM62 matrix with expected value 10^-5. Jaccard clustering was then performed twice, once to form within-species clusters of paralogous genes and a second time to derive a set of multi-species COGs. In the first case, we used an 80% identity cutoff over a minimum of 70% of the length of the smallest protein and a link score of 0.6. In the second, we set the identity cutoff at 50% over >70% length, and Jaccard coefficient cutoff of 0 for edge pruning. In order to minimize ambiguous amino acid homology assignments that can occur when indels are present, peptides were chosen that have the smallest protein size variation across all species (Additional file 7 for list). Protein sequences within a COG were aligned with ClustalW using default parameters 76 .

Three methods were used to reconstruct the phylogenetic relationships among the sequences: i) maximum parsimony (MP) with exhaustive tree search, and otherwise default parameters as implemented in PAUP v4.0b10 77 , ii) Bayesian inference (BI) implemented in MrBayes v3.1.2 78 , and iii). maximum likelihood (ML) estimation as implemented in RAxML v.7.2.8 79 . Initially, the individual protein datasets were analyzed with maximum parsimony, with branch support assessed with 1000 bootstrap replicates. The 150 datasets were subsequently concatenated and again analyzed using MP with exhaustive tree search and similar bootstrapping procedure. For BI of the concatenated datasets, the amino acid transition matrix was set to a mixture of models with fixed rate matrices (Poisson, Jones, Dayhoff, Mtrev, Mtmam, Wag, Rtrev, Cprev, Vt and Blosum) of equal prior probabilities, and otherwise default parameters. Four runs of MrBayes were conducted, each with 4 chains. Two of those runs ran for 567,000 generations and the other two ran for 2,365,000 generations. Convergence was achieved (Potential Scale Reduction Factor, PSRF=1.00) for all model parameters estimated, including tree length (mean-1.79), the amino acid model (Wag, with posterior probability=1.00), and the tree topology. ML analyses of the concatenated datasets all implemented a gamma model of rate heterogeneity with estimation of the proportion of invariable sites. Six separate amino acid substitution models (JTT, MTREV, WAG, RTREV, CPREV, VT) were utilized. Branch support was assessed with 1000 bootstrap pseudoreplications.

Note: Two additional genome sequences of interest were published while this article was in review. T. orientalis does not transform it’s host cell and has been compared to T. parva and T. annulata to identify putative mediators of leukocyte transformation 80 . B. microti phylogeny was analyzed using 316 single copy genes and found to clade separately to T. parva, T. annulata and B. bovis 81 .