A total of 530 blood samples were collected from cattle in four Sudanese regions using FTA™ cards (Whatman Biosciences, United Kingdom). The four regions were north (n = 69), central (n = 195), east (n = 158), and west (n = 108). North and central regions were designated endemic region, while east and west were designated new extension regions. Information on the sampling locations, whether from endemic or new extension regions, total number examined and T. annulata positive samples by PCR are provided in Table 1.

Extraction of DNA from cattle blood samples and Atbara vaccine strain was carried out using the PureLink™ Genomic DNA Mini extraction kit (Invitrogen, Germany). For diagnosis of T. annulata, the primer used was SSU rRNA gene 989 5′AGT​TTC​TGA​CCT​ATC​AG3′ and the reverse primer was 1,347 5′TGCACAGACCCCAGAG G 3′ giving an amplicon of 370 bp (Allsopp et al., 1993; Taha et al., 2013).

In this study, the primers used were designed by Weir et al. (2007). For detection in capillary electrophoresis, the forward primer was labeled with standard labeling dyes at the 5′ end (Supplementary Table S1). The PCR amplification was carried out as described in Salih et al. (2018). For negative control, the nuclease free water was used, while DNA extracted from a schizont-infected lymphocyte culture derived from T. annulata Ankara strain was used as a positive control.

The ABI 3730 Genetic Analyzer (Applied Biosystems-USA) was used to analyze the PCR amplicons at the BecA-ILRI Hub, SegoliP sequencing unit, Nairobi, Kenya. For size fractionation, the Gene Scan 500 LIZ internal lane size standard (Applied Biosystems-USA) was employed. The Gene Mapper tool (Applied Biosystems-USA) was used to score the results, which allowed for the resolution of 1 base pair (bp) changes with many products from a single PCR reaction. The predominant allele was determined as the one with the biggest area under the curve, and amplicons with highest peak height were scored. Allelobin software (Idury and Cardon, 1997) was used to re-sized all Gene Mapper data based on consensus sequence repeats of each marker (Table 2). The inaugural form of file, designated multi locus genotype (MLG) consisted of genotypes created from only the predominant allele at each locus (Weir et al., 2007). The allelic profile dataset, on the other hand, contained genotypic profiles derived from all alleles observed at each locus (where minor peaks were greater than 33% the height of the predominant allele present). The MLG file was used to determine population genetic diversity and structure, while the allelic profile file was utilized to calculate the multiplicity of infection (MOI), as well as to rule out linkage disequilibrium as a null hypothesis.

Arlequin v. 3.5 http://cmpg.unibe.ch/software/arlequin 35/ (Excoffier, and Lischer, 2010) was used to calculate the expected heterozygosity, as Theileria is haploid and heterozygosity cannot be observed directly. To investigate the genetic relationships between population, principal component analysis (PCA) was calculated in GenAlEx6.5 (Peakall and Smouse, 2006; 2012). Analysis of molecular variance (AMOVA) was performed using ARLEQUIN to test for hierarchical population structure. Nei’s genetic distance (D) (Nei, 1978) was calculate between each group of samples from different populations and the vaccine strain using the genetic data analysis tool (GDA) (http://lewis.eeb.uconn.edu/lewishome/gda.html).

The standardized index of association (I ^S _A ) between each group of samples was estimated using the LIAN 3.7 program, as well as, the degree of linkage disequilibrium (LD) within and between populations (Haubold and Hudson, 2000). After each population was studied separately, the samples were pooled and processed as a single dataset.

STRUCTURE 2.3.4 (http://pritchardlab.stanford.edu/structure.html) was used to investigate population structure employing Bayesian clustering analysis with sample sites as a basis and the admixture scenario with linked allele polymorphism (Pritchard et al., 2000; Evanno et al., 2005). Initial runs of one million steps were used to investigate the datasets (burn-in of 20%). For every value of K scale from one (considering all are T. annulata) to five (assuming all the five populations are genetically distinct), triplicates were performed. To identify which K produced the greatest representation of the data, STRUCTUREHARVESTER 0.6.1 (Earl and von Holdt, 2012) was employed. CLUMPP 1.1 (Jakobsson and Rosenberg, 2007) and DISTRUCT 1.1.2 (Rosenberg, 2004) were used to parse and format the data in order to assess the STRUCTURE output. CLUMPP 1.1 aligns cluster assignments across duplicate analyses, while DISTRUCT 1.1.2 assists with visual representation.

MOI was considered as “existence of numerous genotypes per isolate” when more than one allele was detected at a locus and the smaller peaks were exceed 33% of the height of the predominant allele expressed (Weir et al., 2011; Salih et al., 2018). The mean number of alleles across all nine loci was determined for every sample, and this number was used to indicate the multiplicity of infection in that sample. The multiplicity of infection for each population was calculated by taking the overall mean for each sample’s index value.

T. annulata DNA was tested in 530 cattle blood samples. The SSU rRNA PCR assay verified 246 (46.4%) samples positive for T. annulata DNA which were subjected to genotyping in addition to the vaccine strain (Table 1). Distribution of the positive samples were as follow, endemic regions n = 156 (North n = 36, Central n = 120) and new extension regions n = 90 (East n = 47, West n = 43) (Figure 1).

In all of the samples, each marker was highly polymorphic. The polymorphic information content (PIC) of marker TS8 had the highest (0.87), whereas TS9 had the lowest (0.36) (Table 2). This finding argued in favor that these markers could be effective in determining linkage disequilibrium analysis in T. annulata populations. The existence of more than one allele at one or more loci confirmed the presence of several genotypes in the samples. For each marker, the number of alleles identified varied from four in TS9 to 22 for TS8 with the mean of 12.44 per marker (Table 2). The dominant allele frequencies varied from 0.25 (TS25) to 0.65 (TS9), with an average of 0.46 (Table 2).

Principal components analysis (PCA) revealed that there is no clustering according to geographical origin (Figure 2). Two sub-structures with a mix of all four populations in both clusters and the vaccine stain being aligned with left-lower cluster were demonstrated, indicating that the parasite populations are rather distinct, with considerable genetic mixing and gene flow between parasites in the four distinct geographical populations investigated.

The allelic profile data set was examined to see if the T. annulata populations observed in Sudan were in linkage equilibrium or disequilibrium. When all the four sub-populations were analysed together (as a single population), the (I ^S _A ) was positive and greater than zero and the pairwise variance (V _D) was more than the 95% critical value (L) suggesting that the merged populations are in linkage equilibrium (LE) (Table 3). The analysis was performed for each population individually to assess for geographic sub-structuring, and three of the populations central, east and west, were shown to be in linkage disequilibrium (LD) (Table 3). Only 11% of the genetic variation was explained by variations between populations, which account for a considerable portion of the genetic diversity (89%) detected within populations (Figure 3).

Estimating Nei’s genetic distance (D) between each of the four regionally sampled populations as well as between them and the vaccine strain, was used to evaluate genetic differentiation between the four populations (Table 4). The genetic differentiation between central and east populations (D = 0.82) was greater than that observed between the east and west populations (D = 0.64). The population with the lowest genetic distance from the vaccine genotype was west (D = 0.39), while the most genetically similar was north (D = 0.62) (Table 4).

Based on the Evanno et al. delta K technique, the STRUCTURE results imply that K = 3 is the optimal number of genetic groups to define the genotypes of Sudanese T. annulata populations as well as in T. annulata vaccine strain (Figure 4). The three clusters are designated as gene pool 1, 2 and 3 respectively. Gene pool 1 (purple colour) prevailed in central and east, while gene pool 2 (blue colour) were most prevalent in north and vaccine, and pool 3 (yellow colour) predominated in west (Figure 4). In the vaccine strain, gene pool 2 appears to be more common than gene pool 1.

Multiple genotypes were observed in T. annulata populations from the four geographic regions, with multiple alleles being found at one or more loci. The mean number of alleles for the nine loci was determined for each sample, to obtain an index value that denoted multiplicity of infection. Table 5 summarizes the multiplicity of infection for each population and across all four populations analyzed. North and central populations had high mean of values of 1.75 and 1.59, respectively, while east and west had values of 1.33 and 1.03, respectively.