Targeted Next-Generation Sequencing and Informatics as an Effective Tool to Establish the Composition of Bovine Piroplasm Populations in Endemic Regions

Abstract

Protists of the genera Babesia and Theileria (piroplasms) cause some of the most prevalent and debilitating diseases for bovines worldwide. In this study, we established and used a next-generation sequencing-informatic approach to explore the composition of Babesia and Theileria populations in cattle and water buffalo in a country (Pakistan) endemic for these pathogens. We collected individual blood samples from cattle (n = 212) and water buffalo (n = 154), extracted genomic DNAs, PCR-amplified the V4 hypervariable region of 18S small subunit rRNA gene from piroplasms, sequenced amplicons using Illumina technology, and then analysed data using bioinformatic platforms. The results revealed piroplasms in 68.9% (252/366) samples, with overall occurrence being markedly higher in cattle (85.8%) than in water buffaloes (45.5%). Babesia (B.) occultans and Theileria (T.) lestoquardi-like species were recorded for the first time in Pakistan, and, overall, T. annulata was most commonly detected (65.8%) followed by B. bovis (7.1%), B. bigemina (4.4%), and T. orientalis (0.5%), with the genetic variability within B. bovis being pronounced. The occurrence and composition of piroplasm species varied markedly across different agro-ecological zones. The high detection of T. annulata in asymptomatic animals suggested a relatively high level of endemic stability of tropical theileriosis in the bovine population.

Keywords: 18S ribosomal RNA; Babesia; Theileria; bovines; informatics; protist populations; targeted next-generation sequencing.

Figure 1

Map of Pakistan showing the numbers of cattle and buffaloes in which piroplasm species were detected in Jhelum (1), Jhang (2), Layyah (3), Okara (4), Bahawalpur (5), Sukkur (6), and Thatta (7) districts of Punjab and Sindh provinces. Abbreviation: FATA, Federally Administered Tribal Area; KPK, Khyber Pakhtunkhwa.

Figure 2

Phylogenetic relationships of amplicon sequence variants (ASVs) of Babesia species determined in this study (coloured symbols) with reference sequences from previous studies. The relationships were inferred based on separate analyses of 18S rRNA sequence data (aligned over 386 bp) using the Maximum Likelihood (ML) and Neighbor-Joining (NJ) methods. Babesia microti (AY918951) was used as outgroup. Nodal support values are indicated as bootstrap values for the NJ (first) and ML (second) method. The scale bar indicates distance.

Figure 3

Phylogenetic relationships of amplicon sequence variants (ASVs) of Theileria species determined in this study (coloured symbols) with reference sequences from previous studies. The relationships were inferred based on separate analyses of 18S rRNA sequence data (aligned over 386 bp) using the Maximum Likelihood (ML) and Neighbor-Joining (NJ) methods. Theileria ornithorhynchi (KT937390) was used as an outgroup. Nodal support values are indicated as bootstrap values for the NJ (first) and ML (second) method. The scale bar indicates distance.

Figure 4

An UpSetR plot of Babesia and Theileria species detected in the tested bovine population in Pakistan by using PCR-based next-generation (deep) sequencing (NGS). The orange horizontal coordinate columns represent the number of samples which were test-positive for each piroplasm species. The blue vertical columns show the numbers of bovines that are positive for single or coinfections of piroplasm species. The green lines connecting green dots indicate coinfecting piroplasm species.