Genetic Diversity and Potential Paths of Transmission of Mycobacterium bovis in the Amazon: The Discovery of M. bovis Lineage Lb1 Circulating in South America

Bovine tuberculosis (bTB) has yet to be eradicated in Brazil. Herds of cattle and buffalo are important sources of revenue to people living in the banks of the Amazon River basin. A better understanding of Mycobacterium bovis (M. bovis) populational structure and transmission dynamics affecting these animals can significantly contribute in efforts to improve their sanitary status. Herein, we sequenced the whole genome of 22 M. bovis isolates (15 from buffalo and 7 from cattle) from 10 municipalities in the region of the Lower Amazon River Basin in Brazil and performed phylogenomic analysis and Single Nucleotide Polymorphism (SNP)-based transmission inference to evaluate population structure and transmission networks. Additionally, we compared these genomes to others obtained in unrelated studies in the Marajó Island (n = 15) and worldwide (n = 128) to understand strain diversity in the Amazon and to infer M. bovis lineages. Our results show a higher genomic diversity of M. bovis genomes obtained in the Lower Amazon River region when compared to the Marajó Island, while no significant difference was observed between M. bovis genomes obtained from cattle and buffalo (p ≥ 0.05). This high genetic diversity is reflected by the weak phylogenetic clustering of M. bovis from the Lower Amazon River region based on geographic proximity and in the detection of only two putative transmission clusters in the region. One of these clusters is the first description of inter-species transmission between cattle and buffalo in the Amazon, bringing implications to the bTB control program. Surprisingly, two M. bovis lineages were detected in our dataset, namely Lb1 and Lb3, constituting the first description of Lb1 in South America. Most of the strains of this study (13/22) and all 15 strains of the Marajó Island carried no clonal complex marker, suggesting that the recent lineage classification better describe the diversity of M. bovis in the Amazon.

Bovine tuberculosis (bTB) has yet to be eradicated in Brazil. Herds of cattle and buffalo are important sources of revenue to people living in the banks of the Amazon River basin. A better understanding of Mycobacterium bovis (M. bovis) populational structure and transmission dynamics affecting these animals can significantly contribute in efforts to improve their sanitary status. Herein, we sequenced the whole genome of 22 M. bovis isolates (15 from buffalo and 7 from cattle) from 10 municipalities in the region of the Lower Amazon River Basin in Brazil and performed phylogenomic analysis and Single Nucleotide Polymorphism (SNP)-based transmission inference to evaluate population structure and transmission networks. Additionally, we compared these genomes to others obtained in unrelated studies in the Marajó Island (n = 15) and worldwide (n = 128) to understand strain diversity in the Amazon and to infer M. bovis lineages. Our results show a higher genomic diversity of M. bovis genomes obtained in the Lower Amazon River region when compared to the Marajó Island, while no significant difference was observed between M. bovis genomes obtained from cattle and buffalo (p ≥ 0.05). This high genetic diversity is reflected by the weak phylogenetic clustering of M. bovis from the Lower Amazon River region based on geographic proximity and in the detection of only two putative transmission clusters in the region. One of these clusters is the first description of inter-species transmission between cattle and buffalo in the Amazon, bringing implications to the bTB control program. Surprisingly, two M. bovis lineages were detected in our dataset, namely Lb1 and Lb3, constituting the first description of Lb1 in South America.

INTRODUCTION
Mycobacterium bovis (M. bovis) is a member of the Mycobacterium tuberculosis complex (MTBC) and is the leading causative agent of bovine tuberculosis (bTB), an OIE (World Organization for Animal Health) notifiable disease that affects mainly cattle, buffalo, and other domesticated and wild animals, but can also be transmitted to humans (zoonotic TB) (1,2). bTB is distributed worldwide but has very low prevalence in most industrialized countries and has even been eradicated in few nations. However, the disease remains a major problem in developed countries with wildlife reservoirs that end up transmitting the pathogen to domestic livestock and vice-versa, and in developing countries where inefficient bTB control programs result in high disease endemicity and spread (1,2). In the Brazilian Amazon, few studies aiming to better understand bTB's epidemiology were performed (3)(4)(5)(6). The prevalence within animals in the area ranged from 0.1% (4) to 5.4% (6) and the major risk factors associated with bTB were the introduction of new animals into the herds (4), the buffalo species, herds with more than 100 animals, and the presence of cattle and buffalo in the same farm (6).
MTBC members evolved from a most recent common ancestor with M. canettii, and are characterized as clonal species demonstrating high genomic similarity (7). Currently, the use of whole-genome sequencing (WGS) to understand tuberculous mycobacteria populational structure is widespread and provided the basis for outbreak tracing and phylogenetic analysis resulting in the classification of human-adapted MTBC into 8 lineages, with M. tuberculosis accounting for L1 to L4 and L7-L8, and Mycobacterium africanum comprising of L5 and L6 (8,9). On the other hand, M. bovis has been historically classified by Clonal Complexes (CCs), which are identified by genomic deletions, few Single Nucleotide Polymorphism's (SNPs), and/or spoligotypes patterns (10). Accordingly, four different M. bovis CCs have been described presenting distinct geographical distribution patterns: African 1 and 2 restricted to Africa, European 2 commonly found in the Iberian Peninsula, and European 1 distributed globally (11)(12)(13)(14). With the advent of WGS, recent studies at a global scale provided insights into the population structure and evolution of M. bovis lineages (9), showing that CCs do not represent the whole genomic diversity of the isolates (9,15,16) and suggesting the existence of at least four M. bovis lineages, named Lb1 through Lb4, and three "unknown groups" (9). With the populational structure of M. bovis based on WGS starting to be unveiled, additional studies covering different geographic locations are needed to better comprehend worldwide disease spread and to provide new insights regarding the use of genomes to understand disease transmission at the herd and farm levels.
Molecular epidemiological investigation has proved to be a useful tool for TB control and surveillance, which allows us to better understand the dynamics of disease transmission and precisely identify the infectious agent (17). In addition, the knowledge regarding strain diversity within host species has special contribution in areas under risk of zoonotic TB occurrence, thereby providing new insights in strain distribution that may help establishing strategic measures for TB control and prevention (18,19). In Brazil, few and dispersed molecular epidemiologic studies of M. bovis have been reported (20)(21)(22)(23)(24), and the WGS characterization of the pathogen is just starting to be unraveled (24)(25)(26). While few studies focused on areas of high dairy herd productivity (25,26), a recent study evaluating M. bovis genomes obtained from buffalo and cattle of the Marajó Island, Northern Brazil, was performed and showed the existence of a monophyletic group without CC classification (24). Therefore, the aim of this study was to apply whole-genome and SNP-based phylogenomic analyses to obtain novel information regarding the genetic diversity of M. bovis strains circulating in buffalo and cattle from the region of the Lower Amazon River Basin. We believe that such information will guide policy development and strategies to contain the disease in livestock, and thus reduce the risk associated with transmission to humans.

DNA Extraction
M. bovis isolates were reactivated in Stonebrink media and incubated until positive growth at 37 • C. DNA extractions from colonies suggestive of M. bovis for genomic sequencing were performed according to the protocol of van Embden et al. (27), with modifications. Initially, for inactivation, 2-3 colonies were resuspended in 400 µl TE buffer (10 mM Tris-HCI and 1 mM ethylenediaminetetraacetic acid-EDTA, pH 8.0) and heated at 80 • C for 30 min. Subsequently, 50 µl of lysozyme (10 mg/ml) were added and incubated at 37 • C for 1 h. Then, 75 µl of 10% SDS (sodium dodecyl sulfate) and 10 µl of proteinase K (10 mg/ml) were added and incubated at 65 • C for 10 min. Next, 100 µl of 5M NaCl (sodium chloride) and 100 µl of CTAB (cetyltrimethylammonium bromide) were added, followed by stirring and incubation at 65 • C for 10 min. After that, 750 µl chloroform/isoamyl alcohol (24:1) were added, stirred and centrifuged at 12,000 g for 5 min. The aqueous phase (surface) was transferred to another tube, 450 µl of isopropanol were added and incubated at −20 • C for 30 min and then centrifuged again at 12,000 g for 15 min at room temperature. The supernatant was discarded, and the pellet washed once with 1 ml of ice-cold ethanol (70%) with centrifugation at 12,000 g for 5 min. After drying the tube by evaporation at room temperature, the DNA was resuspended in 20 µl of TE buffer and stored in a freezer at −20 • C. The quality and concentration of the extracted DNAs were evaluated using Nanodrop (Thermo Fisher Scientific). Procedures were performed in a Biosafety Level 3 Laboratory located at the Embrapa Gado de Corte, Campo Grande, Brazil.

Genome Sequencing
WGS was performed at the NGS multi-user platform of Oswaldo Cruz Foundation (FIOCRUZ), Rio de Janeiro, Brazil. Briefly, DNA quantification was performed using the Qubit TM dsDNA HS Assay Kit (Thermo Fisher Scientific, Waltham, USA) and the Agilent High Sensitivity DNA Kit (Agilent, California, USA). WGS of M. bovis isolates was carried out on a HiSeq instrument (Illumina, San Diego, CA) using HiSeq Rapid SBS Kit v2 (200 cycles) chemistry and the Nextera DNA Flex Library preparation kit (Illumina, San Diego, CA) according to the manufacturer's instructions. Sequencing reads were deposited in Sequence Read Archive (BioProject number PRJNA675550), NCBI and accession numbers are described in Supplementary Data Sheet 1.

Genome Quality Assessment and Regions of Difference Identification
Obtained reads were trimmed using a Trimmomatic version 0.38 (28) for adapters and low-quality base removal (sliding window 5:20). Trimmed reads were then evaluated for reads size, per base and read sequence quality, presence of adapters, and GC content using FastQC (29). The GC content had to be around 65% (which is common to mycobacterial genomes) and without multiple peaks (i.e., possible contamination with sequences of different GC content) per quality criteria.
As to confirm that genomes were from M. bovis, reads were mapped against M. tuberculosis H37Rv using Burrows-Wheeler Aligner (bwa-mem) (30) and the positions according to reference genome of RD1 (4,354,000-4,358,331 nt), RD4 (1,696,017-1,708,748 nt), and RD9 (2,330,880-2,332,100 nt) were evaluated for coverage as previously described (31). The genome was considered as M. bovis species when RD1 was absent (i.e. region was intact) and RD4 and RD9 were present (i.e. regions were deleted) (32). Obtained coverage against M. tuberculosis H37Rv was also used as quality criteria, considering 95% as a minimum mapping percentage for a genome to be included in the analysis.

Spoligotyping and Clonal Complexes
Spoligotypes were also investigated in silico using SpoTyping (31). Identified genetic spacers were processed in the M. bovis Spoligotype Database (www.mbovis.org) to retrieve a spoligotype pattern and SB number. The CCs African 1 (Af1) and 2 (Af2) and European 1 (Eu1) and 2 (Eu2) were evaluated as previously described (9). Briefly, the SNP in the guaA gene was investigated in the bam files generated from read mapping against M. tuberculosis H37Rv, by checking the position 3,813,236. For the RDs, the same bam files were used to investigate read depth using samtools depth (33) and GNU parallel 2018 (34) for the following regions: RDEu1 (1,768,074-1,768,878 nt), RDAf1 (665,042-668,394 nt), and RDAf2 (680,337-694,429 nt) also as described previously (9).

Phylogenetic Reconstruction
As to evaluate genetic diversity of M. bovis in the geographic region, 15 quality-approved M. bovis genomes previously sequenced and obtained from cattle and buffalo of the Marajó Island (24) (ENA accession number ERP116404) were used to evaluate the phylogenetic relatedness with the M. bovis genomes sequenced in this study. The Marajó Island is geographically close to the targeted region analyzed herein (Figure 1), together with the Lower Amazon River Basin are the regions of major concentration of bufallo in the country, and these M. bovis genomes (24) were the only strains of the North of Brazil sequenced up until this study. These genomes were qualityassessed and SNPs were obtained as described above. A matrix of concatenated SNPs of the M. bovis genomes was constructed as described (9) and used to estimate a maximum likelihood (ML) phylogeny. RAxML version 8.2.12 (37) was used to construct this phylogenetic tree, selecting the GTRCAT model and autoMRE for best-scoring ML tree and a maximum of 1,000 bootstrap inferences. Genomes of Mycobacterium caprae (M. caprae) (ERR1462591, ERR1462625, ERR1462617, ERR1462581) were also included in the SNP matrix to serve as outgroup.

Minimum Spanning Tree
A pairwise SNP-distance matrix and a minimum spanning tree were constructed using PHYLOViZ (38) with default parameters and using a concatenated SNP matrix of the M. bovis genomes as input.

Phylogenomic Analysis With Lineage Representatives of M. bovis Genomes
Mycobacterium bovis genomes obtained from this study and those from Marajó Island (24) were compared against a collection of 106 genomes representing the four M. bovis lineages (Lb1-Lb4) and unknown groups 1, 2, and 3 from a recent study that evaluated worldwide distribution and evolution of M. bovis genomes (9) (Supplementary Data Sheet 1). In addition, 21 quality-approved M. bovis genomes from France (39), which include representatives of the recently suggested CC called European 3 (40), were also included in the analysis were included to construct the phylogenetic tree using the ML approach as described above.

Pairwise SNP-Comparisons
From the SNP-distance matrix, pairwise distances distributions between M. bovis genomes obtained from cattle vs. buffalo (Amazon dataset) and between M. bovis genomes originating from the Lower Amazon River Basin and Marajó Island were compared using the non-parametric Mann Whitney-test in R software. A result was considered statistically significant when p-value ≤ 0.05.  (Figure 1). The genotypic characterization by spoligotyping revealed 5 distinct profiles ( Table 1). The predominant spoligotype was SB0822, detected in 10 isolates, followed by SB0295, representing 6 isolates; the patterns SB1190 and SB1800 were identified in 2 isolates each; the pattern SB0121 had one representative; and one isolate without described spoligotype pattern. Surprisingly, out of 19 M. bovis isolates previously typed using an experimental technique (20), 12 were discordant. Disagreement between in silico and experimental spoligotyping has been previously described (41,42). Finally, comparing hosts, our results show a higher number of spoligotypes patterns in buffalo when compared to cattle ( Table 1). Given the low sample size, further studies should be conducted to confirm if buffalo have consistently higher diversity of spoligotype patterns when compared to cattle in this region.

Out
Among the CCs, only the CC Eu2 was found among 41% (9/22) of the isolates. The remaining 13 samples (59%) were not identified as belonging to any known CC (i.e., Eu1, Eu2, Af1, or Af2). Within buffalo, 40% (6/15) of the samples were Eu2, 60% (9/15) had no CC marker. Within cattle, 42.8% (3/7) of the samples demonstrated the marker of CC Eu2, and the remaining of the isolates (4/7) were not classified within any described CC. In the generated phylogenetic tree, with M. bovis strains from the Lower Amazon River Basin and Marajó Island, the host classes (buffalo and cattle) are found dispersed among different clades (Figure 2), while M. bovis genomes from Marajó Island clustered together, appearing more closely related, and genomes from different Amazonas municipalities did not follow a clear clustering pattern according to geographic region.
Pairwise-SNP comparisons (Figures 3, 4) and minimum spanning tree (Figure 5) of the M. bovis genomes obtained herein and from the Marajó Island (24) showed a highly diverse genomic dataset, with pairwise SNP-distances varying from 8 to 711 (from 8 to 100 among M. bovis genomes from Marajó Island and from 12 to 711 among M. bovis genomes of the Amazon) (Figure 3). Although the overall distribution was not significantly different ( Figure 4A), M. bovis genomes from the Lower Amazon River region tend to show higher pairwise SNP-distances than M. bovis genomes from the Marajó Island (Figure 3). There was also no significant difference of M. bovis genetic diversity distribution using the dataset from the Lower Amazon River region as a function of host ( Figure 4B; host information from Marajó Island was not available in ENA).
Reflecting this high diversity, based on current M. tuberculosis-based SNP threshold to infer recent transmission links (∼12 SNPs) (44)(45)(46), only three possible active transmission clusters can be observed, one in Marajó Island (host unknown), one connecting the municipalities of Urucará and Parintins (between two buffalo, TB051 and TB050), and one connecting the municipalities of Careiro da Várzea, Novo Céu and Itacoatiara (involving cattle and buffalo, TB54, TB46, and TB41), demonstrating recent transmission between different herds and hosts (Figure 5). Isolates comprising each transmission clusters also had the same spoligotype pattern (SB0822, SB0822, and SB0295, respectively).
Based on the proposed four major global lineages of M. bovis (Lb1, Lb2, Lb3, and Lb4) (7) our setting would be composed by two lineages, Lb1 with two buffalo isolates (TB052 and TB042) and Lb3 with six buffalo and three cattle isolates with the CC marker Eu2 and seven buffalo and four cattle isolates without CC markers (Figure 6 and Table 2). All isolates from Marajó Island were also identified as being from Lb3 without a CC marker. The presence of Lb1, infecting two buffalo, is the first description of this M. bovis lineage in Brazil and requires further investigation into the actual origin of these isolates.
In order to evaluate if the high genetic diversity observed in M. bovis from the Lower Amazon River region was only due to the presence of the two highly divergent Lb1 genomes, we compared the distributions of pairwise-SNP distances between the Lower Amazon River region and Marajó Island for the Lb3 only ( Figure 4C). No significant difference in M. bovis genetic diversity was observed ( Figure 4C). However, the maximum SNP-distance between two genomes observed in the Marajó Island was 100, while in the M. bovis genomes of the Amazon (Lb3 only) continued to be high. There was also no

DISCUSSION
To investigate the clonality and population structure of M. bovis in the study area at first, we relied on the use of spoligotyping in the characterization of 22 M. bovis isolates from 15 buffalo and 7 cattle from June 2016 to October 2017. The spoligotype SB0822 in Brazil was first described in our previous study (20) and we agree with the recent study in Marajó Island in Pará state (24) which seems to demonstrate that SB0822 is the predominant spoligotype in the Amazon region, opposite to a national study with 143 samples from 10 states that found only one occurrence of SB0822 (21). Moreover, SB0822 has been previously described in France (47), cattle in Spain (48) and Portugal (49), and buffalo in Colombia (50), and overall agree with the history of the livestock in the region where the first animals introduced came from Cape Verde (a former Portugal's colony) initially to Marajó Island and from there expanded to the floodplains of the Lower Amazon, on the banks of the Amazon River (51).
Our results also show a higher number of spoligotype patterns of M. bovis in the Lower Amazon River Basin compared to the Marajó Island, which can be explained due to the frequent movement of animals in our area of study and the isolation of the herds in the Marajó Island. It is important to highlight that the diversity of M. bovis found within buffalo in our sampling is unlikely a result of recent introduction of animals from other Brazilian states since there has been no buffalo imported from other states to the Amazon region and this is probably maintained by constant reinfection from reservoir animals. Accordingly, SB0121 (the most prevalent spoligotyping in Brazil) was found in only one sample, which may reflect the low transit of animals from other states to the area of this study.  Our results show overall higher genetic diversity of M. bovis genomes obtained from different municipalities of Lower Amazon River region when compared to the Marajó Island. This is most likely due to the geographic isolation of the island with lower chances of animal importation over time. We also show possible transmission links between buffalo and cattle from different herds but from close geographic proximity. As wildlife reservoirs have not been identified in Brazil thus far, this transmission may have occurred due to infected animal transit and/or introduction into different herds, showing the presence of a two-host system allowing inter-species transmission in the region.
The detection of an inter-species transmission link, the intertwined phylogenetic dispersal of M. bovis obtained from cattle and buffalo, and the absence of significant difference in M. bovis genetic diversity in cattle vs. buffalo suggest that contact rate between different hosts and consequent geographic proximity likely played a more important role in determining the host range of M. bovis in this region than host species, agreeing with recent studies (10, 24). Currently the most accepted hypothesis is that M. bovis is not a specialized pathogen, i.e., can affect several host species irrespective of its genetic makeup. However, we cannot neglect the possibility of differential host susceptibility to different strains or lineages of M. bovis, which may allow one particular strain or lineage to thrive in a specific host. Herein and in the study in Marajó Island (24), for instance, there is only 40% power to determine a possible significant difference in distribution of M. bovis strains between buffalo and cattle. Therefore, studies with a higher sample size comparing the dynamics of M. bovis infection within buffalo and cattle are needed to definitively clarify this research question.
The results from our study support the fact that current CCs cannot represent the whole diversity of M. bovis strains.
Interestingly, this is the first time that is described in Brazil isolates of M. bovis originating from the M. bovis lineage Lb1 (9). In a recent global phylogenomic study of M. bovis, strains of Lb1 were shown to emerge from older nodes than Lb3 in the phylogenetic tree and were detected in Eritrea, Ethiopia, Tanzania, Uganda, Tunisia, France, Spain, Italy, Switzerland, and the United States (9). However, some strains identified in this lineage carry the CC Af2 marker (9) which is found at high frequency in bTB cases from East Africa (12), demonstrating the important ties of this lineage to the African continent. One hypothesis to be looked at is that due to the proximity of these countries with Portugal and its colonies; these strains might have been first introduced in the Amazon region during Brazil colonization.
Results from this study, along with others (9,16,39,40), can be used to refine the understanding of M. bovis lineage Lb1. With current genomic data, this lineage is composed of two main clusters: one made of M. bovis genomes carrying the CC Af2 marker, and the other without known CC markers (9,16). Herein, the Lb1 cluster without CC Af2 marker includes 19 out of the 21 French M. bovis genomes of SB0120 sequenced by Hauer et al. (39) and included in this study. Recently, Branger et al. (52) suggested that this phylogenetic cluster be called CC European 3, also based on their previous work (39). In this study and in others (9,16), this cluster is composed of M. bovis genomes of many spoligotype patterns (e.g., SB0120, SB0134, SB0828, SB0948, SB1517; Supplementary Data Sheet 1) and of BCG vaccine strains (16,39). Although 19 specific SNPs have been provisionally suggested to be specific of Lb1 (9) and 5 SNPs of its non-Af2 cluster [i.e., CC Eu3; cluster I in Hauer et al. (39)], further studies using comprehensive global datasets should be conducted to confirm or identify definitive genetic markers. The remaining two French M. bovis genomes, SRR7851366 and SRR7851376, grouped with genomes of Lb3 and "unknown 3" group, respectively; their disparate position on the phylogenetic tree corroborates the finding of Hauer et al. (39).
The finding of isolates from Marajó Island belonging to Lb3 without CC marker reinforces the previously described (24) existence of a unique M. bovis clade in the island that was likely introduced in a single event. In contrast, our results suggest that M. bovis was introduced into the Lower Amazon River region as three different events, for which the temporal order remains to be evaluated. One introduction is related to the neighbor cities of Parintins and Urucara, with strains Lb1. Another introduction occurred with Lb3 strains without the CC Eu2 marker, probably with the same origin from the Marajó Island. And finally, an additional introduction is observed with Lb3 strains carrying the CC Eu2 marker, likely from cattle imported from other states and spreading to buffalo. In 2019, a study with 90 samples of cattle lesions suggestive of bTB from the states of Goias, Mato Grosso, Mato Grosso do Sul, Minas Gerais, São Paulo, Tocantins, and Pará found that 14.4% (13/90) belonged to the CC Eu1 and 81.1% (73/90) to the CC Eu2, while 4.65% (4/90) were not identified as any of the four known complexes (53). As the isolates without CC markers were not classified into lineages, the data collected seems insufficient to reveal the true epidemiological picture of the bTB in Brazil. Knowing the genetic profile and understanding the transmission routes of M. bovis in the Amazon and elsewhere is essential in order to focus on public health and veterinary resources to contain bTB.
A limitation of this study is the sample size. The small number of M. bovis isolates may not be representative of the whole bacterial and animal populations of the Lower Amazon River Basin. Results must be interpreted in light of this fact, and efforts should continue to isolate and study additional M. bovis strains from the region.
From our study, we can make the following conclusions: (1) The M. bovis CCs classification cannot cover the whole diversity of M. bovis strains present in the Amazon region; (2) The presence of M. bovis strains Lb1 infecting buffalo requires further investigation into the actual origin of these isolates, showing that the true global diversity of M. bovis strains remains to be discovered, likely influenced by cattle trade over history; and (3) The M. bovis classification in lineages by SNP-based phylogenetic analyses seems to better cover the diversity of M. bovis strains present in the Amazon region compared to CC classification.

DATA AVAILABILITY STATEMENT
Sequencing reads were deposited in Sequence Read Archive (Bio Project Number PRJNA675550), NCBI and Accession Numbers are described in Supplementary Data Sheet 1.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because no animal was used in this study.

AUTHOR CONTRIBUTIONS
PC and JBK designed experiments, coordinated the research team, interpreted the data, wrote and revised the manuscript. HT and CS coordinated field sampling. CS coordination of field sampling. RA reviewed the manuscript. CZ designed experiments, assessed reads quality, performed phylogenomic analysis, interpreted the data, wrote and revised the manuscript. TTP performed pairwise SNP comparisons, constructed the minimum spanning tree, and interpreted the data. AG designed experiments, supervised CZ and TTP, interpreted the data, wrote and revised the manuscript. FA and TNP cultured M. bovis, and performed PCR and DNA extraction. AD performed whole genome sequencing. All authors contributed to the article and approved the submitted version.