Transmission Network of Deer-Borne Mycobacterium bovis Infection Revealed by a WGS Approach

Bovine tuberculosis (TB) is a zoonotic disease, mainly caused by Mycobacterium bovis. France was declared officially TB free in 2001, however, the disease persists in livestock and wildlife. Among wild animals, deer are particularly susceptible to bovine TB. Here, a whole genome sequence (WGS) analysis was performed on strains with the same genetic profile—spoligotype SB0121, Multiple Loci VNTR Analysis (MLVA) 6 4 5 3 11 2 5 7—isolated from different types of outbreaks, including from deer or cattle herds, or zoological or hunting parks where the presence of infected deer was a common trait in most of them. The results of the phylogeny based on the SNP calling shows that two sub-clusters co-exist in France, one related to deer bred to be raised as livestock, and the other to hunting parks and zoos. The persistence over almost 30 years of sporadic cases due to strains belonging to these clusters highlights the deficiency in the surveillance of captive wildlife and the need for better monitoring of animals, especially before movement between parks or herds.


Introduction
Bovine tuberculosis (bTB), mainly caused by Mycobacterium bovis, is an important re-emergent zoonotic disease in Europe [1]. Although France has been officially bTB free (OTF) since 2001, the persistence of the disease in livestock and its occurrence in some areas of wildlife are matters of great concern [2]. This ancient disease is no longer considered just a cattle-borne problem, but a concern for multihost communities that include wildlife species such as wild boar (Sus scrofa), red deer (Cervus elaphus), and badgers (Meles meles) [3].
In Europe, red and fallow deer (Dama dama) are frequently reported to be infected by M. bovis or M. caprae [4]. In the wild population, red deer are considered to play an important role in the dissemination of the disease as an important number of infected animals disclose generalized patterns of clinical tuberculosis [5]. Farmed deer are also very susceptible to M. bovis, as a high incidence of animals with generalized infection is commonly observed in deer herds [6]. Gross lesions are generally observed in lymph nodes of the head, thorax and abdomen [4]; most frequently in the retropharyngeal lymph nodes for red deer, and the thoracic region for fallow deer [7].
Interaction among different wildlife species, and between wildlife and human beings, may now be more frequent due to the increasing amount of captive wildlife in zoos or other types of animal parks [8]. The zoonotic transmission of Mycobacterium tuberculosis complex (MTBC) pathogens has been reported worldwide: e.g., an M. bovis infection in zookeepers exposed to a white rhinoceros [9], an M. pinnipedii infection in animal handlers linked to a sea lion colony [10], and an M. tuberculosis infection of staff members exposed to elephants [11,12]. Infrastructure where the public has potential contact with captive wildlife is a serious emerging concern for zoonotic transmission, especially for immunocompromised individuals [8]. The same concern arose after the discovery in an Iranian zoo of an M. bovis infection in free-roaming fallow deer that had been in close contact with the public [13].
MTBC members are excellent microorganisms for conducting molecular epidemiology studies on, by virtue of their low mutation rate nature. In the case of bTB, two genotyping methods are widely used-spoligotyping [14] and Multiple Loci VNTR Analysis (MLVA) typing [15]-both of which establish differentiation between M. bovis strains, taking into account relative genetic stability. The combination of these two techniques allows for a very fine differentiation of strains. In France, 497 different genotypes could be identified between 1978 and 2013, in a collection of M. bovis strains isolated from cattle and wild animals [2,16]. The vast diversity of the molecular profiles obtained by combining these two techniques makes it possible, by comparing the profiles of isolated strains of cattle and wild animals, to determine the origin of the infection of a very large number of outbreaks and to highlight possible inter-species transmission [17]. However, in areas where the disease is still highly prevalent, unique dominant genetic profiles-per-zone are shared by almost all isolates, which makes the reconstitution of a transmission chain impossible [17,18]. Molecular techniques with finer resolutions are necessary in such contexts. Analyses of mycobacterial whole genome sequences (WGS), which detect genomic changes at a very small scale, hold promise as a strain discrimination alternative [19]. This technique has particularly been used to better understand the human-to-human transmission of different mycobacteria, such as M. leprae [20] or M. tuberculosis [21]. WGS data have also been used to reconstruct transmission scenarios of M. bovis infection [22,23].
The aim of this study was to use WGS analysis in order to investigate the transmission network of M. bovis with a particular genotype, mainly identified in infected red deer belonging to French farms, zoos, or hunting or other type of parks.

Ethical Statement
In France, bTB is a notifiable disease with an eradication program in cattle and surveillance in free range wildlife. The official methods for the diagnosis of this disease are culture, PCR and histopathology. Therefore, all the samples included in this study are issued from animals analyzed within an official context. No purposeful killing of animals was performed for this study. All samplings were in complete agreement with national and European regulations. No ethical approval was necessary.

Bacterial Strains
This study included 24 field M. bovis isolates belonging to the French bTB National Reference Laboratory (NRL, Anses) collection ( Table 1). The source population was animals whose samples were analyzed because (i) they presented non-negative skin tests and/or γ-interferon tests results and were culled for diagnostic purposes; (ii) they were culled as a result of the total or partial slaughter of their herd of origin after confirmation of the herd's infection, or (iii) they presented macroscopic bTB-like lesions at routine abattoir inspection. Bacterial culture was performed following the protocol established by the French NRL (NF U 47-104) for the isolation of M. bovis. The identity of the isolated Mycobacterium tuberculosis complex colonies was confirmed by DNA amplification as described by Hénault et al. [24], and M. bovis was confirmed by spoligotyping and VNTR typing as described below. For livestock breakdowns, a herd-based criterion was used for strain choice: at least 1 strain/outbreak was selected. Spoligotyping was performed as described by Zhang et al. [14], using TB-SPOL kits purchased from Beamedex ® (Beamedex ® SAS, Orsay, France) on BioPLex200/Luminex 200 ® as described by Hauer et al. [2]. Spoligotypes have been named according to an agreed international convention (www.mbovis.org). A Multi-Locus Variable numbers tandem repeats Analysis (MLVA) profile identification [25] was performed by Genoscreen (Lille, France), using PCR amplification targeting genetic loci including mycobacterial MIRU-VNTR. Analysis was based on 8 loci, ETR A (VNTR2165), ETR B (VNTR2461), ETR C (VNTR577), ETR D (MIRU4 or VNTR580), QUB 11a (VNTR2163a), QUB 11b (VNTR2163b), QUB 26 (VNTR4052) and QUB 3232 (VNTR3232), chosen in the framework of a European consortium based on their degree of polymorphism and their ability to discriminate local strains [26].

Whole Genome Sequencing and Variant Calling
Whole genome paired end (2 × 250 bp) sequencing was performed on Rapid Run HiSeq 2500 (Illumina) using NextEra XT libraries (Illumina, San Diego, CA, United States), and indexed following the manufacturer's recommendations (Illumina, San Diego, CA, United States) by Genoscreen (Lille, France). The theoretical coverage was expected to be higher than 50×. Raw sequences were aligned to the reference genome AF2122/97 (Genbank accession NC_002945.4) using BioNumerics v7.6 (Applied Maths). A wgSNP analysis was performed with the strict SNP filtering template, i.e., removing positions with at least one unreliable base, ambiguous base, or gap; removing nondiscriminatory positions; with a minimum of 10× coverage; with a minimum distance of 12 nucleotides between SNPs. The evolutionary history of this dataset was inferred on MEGA7 software [27], using the Maximum Likelihood method based on the General Time Reversible (GTR) model (model selected based on jmodeltest v2.1.10 [28]). A bootstrap analysis was performed with 500 replicates.
WGS data were submitted to the EBI/EMBL. The accession Numbers are listed in Table 1.
These isolates were analyzed by WGS. In total, we identified 770 SNPs compared to the M. bovis AF2122/97 reference genome (supplementary Table S1). The phylogenetic analysis highlighted different clusters, all enrooted by the wild boar strain isolated in 1979 (Figures 1 and 2). A distance matrix of SNPs was generated using MEGA (Table 2) and the number of SNPs between each strain was determined. A first cluster supported by a bootstrap value of 100% is defined by 47 SNPs. This cluster is sub-divided into two sub-clusters: the first one, which is phylogenetically well supported (cluster 1.1, 98%), containing strains isolated from the cattle or deer of livestock herds; the second one, which is less robust (cluster 1.2, 67%), contains strains isolated from zoo animals or belonging to hunting parks ( Figure 2). A low number of SNPs defined these sub-clusters-four and one respectively. Cluster 2 groups all the strains with spoligotype F70/SB0295, isolated from one French cow and from foreign-born animals, and is defined by 24 SNPs. Cluster 3, defined by 51 SNPs, contains strains isolated from foreign-born animals with spoligotype GB54/SB0121.

Discussion and Conclusions
This investigation, based on genomic data, allowed us to reveal the circulation of closely related M. bovis strains that seem to be confined to deer farms and animal parks. Although these types of strains have been regularly found since the 1990s, no geographical link seems to exist between most herd cases. WGS analysis showed a main cluster in France, divided into two very close sub-clusters (with just five SNP differences)-cluster 1.1 related to cattle and deer herds reared as livestock, and cluster 1.2 related to recreational or hunting animal parks. This analysis confirms the infection by the same strain of a bison and a wapiti from the same park. It also shows a link between the strain isolated recently in a fallow deer and two wild boar strains in 2007 and 2012 [29]. An epidemiological link exists between these three cases: the wild boars belonged to the same hunting park and the fallow deer, found infected but in another animal park, was originally issued from the same hunting park as that of the wild boar cases. Altogether, the common point in all these French cases was the presence of deer in the parks/herds or close to cattle herds. This genotype has been circulating in France for the last 30 years, however such strains have been isolated sporadically which clearly shows that there is a lack of appropriate surveillance in animal parks. Given that primates are not maintenance hosts of M. bovis, the case described in the monkey-acting as a sentinel of the infection-could be the result of the under-detection of another wild ungulate case in the park. The low number of screening tests, such as the skin test, conducted on deer is the result of several factors, such as the difficulty to perform the test on these wild animals, the potential cross-reaction with environmental mycobacteria, or the inherent technical variability of the test [30].
Our data highlight the limit of classical genotyping techniques that can be overcome by the use of whole genome SNP analysis. Strains that shared the same genotype (spoligotype-VNTR) from French and foreign-born animals were genetically different and very distant. The 16 French strains gathered in cluster 1 seem to be closer (between 0 and 46 SNPs) than the strains of cluster 2 and 3. Indeed, some of them have very few genetic differences. Two strains isolated from a bison and a wapiti from the same animal park are identical. The WGS analysis allows us to confirm a link between two strains isolated from one wild boar in 2007, another in 2012 and a strain isolated from a fallow deer in 2017. The total number of SNPs for cluster 1 is 107, 50 for cluster 1.1, and 57 for cluster 1.2. On the other hand, the number of differences in the remaining clusters is higher with fewer strains, i.e., 163 SNPs for cluster 2 and 199 SNPs for cluster 3. Table 2. Estimates of Evolutionary Divergence between Sequences. The number of base differences between sequences are shown.

Discussion and Conclusions
This investigation, based on genomic data, allowed us to reveal the circulation of closely related M. bovis strains that seem to be confined to deer farms and animal parks. Although these types of strains have been regularly found since the 1990s, no geographical link seems to exist between most herd cases. WGS analysis showed a main cluster in France, divided into two very close sub-clusters (with just five SNP differences)-cluster 1.1 related to cattle and deer herds reared as livestock, and cluster 1.2 related to recreational or hunting animal parks. This analysis confirms the infection by the same strain of a bison and a wapiti from the same park. It also shows a link between the strain isolated recently in a fallow deer and two wild boar strains in 2007 and 2012 [29]. An epidemiological link exists between these three cases: the wild boars belonged to the same hunting park and the fallow deer, found infected but in another animal park, was originally issued from the same hunting park as that of the wild boar cases. Altogether, the common point in all these French cases was the presence of deer in the parks/herds or close to cattle herds. This genotype has been circulating in France for the last 30 years, however such strains have been isolated sporadically which clearly shows that there is a lack of appropriate surveillance in animal parks. Given that primates are not maintenance hosts of M. bovis, the case described in the monkey-acting as a sentinel of the infection-could be the result of the under-detection of another wild ungulate case in the park. The low number of screening tests, such as the skin test, conducted on deer is the result of several factors, such as the difficulty to perform the test on these wild animals, the potential cross-reaction with environmental mycobacteria, or the inherent technical variability of the test [30].
Our data highlight the limit of classical genotyping techniques that can be overcome by the use of whole genome SNP analysis. Strains that shared the same genotype (spoligotype-VNTR) from French and foreign-born animals were genetically different and very distant.
The use of WGS revealed an unsuspected transmission network between hunting parks and zoos that needs to be further explored. These results reinforce the suspicion of an infectious deer-breeding network that could be at the origin of the bacillus dissemination. Indeed, this transmission could have occurred through inter-breeding animal exchanges. Our results raise awareness about the infectious risk of wild animals in captivity and the consequences of trading non-appropriately-tested animals. P22 ELISA, a currently available, highly specific and sensitive new test, could contribute to a more accurate diagnosis of TB in this species [30].
This study raises the problem of the zoonotic transmission of the disease, especially for cases in animal parks where close contact between free-roaming animals and the public is allowed. The fact that deer often develop generalized TB, and thus that they probably excrete the bacillus at significant rates, increases the risk to humans. This study therefore shows the need for the strengthened monitoring of these animals with regular ante-mortem tests (as for livestock) before allowing animal movement between parks.