Population structure and transmission of Mycobacterium bovis in Ethiopia.

Bovine tuberculosis (bTB) is endemic in cattle in Ethiopia, a country that hosts the largest national cattle herd in Africa. The intensive dairy sector, most of which is peri-urban, has the highest prevalence of disease. Previous studies in Ethiopia have demonstrated that the main cause is Mycobacterium bovis, which has been investigated using conventional molecular tools including deletion typing, spoligotyping and Mycobacterial interspersed repetitive unit-variable number tandem repeat (MIRU-VNTR). Here we use whole-genome sequencing to examine the population structure of M. bovis in Ethiopia. A total of 134 M. bovis isolates were sequenced including 128 genomes from 85 mainly dairy cattle and six genomes isolated from humans, originating from 12 study sites across Ethiopia. These genomes provided a good representation of the previously described population structure of M. bovis, based on spoligotyping and demonstrated that the population is dominated by the clonal complexes African 2 (Af2) and European 3 (Eu3). A range of within-host diversity was observed amongst the isolates and evidence was found for both short- and long-distance transmission. Detailed analysis of available genomes from the Eu3 clonal complex combined with previously published genomes revealed two distinct introductions of this clonal complex into Ethiopia between 1950 and 1987, likely from Europe. This work is important to help better understand bTB transmission in cattle in Ethiopia and can potentially inform national strategies for bTB control in Ethiopia and beyond.

Files S1 and S2 (available in the online version of this article)

Raw sequencing data for all sequenced isolates have been deposited in the European Nucleotide Archive (ENA) under project PRJEB32192. Accession numbers are detailed in File S1.

The R code used to analyse the data in this manuscript is available on Github: https://github.com/avantonder/Ethiopia2020.

Bovine tuberculosis (bTB) is endemic in cattle in Ethiopia, a country that hosts the largest national cattle herd in Africa. Tracing the source of infections of livestock is an important tool for understanding the epidemiology of bTB and defining control/eradication strategies. This study performed the first WGS-based analysis of isolates from Ethiopia and attempted to include isolates from multiple herds in different parts of the country. Our study highlighted the diversity of in Ethiopia, identified within- and between-herd transmission often over long distances and showed that there had been at least two introductions of clonal complex European 3 (Eu3) into Ethiopia in the early 1970s. This work is important as it helps to better understand bTB transmission in cattle in Ethiopia and will potentially inform national strategies for bTB control in Ethiopia and beyond.

Introduction

is one of several highly related subspecies of the complex (MTBC) which causes tuberculosis (TB) in humans and a range of domesticated and wild animals [1]. While the main human pathogen of this complex, Mycobacterium tuberculosis sensu stricto, causes approximately ten million new TB cases yearly [2], is the main causative agent of TB in cattle, also known as bovine TB (bTB). The actual global prevalence of bTB in cattle is not known but recent estimates suggest that around 7 % of all cattle around the world are likely affected, having a significant impact on their productivity [3]. is also capable of being transmitted to humans and it is believed that at least 1.5 % of all human TB cases are due to zoonotic TB [2].

Bovine TB is present across all continents and many countries in Europe, Australasia and the Americas have been able to eliminate or control bovine TB in their national herds, often through costly test-and-slaughter programmes [4, 5]. In countries without adequate resources to control the disease, bTB is frequently endemic [2]. This is the case for Ethiopia which has the fifth largest national cattle herd in the world with over 60 million animals, of which the vast majority are local zebu breeds reared in extensive farming systems [6]. bTB is endemic at a very low level in these extensive husbandry systems but thrives in the intensive dairy sector, with mainly high milk-yield Holstein-Friesian (H-F) or H-F/zebu cross-bred dairy cattle, where the prevalence is high [7]. This is particularly true in the well-established dairy belt in central Ethiopia where several studies have recorded over 25 % bTB prevalence in animals [8, 9] (Almaw et. al. under review). With the current increase in urbanization in Ethiopia, demand for milk is increasing around urban centres [10] and consequently, the intensive dairy sector is expanding and emerging in new centres across the country. Such expansion is associated with trading of dairy cattle between herds and regions, leading to increased risk of disease transmission in this sector.

To address questions around transmission and control of bTB in the highly affected dairy sector, a project named Ethiopia Control of Bovine TB Strategies (ETHICOBOTS) was funded by the Zoonoses and Emerging Livestock Systems programme (UK). One aim of this multi-disciplinary project was to investigate the population structure of in the Ethiopia dairy sector and to explore epidemiological links. Previous studies have used strains from across Ethiopia to establish a basic picture of the population through conventional methods including spoligotyping and Mycobacterial interspersed repetitive unit-variable number tandem repeat (MIRU-VNTR) [8, 11–19]; this contributed to the definition of the African 2 clonal complex of (Af2) confined to Ethiopia and East Africa [20]. However, these techniques have their limitations: both MIRU-VNTR and spoligotyping have less discriminatory power and identical spoligotype patterns can be found in phylogenetically unrelated strains [21-24]. Whole-genome sequencing (WGS) and SNP-based analyses can now provide more detailed and discriminating phylogenetic analyses [25, 26]. Also, due to a lack of recombination in the MTBC, SNPs exhibit very low degrees of homoplasy [22, 24].

A small number of isolates (n=14) from neighbouring Eritrea have previously been analysed using WGS and established the presence of an undescribed clonal complex since named Eu3 [27, 28]. In this study we performed genomic analyses of 134 Ethiopian isolates collected from dairy cattle, humans and a dromedary from several study sites across the country, giving the first WGS-based picture of the population structure of in Ethiopia. We further analysed a collection of 107 Eu3 genomes, including 40 from Ethiopia, to date the most recent common ancestor (MRCA) of this lineage and to estimate the date of introduction into Ethiopia.

Methods

Isolate selection, culturing, and preparation of genomic DNA

Metadata on the host species of the 134 . isolates selected for this study are summarized in File S1. The vast majority (n=121) originated from 76 dairy cattle (12 H-F, 8 Zebu and 56 H-F/Zebu cross-bred dairy cattle), collected from known dairy farms in Mekele, Gondar and Alage (near Hawassa) as well as from the large dairy belt in central Ethiopia (Addis Ababa, Holeta, Sebeta, Bishoftu, Sendafa and Sululta) while six isolates were collected from zebu cattle and one from a dromedary slaughtered at Gondar, Addis Ababa, Butajira, Negele and Filtu abattoirs, respectively (Fig. 1a, File S1). The latter six, the dromedary sample, plus 29 isolates from dairy cattle from central Ethiopia were collected in previous studies [8, 11, 14] and frozen stocks (archived at −80 °C) were used for re-culturing and subsequent genomic DNA preparation. Ninety-two of the isolates from dairy cattle were collected as part of the ETHICOBOTS project. Eighty-six of these isolates originated from tuberculin reactor animals that were slaughtered and tissue samples of suspected TB lesions were processed and cultured for mycobacteria according to a previously described protocol [29]. The remaining six isolates originated from milk collected from tuberculin reactors. These milk samples were processed based on published methodology [30-32] and cultured as described for the tissue samples. The six isolates of human origin were isolated from pulmonary TB cases (sputum) or from TB lymphadenitis cases (fine-needle aspirates), from study sites overlapping with the cattle samples (File S1) [33, 34] (H. Taye, in preparation). All isolates included in this study were isolated between 2006 and 2018 (File S1). Genomic DNA of each isolate was extracted from heat-inactivated cells using the DNeasy Blood and Tissue Kit (Qiagen). Before submission for WGS, all genomic DNA samples were confirmed as by RD4 deletion typing [11].

Fig. 1.

(a) Map of Ethiopia showing the location of isolation for each sequenced isolate coloured by the host .(b) Maximum-likelihood phylogenetic tree of 134 Ethiopian genomes rooted using AF2122/97.

Whole genome sequencing and sequence analysis

WGS was performed in four rounds using the Illumina HiSeq 2500, MiSeq, NextSeq 500, and HiSeqX platforms to produce paired-end reads of between 50 and 150 base-pairs in length. Raw sequencing reads were deposited at the European Nucleotide Archive (ENA) under project PRJEB32192; all accessions used in this project are listed in File S1. FastQC v0.11.9 [35] was used to generate basic quality control metrics for the raw sequence data. Sequence reads were classified using Kraken v0.10.6 [36] and the abundance of the classification was refined to a single level using Bracken v1.0 [37]. Samples with less than 70 % of all reads assigned to the genus were excluded. In silico spoligotyping was performed using SpoTyping v2.1 [38] and the binary spoligotype representations were queried against the spoligotype database (www.mbovis.org) to extract spoligotypes named in the format of SBXXXX. Clonal complexes were assigned to samples using RD-analyzer v1.0 [39] with samples not identified as belonging to previously described clonal complexes [European 1 (Eu1), European 2 (Eu2), African 1 (Af1), and African 2 (Af2)] designated as ‘Other’ [20, 40–42]. Further assignment to clonal complex was based on the phylogenetic lineages recently identified by Loiseau et al. [43].

Sequence reads were mapped to the AF2122/97 reference genome (NC0002945) using BWA v0.7.17 (Burrow-Wheeler Aligner) (minimum and maximum insert sizes of 50 and 1000, respectively) [44]. SNPs were called using SAMtools v1.2 mpileup and BCFtools v1.2 (minimum base call quality of 50 and minimum root squared mapping quality of 30) as previously described [44, 45]. Samples with reads mapping to less than 90 % of the AF2122/97 reference were excluded. Genomic regions consisting of GC-rich sequences such as PPE proteins and PE-PGRS repeats were masked in the resulting alignment using previously published coordinates [46]. Variable sites were extracted from the masked alignment using snp-sites v2.5.1 [47]. Maximum-likelihood phylogenetic trees were constructed using IQ-tree v1.6.5 [48] (constant sites added to alignment, extended model selection and 1000 bootstraps) and the resulting trees were annotated and rooted in iTOL [49]. Pairwise SNP distances for all genomes were calculated using pairsnp (https://github.com/gtonkinhill/pairsnp).

To provide wider context for the isolates sequenced in this study, a global representative dataset (n=352; File S1) was assembled. Genomes were selected based on the following criteria: they were previously published [27, 43, 50–54] and represented the currently understood diversity. All previously characterized clonal complexes (Af1, Af2, Eu1, Eu2 and Eu3) and uncharacterized lineages (Unknown1 to Unknown8) from Loiseau et al. [43] were included. Due to the large number of available sequences for clonal complexes Eu1 and Eu2, a random selection of 100 genomes was chosen for each. Sequence data was downloaded from the ENA and trimmed using Trimmomatic v0.33 [55]. Sample QC, spoligotype assignment, mapping and phylogenetic tree construction were performed as above. The tree was rooted with a isolate.

Spatial analysis

A map of the geographical locations of isolate collection (latitude and longitude) was constructed using R and the ggmap library [56, 57]. The spatial position of genomes from each clonal complex was plotted and a convex hull (the smallest polygon incorporating a given set of points), calculated using the getConvexHull function from the R library contoureR [58], was drawn around genomes from each clonal complex. Pairwise geographic distances between each isolate in kilometres (km) were calculated using the distHaversine function from the R library geosphere [59]. The associations between genetic and spatial distance for clonal complexes Af2 and Eu3 were tested with a Mantel test (1000 permutations to assess significance) implemented using the R library ade4 [60].

Putative transmission clusters

Two different maximum inter-isolate pairwise SNP thresholds were used to define transmission clusters: 15 SNPs to capture older transmission events and five SNPs to capture recent transmission events. The transmission clusters were calculated using the R library iGRAPH [61] and plotted and annotated using the R library ggraph [62].

Molecular dating of Eu3 clonal complex

beast v1.8.4 [63] was run on a SNP alignment of all Ethiopian (n=40) and published (n=67) Eu3 genomes (File S1), using tip sampling dates for calibration. Three runs of 108 Markov chain Monte Carlo (MCMC) iterations were performed using a HKY substitution model, strict or constant molecular clock and constant or exponential population size and growth (12 separate runs). The performance of each model was assessed through the comparison of posterior marginal likelihood estimates (MLE) and the model with the highest Bayes factor [64] (strict clock/constant population size) was selected (Table S1). These three MCMC runs were combined using LogCombiner v1.8.4 (10 % burnin) and convergence was assessed [posterior effective sample size (ESS) > 200 for each parameter] using Tracer v1.6. A maximum clade creditability tree summarizing the posterior sample of trees in the combined MCMC runs was produced using TreeAnnotator v1.8.4. The resulting tree was annotated using ggtree [65]. The R library TIPDATINGBEAST [66] was used to confirm a temporal signal in the dataset; briefly, the dates for each sample were resampled to generate 20 new datasets with randomly assigned dates. beast was then run on these new datasets using the same strict constant priors.

Results

Study sites and population

A total of 134 . isolates were successfully sequenced and passed QC; these included 127 isolates collected from 85 cattle, six isolates from six human TB patients and a single isolate from a dromedary. The origin of isolation and WGS-based genotypes of these isolates are listed in Table 1 (the complete set of metadata can be found in File S1). The majority of the isolates were collected from the large dairy belt in central Ethiopia, in the surroundings of Addis Ababa (Fig. 1a). In silico spoligotype assignment revealed a total of 22 different spoligotypes with the most prevalent types being SB0134 (n=36), SB1176 (n=36) and SB0133 (n=28). Three different clonal complexes were observed in the dataset (Fig. 1b): African 2 (Af2; n=92), European 3 (Eu3; n=40) and Unknown8 (n=2).

Table 1.

Metadata of 134 sequenced isolates included in this study, assigned to the African 2, the European 3, and the Unknown8 clonal complexes

Clonal complex	Samples (n)	Spoligotype (n)	Study site (n)	Host (n)
African 2	92	SB1176 (36), SB0133 (28), SB0912 (8), SB0303 (4), Unknown (3), SB1477 (2), Singletons (11)	Sululta (26), Mekele (15), Bishoftu (11), Addis Ababa (11), Holeta (10), Sebeta (6), Alage (6), Negele (3), Others (4)	Bovine (87), human (4), dromedary (1)
European 3	40	SB0134 (35), Singletons (5)	Mekele (14), Addis Ababa (6), Sebeta (5), Sululta (4), Gondar (4), Holeta (3), Butajira (2), Sendafa (2)	Bovine (40)
Unknown8	2	SB1476	Mekele (1), Bishoftu (1)	Human (2)

Comparison of sequenced genotypes with previously published genotypes

Spoligotype information for previously published isolates from animals, mainly cattle, was collated and the likely clonal complex inferred [67] (Table 2). The most frequent Af2 and Eu3 spoligotypes (SB1176, SB0133, SB0912, SB1477 and SB0134) previously identified in Ethiopia were also represented amongst the spoligotypes of the sequenced isolates in this study at similar frequencies. SB1476 (Unknown8), previously found exclusively in Ethiopian cattle, was found in two of the six isolates collected from humans.

Table 2.

Genotyping data and inferred clonal complex of previously published isolates from Ethiopia

Clonal complex	Samples (n)	Spoligotype (n)	Study site (n)	Host (n)	Reference (n)
African 2	190	SB1176 (104), SB0133 (50), SB1477 (7), SB0912 (6), SB1490 (4), SB0933 (3), SB1491 (2), SB1942 (2), SB1983 (2), Singletons (10)	Addis Ababa (64), Holeta (41), Kombolcha (26), Negele (24), Jinka (12), Adama (7), Woldiya (4), Gondar (3), Fiche (3), Ghimbi (2), Butajira (2), Melga-Wondo (1), Hawassa (1)	Bovine (188), camel (2)	Berg 2009 [11](48), Ameni 2007 [13](41), Biffa 2010 [15](27), Gumi 2012 [31](24), Firdessa 2012 [8](19), Ameni 2010 [12](13), Tsegaye 2010 [18](7), Mekibeb 2013 [17](6), Ameni 2013 [14](3), Mamo 2011 [19](2)
European 3	24	SB0134 (19), SB1517 (3), Singletons (2)	Addis Ababa (21), Yabello (2), Gondar (1)	Bovine (24)	Firdessa 2012 [8](12), Biffa 2010 [15](6), Berg 2009 [11](5), Tsegaye 2010 [18](1)
Unknown8	18	SB1476	Ghimbi (8), Gondar (7), Butajira (2), Jinka (1)	Bovine (18)	Berg 2009 [11](18)

Phylogenetic relationships of genomes

The phylogenetic tree of 134 Ethiopian genomes, rooted with the Eu1 reference AF2122/97, shows a distinct phylogenetic structure with the three included clonal complexes clearly segregating in the phylogeny (Fig. 1b). The Unknown8 clonal complex is basal followed by a split leading to Eu3 and Af2. Three distinct sub-lineages were observed for Af2 and each was associated with a predominant spoligotype (SB0303, SB0133 and SB1176).

Spatial distribution of clonal complexes

Fig. 2a shows a spatial distribution of the 134 isolates overlaid with the smallest polygon incorporating the isolates from each clonal complex. There appears to be a distinct geographical distribution with Eu3 isolates more likely to be distributed to the north and west of the country whilst Af2 isolates have a more southern and eastern distribution. There is a significant overlap between these two clonal complexes in central Ethiopia around Addis Ababa, likely reflecting a greater density of infected cattle from across the country.

Fig. 2.

(a) Spatial analysis of distribution of isolates coloured by clonal complex. Each polygon represents the minimum convex polygon of the sampled locations of the isolates from each clonal complex. (b) Scatterplot of SNP distance against geographic distance for all pairs of genomes. (c) Histogram of all pairwise SNP distances. (d) Boxplot of all pairwise SNP distances separated by clonal complex. (e) Boxplot of all pairwise SNP distances separated by within- and between-animal. (f) Histogram of within-animal pairwise SNP distance. The insert shows the 0–75 SNP range zoomed in.

A comparison of pairwise genetic distance (SNPs) and geographic distance (kilometres) for all pairs of genomes is shown in Fig. 2b. The results of the Mantel tests for Eu3 and Af2 (abattoir isolates removed) showed that there was no association between genetic and spatial distance (Eu3 observation=0.10; simulated P-value: 0.001; Af2 observation=0.03; simulated P-value: 0.001). Instead, the pattern observed reflects the large geographic distribution of these genetic lineages in Ethiopia.

Genetic diversity and putative transmission

The distribution of pairwise SNP distances within the dataset was multimodal reflecting the population structure of the dataset with distinct modes observed for the three clonal complexes as well as the distinct sub-structure of Af2 (Fig. 2c). Categorizing the pairwise SNP distances by clonal complex showed that there is considerably more diversity within Af2 compared to Eu3 and Unknown8 with the maximum pairwise SNP distance within Af2 being 418 SNPs (median 358 SNPs), compared to 82 SNPs (median 36 SNPs) and 11 SNPs (median 11 SNPs) for Eu3 and Unknown8, respectively (Fig. 2d). The within-host pairwise SNP diversity had a median of ten SNPs (range 1–647 SNPs) whilst the between-host pairwise SNP diversity had a median of 380 SNPs (range 2–733 SNPs; Fig. 2e). Further examination of the within-host diversity revealed three distinct peaks within the distribution of pairwise SNP distances with the majority of within-host isolates being within 1 and 30 SNPs of each other (Fig. 2f). The maximum within-host pairwise SNP distance of 647 SNPs was observed in an animal (ser-002) with three isolates, two from Af2 and one from Eu3.

A total of 17 putative transmission clusters were defined using the conservative pairwise SNP threshold of 15 SNPs. The clusters varied in size between two and 19 isolates, and there were 35 isolates that were not assigned to a cluster (singletons; n=35; Fig. 3b). Extensive within- and between-herd transmission was observed with both short- and long-distance transmission occurring. Examples of isolates from the same animal being in different transmission clusters were observed; for instance, isolates from animal ALB720 were found in two different transmission clusters along with isolates from animal 4589 (Fig. 3b). These animals were from herds N (Alage) and Y (Mekele), which are approximately 662 km apart. The largest transmission cluster, containing 19 isolates from 16 animals, showed clear epidemiological links between the Sululta abattoir (SA) and herds from around Sebeta (A, O, P, Q) and Bishoftu (S, T), whilst the second largest transmission cluster, containing 14 isolates from 11 animals showed links to abattoirs in and around Addis Ababa (AA, BA, SA) and herds from as far away as Gondar (U) and Mekele (X, Y).

Fig. 3.

(a) Map of Ethiopia showing the locations of the herds, abattoirs and hospitals from which the isolates were sourced (letter coding from File S1). The region around Addis Ababa was magnified in the insert. (b) Putative transmission clusters defined using a pairwise SNP threshold of 15 SNPs. (c) Putative transmission clusters defined using a pairwise SNP threshold of five SNPs. Nodes are coloured by animal and labelled with the herd, abattoir or hospital of isolation. Edges coloured in blue represent within-host links whilst edges coloured in red represent between animal links. For simplicity, clusters where n <2 are not shown.

When using the more stringent threshold of five pairwise SNPs, ten putative transmission clusters comprising between two and six isolates were defined; 104 isolates were not assigned to a cluster (Fig. 3c). Transmission clusters with both within- and between-herd transmission were defined. Two of the transmission clusters showed evidence of long-distance transmission (>230 km apart) between herds in Gondar (U, V) and Mekele (Y, Z).

Global and European 3 specific phylogenetics

Fig. 4a shows a maximum-likelihood phylogeny of 485 . genomes from 22 countries (including the genomes from this study), collected from 20 different host species between 1983 and 2018. All previously described clonal complexes along with one undescribed clonal complex (labelled as ‘Other’) are represented in the collection. The structure of the phylogeny is consistent with recently published work [43] and shows that pyrazinamide (PZA)-susceptible from Malawi are basal to the other clonal complexes, which are all PZA-resistant. The first of the next two ancestral splits in the phylogeny leads to the sister clades Af2 and Eu3, where the majority (n=132) of the Ethiopian genomes are contained, and the clade containing all other clonal complexes (Other, Af1, Eu1, Eu2 and Unknown3-8). The remaining two Ethiopian genomes (SB1476/Unknown8) are the outgroup to the globally distributed Eu1 clonal complex (Fig. 4a).

Fig. 4.

(a) Global maximum-likelihood phylogeny of 485 genomes rooted using . (b) Time‐calibrated maximum clade credibility tree of European 3 genomes.

The 40 Ethiopian Eu3 genomes were combined with 67 other Eu3 genomes from ten countries and the resulting phylogeny was dated using beast (Fig. 4b). The median date of the MRCA for this clade was estimated to be 1674 [95 % CI: 1557–1769]. The dated phylogeny shows two potential Eu3 introduction events into Ethiopia (Eritrea was part of Ethiopia until 1993) around the same time period: 1975 [95 % CI: 1959–1987] and 1974 [95 % CI: 1959–1984] (Fig. 4b). The Eu3 phylogeny also contains nine BCG vaccine isolates collected in Malawi and the UK and these have a MRCA dating back to 1950 [95 % CI: 1926–1968]. The results of the dated tip randomization analysis are shown in Fig. S1. The substitution rate in the observed data did not overlap with the estimated substitution rates in the randomized datasets showing that the temporal signal observed was not obtained by chance.

Discussion

This is the first study to use WGS to examine the population structure of in Ethiopia. A total of 134 isolates from 12 study sites across Ethiopia were successfully sequenced and three distinct clonal complexes, which appeared to be mostly geographically segregated, and 22 different spoligotypes were observed in the dataset. The in-depth genome analysis showed there was no association between geographic and genetic distance and differing levels of genetic diversity were observed amongst the three clonal complexes. A range of within-host diversity was also observed with isolates from the same animal often being in different transmission clusters, as well as evidence of both short- and long-distance transmission with isolates in transmission clusters collected hundreds of kilometres apart. Molecular dating of a collection of Eu3 genomes which included the Ethiopian genomes in this study estimated the date of introduction of different Eu3 sub-lineages into Ethiopia and Eritrea during the same time period around the early 1970s (CI: 1958–1987).

To get a better understanding of how well the genomes analysed in this study represented the known population structure of in Ethiopia, we collated spoligotype information for the vast majority of Ethiopian isolates from animals (mainly cattle) published to date [67]. When comparing the spoligotypes in the current study (Table 1) with those previously published (Table 2), it was clear that the most frequent spoligotypes of Af2 and Eu3 identified in Ethiopia (SB1176, SB0133, SB0912, SB1477 and SB0134) were also represented in our dataset with their prevalence observed at a similar frequency. Spoligotype SB1476 has also been frequently observed in cattle in previous studies [11] but we provide the first evidence of this spoligotype in humans, possibly as a result of zoonotic transmission. On this basis, the isolates sequenced in this study appear to be representative of the population structure of in Ethiopia.

Examination of the geographical distribution of the clonal complexes, in particular Af2 and Eu3, contained within the dataset showed a potential geographic divide between their distributions with an overlap in the herds around Addis Ababa. The number of isolates from outside this region is small but given the expected movement, often long-distance, of animals between different herds, the geographic separation of clonal complexes is still distinct. The reason for this pattern is unclear but may be due to historical isolation of different parts of the country due to geographic distance; it is interesting that Af2, the indigenous lineage in Ethiopia, was not found in the north west of the country but this may simply be due to a lack of sampling. We could find no evidence of a relationship between genetic and geographic distance for the different clonal complexes implying that the population in Ethiopia is well-mixed and maintained by both short- and long-distance cattle movements. There was a large difference observed in genetic diversity between clonal complexes Af2 and Eu3; this is likely due to the long-term endemic nature of Af2 in Ethiopia which has allowed for significant genetic divergence over time with the emergence of clear sub-lineages within this clonal complex observed. Conversely, the comparatively recent introduction of Eu3 into Ethiopia (see below), has not allowed enough time for more genetic diversity to emerge. As expected, considerable genetic diversity was observed between isolates from different animals with some isolates, from different clonal complexes, being as much as 733 SNPs apart. Less expected was the range of within-host diversity observed. The majority of isolates from the same animal were up to 30 SNPs (median=10 SNPs) apart from each other; given the previously estimated mutation rate of of 0.15–0.53 SNPs per genome per year [68, 69] it is likely infections are being maintained over a long period of time. Less closely related isolates from the same animal were also found showing that multiple infections by different strains was taking place amongst our samples.

With bTB being endemic in Ethiopia and the intensive dairy sector highly affected, there is considerable interest from stakeholders to explore control strategies for this disease. One of the aims of the ETHICOBOTS project was to use WGS to identify potential transmission of bTB in Ethiopian dairy cattle. There is currently no clear consensus as to what the most appropriate SNP threshold use for defining transmission clusters in either bTB or M. tuberculosis sensu stricto with thresholds of 5, 12 and 15 SNPs having been previously applied [70-72]. Thus, we defined transmission clusters both using a conservative threshold of 15 SNPs (Fig. 3b), which would allow for both the identification of older transmission events and account for varying rates of mutation within the animals sampled, and a stricter threshold of five SNPs (Fig. 3c), which would identify more recent transmission. Using either threshold, we could find clear evidence of transmission occurring between animals from the same herd as well as animals from different herds and from animals sampled at local abattoirs. Using the more conservative threshold of 15 SNPs, we found evidence of isolates from the same animal in different transmission clusters, showing that animals are likely being re-infected with different strains. There was also strong evidence of both short- and long-distance transmission with isolates in some transmission clusters being hundreds of kilometres apart. Despite the sparse sampling, the fact that we were able to identify various types of transmission events even with a strict threshold of five SNPs suggests that there is considerable ongoing transmission amongst cattle in Ethiopia. These findings are not surprising as very few farms in Ethiopia have the ability to control bTB in their herds, e.g. by test-and-slaughter based on the tuberculin skin-test and there is considerable long-distance trade between affected herds. Instead, this uncontrolled chronic disease results in cattle infected with being likely to survive for a long enough time to be subject to exposure to multiple bTB strains and/or cattle trade (short- or long-distance), increasing the risk of disease transmission and reinfection.

Depending on the SNP threshold used, there is considerable heterogeneity with respect to the size of defined transmission clusters and contribution from individual farms. With the more conservative threshold of 15 SNPs, genomes isolated from Farm Y (Fig. 3a) were part of three clusters suggesting epidemiological links both locally and nationally. Although we must be cautious in interpreting patterns within such a sparse sample and higher SNP threshold, this pattern is consistent with so-called ‘super-spreading’ behaviour where a small number of herds contribute disproportionately to transmission. The potential for such super-spreading, with an opportunity to target herds for control, was evident in our movement-based network analysis [73].

Additionally, we wanted to examine the possibility of zoonotic transmission between humans and cattle in Ethiopia. The minimum pairwise SNP distance between any pairs of human and animal isolates was 41 SNPs (median 658 SNPs) providing no evidence of recent zoonotic transmission from the cattle or herds that were sampled (potential transmission events would be highlighted by small pairwise SNP distances). This was unsurprising, given there were only six isolates from humans included in this study. However, two of the human isolates (of type SB1476) were only 11 SNPs apart from each other suggesting a potential epidemiological link; these isolates were from Bishoftu and Mekele, approximately 581 km apart, suggestive of long-distance travel of infected individuals or animals. Overall, these data provide clear evidence of considerable transmission between cattle but denser sampling will be required to establish specific evidence of zoonotic transmission. Denser, more targeted sampling would also allow for more sophisticated transmission analyses using tools such as TransPhylo [74], which account for within-host diversity and can estimate unsampled cases within transmission networks.

Recent work has suggested an East African origin for [43]. Given its position in the global phylogeny and confined geographical distribution [20], Af2 has likely been circulating in East Africa for hundreds of years; unfortunately, the lack of genomes from other parts of the world to provide context to the Ethiopian Af2 genomes, due to the small number sequenced, means confirming this in Ethiopia is currently not possible. However, we can show that there is a clear phylogenetic structure in this clonal complex with the Ethiopian genomes dividing into three distinct lineages. The position of Eu3 as the sister group to Af2 also implies a likely East African origin; however, based on the samples available for this study, the basal position of French isolates (Fig. 4b) suggests that the currently circulating Eu3 lineage in Ethiopia may in fact be European in origin. One possible explanation for this is that the ancestor of Eu3 was brought to Europe a few hundred years ago and that the modern Ethiopian and Eritrean genomes are descendants of that population, not of the ancestral population that may or may not still be circulating somewhere in East Africa. What is reasonably clear, given the long branches and subsequent expansion, is that there were two introductions of Eu3 sub-lineages, consistent with the study that first analysed the Eritrean WGS data [27], into Ethiopia between 1958 and 1987 (which also included Eritrea during that time) with the median estimate for introduction being in the early 1970s. Given the very similar dates, it is possible that these were part of the same series of cattle imports. In terms of likely origin of these imports, France should be viewed as a proxy for the diversity seen in mainland Europe (these samples were chosen to represent the diversity seen in France [50]), so the actual origin may be elsewhere in Europe. There are historical records of the first dairy cattle being imported into Ethiopia around 1950 as part of the United Nations Relief and Rehabilitation Administration (UNRRA) [75] with further subsequent imports from Kenya in 1959 [76]. Several livestock and dairy development projects that took place in the 1950–70s and funded by Sweden and the World Bank may have brought in dairy cattle of exotic breeds from overseas [3].

The other clonal complex found in Ethiopia, Unknown8, represented by two human genomes with the spoligotype SB1476, has thus far only been found in Ethiopia [11]. The position of this lineage in the global tree, and the hypothesized East African origin of [43], suggests that this clonal complex may be the ancestor of Eu1, the most prevalent and geographically distributed lineage known to date. However, further work would need to be done to confirm this through the collection of larger numbers of isolates with spoligotype SB1476.

This study performed the first WGS-based analysis of isolates from Ethiopia and attempted to include isolates from multiple herds in different parts of the country. There is considerable genetic diversity amongst in Ethiopia with multiple clonal complexes circulating and that they were likely to have been introduced in the country at different time-points. This work is important as it helps to better understand bTB transmission in cattle in Ethiopia and will potentially inform national strategies for bTB control in Ethiopia and beyond.

Click here for additional data file.

Click here for additional data file.