Phylogenomic analysis of Mycoplasma bovis from Belgian veal, dairy and beef herds

M. bovis is one of the leading causes of respiratory disease and antimicrobial use in cattle. The pathogen is widespread in different cattle industries worldwide, but highest prevalence is found in the veal industry. Knowledge on M. bovis strain distribution over the dairy, beef and veal industries is crucial for the design of effective control and prevention programs, but currently undocumented. Therefore, the present study evaluated the molecular epidemiology and genetic relatedness of M. bovis isolates obtained from Belgian beef, dairy and veal farms, and how these relate to M. bovis strains obtained worldwide. Full genomes of one hundred Belgian M. bovis isolates collected over a 5-year period (2014–2019), obtained from 27 dairy, 38 beef and 29 veal farms, were sequenced by long-read nanopore sequencing. Consensus sequences were used to generate a phylogenetic tree in order to associate genetic clusters with cattle sector, geographical area and year of isolation. The phylogenetic analysis of the Belgian M. bovis isolates resulted in 5 major clusters and 1 outlier. No sector-specific M. bovis clustering was identified. On a world scale, Belgian isolates clustered with Israeli, European and American strains. Different M. bovis clusters circulated for at least 1.5 consecutive years throughout the country, affecting all observed industries. Therefore, the high prevalence in the veal industry is more likely the consequence of frequent purchase from the dairy and beef industry, than that a reservoir of veal specific strains on farm would exist. These results emphasize the importance of biosecurity in M. bovis control and prevention.


Introduction
Mycoplasma bovis (M. bovis) causes mostly pneumonia, arthritis, otitis in calves and mastitis in adult cattle [1,2] resulting in high antimicrobial use (AMU) and enormous economic losses in cattle farming sectors worldwide [2][3][4]. In Belgium, 100% of the veal farms are seropositive for M. bovis [5,6], whereas M. bovis is involved in 33% of acute pneumonia outbreaks in beef and dairy farms [7]. Treatment of M. bovis is frequently unsatisfactory, probably due to a combination of intrinsic and acquired antimicrobial resistance, immuno-evasive properties of the pathogen and failure of the animal to generate an effective immune response [8,9]. Together with the absence of an effective commercially available vaccine, the control of M. bovis is particularly challenging.
A contemporary fear is that the veal sector, currently combining a high AMU and a farm level M. bovis prevalence of 100%, is a reservoir for multi-resistant sectorspecific M. bovis strains [6,10,11]. Currently, there is insufficient knowledge about the epidemiology of circulating M. bovis strains to answer this question. Several epidemiological studies observed clonal emergence and identified dominant lineages of the M. bovis bacterium, based on antimicrobial resistance patterns and different strain typing methods [12][13][14][15]. In the past, different approaches were used to subtype M. bovis strains, including random amplification of polymorphic DNA (RAPD), arbitrarily primed PCR (AP-PCR), amplification fragment length polymorphism (AFLP), pulsed-field gel electrophoresis (PFGE), multiple-locus variable-number tandem repeat (MLVA), and multi-locus sequence typing (MLST) [13,14,[16][17][18]. Unfortunately, results from these typing methods are difficult to compare and only focus on a small fraction of the genomic information, resulting in a limited insight in the genetics and incongruence among studies [17,19,20]. Therefore, whole genome sequencing (WGS) could be a great opportunity, considering its highly discriminative capacity and reproducibility compared to older typing methods [21,22].
Several studies already investigated whether specific M. bovis strains were associated with affected organs, such as udder, respiratory tract or joints [14,23,24], geographical location [23,24] or health status [23,25]. Only one study determined epidemiology based on AP-PCR in three farms from three different husbandry conditions (dairy calf ranch, feedlot and closed beef herd), presuming that management factors could influence the distribution of M. bovis [16]. These husbandry conditions are not comparable with the three main sectors in Europe. In Europe, a lot of short-distance movements of cattle between farms is seen, and the veal industry is an important side market of the dairy and beef industry [26,27]. It is currently not clear whether sector-specific M. bovis strains are present and what their genetic relation is to previously sequenced M. bovis isolates. Therefore, the present study first evaluated the molecular epidemiology and genetic relatedness of M. bovis isolates obtained from Belgian beef, dairy and veal farms. Furthermore, it studied the relationship of these isolates to M. bovis strains from other countries.

Library preparation and MinION long-read sequencing
Quality-checked native M. bovis DNA was immediately used for library preparation using the Rapid Barcoding Sequencing Kit (SQK-RBK004; Oxford Nanopore Technologies (ONT)), following manufacturers' instructions. For each run, ten field strains, one positive control (PG45) and one negative control (sterile broth) were multiplexed (400 ng DNA per sample). A new R9.4.1 Flow cell (ONT) was used for a 48 h sequencing run on Min-ION device (ONT). Raw fast5 read files were collected using MinKnow v.3.6.5.
Overall consensus assembly accuracies were verified by comparing total Single Nucleotide Polymorphisms (SNPs) using the CSI phylogeny package (v1.4, Center for Genomic Epidemiology, Denmark; [33]) as compared to the M. bovis PG45 type strain (NC_014760.1) reference sequences. To validate the use of long-read sequencing, SNPs of ten independent M. bovis PG45 assemblies were compared to those in a single MiSeq experimental dataset. All M. bovis consensus genomes are available for download on the NCBI GenBank database under the BioProject PRJNA639688 and accession numbers (SAMN15246515-SAMN1524662). Sequencing summaries can be found in Additional file 1.

Phylogenetic analysis
Phylogenetic analysis was performed on all newly generated consensus sequences alone or in combination with 250 previously published M. bovis sequences using the FastTree-based CSI Phylogeny v1.4 (see Additional file 2). All analyses included the M. bovis PG45 type strain (NC_014760.1) as reference and M. agalactiae PG2 (NC_009497.1) as outgroup. Resulting Newick files were visualized with MEGA-X software [34].

Cluster and strain determination
Due to the lack of relatedness criteria for SNP typing schemes of M. bovis and the need to establish these per organism and experimental design [35], clusters were defined by visual inspection of the phylogenetic tree and by taking into account bootstrap support. In addition, the matrix of pairwise SNP counts was extracted from CSI Phylogeny for further inspection. Mean SNP differences were calculated between within-cluster isolates, and outliers were defined with the 1.5xIQR rule, using the Outlier Calculator (https ://miniw ebtoo l.com/outli er-calcu lator /).

Geographical distribution
Esri ® ArcMap ™ (version 10.7.1) software was used to visualize the geographical distribution and density of M. bovis isolates over Belgium. Herd size was based on the national Identification and Registration system, containing, on the first of January 2017, a total of 23 995 cattle herds in Belgium (23 733 conventional herds; 262 veal), and a total of 2 517 850 cattle. The spatial distribution of the Belgian cattle (both cattle and veal calves) was displayed using kernel smoothing. Coordinates of Belgian cattle herds were converted into a continuous raster using the kernel density estimation function weighted by number of cattle (Spatial Analyst, ArcMAP X, ESRI, Redlands, CA, USA).

Phylogenetic analysis of Belgian isolates
A median sequencing depth of 618X (range: 32X-2689X) was obtained from long reads with an average N 50 read length of 5706 ± 1514 bp for all Belgian M. bovis isolates. First, implementation of long-read sequencing of M. bovis genomes was validated by comparing total SNPs from ten independently sequenced M. bovis PG45 sequences and a single M. bovis PG45 MiSeq dataset, showing 53 ± 3 SNPs and 27 SNPs difference, respectively, compared with the 1 003 404 bp of the M. bovis PG45 reference genome (NC_014760.1). The observed average SNP difference of 0.005% for the long-read sequencing approach was considered acceptable to allow meaningful phylogenetic analyses. In addition, control strain M. bovis PG45 results were mutually compared over all ten runs, showing a mean SNP difference of 20 (range 8-30, standard deviation 4.6), which was also demonstrated an acceptable inter-experimental variation.
Taking into account all Belgian M. bovis strains, and also including the outgroup M. agalactiae, 51.4% of the M. bovis genome or 515 324 nucleotide positions were used for phylogeny. The minimum and maximum SNP differences among Belgian M. bovis isolates were 33 and 3775, respectively.
Between and within clusters, no association could be observed for the different cattle sectors or year of isolation ( Figure 1). Two different isolates from the same herd (veal) and same sampling period (Mb49 and Mb50) did not cluster together (II and V). All clusters persisted in Belgium for at least 1.5 consecutive years throughout the country. M. bovis strains isolated from the middle ear (n = 3) were clustered within cluster IV, while those obtained from milk, joint and other samples were scattered over different clusters (Figure 1). Finally, no clear association between geographic location of sampled farms and M. bovis clusters was observed, as shown in

Discussion
In this study, one-hundred M. bovis isolates from different Belgian cattle sectors (beef, dairy or veal) were phylogenetically compared to investigate whether sector-specific strains exist and whether such strains are related to M. bovis strains previously isolated and sequenced worldwide.
In this study, we chose to apply the ONT long-read secquencing approach [36], because no default WGS approaches are defined for Mycoplasma spp. and short-read sequencing biases have been described for genomes with highly repetitive regions.
WGS approaches have become more attractive over the last years, as the cost for next-generation sequencing has significantly reduced. Single Nucleotide Polymorphism (SNP) analysis using Illumina short read data from M. bovis isolates already showed to be an effective way for M. bovis genotyping [24,25]. Longread nanopore sequencing (Oxford Nanopore Technologies) is known to create much faster results, and was recently applied in veterinary medicine as well [37,38]. Yet, lower single read accuracies are currently obtained with ONT in comparison to Illumina [39]. Therefore, the implementation of long-read sequencing to generate M. bovis genome assemblies was verified, showing only on average 53 SNPs difference of the long-read approach with the publically available M. bovis PG45 genome, representing an acceptable error rate of 0.005%. As a result, the authors believe that nanopore sequencing is a highly accessible tool, allowing practical use outside academia in routine diagnostics and real time surveillance.
From the study results, several interesting observations were made. First of all, the obtained M. bovis isolates belonged to at least five different M. bovis clusters, of which three dominant clusters were identified. This is in agreement with an Israeli study based on WGS-SNP, where six clusters were observed of which three were dominant. Remarkably, one cluster contained more than 50% of the isolates in that study [25]. Several other studies also showed one or two dominant lineages, although different typing methods were used [13,15,18,24,40]. In contrast to Aebi et al. [41], where mostly herd-specific M. bovis isolates were seen in Switzerland, we observed close relatedness of M. bovis isolates over the different herds. This might be a result of more frequent purchasing cattle from different origins and transportation in Belgium, because 40% of cattle is transported at least once (and up to eight times) over a 5-year lifespan in Belgium [7,16,17,26,42]. As such, the higher heterogeneity observed in cluster I, IV and V compared with clusters II and III, may be caused by different rates of genetic drift between clonal lines [17].
Secondly, no sector-specific strains or clusters were identified in the present study, which does not entirely come as a surprise. In Belgium, veal calves are purchased from both dairy and beef farms, and fattened and slaughtered in specialized veal farms and slaughterhouses, respectively [27]. Also, approximately 15-20% of the farms in Belgium are mixed farms. Therefore, contact among different cattle sectors is intense. In addition, herd visitors and artificial insemination might play a role in spreading of M. bovis or introducing new strains on farms [43][44][45].
Thirdly, when we take a closer look at how the different clusters have spread over Belgium, no clear association with location was observed. This was also concluded in studies performed in the UK and USA [17,23,42]. Nevertheless, in the provinces of Antwerp and Western Flanders seem to be hotspots for M. bovis outbreaks. Besides a high number of local transports, Antwerp and the Flanders area are also the main gates for cattle import, which makes these areas predisposed for the introduction of new M. bovis strains [26]. Unfortunately, M. bovis genomes were not available for isolates obtained from the top import countries for Belgium, which are Germany and the Netherlands.
Although this study was not designed to draw definitive conclusions about year of isolation and affected organs, some preliminary observations can be made. For example, no association between M. bovis strain and year of isolation was observed. On the other hand, we saw that representatives of all clusters persisted for at least 1.5 consecutive years on Belgian territory. The persistence of strains within a country or herd has been described before [16,24,41]. Furthermore, shifts between dominant lineages from older to new strains have been reported before as well [13,15,46,47]. Also, we did not observe an association between cluster and affected organ, which is in line with previous studies [14,23,24]. Yet, it was remarkable that all isolates obtained from the middle ear were clustered, which could suggest the middle ear as possible predilection site for certain M. bovis strains. However no definitive conclusions can be drawn as there were only few isolates obtained from this isolation site in the present study. In addition, we isolated different M. bovis strains (Mb49, Mb50) from two veal calves on the same farm at the same time. The observation of two different strains in one herd or even one animal has been described before [14,16,42,48], in contrast to Arcangioli et al. [46], who isolated only one identical dominant profile in the same feedlot.
Finally, it was evident that Belgian isolates were mostly related to European and Israeli M. bovis isolates, even though only a few genomes of European M. bovis isolates have been published in the NCBI database. This seems plausible, as Belgian farmers mostly purchase cattle from European farms, while Israel also partly imports cattle from Eastern-Europe. The fact that Israeli isolates are often related to Chinese and Australian strains, is also due to import of cattle, as outlined in detail elsewhere [25,40]. Some of the Belgian outlier strains were related to American isolates, which might be explained by the fact that M. bovis was first isolated in the USA and outbreaks in Europe were only seen years later. So, we can only speculate whether these outlier strains could have been imported or evolved geographically distinct from each other. Clusters of the Belgian isolates were not clustered exactly in the same way in the Belgian vs. the worldwide phylogenetic tree. A possible explanation could be the loss of overall coverage between the construction of both phylogenetic trees (51% for the Belgian to 40% worldwide). This might be due to (1) heterogeneity among isolates worldwide and/or (2) the use of genomes obtained by different laboratories, using different sequencing protocols, as the quality can be influenced by strain maintenance, DNA extraction, library preparation, sequencing, and the bioinformatics analysis [49].
In conclusion, multiple M. bovis clusters were circulating in Belgium in 2014-2019, and were persisting for several years. Neither the veal industry, nor any other cattle industry could be identified as source of strain persistence. Connections between dairy, beef and veal industry are intense and M. bovis appears to easily spread among these sectors. The M. bovis issues in the veal industry seem more likely the consequences of strain import from dairy and beef, rather than persistence of a limited number of veal specific strains. This information can contribute to better control and prevention of M. bovis infections by improved biosecurity.