Molecular Detection and Phylogenetic Analyses of Babesia spp. and Theileria spp. In Livestock in Bangladesh

1. Introduction

Piroplasms, Babesia spp. and Theileria spp., are tick‒borne protozoan parasites causing piroplasmosis in livestock and wild animals worldwide [1]. The species of the genus Babesia infect red blood cells (RBC), while Theileria spp. infect both RBC and white blood cells (WBC). However, the transmission of piroplasms occurs by different tick species of the family Ixodidae [2]. The economic impact on the livestock industry by piroplasmosis may be due to loss of production, reduced working efficiency, cost of treatment and prevention, and morbidity and mortality [3,4].

Bovine clinical babesiosis, caused by B. bigemina, B. bovis, B. naoakii and B. divergens, leads to intravascular haemolytic anaemia [5,6] in cattle and small ruminants. Among these Babesia spp., B. bovis is the most pathogenic species and cause neurological and respiratory disorders in animals, leading to death. Whereas several other Babesia spp., B. ovata, B. major and B. occultans, are known to be less pathogenic and cause subclinical infections. In contrast, B. ovis is the most pathogenic for sheep and goats in tropical and subtropical areas and is characterized by fever, anaemia, icterus, hemoglobinuria and death [7]. Babesiosis is mainly transmitted by Rhipicephalus (Boophilus) microplus.

On the other hand, clinical theileriosis is caused by T. parva and T. annulata, which are transmitted by Hyalomma spp. Theileria parva causes East Coast Fever (ECF) in eastern, central and southern Africa, whereas T. annulata causes tropical (Mediterranean) theileriosis in North Africa, southern Europe and Asia [8]. Moreover, T. mutans, T. tautoragi and T. velifera are also reported to cause bovine theileriosis in Africa [9,10]. Theileria annulata causes lymphoproliferative disease [11] manifested by fever, inappetence, lymphadenopathy, icterus, tachycardia, dyspnoea and weakness [12]. Mortality due to tropical theileriosis may range from 90% in the newly introduced exotic breed to nearly 5% in indigenous breeds [13]. Therefore, the disease threatens livestock production and improvement in developing countries [11]. On the contrary, T. orientalis, T. sergenti and T. buffali are regarded as non‒lymphoproliferative benign theileriosis [14], infecting a wide variety of hosts globally [15]. Although T. orientalis is known to be benign, disease outbreak has been reported from several countries [16,17,18,19,20].

Bangladesh is an agriculture-based developing country in South Asia. It lies between 23° 41' 39.52" N and 90° 20' 39.67" E. Livestock being an integral part of agriculture, contributes 1.90% to gross domestic products (GDP) and 16.52% to agricultural products [21]. The country has 24.7 million cattle head, 1.5 million buffaloes, 26.7 million goats and 3.7 million sheep. The livestock sector provides 20% full-time and 50% partial employment to the rural people. Along with the fast-growing economy of the country, the demand for products is also increasing. Unfortunately, tick‒borne diseases (TBDs), particularly piroplasmosis, pose a significant threat to animal upgradation programs. The tick species previously reported from Bangladesh are R. microplus, R. sanguineus, Hyalomma anatolicum, Haemaphysalis bispinosa, and Amblyomma testudinarium [22,23,24]. In contrast, the occurrence of TBPs reported so far are bovine babesiosis, theileriosis and anaplasmosis from some selected areas. Most of these reports are mainly based on thin blood smears [25,26,27,28,29], while a few are based on molecular tools from selected areas [30,31]. Therefore, a complete scenario on piroplasms in different livestock species from all over the country is deficient. In this study, blood samples from three animal species in five different locations were analysed to screen piroplasms.

2. Materials and Methods

2.1. Ethics Statement

Approval for animal sampling was obtained from the Department of Livestock Services (DLS), Ministry of Fisheries and Livestock, Government of the People’s Republic of Bangladesh. Verbal consent was attained from the animals' owners through providing them with detailed study objectives. Blood samples were collected by registered veterinarians through proper restraining to avoid any injury to the animals. Moreover, the ethical guidelines for the use of animal samples approved by Obihiro University of Agriculture and Veterinary Medicine, Obihiro, Hokkaido 080‒8555, Japan (Animal experiment approval ID number: 22‒23) were followed during sample collection.

2.2. Study Sites and Sample Collection

From June, 2021 to March, 2022, blood samples were collected from cattle (Bos indicus), gayal (Bos frontalis) and goats (Capra hircus) in five different areas of Bangladesh, namely Jhenaidah, Bogura, Sirajganj, Bandarban districts, and Naikhongchari upazila (sub‒district) (Figure 1). Jhenaidah district was selected because of its border with India in the West, Bogura was chosen as a representative of the northern districts, Sirajganj district was selected for its unique livestock rearing practices called Bathan, acres of fallow grassland where animals are housed and maintained in the dry season (December‒June), Bandarban district and Naikhongchari sub‒districts were selected as representatives of hill districts. Although Naikhongchari is a sub‒district of Bandarban, it was considered separately because it is bordered by Myanmar in the southeast and a coastal district, Cox`s Bazar, in the West. Blood samples were collected randomly from 276 apparently healthy animals, namely cattle (Bos indicus; 174), Gayal (Bos frontalis; 9) and goats (Capra hircus; 93) in Jhenaidah (29), Bogura (14), Sirajganj (107), Bandarban (105) and Naikhongchari (21). Approximately 2 ml of blood from the jugular vein of each animal was collected in EDTA coated vacutainer tube (BD Bioscience, NJ, USA). All three species of sampled animals were categorized into two groups, namely young (<2 years old) and adult (≥2 years old). The age of the animals was confirmed through dentition of the animals and farmer’s record books. The collected blood samples were stored in a cool box in the field and refrigerator (4⁰C) in the laboratory.

2.3. Dried Blood Spot Preparation on FTA^TM Elute Micro Card and DNA Elution

A total of 30 μl blood was withdrawn from the collection tube and dispensed in a concentric circular motion on one of the circles of Whatman FTA^TM elute micro card (GE Healthcare, Buckinghamshire, UK). Four circles of one FTA^TM elute micro card were loaded with blood samples from four different animals. The loaded cards were allowed to dry thoroughly for at least three hours and stored at room temperature until DNA elution. Genomic DNA was eluted from the loaded circles of FTA^TM elute micro card following the manufacturer’s guidelines. In brief, the cards were placed on a cutting mat (2.5″ X 3.0″), and a 3 mm disc was punched out from each of the circles using a UniCore punch kit (QIAGEN, Hilden, Germany) and placed in 1.5 ml sterile microcentrifuge tube. The disc was then rinsed in 500 μl sterile water by vortexing three times for five seconds, and the water was removed using a sterile pipette. The washed disc was then centrifuged at 11000 rpm for 30 seconds and excess water was removed. After adding 50 μl sterile water, the tubes containing washed discs were placed in a heating block at 98⁰C for 30 minutes. Following heat treatment, the samples were vortexed and centrifuged at 11000 rpm for 30 seconds to separate the disc from the eluate. Finally, the discs were removed from the tube, and the eluted DNA was stored at ‒30⁰C.

2.4. Molecular Detection of Piroplasms

Genomic DNA was screened for Babesia spp. (B. bigemina, B. bovis, B. naoakii and B. ovis) and Theileria spp. (T. annulata and T. orientalis) by species‒specific polymerase chain reaction (PCR) assays (Table 1). Notably, B. ovis was only screened from goat blood samples. DNA fragments of the target genes were amplified by PCR in a total volume of 10 μl, containing 0.05 μl One Taq DNA Polymerase (New England Biolabs, Ipswich, UK), 0.2 μl of deoxynucleotide triphosphate mix (dNTPs; 10mM), 0.2 μl of each primer (10 μM), 2.0 μl of 5X One Taq standard buffer, 1 μl template DNA and 6.35 µl of UltraPure™ DNase/RNase-free distilled water (Invitrogen, MA, USA). The thermal cycle conditions for the PCR essays in this study were adopted from the previous studies [7,11,32,33,34] The amplified PCR products were then subjected to electrophoresis in 1.5% agarose gel, stained with ethidium bromide, and visualized under UV-transilluminator.

2.5. Sequencing of the PCR-Positive Samples

For sequencing, 2-7 PCR-positive samples from each pathogen were randomly selected. The PCR amplicons were extracted from the gel and purified by using Nucleo Spin® Gel and PCR Clean-up kit (Macherey Nagel, Düren, Germany). The concentration of the extracted products was measured by NanoDrop 2000 spectrophotometer (ThermoFisher Scientific, Waltham, MA, USA). Sequence reactions were performed using BigDye™ Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, California, USA) and ABI Prism 3100 Genetic Analyzer (Applied Biosystems). The resultant sequences were trimmed and assembled using the CodonCode Aligner version 9 (CodonCode Corporation, MA, USA) to get consensus sequences. GenBank BLASTn analysis was conducted to confirm the identity to the sequences already registered in the GenBank database. Shared percent identities among the sequences of each pathogen were calculated through EMBL-EBI Clustal Omega multiple sequence alignment.

2.6. Phylogenetic Analyses

Phylogenetic analyses for the sequences of B. bigemina (BbigRAP-1a), B. bovis (BboSBP-4), B. naoakii (AMA-1), B. ovis (ssu rRNA), T. annulata (Tams-1) and T. orientalis (MPSP) were conducted by maximum likelihood (ML) method using MEGA XI [35]. The ML method was used to select the best nucleotide substitution model based on Bayesian information criterion for ML analyses, and the phylogeny test was carried out by bootstrap method with 1000 replications.

2.7. GenBank Accession Numbers

Accession numbers for the sequences generated in this study were obtained through depositing in the GenBank database of the National Center for Biotechnology Information (NCBI), using BankIt for the genomic DNA sequences and the ribosomal RNA submission portal (submit.ncbi.nlm.nih.gov/subs/genbank/) for the ribosomal RNA sequences. Assigned accession numbers for the sequences generated in this study are OQ162124‒ OQ162130 (B. bigemina), OQ144958 and OQ144959 (B. bovis), OQ148404 and OQ148405 (B. naoakii), OQ130581 and OQ130582 (B. ovis), OQ162131‒ OQ162136 (T. annulata) and OQ144960‒ OQ144964 (T. orientalis).

2.8. Statistical Analyses

The data comprised of locations (Jhenaidah, Bogura, Sirajganj, Bandarban and Naikhongchari) and animal parameters (age, gender and species) were considered independent variables. Data for low detection rates were excluded from the analyses. Pearson’s chi-square and Fisher’s exact test were used to assess the association of TBP detection rates in different study locations, age, sex and species of the animals in GraphPad Prism 8 (GraphPad Software, CA, USA). A p-value was considered significant when it was <0.05.

3. Results

3.1. Overall Prevalence

Among 276 animals screened for Babesia spp. and Theileria spp. by PCR assays, 203 animals (73.55%) were found positive for at least one of the six piroplasms, while 73 animals (26.45%) were not infected by any of the piroplasms. The piroplasms detected in this study are B. bigemina (49.28%), B. bovis (0.72%), B. naoakii (1.09%), B. ovis (32.26%), T. annulata (6.52%) and T. orientalis (46.01%) (Table 2). A significant difference in the prevalence of B. bigemina and T. orientalis was detected among different study locations (Table 2). A geographic separation was observed in the prevalence of B. bovis (Siranganj and Bandarban districts), B. naoakii (Sirajganj district) and T. annulata (Badarban district and Naikhongchari sub‒district). The prevalence of B. bigemina (54.60%) and T. orientalis (67.24%) was significantly higher in the cattle than that in gayals and goats (Table 3). In respect to the age groups, a significantly higher prevalence of T. orientalis was observed in adult goats (≥2 years old) (Table 4). On the other hand, no significant difference in the prevalence was observed between the sexes of animals (Table 5).

3.2. Co‒Infections with Different Piroplasms

The animals infected with two or more pathogens were considered co‒infections. Among 276 animals, 109 animals (39.49%) were found to be infected with two or more piroplasms (Table S1). The number of piroplasms in co‒infections ranged from double to triple. Infections with double pathogens contributed the majority of co‒infections (104/109; 95.41%), whereas triple infections were observed only in 4.59% animals (5/109). The most recurrent combinations of co‒infections were with B. bigemina and T. orientalis (79/109; 72.48%), followed by B. bigemina and B. ovis (12/109; 11.0%) and T. annulata and T. orientalis (7/109; 6.42%). Among the study locations, the highest prevalence of co‒infections was observed in Bandarban (51/105; 48.57%), followed by Naikhongchari (8/21; 38.10%), Sirajganj (40/107; 37.38%), Bogura (4/14; 28.57%) and Jhenaidah (6/29; 20.69%). No significant difference in the prevalence of co‒infections was observed among the study locations.

3.3. Gene Sequence Analyses

Identity for the sequences of each detected pathogen was compared among themselves and with the reference sequences in the Genbank. The percent nucleotide identities of seven BbigRAP‒1a (B. bigemina) sequences (OQ162124‒ OQ162130) in this study ranged from 99.76‒100%. Additionally, these sequences shared 99.75‒100% identity values with MH790974 (Bangladesh) and MG210824 (Tanzania). The shared nucleotide percent identity value between the two isolates of B. bovis (BboSBP‒4) generated in this study (OQ144958 and OQ144959) was 99.40%. These two isolates of B. bovis shared 95.53‒95.96% identity to ON012668 (Kenya) and ON012677 (Australia). In case of B. naoakii (AMA‒1), two sequences (OQ148404 and OQ148405) generated in this study shared the percent identity of 99.73% between them, while these two sequences shared 99.19‒99.46% identity to the isolate from Sri Lanka (LC385894). Furthermore, the ssu rRNA sequences of B. ovis (OQ130581 and OQ130582) displayed a percent identity of 96.58% with each other. Interestingly, they showed a higher identity of 99.45‒99.64% and 99.41‒99.61% with OP003548 (Philippines) and KF723611 (Tunisia), respectively. All six sequences of T. annulata (Tams‒1) in this study were identical. However, these sequences had 99.56% identity to the sequences from Turkey (AF214914) and Egypt (MZ197898, MN251046 and MN251047). On the contrary, the percent identity among the sequences (OQ144960‒ OQ144964) of T. orientalis ranged from 82.22‒99.61%. Nevertheless, these sequences had an identity of 99.87% with LC438477 (Sri Lanka) and 100% with AB560818 (Viet Nam).

3.4. Phylogenetic Analyses

The phylogenetic trees of B. bigemina (BbigRAP-1a), B. bovis (BboSBP-4), B. naoakii (AMA-1), B. ovis (ssu rRNA), T. annulata (Tams-1), T. orientalis (MPSP) were constructed along with the respective reference sequences from the Genbank database. In the phylogram, all the sequences of B. bigemina (BbigRAP-1a) clustered together in the same clade, and displayed a close relationship to those reported previously from Bangladesh (MH790974), Tanzania (MG210824), South Africa (MK481015), Benin (KX685380) and Burkina Faso (OK323209) (Figure 2). The phylogram constructed from BboSBP-4 gene sequences of B. bovis showed that the isolates from Bangladesh were clustered in one clade along with those from Australia (ON012677) and Kenya (ON012668) (Figure 3). Moreover, B. naoakii sequences (AMA-1) in this study were included in the same clade and had a close phylogenetic relationship with LC385894 (Sri Lanka) and LC486015 (Viet Nam) (Figure 4). Similarly, B. ovis (ssu rRNA) isolates from Bangladesh formed a monophyletic clade with those from Tunisia (KF23611 and KP670199), Philippines (OP003548), Iran (KY581550 and KY581551), Iraq (MN309742 and MN560046), Turkey (KY867435 and MG569902), Bosnia and Herzegovina (OM758310), Spain (AY533146) and Portugal (KJ829366) (Figure 5). However, all T. annulata isolates formed a monophyletic clade which was sister to the clade formed by isolates from Egypt (MN251046, MN251047 and MZ197898) and Turkey (AF214914) (Figure 6). On the contrary, the T orientalis (MPSP) isolates from Bangladesh were distributed in two different clades in the phylogram, representing two different genotypes, type 5 and 7 (Figure 7).

4. Discussion

This study was conducted in wide geographic areas, eight agro‒ecological zones (AEZs) of Bangladesh [36], including two predominant livestock species (B. indicus and C. hircus) and one nearly endangered species (B. frontalis) [37]. Moreover, this is the first study on the detection of tick‒borne pathogens from livestock in Jhenaidah, Bogura and Bandarban. In addition, to the best of our knowledge, this study provides the first molecular data on the piroplasms in goats and gayals from Bangladesh.

In this study, a very high prevalence (73.55%) of piroplasms was observed in livestock from Bangladesh. Among the piroplasms, B. bigemina, B. bovis, B. naoakii, B. ovis, T. annulata and T. orientalis were found to infect the animals screened. Moreover, a high rate of co‒infections with two or more pathogens was also detected in this study. Among the four Babesia spp. screened, the prevalence of B. bigemina was the highest (49.27%), followed by B. ovis (32.26%), B. naoakii (1.09%) and B. bovis (0.72%). Although the prevalence of B. bovis is supported by the previous report [30] (0.5%), the prevalence of B. bigemina and B. naoakii was much lower in the previous studies [30,31]. The lower prevalence in the previous studies might be due to the inclusion of samples from limited areas and the selection of crossbred cattle [31]. Crossbred cattle are mainly managed in an intensive system (except Sirajganj district) and therefore do not get exposure to ticks in the pasture. A relatively high prevalence of B. ovis (32.26%) was observed in goats. Although Black Bengal goats are famous worldwide for their meat and skin, they are always neglected to be monitored for TBPs. In this study, the existence of B. ovis in goats in Bangladesh is going to be confirmed by sequencing for the first time. As B. ovis has been regarded as the most pathogenic for sheep [7], a detailed study on the pathogenicity of this piroplasm in Bangladesh is urgently needed. Among different study locations, the prevalence of B. bigemina was significantly higher in Sirajganj (61.68%) and Bandarban (48.57%). Moreover, B. bovis was only found in Sirajganj and Bandarban, while B. naoakii was only detected in Sirajganj. In addition, the detection rate of B. ovis was also higher (37.09%) in Sirajganj. In relation to the livestock species, the prevalence of B. bigemina was significantly higher in cattle (54.60%). The higher prevalence of B. bigemina in cattle might be due to the rearing system, semi‒intensive in Bandarban and extensive in Sirajganj in the dry season, which exposes them to vector tick, R. (Boophilus) microplus. The phylogenetic analyses suggest that each of Babesia spp. formed a monophyletic clade with other reference sequences in the respective phylogram, indicating that a single genotype for each species of Babesia is circulating in Bangladesh.

Two Theileria spp., T. annulata and T orientalis, were screened from livestock in this study. The prevalence of T. annulata in this study was supported by the previous reports [25,31]. Among different study areas, T. annulata was only detected in the hill tracts (Bandarban and Naikhongchari). Hilly areas are rich in vegetation where the vector ticks, Haemaphysalis spp., are abundant. In addition, semi‒intensive rearing practices of livestock in the hills allow them to be exposed to vectors more frequently than in other parts of the country. On the other hand, T. orientalis was found as the third most abundant (46.01%) species, which is in line with the findings of the previous study [30]. The prevalence of T. orientalis is significantly higher in Bandarban (68.57%) and Sirajganj (29.90%), which might be due to the higher exposure of livestock to vector, H. bispinosa (semi‒intensive in Bandarban and extensive system in dry seasons in Sirajganj). The detection rate of T. orientalis was significantly higher in bovine, cattle (67.24%) and gayals (44.44%), than in caprine (7.53%). The lower prevalence of T. orientalis in goats could be due to their intensive management system. Moreover, a significantly higher prevalence of T. orientalis was observed in goats over two years old. The lower prevalence in young goats might be associated with the acquired immunity from the mother and/or less exposure to tick vectors. In the phylogram, the isolates T. annulata from Bangladesh were identical and formed a monophyletic clade with the isolates from Egypt and Turkey. On the contrary, the isolates of T. orientalis were distributed into two different clades, corresponding to Type 5 and 7, in the phylogram constructed from MPSP gene. Among 11 different genotypes of T. orientalis, Type 2 [20] and Type 7 [17] are reported to be associated with clinical cases. Nevertheless, two genotypes, including Type 7, were detected among the five sequences generated in this study. Therefore, a considerable number of sequences from all the study areas may provide precise data on the circulating genotypes of T. orientalis in Bangladesh.

A very high prevalence of co‒infections (39.49%) was observed among the animals examined, which is supported by previous reports [30]. The high prevalence of co‒infections with different piroplasms is associated with the distribution of tick vectors in the country. Ripicephalus (Boophilus) microplus and H. bispinosa are the most abundant tick species throughout Bangladesh (unpublished). These tick species are the predominant vectors for piroplasms in Bangladesh. Therefore, the prevalence of multiple infections is high in livestock in the country.

In general, almost all of the piroplasms screened in this study were highly prevalent in Sirajganj and Bandarban. The high prevalence in Sirajganj might be due to the Bathan system. On the other hand, Bandarban is a hilly area where the humidity, rainfall and vegetations are higher than in other parts of the country, providing suitable conditions for the breeding of different tick species. In the hills, animals are allowed to graze freely throughout the day, further allowing them to be infested by ticks. A geographic separation was observed in B. bovis (Sirajganj and Bandarban), B. naoakii (only in Sirajganj) and T. annulata (Bandarban and Naikhongchari. Although the exact reason behind the geographic separation of piroplasms is still unknown, proper distribution mapping may aid in preventing and controlling these piroplasms in Bangladesh.

5. Conclusions

This study provides the molecular prevalence of piroplasms in goats, the second most important livestock species, and gayals (besides cattle) from Bangladesh for the first time. Moreover, this is the first report on piroplasms from livestock in the studied areas, except Sirajganj and Naikhongchari. In addition, a previously reported piroplasm, B ovis, is going to be confirmed by sequencing here for the first time from Bangladesh. Therefore, our study focuses on the necessity for further investigation on the piroplasms in Bangladesh. However, the inclusion of more samples from other geographies may provide a better picture of the piroplasms of livestock in the country for formulating future control and prevention strategies.