The Program for Monitoring and Control of Bovine Tuberculosis in Bulgaria includes annual, one-time intradermal tuberculin testing with bovine PPD tuberculin of all cattle over 42 days of age. Doubtfully PPD tuberculin-reacted animals undergo differential tuberculin testing with bovine and avian PPD tuberculin on day 42 after the first tuberculin test. Animals positively reacted to regular and differential tuberculinization are sent for sanitary slaughter to isolation slaughterhouses. In herds with detected tuberculosis, all cattle over 42 days of age are tuberculinized every 2 months until negative results are obtained in the intradermal and laboratory tests. The animals are tested again after 6 months and in case of repeated negative results, the herd is considered free of tuberculosis. In 2015-2020, 66 TB-affected farms with 851 TB-positive animals were identified in 11 provinces across the country. Lymph nodes from slaughtered cattle were sent to the National Reference Laboratory for bacteriological and PCR тesting and confirmation of results. All 851 cases were confirmed by both PCR and bacteriology, as requested by National bTB Control Program implemented since 2015.
Position of Bulgarian M. bovis/M. caprae within global WGS-based phylogeny
The present research included a representative sample of 40 M. bovis/M. caprae isolates from lymph nodes of TB-positive cattle available in the strain collection in National Diagnostic and Research Veterinary Medical Institute “Prof. Dr. G. Pavlov”, Sofia, Bulgaria. These isolates represent the areas with the most affected farms and those with isolated cases of the disease that represent 60.6% of the infected herds in 2015-2020.
Of the 40 DNA samples from bacterial isolates, 37 were found of sufficient quality and concentration to be suitable for WGS. Of them, 3 isolates had low-quality sequencing reads and were excluded from further analysis and 34 isolates were successfully SNP genotyped (Additional file 1 with Table S1). Phylogenetic analysis was carried out on the 34 Bulgarian isolates along with a selection of genomes from a global dataset of the M. bovis WGS data compiled and curated at the National Veterinary Services Laboratories, United States Department of Agriculture, Ames, IA, USA, including publicly available genomes in the NCBI Sequence Read Archive (Additional file 2 with Table S2). High-quality SNPs were identified that clustered the isolates into smaller more manageable groups. Groups were created based on the number of isolates as well as the evolutionary distance. We identified the most recent common ancestor (MRCA) in the groups. Then, the most closely related isolate was selected based on having the fewest number of SNPs differences, and a pairwise comparison of SNPs to the MRCA was recorded.
In addition, we performed in silico spoligotyping of the analyzed strains and assignment of the international code number according to Mbovis.org database (Table S1). In total, three spoligotypes were identified and included 28 M. caprae isolates (SB0418), 5 M. bovis isolates (SB0120) and 1 M. bovis isolate (SB0339). Some of the studied isolates were previously subjected to classical macroarray-based spoligotyping [23]. Both in silico and experimental spoligotypes were available for 24 isolates and showed 100% concordance.
The Bulgarian M. bovis/M. caprae isolates fall into 3 major phylogenetic groups identified in this work as Mbovis-13, Mbovis-37 and Mcaprae-27B (Fig. 1). The groups were different in size and phylogenetic relationships. The Mbovis-37 group included 5 isolates from Bulgaria. These isolates share a distant common ancestor with isolates from Poland (Fig. 2A). The most recent common ancestor to these isolates has only accumulated 2 SNP (red bracket) since sharing a common ancestor with the isolate from Poland (green arrow). Furthermore, the Mbovis-37 isolates from Bulgaria are very closely related, all being ≤11 SNP from a shared common ancestor (Fig. 2B). They differ by 22-29 SNPs from their closest known relative (023056-351-19) an isolate from Poland but are distinct from all other isolates by (> 200) SNPs. Mbovis-37 group does not belong to any of the clonal complexes currently defined for M. bovis (Eu1, Eu2, Af1, Af2, and BCG-like).
Only one Bulgarian isolate belongs to the Mbovis-13 group (Fig. 1). Based on the global SNP framework, this group corresponds to the European Complex 2, thought to have evolved on the Iberian Peninsula and distributed throughout western continental Europe and some regions of northern Africa [18, 24, 25]. The Mbovis-13 isolate shares its most recent common ancestor with two isolates from humans in Germany (red arrow); altogether this German/Bulgarian group shares a common ancestor with an isolate from France (purple arrow) (Fig. 3).
Finally, the largest group of 28 M. caprae isolates clustered in the Mcaprae-27B branch (Fig. 4). This branch includes isolates from Poland, Spain, Germany, the Republic of Congo and one of unknown origin (labeled as Rhine in Fig. 4). The isolates from Bulgaria share the most recent common ancestors with isolates from Spain (pink arrows). Almost all the M. caprae isolates from Bulgaria form a unique clade of closely genetically related isolates, sharing a common ancestor within ≤15 SNP. Three distinct subclusters were observed within this branch that we have labeled as A, B, and C for convenience of discussion of their geographic specificity at the intra-country level (Fig. 5).
We additionally compared our WGS results to the recently proposed new classifications to define M. bovis phylogenetic lineages and sublineages [18, 26, 27]. Zwyer et al. [27] proposed to divide M. bovis, M. caprae and M. orygis in three main phylogenetic lineages named La1 (further subdivided into La1.1-La1.8), La2, and La3. We assigned the Bulgarian clusters to these lineages as follows: Mbovis-37 to La1.7; Mcaprae-27B to La2; Mbovis13 to La1.7.1. Zimpel et al. [18] proposed four distinct global lineages of M. bovis labeled as Lb (Lb1, Lb2, Lb3, and Lb4). Based on their classification, Bulgarian isolates of the Mbovis-13 and Mbovis-37 clusters were assigned to Lb3.
A subset of sequences representing each of the clonal complexes described by Loiseau et al. [26] was retrieved from NCBI and a phylogenetic analysis was performed adding the Bulgarian isolates to determine their location in the phylogeny. The Bulgarian sequences fell in European Complex 2, other unknowns, and M. caprae. Reference genomes M. bovis AF2122/97, M. bovis BCG, M. bovis AN5, and M. bovis Ravenal were also included for additional perspective (Additional file 3 with Fig. S1).
The geographic diversity of the studied isolates is illustrated in the map in Fig. 5A. Our previous study provided the first insight into phylogeography of M. bovis/caprae in Bulgaria based on spoligotyping. In spite of the well-known limitation of that method, a certain gradient was observed at the within-country level. Classical M. bovis types predominated in the Northeast, while M. caprae was prevalent in Central and Southwestern Bulgaria. By spoligotyping, these M. caprae isolates belonged mostly to the Central/Eastern European cluster [23].
In the present study, higher resolution WGS/SNP typing confirmed these findings. The identified M. bovis isolates were from Northeastern Bulgaria. The five isolates from the Mbovis37 group (Fig. 2A) did not belong to any well-defined cluster and were found in the districts of Razgrad and Shumen. The most genetically distinct isolate from the Mbovis13 group (Fig. 3) was from the Dobrich district.
The most numerous sample of isolates belong to the Mcaprae-27B group (Fig. 4A) and its three subclusters A, B and C predominate in the central and southwestern part of the country (Fig. 5). M. caprae isolates of all three subclusters have been identified in the districts of Pazardjik and Plovdiv. In addition, subcluster A includes an isolate from Kardjali, subcluster B – two isolates from Pernik, and subcluster C – an isolate from Stara Zagora district (Fig. 5). These three subclusters are also clearly distinguished on the minimum spanning tree of the Mcaprae-27B group (Fig. 6).
The identified phylogenetic relationships of the studied isolates from cattle in Bulgaria and other countries observed in Figs. 2 and 3 may be partly explained by livestock trade contacts with regard to livestock commercial chains. Concerning this particular study, data on the live cattle export from and import to Bulgaria from available reports are helpful but changing trade directions through years should be taken into consideration. For example, the recent annual report on livestock 2018-2019 in Bulgaria [28] mentioned the Czech Republic, Germany, and Hungary as sources of live cattle imports to Bulgaria, but not Spain nor Poland.
On the other hand, the very long branch (144 SNPs) that separates Bulgarian Mbovis-37 isolates from most of those from Poland, poses a different scenario (Fig. 2a). Given such a large genetic distance, it is difficult to decipher the source of these isolates.
The Polish isolate 20-023056-351-19_2019_cattle_PL was the closest relative to Bulgarian Mbovis-37 with the distance of 22-29 SNPs. This strain came from Płońsk county in the Masovian Voivodeship in central Poland and it was the only isolate from this location in the national database. Unfortunately, the epizootic investigation failed to identify the source of the infection. This individual in the past was in several different farms where tuberculosis was not recorded. However, there are suspicions that the infection occurred in the Płock county, where the tested cow lived in 2014. In 2014, many bTB cases were recorded in this county. There is a high probability that M. bovis transmission may have occurred between neighboring herds. A distance of 22-29 SNPs separating this Polish isolate from the closest Bulgarian isolate seems too long to speculate about recent links. However, when this distance is translated into > 30 years to their divergence this is reminiscent about the time when two countries were economically united within the Soviet bloc. While the well-known Warsaw Treaty was a military union, there was also The Council for Mutual Economic Assistance focused on economic development and exchange that existed from 1949 to 1991 (https://en.wikipedia.org/wiki/Comecon).
The Bulgarian M. caprae are also quite distant from their closest known relatives from Spain in the Mcaprae-27B group, with at least a 150 SNP genetic distance (Fig. 4A). Similar to the above reasoning, they are distantly related by WGS analysis. It is possible that more closely related intermediate strains may be extinct because of slaughter.
Genetic diversity of Bulgarian M. caprae/M. bovis and insight into transmission
Previously, it was suggested for human M. tuberculosis studies that epidemiologically linked isolates consistent with recent transmission differ by five or fewer SNPs. In contrast, isolates differing by more than 12 SNPs are not directly linked, while pairs differing by six to 12 SNPs are indeterminate [29]. The same study calculated the mutation rate within M. tuberculosis genome to vary from 0.3 to 0.7 SNP/genome/year while 0.5 is generally accepted as a mean value. Previous M. bovis studies estimated mutation rate in different settings and with wildlife hosts from 0.2 SNP/genome/year (95% HPD: 0.1 – 0.3) [30], to 0.37 (95% HPD: 0.24–0.51) [31] and 0.53 (95% HPD: 0.22-0.94) [32]. In view of these heterogeneous estimations, we believe that the mean 0.5 mutation rate used in our study is acceptable.
In view of these estimations, we interpreted minimum spanning trees (MST) built for two groups of the Bulgarian isolates, Mcaprae-27B (28 isolates, Fig. 6) and Mbovis-37 (5 isolates, Fig. 7). Taken together, 12 Mcaprae-27B and 3 Mbovis-37 isolates were separated by 1 to 6 SNPs within four clusters that could correlate with the relatively recent transmission. Other isolates were more divergent and could not be assigned to the shared clusters.
The MST of the largest group M. caprae-27B shows a high diversity (Fig. 6). Three isolates with subcluster B (upper part of the tree) differ by 1-3 SNPs: two isolates from Pernik and one from Pazardjik, both locations are neighboring compared to other locations in this study. Also, in the central part of the tree, we observe two isolates that differ by 4 SNPs, from neighboring provinces of Plovdiv (parental isolate, 2015) and Kardjali (distal derived isolate, 2019). These two isolates are separated by 4 years and their location on the tree (the most recent isolate being the distal one) and their genetic distance of 4 SNPs roughly correlate with the expected mutation rate of 0.3-0.7 SNPs/genome/year. Interestingly, all three M. caprae isolates from Kardjali were very genetically distant from one another. The same is true for isolates from Stara Zagora. This may be due to the undersampling of these locations, however, even the largest sample of isolates from Plovdiv showed significant genetic divergence.
The MST tree of the 5 M. bovis isolates (Mbovis-37 group) revealed one isolate in the central position and two of the derived isolates differing from it in 6 SNP (Fig. 7). This may indicate a recent transmission however these isolates were recovered in 2015-2016 and no further isolates were identified. Two other M. bovis isolates from this group differed from the central node by 10 or more SNPs and were therefore not likely to have undergone recent transmission.
A closer look at the year of isolation reveals that in some instances more distal and derived isolates were isolated earlier than preceding “parental” isolates (Figs. 6 and 7). This discrepancy may suggest that true parental isolates were not effectively isolated because the infected animals were not identified. The discrepancy could not be due to the under-testing because all herds are regularly tuberculin-tested in Bulgaria. Another plausible explanation for this could be that these animals had a latent form of bTB that was not detected during surveillance [33, 34]. The role of wildlife in the transmission of bTB in Bulgaria should be also taken into consideration given the published reports about M. bovis confirmed in badger and wild pigs [12,13,14]. However, these reports date back to 1995-2003, and no papers on bTB in wildlife in Bulgaria were published recently.
Interestingly, M. caprae isolates in this study were all from the southern and southwestern parts of Bulgaria. M. caprae was initially isolated from goats and is noteworthy that sizeable goat populations are present in the neighboring region of Northern Greece. While data on M. caprae zoonotic TB in Greece is limited, dairy goat farm outbreaks and goat to human transmission in northern Greece have been documented since 2005 ([35] and references therein). In Bulgaria, goat breeding is not much developed. Usually, private farmers keep from one to several animals, which are not subject to monitoring and control for the presence of tuberculosis. The above recent study in Northern Greece described M. caprae transmission from goat to the goat breeder [35]. The MIRU-VNTR phylogenetic analysis together with reference isolates from MIRU-VNTRplus.org placed that Greek isolate within the branch dominated by SB0418 spoligotype. Noteworthy, all but one M. caprae isolates in this Bulgarian study also presented this spoligotype. In this view, a link between neighboring areas of northern Greece and southwestern Bulgaria appears plausible but should be addressed by WGS analysis.
