DNA from laboratory cultured strains of pathogens or verified single-pathogen infections were used as single-species control samples to evaluate the specificity of primer pairs, and this panel included Anaplasma phagocytophilum, Ehrlichia ruminantium, Babesia bigemina, Babesia bovis, Theileria annulata, Theileria mutans, Theileria parva and Theileria taurotragi (for more details see Table 1).

Samples originally collected from the Infectious Diseases of East African Livestock (IDEAL) calf cohort study conducted in Kenya from 2007 to 2009 were used (8). In this cohort study, calves were followed from birth up to 1 year old. These calves did not receive any preventive vaccines or treatments during this period in order to measure true cumulative pathogen exposure. For the purposes of tool validation, we examined samples from a subset of 31 calves from the IDEAL study in three categories: (1) nine animals with clinical episodes that survived, (2) nine animals without clinical episodes that survived and (3) 13 animals not categorised. In addition, a subset of eight of these IDEAL samples, underwent sequencing across different runs to assess reproducibility of results. These particular samples will be henceforth referred to as “repeated samples.” A total of 279 Kenyan cattle blood samples from the above-mentioned animals were tested in this study (Supplementary Table S1). These subsets were selected as they can provide insight into tool validation with respect to calves that are more likely to be infected (‘clinical disease’ cohort), as well as the degree of infection/co-infection occurring in calves that did not undergo clinical episodes (‘no clinical episodes’), and animals that were previously diagnosed to be infected with African trypanosomes, to assess diversity in animals infected with a pathogen known to cause immunosuppression (46) and therefore potentially predispose to multiplicity of infection (‘trypanosomes’).

DNA was extracted from 100 μL of blood originally collected from the calves in the IDEAL study and stored in the ILRI liquid nitrogen biobank, using the Qiagen DNA blood and tissue kit according to the manufacturer’s recommendation, with the modification of incubating blood with lysis buffer for 30 min to allow full lysis of the cattle blood. A control extraction was also carried out using distilled water as a substrate, and this was subsequently used as a negative control to test for any cross-contamination that may have occurred during DNA extraction.

For the amplification of parasites from the Theileria and Babesia genera, the V4 hypervariable fragment (378–424 bp) of the 18S rRNA gene was targeted. Primers were adapted from the primers previously described by Bishop et al. (47); RLB-F (28) and RLB-R2 (48) (Table 2).

For Anaplasma sp. and Ehrlichia sp., primers were designed in the homologous regions of the variable regions of 16S rDNA based on the 23 Anaplasma sp. and Ehrlichia sp. sequences obtained from NCBI database (Supplementary Table S2) and their specificity was confirmed by blasting the primer sequence against the NCBI database. Different combinations of primers, predicted to give an amplicon size of between 100 and 600 bp, were tested using a panel of control samples with known pathogen content, to validate their ability to amplify samples with Anaplasma sp. and Ehrlichia sp., but not other haemopathogens. Based on the validation results, the AE-F4 and AE-R3 primer pair was selected (Table 2).

In order to generate amplicons with unique identifiers ready for sequencing, two rounds of PCR were performed. The first-round targeted amplification of the genus specific sequences of the 16S/18S rDNA locus using primers which included adapter sequences to facilitate the second round of PCR. The second round of amplification was performed to anneal the Nextera II multiplex identifier tags (MID) and sequencing primers (Supplementary Table S3).

The two sets of genus specific PCRs; Theileria/Babesia (referred to hereafter as “ThBa”) and Anaplasma/Ehrlichia (referred to hereafter as “AnEh”) were conducted in separate 20 μL reactions, each containing 4 μL Phusion high-fidelity PCR master mix (New England Biolabs), 0.2 μL HF Phusion Taq Polymerase, 10 mM dNTPs, 0.6 μL DMSO, 10 μM of each forward and reverse primer, and 2 μL of template DNA. Cycling conditions were as follows: Theileria/Babesia: 98°C (30 s), 25 cycles of 98°C (10 s), 60°C (30 s) and 72°C (30 s), with 10 min at 72°C; Anaplasma/Ehrlichia: 98°C (30 s), 25 cycles of 98°C (10 s), 64°C (20 s), 72°C (20 s) with 10 min at 72°C. For initial assessment of amplification, some of the amplified products were subjected to electrophoresis in a 1.5% agarose gel in TAE buffer and visualised under ultraviolet light. PCR products were then purified using AMPure XP magnetic beads (Beckman Coulter) according to manufacturer’s instructions.

The second round PCR was performed by using combinations of 16 forward (N502 to N511, N513 to N522) and 12 reverse barcoded primers (N701 to N715) (Supplementary Table S3) so that each individual sample in a pool was assigned a unique barcode combination. PCRs were performed using the same reaction mix as described for the first round PCR, except for the addition of 10 μM barcoded forward and reverse primers, and 1 μL of first round PCR product as a template. Cycling conditions were as follows: Theileria/Babesia: 98°C (30 s), 10 cycles of 98°C (10 s), 64°C (30 s), 72°C (30 s), and 10 min at 72°C. Anaplasma/Ehrlichia: 98°C (30 s), 10 cycles of 98°C (10 s), 64°C (20 s), 72°C (20 s), and 5 min at 72°C. 10 μL of each sample was pooled to create a master sequencing library pool. Then, 100 μL from each pool were loaded onto a 1.5% agarose gel. The products were then purified using by Qiagen QIAquick gel extraction kit (Qiagen) and subsequently with AMPure XP magnetic beads (1X), following the manufacturer’s instructions (Beckman Coulter). Each pool was then quantified using Qubit (Invitrogen) according to the manufacturer’s instructions. Finally, 70 μL of purified pool, at a concentration of 50 nM for AnEh and 22 nM for ThBa, was submitted for sequencing on an Illumina MiSeq platform using a 500-cycle paired-end reagent kit (MiSeq Reagent Kits v2, 2 × 250 bp paired-end reads) with 10% PhiX Control v3 (Supplementary Figure S2). Each of the pools were sequenced in a separate lane.

The first round PCR product of all samples was assessed for positivity on an agarose gel. Based on the results observed on the gel, we proceeded to create two different libraries for each infection. The first library consisted of samples that showed the expected positive band on the gel, and included selected positive and negative control samples. The second library comprised samples that did not show any visible bands by gel electrophoresis in the first PCR. These samples were pooled separately, also including selected positive and negative samples, and eight repeated samples from the first library to assess the haemobiome tool’s reproducibility.

In order to compare the haemabiome results with data derived from well characterised species-specific assays, specific primers targeting T. parva sporozoite microneme-rhoptry surface antigen (p104) and the A. bovis 16S rRNA gene were amplified by a two-step semi-nested PCR using forward and reverse primers as previously described (49, 50). For T. parva p104, the first round PCR was performed using 20 μL reactions containing 4 μL Phusion high-fidelity PCR master mix (New England Biolabs), 0.2 μL HF Phusion Taq Polymerase, 10 mM dNTPs, 10 μM of each forward and reverse primer, 11.8 μL ddH20 and 1 μL of template. p104 first round PCR cycling conditions were 98°C (30 s), [98°C (10 s), 62°C (30 s), 72°C (20 s min)] for 30 cycles, and 2 min at 72°C. For second round PCR, 1 μL of first PCR product as template with the same conditions, except the annealing temperature which was 66°C.

For the A. bovis 16S gene, PCR was performed using 20 μL reactions containing 10 μL Q5 High-Fidelity Reaction Buffer (New England Biolabs), 10 μM of each forward and reverse primer, 7 μL ddH20 and 1 μL of template. The first round PCR cycling conditions were 98°C (1 min), [98°C (1 min), 56°C (1 min), 72°C (1:30 min)] for 35 cycles and 5 min at 72°C.

For the second round PCR, 1 μL of first PCR product was used as template with PCR conditions of 98°C (30 s), [98°C (10 s), 62°C (30 s), 72°C (30 s)] for 30 cycles and 2 min at 72°C.

Amplicon libraries were submitted to Edinburgh Genomics Core facility to generate 250 bp paired end sequencing reads on Illumina MiSeq v2 platform. The raw sequencing data has been submitted to the European Nucleotide Archive and is available under the project accession number PRJEB79313. Sequencing reads were deconvoluted into data from individual samples based on the barcoding combinations. To ensure data quality, the raw sequences were evaluated using Fastqc and poor quality reads with a Phred score below 28 were removed using Sickle (51). The resulting high quality paired end reads were merged using Flash (52) to produce extended amplicon sequences. These merged sequences were then separated based on the presence of primers used for AnEh and ThBa. Within data from each sample, sequences that were 100% identical were grouped into clusters. However, clusters were excluded from further analysis if they met the following criteria: (1) they represented a single copy sequence read, (2) they were predicted to be formed as a PCR chimera from two other more frequent cluster sequences within the same sample, or (3) if they differed by only one nucleotide from a cluster sequence present in the same sample with higher read counts and a fold change of 3, as these clusters may have arisen from sequencing or PCR errors and could not confidently be defined as distinct. The remaining clusters were compared to the SILVA database v138.1 (18S/16S rRNA) using BLAST (53). Following the data from the control samples, specific criteria were established to define positive infections in each sample. These measures comprised the following rules: (a) establishing read frequency thresholds of 1,000 for AnEh and 500 for ThBa mapped clusters (Figures 1A,B, 2A,B), (b) setting percentage identity thresholds of 99 and 97% for AnEh and ThBa, respectively, and (c) ensuring sequence lengths were within the expected ranges of 160–162 bp for Anaplasma/Ehrlichia and 330–378 bp for Theileria/Babesia, respectively (Supplementary Figure S3).

The obtained partial 16S/18S sequences were aligned with reference sequences from Genbank and existing literature listed in Supplementary Tables S2, S4, using Clustal W with Bioedit software version 7.2.5 (54). Phylogenetic trees were constructed using trimmed partial 16/18S gene sequences from both the IDEAL sample data and aforementioned reference sequences.

To determine the best available evolutionary models for Bayesian inference (BI), we used JModelTest version 0.1.1.40 (55) and selected models based on Akaike Information Criterion (AIC). The Bayesian analysis was performed using MrBayes 3.2.6_168 (56) through the phylogeny web service.¹ The Bayesian analysis processed two simultaneous independent runs with four chains each proceeding for 1,000,000 generations. Trees were sampled every 250 generations. The first 25% of iterations were excluded as burn-in. Resulting trees were combined to create a majority-rule consensus tree, from which posterior probabilities were obtained. A Paracoccus sp. (GenBank accession KP003988) and the Hemalivia stellate (GenBank accession KP881349) sequences were used as outgroups for constructing the AnEh and ThBa trees, respectively, following the study by Chiuya et al. (57).

Species-specific nested PCR results were taken as the “gold standard” and compared with the Miseq results. Using the constructed gold standard, sensitivity and specificity were estimated using the “medcalc” software version 19.2.6 (58). Agreement between the different diagnostic tests (species-specific nested PCR and Miseq) assessing the presence of A. bovis and T. parva was calculated. A kappa measure of agreement test was performed to compare the performance of the two tests; κ-value < 0 indicates no agreement beyond chance. A κ-value between 0.41 and 0.60 indicates a moderate level of agreement while a κ-value between 0.81 and 0.99 indicates almost perfect agreement (59). Other statistical analyses were carried out using R version 3.5.1.

To validate the specificity of the Haemabiome tool, we performed the two PCRs using the primer sets described in Table 2 on the set of negative control samples (water negative controls processed using the DNA extraction protocol). We used these results to set the selected filtering thresholds for each pathogen genus (AnEh <1,000 reads and ThBa <500 reads, respectively). Although there are small numbers of reads present in the negative control samples, these are generally well below thresholds and can be confidently assigned as ‘false positives.’

Moreover, for the AnEh primers several bacterial genera were present in some negative controls, which were previously reported as contaminants in DNA extraction kits, PCR and other laboratory reagents. This included Acinetobacter, Actinomyces, Bradyrhizobium, Corynebacterium, Cutibacterium, Delftia, Moraxella, Ralstonia, Stenotrophomonas, Streptococcus, and Xanthomonas (60, 61). Although detected with very few reads (<100), other bacterial genera such as Gordonia and Nocardioides could also be considered potential contaminants of the 16S rDNA NGS process (Figure 1A) for the AnEh primers. However, some negative control samples from the AnEh PCRs did generate read numbers that approached the threshold of 1,000 reads (Figure 1A). When checked, the sequences in these reads aligned with “uncultured bacterium,” indicating clear false positives. These data provided confidence that the tool was not generating false positives above the defined threshold (Figures 1A,B).

Sequencing replicates were carried out using the species-specific positive control DNA (Table 1) to validate the consistency, accuracy and reliability of deep amplicon sequencing. Read numbers for each replicate set are shown in Figure 2. The results showed that the Haemabiome tool successfully amplified and detected all positive control samples, albeit with varying read numbers between different sequencing runs (Figure 2).

To assess reproducibility, eight field samples were subjected to repeat testing in different sequencing runs for both AnEh and ThBa primer sets (Figure 3). In each distinct run, the same pathogen species were consistently detected within each tested sample, and the number of reads linked to each species showed similar relative patterns between repeats. This consistent detection and comparable pattern of read counts indicated that the Haemabiome tool successfully identified the presence of the same pathogen species in both runs in similar proportions, confirming reliability and accuracy of the approach for each pathogen genus across independent sequencing runs.

Finally, as part of the validation we utilised known species-specific primers for the selected individual pathogens, to confirm the Miseq data generated using the Haemabiome tool and the IDEAL samples. For the amplification of A. bovis a nested PCR targeting the species-specific16S gene was used (50). For T. parva we used a species-specific nested PCR targeting the T. parva p104 gene (49).

A total of 279 samples from the 31 animals were screened. Of these, 179 (64.8%) were positive for A. bovis by species-specific nested PCR, compared to 138 (49.4%) samples testing positive using the Haemabiome tool (Table 3). The calculated κ-value, which measures the agreement between the two tests, for A. bovis detection, indicated a substantial level of agreement (κ = 0.65) (Table 4).

A total of 83 (29.7%) samples were positive by the p104 nested PCR for T. parva, while only 23 (8.2%) samples were positive on the Miseq assay (Table 3). The calculated kappa value of 0.26 indicated only slight or weak agreement between the PCR and Miseq tests for the presence of T. parva infections (Table 4).

Using the species-specific PCR results as the gold standard we estimated the sensitivity and specificity of the Haemabiome tool (Table 4). The sensitivity and specificity of the Haemabiome tool were estimated as 77.9 and 96%, respectively for A. bovis, resulting in an overall accuracy of 82.4%. On the other hand, the Haemabiome tool had an estimated sensitivity of 22.9% and specificity of 98% for T. parva, resulting in an accuracy of 81.7%.

This suggests that the Haemabiome tool may be significantly less sensitive than the nested PCRs, which is to be expected given the high sensitivity of the latter tests.

Figure 4 shows a comparison of the Haemabiome tool with species-specific PCR per calf and at different sampling times to confirm if the same sample is positive or negative for both approaches. This figure clearly demonstrates that samples detected as positive by the gold standard method were almost all positive by Miseq (Figures 4A,B). Additionally, the results obtained from the field samples illustrate the power of the Haemabiome tool in accurately identifying the multiplicity of infections.

Initially all 279 samples were assessed by gel electrophoresis for presence of visible bands after the first round of PCR. Out of these, 108 displayed the expected band on the gel for AnEh, of which 89 were also positive for on the Miseq sequencing results, while the remaining 19 samples were negative. By comparison, 100 of the 171 AnEh samples that did not exhibit the anticipated band on gel electrophoresis, were positive for AnEh in the Miseq data, despite the lack of visible band in gel electrophoresis. Out of the 279 samples tested using the ThBa PCR 156 were positive based on gel electrophoresis, with 146 of those samples confirmed positive on Miseq sequencing. However, among the 123 samples classified as negative by gel electrophoresis, 73 were identified as positive for ThBa infection using the Haemabiome tool (Supplementary Table S5).

Among the detected Anaplasma species, A. bovis was found in 49.5% (138/279) of samples, and A. platys in 42.6% (119/279), making them the most abundant species in this sample set. Uncultured Anaplasma sp. Saso was the third most detected species (25.1%; 70/279) followed by Candidatus Anaplasma boleense, A. phagocytophilum (1.4%; 4/279) and Uncultured Anaplasma sp. and Anaplasma sp. clone ZJ06 (1%, 3/279). In terms of Ehrlichia, E. minasensis (or E. canis) was the most abundant species, detected in 6.4%, (18/279) of samples, followed by E. ruminantium at 2.5% (7/279). Overall, 67.7% (189/279) of samples were positive for any Anaplasma/Ehrlichia infection (Table 5).

The sequences from the 138 Miseq positive samples within the Anaplasma genus resolved into 7 distinct clades (Figure 5). The first Anaplasma clade consisted of all A. phagocytophilum, A. phagocytophilum related and Candidatus Anaplasma boleense sequences. This included three sequences from the present study that shared 100% similarity with a previously published Candidatus Anaplasma boleense isolate (GenBank accession, KU586025) and four sequences with previously published A. phagocytophilum sequences (i.e., GenBank accession, CP006617, U02521).

The second clade is represented by 70 sequences obtained from samples that share 100% nucleotide identity with Uncultured Anaplasma sp. Saso strain (GenBank accession, KY924885), including sequences previously isolated from Kenya by Okal et al. (62). The third clade comprised a number of A. bovis isolates with 138 (i.e., all those positive on the AnEh PCR) sequences from this study, with 100% identity to a published A. bovis from Okal et al. (62) in Kenya (GenBank accession, U03775). A small clade with three identical study sequences were also 100% identical to a sequence of Uncultured Anaplasma sp. clone ZJ06 strain (GenBank accession, JN862825) isolated from a cow in China.

A fifth clade consisted of 119 study sequences, which had previously been confirmed as A. platys or A. platys-like sequences. This group included isolates from the study conducted by Okal et al. (62) (GenBank accession, MW019850 and MW019879), as well as a Candidatus Anaplasma camelii sequence (GenBank accession, MH936009) derived from a Kenyan camel. Obtained sequences shared 100% nucleotide identity and 99.4% identity with other A. platys sequences (GenBank accession, KX447502 and EF139459), respectively.

A sixth clade was composed of Anaplasma species that were not detected in our field samples, including A. marginale, A. centrale, A. ovis and other uncultured Anaplasma species. A seventh clade comprised previously published E. ruminantium sequences, along with 7 sequences from present study. These sequences grouped into the clade with 100% sequence similarity with the counterparts in GenBank (CR767821, CR925677, X61659). Finally, the last clade encompassed the E. canis/E. minasensis group, which included 18 study sequences. These sequences demonstrated 100% sequence identity with counterparts in the database from different hosts (GenBank accession, AY394465, CP025749, EF139458, and EU106856 for E. canis and OQ136683, MT163432 and NR148800 for E. minasensis), and >99% of identity with the other published sequences (Figure 5).

T. mutans was the most abundant species and detected in 68.1% (190/279) of samples, followed by Theileria sp. strain MSD in 50.9% (142/279) of samples. T. velifera was the third most detected species at 35.5% (99/279). Finally, T. parva and T. taurotragi were detected in 8.2% (23/279) and 7.2% (20/279) of samples, respectively. In contrast to the high detection rate of Theileria species, B. bigemina was found in only 1.1% (3/279) and B. bovis in only 0.3% (1/279) of samples. Overall, 78.5% (219/279) of the samples were identified as hosting one or more Theileria of Babesia species using the Miseq sequencing approach (Table 5).

The Theileria and Babesia sequences could be grouped into five clades, and the phylogenetic tree is shown in Figure 6. The first clade contained only an unique sequence from this study, detected in a single sample, which clustered with a B. bovis published sequence (GenBank accession, L19077 and KF928959) with >99.1% nucleotide identity. The second clade included one sequence identified in three samples this study, which had 100% nucleotide similarity with a B. bigemina sequence isolated from zebu cattle in Uganda (GenBank accession, KU206291) as reported by Byaruhanga et al. (63).

A third clade encompassed one sequence which was detected in 190 study samples, which grouped together with 100% nucleotide identity with previously published T. mutans sequences (GenBank accession, MN853559 and MN853556); >99.7% identity was observed with other published sequences, and another sequence which was detected in 142 study samples grouped with Theileria sp. strain MSD (99.7% identity, GenBank accession, AF078816). A fourth clade contained three distinct sequences, detected in 20 study samples. The first sequence was detected in four samples and showed 99.69% nucleotide identity with published T. taurotragi sequences (GenBank accessions MN744953 and L19082). The second sequence was detected in 14 samples and showed 100 and 99.73% nucleotide identity with published T. taurotragi sequences (GenBank accessions MN744953 and L19082, respectively). The third sequence was found in two samples and showed 99.38 and 99.07% nucleotide identity with published T. taurotragi sequences (GenBank accessions L19082 and MT814750, respectively).

Furthermore, a T. parva sequence was detected in 23 study samples, sharing 99.7% nucleotide identity with the published sequences (GenBank accession, MG952921 and HQ895985). The final clade contained a T. velifera sequence, detected in 99 study samples, forming a distinct group with 100% identity with counterparts in the database (Genbank accession, AF097993, JN572702, MN853561 and MN853565). Theileria species such as T. annulata, T. buffeli, T. orientalis and T. sinensis presented in a separate group, although no sequences from this study aligned with these species (Figure 6).