Molecular Phylogeny of Theileria equi in Horses of Khuzestan Province Based on 18S Ribosomal RNA

1. Introduction
Equine piroplasmosis is a tick-borne disease of equids caused by two protozoan parasites, Babesia caballi and Theileria equi (formerly Babesia equi), which are transmitted primarily by ixodid ticks. Clinically, the disease manifests with fever, inappetence, anemia, icterus, hemoglobinuria, and elevated respiratory and heart rates; in severe cases, mortality may occur. Equine piroplasmosis is endemic to most tropical and subtropical regions, posing a significant threat to equine health and international trade. According to the World Organization for Animal Health (WOAH), formerly known as the Office International des Epizooties (OIE), endemic areas include Africa, Asia, Central and South America, Cuba, Southern Europe, and the Middle East [1].
The taxonomic classification of the protozoan parasites responsible for equine piroplasmosis has been debated since their initial discovery, with the status of T. equi remaining particularly controversial [2]. Currently, both B. caballi and T. equi are placed within the phylum Apicomplexa, which also includes other hemoprotozoa such as Plasmodium and Theileria [3]. While B. caballi is widely accepted as a typical member of the genus Babesia, T. equi exhibits morphological and developmental features that align more closely with members of the genus Theileria [1]. Molecular phylogenetic analyses have revealed that T. equi shares characteristics with both Babesia and Theileria, suggesting that it does not fit neatly into either genus. Several studies have suggested making a new genus to accommodate its unique position in family tree, which lies intermediate between the two lineages [4-6]. Nuclear rRNA genes are widely regarded as suitable molecular targets for species-level identification [7, 8]. However, sequence heterogeneity within the 18S rRNA gene of T. equi has been documented in various countries [9-11], raising questions about its taxonomic stability. Although Iran is located within a known endemic region for equine piroplasmosis, and infections caused by T. equi and B. caballi have been reported in equids from multiple provinces [12-14], no molecular characterization or phylogenetic analysis of these parasites has yet been conducted on Iranian isolates.
The purpose of this study was to explore the molecular characteristics and phylogenetic placement of T. equi isolates collected from horses in Khuzestan Province, Iran.

2. Materials and Methods
2.1. Sample collection 
EDTA-treated blood specimens were collected from five pure-bred Iranian Arab horses located in Ahvaz, Khuzestan Province, Iran. Blood was drawn from the jugular vein of each animal using sterile syringes and immediately transferred into vacutainer tubes containing ethylenediaminetetraacetic acid (EDTA) as an anticoagulant. All specimens were handled under standardized conditions to preserve DNA integrity for downstream molecular analyses in the Molecular Biology Laboratory of the Faculty of Veterinary Medicine at Shahid Chamran University of Ahvaz.

2.2. DNA extraction
Genomic DNA was isolated from 50 mL of EDTA-treated whole blood samples using the Rapid Genomic DNA Isolation Kit (MBST, Tehran, Iran). The concentration and purity of the extracted DNA were evaluated by measuring the A260/A280 absorbance ratio using a NanoDrop microvolume spectrophotometers. Purified DNA was aliquoted and stored at -20 °C until further molecular analyses.

2.3. Polymerase chain reaction (PCR) amplification and sequencing 
To detect T. equi, primers targeting the 18S rRNA gene, Bec-UF2 (5′-TCGAAGACGATCAGATACCGTCG) and Equi-R (5′-TGCCTTAAACTTCCTTGCGAT) were used as described [15]. PCR amplification was carried out in 25 μL reaction volumes, including 50 mM KCl, 10 mM Tris-HCl (pH 8.3), 0.6 μM of each primer, 0.4 mM of each dNTP, 2 mM MgCl₂, 100 ng of genomic DNA, and 2.5 U of Taq DNA polymerase. The Thermal cycling process was set up as follows: initial denaturation at 94 °C for 3 minutes, followed by 30 cycles of denaturation at 94 °C for 45 seconds, annealing at 60 °C for 50 seconds, and extension at 72 °C for 60 seconds, with a final extension at 72 °C for 5 minutes. The DNA fragments run on a 1% agarose gel stained with Safe Stain (Sinaclon, Tehran, Iran) and checked under UV illumination. PCR products showing positive amplification for T. equi were excised from the gel and purified using the GF-1 Gel DNA Recovery Kit (Vivantis, Malaysia). The Purified DNA fragments were then sequenced using dideoxy chain termination method on an Applied Biosystems 373 automated DNA sequencer.

2.4. DNA sequence analysis
The amplified DNA fragments were subjected to sequencing on both strands using a dideoxy termination method and run on an Applied Biosystems 373 DNA sequencer. DNA sequence comparisons were performed using the blastn algorithm [16] through the NCBI GenBank database to identify species based on nucleotide homology. Homologous sequences were retrieved from the database for comparative analysis. Multiple sequence alignment of all retrieved and experimental sequences was carried out using the CLUSTAL_W program [17] hosted by the European Bioinformatics Institute. Phylogenetic analysis and genetic distance estimation were conducted by the neighbor-joining method with 1,000 bootstrap replicates, performed by MEGA11 software [18]. The fractional GC content of the nucleotide sequences was calculated using the Sequence Manipulation Suite. The secondary structures of rRNA were visualized via the RNAfold web server, as described [19].

3. Results
3.1. PCR amplification and sequence analysis
The quantity of nucleic acid extracted from the blood samples tested was between 8 and 20 ng/μL. 
PCR amplification was carried out to isolate an 18S rRNA gene fragment (Figure 1).

Specifically, TeKh1, TeKh3 and TeKh4 were successfully amplified and sequenced in both directions. The final sequences were then determined using overlapping regions prior to being run against the database. 
Multiple alignment showed that these three sequences (TeKh1, TeKh3 and TeKh4) were identical. This fragment was 435 bp and named TeKh. The nucleotide composition of the amplified sequence was as follows: 111 A (25.52%), 84 C (19.31%), 125 G (28.74%) and 115 T (26.44%). Thus, the percentage of purine nucleotides (54.25%) in the fragment was slightly higher than the pyrimidine nucleotides (45.75%) and its GC count was 48%. In order to compare the amplified nucleotide sequence with the sequences in the GenBank database, a comparison was performed from the blastn program with the highly similar sequences (megablast) program. 
Taxonomic analysis of TeKH found 108 matches in the order Piroplasmida, which is part of the phylum Apicomplexa. Of these, 104 hits were related to the family Theileriidae, and 86 were linked to the genus Theileria. Out of these 86, 74 were identified as T. equi. The results showed that the sequence of TeKh is 100% identical to the only available sequence T. equi strain Fars 97 small subunit ribosomal RNA gene (MK615933), which was published in the GenBank database from Fars Province of Iran in March 2019 as a direct submission. 

3.2. Multiple alignment of TeKh
After removing the primers and trimming the sequence TeKh, the 391 bp gene fragment of TeKh was used for further analysis. In order to analyze the amplified sequence of T. equi from Khuzestan in a larger context, multiple alignments were conducted based on the TeKh sequence with the retrieved sequences, including the only Theileria sequence from Iran (strain Fars 97 small subunit ribosomal RNA; MK615933), using Clastal_W software. The amplified nucleotide sequence of TeKH in the aligned regions was 98.71% and 99.74% similar to sequences T. equi (OM475523l; Russia) and T. equi (OL638195; Turkmenistan), respectively. Among all the identical regions, only two differences were observed. A single G to A substitution found at position 234 of T. equi from Turkmenistan (OL638195), which is the only non-conservative substitution within this stretch of nucleotide sequence. The Russian sequence (OM475523l) lacked a T at nucleotide position 361. TeKh was identical to T. equi sequence isolated from the Fars province in Iran. 

3.3. Secondary structure prediction of TeKh
The secondary structure of the 18S rRNA sequence from TeKh was predicted using the RNAfold algorithm (Figure 3).

The analysis yielded a thermodynamic free energy of (-121.98) kcal/mol, indicating that the RNA can adopt a highly stable folding conformation. However, the most stable configuration, which is referred to as the minimum free energy (MFE) structure, was rarely populated, exhibiting a zero percent frequency in the predicted ensemble. This observation suggests that the RNA molecule possesses substantial structural flexibility, with multiple energetically favorable folding alternatives. This was further supported by a high ensemble diversity score of 156.57, reflecting a broad range of possible secondary structures within the thermodynamic ensemble.

3.4. Phylogenetic analysis and genetic distance
The phylogeny tree of TeKh, isolated from Khuzestan, was created using the nucleotide sequence of the 18S rRNA gene. This was done using two methods: the Neighbor-Joining method and maximum likelihood. Both methods produced very similar trees with strong support values. The phylogeny tree showed that all T. equi species in cluster A were grouped together. However, T. equi from Turkmenistan (OL638195) and Russia (OM475523) were in separate groups, with bootstrap values of 72% and 44%, respectively. The sequence from Khuzestan (TeKh) and the only sequence from Iran, which was isolated from Fars Province (MK615933), were both in cluster A and were identical. Cluster B was separated from cluster A with a moderate support (69%). More genetic diversity was observed in cluster B. Group B included T. bicornis (AF499604) from South Africa and T. cervi (KT959227) from China. Sequences from Switzerland, which was identified as Babesia equi (KM046920), and from Kenya, which was identified as Theileria sp (KF597078), were also included. The support values for these sequences were 79% and 42%, respectively. The tree was built using the 18S rRNA gene sequence of Candida odintsovae (AB054570) as the out group. 
The genetic distance between T. equi from Khuzestan (TeKh) and other T. equi strains was calculated using MEGA11 (Table 1).

The genetic diversity of TeKh and other T. equi strains varied from 0% (in the T. equi Fars 97 small subunit ribosomal RNA gene; MK615933) to 0.03% (in the T. equi isolate UDSSR 91 small subunit ribosomal RNA gene; OL638195). The sequence of TeKH isolated from Khuzestan province was identical to the sequence from Fars Province in Iran. The genetic diversity between TeKH and other sequences ranged from 2.4% to 4.2%. As expected, the genetic distances between TeKh and the out group was 11.1%.

4. Discussion
Phylogenetic analysis of Theileria, based on the 18S rRNA gene, revealed a high degree of evolutionary relatedness among members of the T. equi species. In the resulting tree, the genera Theileria and Babesia are clearly delineated into two distinct clusters, reflecting their taxonomic separation. Cluster A comprises T. equi isolates recovered from ungulates, all exhibiting an exceptionally low genetic divergence of approximately 0.03%, underscoring their close evolutionary relationship. This cluster includes samples from Saudi Arabia, Turkey, Brazil, India, Paraguay, Russia, Cuba, China, and the USA, all of which demonstrate minimal sequence variability. The high sequence conservation of the 18S rRNA gene in T. equi across diverse regions reflects its essential role in eukaryotic evolution. As part of a multi-copy gene family maintained through concerted evolution, its stability ensures consistent ribosome function,which is evolutionarily advantageous, ensuring reliable ribosome function and cellular homeostasis [20]. This evolutionary constraint is mirrored in the RNA’s folding behavior, where thermodynamic analysis revealed a highly stable secondary structure (−121.98 kcal/mol) alongside a high ensemble diversity, suggesting both resilience and structural flexibility. Collectively, these features highlight the gene’s functional importance and evolutionary persistence.
Phylogenetic analysis revealed that the T. equi sequence from Khuzestan clustered closely with the only available T. equi entry from Fars Province in the GenBank database, indicating high genetic similarity. The negligible divergence between these Iranian isolates suggests that animal movement between regions may contribute to the shared genetic profile, reflecting limited intraspecific variation within local equine populations. Cluster B was distinctly separated from cluster A, supported by a 99% confidence level using the neighbor-joining method. While cluster A comprised tightly grouped T. equi 18S rRNA sequences with minimal divergence, cluster B encompassed a broader diversity of species, including B. equi and B. bennetti. Notably, a unique Theileria sp. (KF597078) isolated from a waterbuck in a region endemic for T. parva was also placed in cluster B. Although there is no direct evidence implicating waterbuck in the transmission of T. parva, the detection of novel Theileria genotypes in this cattle-like species may have implications for wildlife-associated parasite reservoirs [21, 22].The inclusion of T. cervi (KT959227), responsible for a piroplasmosis outbreak in Sika deer (Cervus nippon) in China, further highlights the complexity within cluster B. That outbreak also marked the first detection of T. annulata and B. motasi, parasites typically affecting cattle and sheep, within Sika deer populations [23]. Additionally, a lone sequence of Cytauxzoon sp. (KT361082), identified in a domestic cat from Kenya, was assigned to cluster B. Although morphologically similar to Theileria, Cytauxzoon is distinguished by its ability to replicate in macrophages, whereas *Theileria* species typically replicate in lymphocytes [24].

​​​​​​​5. Conclusion
In conclusion, the T. equi sequences isolated from Khuzestan and Fars provinces in Iran showed no detectable genetic differences, highlighting strong intraspecific similarity. Although the sample size in this phylogenetic analysis was limited, the consistently high sequence identity among T. equi isolates from various geographic regions suggests that further investigation is warranted before the 18S rRNA gene can be confidently employed as a species-specific molecular marker for diagnosing T. equi in equine hosts.

Acknowledgements
This study was supported by the Vice President for the Research Affairs Office at the Shahid Chamran University of Ahvaz, Ahvaz, Iran.

Compliance with ethical guidelines
There were no ethical considerations to be considered in this research.

Funding
This research did not receive any grant from funding agencies in the public, commercial, or non-profit sectors.

Authors' contributions
Conceptualization, study design, data acquisition, experiments, data interpretation, and writing the original draft: All authors; Statistical analysis, review and editing: Abbas Jolodar. 

Conflict of interest
The authors declared no conflict of interest.

Data availability 
The data that support the findings of this study are available from the corresponding author upon reasonable request.

References

1. Mehlhorn H, Schein E. Redescription of Babesia equi Laveran, 1901 as Theileria equi Mehlhorn, Schein 1998. Parasitol Res. 1998; 84(6):467-75. [DOI:10.1007/s004360050431] [PMID]
2. Kappmeyer LS, Thiagarajan M, Herndon DR, Ramsay JD, Caler E, Djikeng A. et al. Comparative genomic analysis and phylogenetic position of Theileria equi. BMC Genomics. 2012; 13:603. [DOI:10.1186/1471-2164-13-603] [PMID]
3. Wise LN, Kappmeyer LS, Mealey RH, Knowles DP. Review of equine piroplasmosis. J Vet Intern Med. 2013; 27(6):1334-46. [DOI:10.1111/jvim.12168] [PMID]
4. Allsopp MT, Allsopp BA. Molecular sequence evidence for the reclassification of some Babesia species. Ann N Y Acad Sci. 2006; 1081:509-17. [DOI:10.1196/annals.1373.076] [PMID]
5. Schnittger L, Rodriguez AE, Florin-Christensen M, Morrison DA. Babesia: A world emerging. Infect Genet Evol. 2012; 12(8):1788-809. [DOI:10.1016/j.meegid.2012.07.004] [PMID]
6. Uilenberg G. Babesia--a historical overview. Vet Parasitol. 2006; 138(1-2):3-10. [DOI:10.1016/j.vetpar.2006.01.035] [PMID]
7. Chae JS, Allsopp BA, Waghela SD, Park JH, Kakuda T, Sugimoto C, et al. A study of the systematics of Theileria spp. based upon small-subunit ribosomal RNA gene sequences. Parasitol Res. 1999; 85(11):877-83. [DOI:10.1007/s004360050651] [PMID]
8. Katzer F, McKellar S, Kirvar E, Shiels B. Phylogenetic analysis of Theileria and Babesia equi in relation to the establishment of parasite populations within novel host species and the development of diagnostic tests. Mol Biochem Parasitol. 1998; 95(1):33-44. [DOI:10.1016/S0166-6851(98)00085-1] [PMID]
9. Bhoora R, Franssen L, Oosthuizen MC, Guthrie AJ, Zweygarth E, Penzhorn BL, et al. Sequence heterogeneity in the 18S rRNA gene within Theileria equi and Babesia caballi from horses in South Africa. Vet Parasitol. 2009; 159(2):112-20. [DOI:10.1016/j.vetpar.2008.10.004] [PMID]
10. Kouam MK, Kantzoura V, Masuoka PM, Gajadhar AA, Theodoropoulos G. Genetic diversity of equine piroplasms in Greece with a note on speciation within Theileria genotypes (T. equi and T. equi-like). Infect Genet Evol. 2010; 10(7):963-8. [DOI:10.1016/j.meegid.2010.06.008] [PMID]
11. Liu Q, Meli ML, Zhang Y, Meili T, Stirn M, Riond B, et al. Sequence heterogeneity in the 18S rRNA gene in Theileria equi from horses presented in Switzerland. Vet Parasitol. 2016; 221:24-9. [DOI:10.1016/j.vetpar.2016.03.003] [PMID]
12. Abedi V, Razmi G, Seifi H, Naghibi A. Molecular and serological detection of Theileria equi and Babesia caballi infection in horses and ixodid ticks in Iran. Tick Borne Dis. 2014; 5(3):239-44. [DOI:10.1016/j.ttbdis.2013.11.008] [PMID]
13. Habibi G, Esmaeilnia K, Hablolvarid MH, Afshari A, Zamen M, Bozorgi S. Microscopic and Molecular Detection of Theileria (Babesia) Equi Infection in Equids of Kurdistan Province, Iran. Iran J Parasitol. 2016; 11(1):86-90. [PMID]
14. Seifi HA, Mohri M, Sardari K. A mixed infection of babesia equi and babesia caballi in a racing colt: A report from Iran. J Equine Vet Sci. 2000; 20:858-60. [DOI:10.1016/S0737-0806(00)80117-3]
15. Alhassan A, Pumidonming W, Okamura M, Hirata, H, Battsetseg B, Fujisaki K, et al. Development of a single-round and multiplex PCR method for the simultaneous detection of Babesia caballi and Babesia equi in horse blood. Vet Parasitol. 2005; 129(1-2):43-9. [DOI:10.1016/j.vetpar.2004.12.018] [PMID]
16. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990; 215(3):403-10. [DOI:10.1016/S0022-2836(05)80360-2] [PMID]
17. Thompson JD, Higgins DG, Gibson TJ. CLUSTAL W: Improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 1994; 22(22):4673-80. [DOI:10.1093/nar/22.22.4673] [PMID]
18. Tamura K, Stecher G, Kumar S. MEGA 11: Molecular Evolutionary Genetics Analysis Version 11. Mol Biol Evol. 2021; 38(7):3022-7. [DOI:10.1093/molbev/msab120] [PMID]
19. Lorenz R, Bernhart SH, Höner Zu Siederdissen C, Tafer H, Flamm C, et al. ViennaRNA Package 2.0. Algorithms Mol Biol. 2011; 6:26. [DOI:10.1186/1748-7188-6-26] [PMID]
20. Eickbush TH, Eickbush DG. Finely orchestrated movements: evolution of the ribosomal RNA genes. Genetics. 2007; 175(2):477-85. [DOI:10.1534/genetics.107.071399] [PMID]
21. Githaka N, Konnai S, Bishop R, Odongo D, Lekolool I, Kariuki E, et al. Identification and sequence characterization of novel Theileria genotypes from the waterbuck (Kobus defassa) in a Theileria parva-endemic area in Kenya. Vet Parasitol. 2014; 202(3-4):180-93. [DOI:10.1016/j.vetpar.2014.02.056] [PMID]
22. Raoofi A, Fatemi M, Bokaie S, Zeighami A. Comparison of Tolerance to Theileriosis in Different Breed of Cattle by Evaluation of Clinical Signs and Response to Treatment. Iran J VetMed. 2020; 14(3):231-7. [Link]
23. Liu J, Yang J, Guan G, Liu A, Wang B, Luo J, et al. Molecular detection and identification of piroplasms in sika deer (Cervus nippon) from Jilin Province, China. Parasit Vectors. 2016; 9:156. [DOI:10.1186/s13071-016-1435-3] [PMID]
24. Reichard MV, Van Den Bussche RA, Meinkoth JH, Hoover JP, Kocan AA. A new species of Cytauxzoon from Pallas’ cats caught in Mongolia and comments on the systematics and taxonomy of piroplasmids. J Parasitol. 2005; 91(2):420-6. [DOI:10.1645/GE-384R] [PMID]