In this study, we downloaded a total of 700 Mycobacterium bovis whole-genomes from the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA). These isolates originated from three previous studies focusing on M. bovis evolutionary dynamics and cross-species transmission in NZ (PRJNA363037; Crispell et al., 2017), UK—Woodchester Park (PRJNA523164; Crispell et al., 2019), and in the USA—Michigan (PRJNA251692; Salvador et al., 2019), respectively. Metadata associated with these isolates included geographic location and host species from which M. bovis samples were extracted (11 in total): six in NZ [stoat, Mustela erminea; porcine, Sus scrofa; cervine (various NZ cervine species unknown); cattle, Bos taurus; brushtail possum, Trichosurus vulpecula; and ferret, Mustela furo], two in the UK (cattle, Bos taurus and the Eurasian badger, Meles meles), and three in the USA (white-tailed deer, Odocoileus virginianus; elk, Cervus canadensis; and cattle, Bos taurus). Although there was only one M. bovis isolate from a stoat host, it was kept in the dataset to assist in further analyses of geographic and subpopulation structure. Bioinformatic statistics associated with the isolates are shown in Supplementary Table 1.

To improve the quality from the paired-end sequences, we used fastp v0.22 (Chen et al., 2018) to remove adapter sequences and low-quality ends. We used a sliding window approach (size 4 bp) to trim reads with window average quality below 30. Additionally, we filtered out all reads less than 15 bp. Once all the reads were processed, we mapped them against the M. bovis reference genome AF2122/97 (NC_002945.4, Genbank accession code PRJNA89) using the Burrows-Wheeler Aligner (BWA) software (Li and Durbin, 2009). We removed duplicated reads using Picard v2.0.1 (Broadinstitute/Picard, 2021/2014) to limit the impact of non-unique or erroneous reads on the SNP calling process. We performed Variant calling using Freebayes v1.3.5 (Garrison and Marth, 2012), and we kept detected SNPs for the downstream analysis if they possessed (i) a phred-quality score (QUAL) above 20, (ii) a mapping quality (MQ) above 59, (iii) did not fall within gene regions that coded for pe, pe/PGRS and ppe genes (a family of redundant sequence genes that are difficult to work with in-silico; Sampson, 2011; Delogu et al., 2017), and (iv) were more than 500 bp from insertion or deletion mutation (INDEL) regions (this is an extra measure to reduce impacts in the alignment that can lead to false positive detections). Furthermore, the isolates used in this study were sequenced either on a HiSeq or MiseqIllumina platform, and we compared specific bioinformatic metrics (read depth, number of SNPs, SNP depth, and bp depth) between isolates to check if the choice of sequencing platform played a role in SNP detection and quality. A SNP alignment was created using BCFtools v1.13 (Li, 2011) ‘consensus’ command and snp-sites v2.5.1 (Page et al., 2016) to integrate the detected SNPs into a copy of the M. bovis reference genome.

Evolutionary relationships among M. bovis isolates were investigated using the Maximum Likelihood software IQtree v2.1.4 (Minh et al., 2020) with the M. bovis SNP alignment used as input. IQTree uses the ModelFinder approach to determine the nucleotide substitution model that best fits the data based on the Bayesian Information Criteria (BIC) comparison between substitution models (Kalyaanamoorthy et al., 2017). After inferring the substitution model, IQtree ran 1,000 bootstrap iterations in order to provide a measure of internal node support ranging from 0 to 100 (Hoang et al., 2018).

To determine if the existence of potential regions with elevated densities of base substitutions in the M. bovis genomes has any effect on the resolution of their phylogeny and nodal support, we compared a phylogeny produced with the original SNP alignment data with the one generated by the Gubbins software (Croucher et al., 2015). Each phylogeny was assessed based on the proportion of highly supported internal nodes (nodal support). After choosing the phylogeny with the best nodal support, additional comparisons between phylogenies were made based on different rooting strategies such as using the M. bovis reference (AF2122/97, NC_002945.4, Genbank accession code PRJNA89) as an outgroup, Minimal Ancestor Deviation (Tria et al., 2017), and Midpoint rooting (Kinene et al., 2016). The best rooting strategy was based on comparisons of Robinson Foulds distance between the different rooted phylogenies, which is a popular metric to ascertain the number of operations required to convert the topology of one tree to another (Pattengale et al., 2007). If this value was large, then the two trees being compared are highly dissimilar. After assessing which phylogeny was more representative of the evolutionary history between isolates, that phylogeny was used to simulate an M. bovis alignment that will be used to generate null distributions of SNP density and selective sweep regions (sections below).

In order to determine M. bovis population structure, we used fastBAPS v1.0.4 (Tonkin-Hill et al., 2019), a model-based clustering approach, to determine clusters of related sequences present in the genetic alignment. This software utilizes a Bayesian hierarchical clustering approach that handles large alignments and quickly discerns sub-populations. As input, fastBAPS received the SNP only alignment and utilized the optimized symmetric prior with the assumption that all allele frequencies at a SNP site are equivalent (making it analogous to a uniform-like prior). The sub-population clusters were extracted directly from the output of fastBAPS.

To investigate if certain regions in the M. bovis genome have significantly higher values of either SNP density or selective sweep sites, a probabilistic hypothesis test was used to find highly significant genomic regions. For that, we produced a simulated alignment from the original alignment with the tool Alisim using the same alignment length, nucleotide substitution model, and phylogenetic tree topology (created from the previously computed Maximum Likelihood tree; Ly-Trong et al., 2022). Alisim next created a random sequence based on the previous specifications and then simulated nucleotide substitutions that were independently added while also conforming to the phylogenetic topology, resulting in a simulated alignment output. From the simulated alignment, we performed two separate sliding window analyses for calculating SNP dense regions and Selective Sweep regions, respectively. We captured the number of SNPs for the SNP density analysis and the support for selective sweep events in the selective sweep analysis. The data generated from the two separate sliding window analyses were then used to create two null distributions that were used to perform hypothesis testing for the SNP dense and selective sweep regions in the original data. Significance was determined by performing a hypothesis test on each window, where the probability of calculating a certain metric in the sliding window from the original data was compared to the probability of seeing the same metric in the null distribution (Figure 1A).

To find SNP dense regions in the genome alignment, we implemented a SNP counting approach that used a sliding window to find regions with significantly higher SNP counts than what would be expected by chance (Figure 1B). A combination of sliding window size (w) and step size (s) values, where w = {100, 500, 1,000, 2,500, 5,000}, s = {50, 100, 300, 500, 1,000}, and w > s, was used to determine the best combination of parameters to the highest number of SNPs in windows that were found to be significant (Tajima, 1991). Once the optimal window and step lengths were inferred, we used those lengths for each region in the genome to identify high SNP counts. Significance was determined by performing a hypothesis test on each window, where the probability of finding a certain number of SNPs in a window was based on the simulated alignment null distribution. The higher the number of SNPs present in a window, the lower the probability that this window evolved due to conventional means. Our implemented approach was simpler than other tools that are meant to find regions with high SNP density. For example, pre-existing software tools, such as Gubbins, detect regions with high SNP counts by first allowing a sliding window to take lengths between 0.1 to 10 kb to identify regions with at least 10 SNPs. Gubbins methodology also incorporates additional hypothesis tests to find both regions with high SNP counts and putative regions of homologous recombination. Our implementation of high SNP regions ensured that we were extracting all of the genomic regions with high SNP counts and not only regions that had high SNP counts and putative support for homologous recombination (which would make us miss some of the high SNP counts regions).

The null distribution should have had characteristics of a Poisson distribution since we were counting the number of SNPs in a certain window size. Therefore, Poisson distribution parameters for this null distribution were estimated using the R package fitdistrplus based on a maximum likelihood method (R Core Team, 2022; Delignette-Muller and Dutang, 2015). For each window that was analyzed, we performed a statistical hypothesis test to determine the probability that the number of SNPs in a window came from the null distribution. Since multiple tests were required for the alignment, to reduce the number of false-positive regions, we used a Bonferroni corrected p-value (Dunn, 1961). If the probability of the number of SNPs based on the null distribution was lower than the Bonferroni corrected p-value {p-value/(the number of tests)}, then that window was recorded as having a significantly higher number of SNPs than expected by chance. The combination of w and s that produced the highest number of significant results was used for subsequent analyses. We merged the overlapping windows to highlight regions that had high numbers of SNPs throughout using BEDtools (Quinlan and Hall, 2010). These discrete regions were then labelled as SNP dense regions (SDRs).

To identify the sites that were the center of the selective sweep regions, we used the tool OmegaPlus v3.03 (Alachiotis et al., 2012) to determine the value omega, which summarizes the extent of linkage disequilibrium (LD) upstream and downstream of a pre-determined genomic position in the M. bovis genome alignment (Figure 1C). The higher the omega value, the higher the possibility that the genomic position is within the proximity of a selective sweep. We computed the omega values for 5,000 equidistant positions within the simulated alignment. Similar to the protocol employed for detecting SNP dense regions, we first characterized the null distribution of omega values as a gamma distribution since the values started at zero, were continuous, and theoretically could have been unbounded. After filtering out the positions that had the value of zero, we again used the fitdistrplus package to estimate the gamma distribution parameters based on the maximum likelihood method. For each window that was analyzed, we performed a statistical hypothesis test with Bonferroni corrected p-values to determine the probability that the evidence for a selective sweep coming from that region was from the null distribution. To determine the appropriate maximum window size, we computed the number of significant results (n) using multiple window sizes w = {150, 300, 500, 750, 1,000, 1,500, 3,000, 5,000, and 10,000}, and divided the number of significant results by w. Since the calculation of the omega statistic consistently increases as w increases, dividing the number of successful regions by w provided a way to find window sizes that did not add extra significant results to the analysis. The maximum value of w was determined by finding the window size that maximized the calculation (number of significant windows)/(window size length). We merged the overlapping windows to highlight genomic regions that possessed high evidence for selective sweep events using BEDtools. These discrete regions were then labelled as selective sweep regions (SSRs).

After finding the genomic regions with statistical support for increased SNP occurrence and selective sweeps, we measured the impact that evolution in these regions has on being able to differentiate isolates based on their ecological grouping. Random Forest models supported by the random Forest package, were used to perform classification of isolates based on evolutionary metrics to elucidate if SDRs and/or SSRs were important individual predictors in correctly classifying the M. bovis data (Figure 2; Breiman, 2001). Specifically, we used the following data as predictors in our models: (i) SDR number of SNPs, (ii) SDR number of INDELs, and (iii) SSR number of SNPs. To measure how the genomic regions we found compare to general characteristics of the M. bovis genomes, we also included other genomic metrics of the isolates in our analysis. GC content was calculated directly from the full M. bovis alignment, whereas number of non-coding SNPs, number of coding SNPs, number of non-coding INDELs, and number of coding INDELs were extracted directly from each isolate’s variant calling format (VCF) file. Prior to model fitting, predictors with a high amount of correlation with other predictors were removed in order to limit redundant data using the R caret package (see Supplementary File 3; Kuhn, 2008). Model performance was calculated based on the overall accuracy in predicting the correct M. bovis isolate membership for particular groupings. A predictor’s importance was measured through its mean decrease in accuracy (MDA), with higher values indicating that models which did not incorporate the predictor suffered by the indicated number of accuracy points. For each random forest analysis, we ordered predictors based on their MDA and recorded the top 20 most important ones. Certain analyses called for predictors to be compared across different models, such as comparing predictor importance across ecological scales, global versus local cattle hosts, and wildlife reservoir versus livestock hosts. For these analyses, Venn diagrams were constructed to observe which predictors were shared or unique amongst the models. Additionally, predictors that saw sharp increases or decreases in importance when compared across models were recorded.

After key genomic regions were identified from each random forest model, we used BEDtools to identify sets of genes that were present within important SDRs and SSRs. In particular, we highlight genes that have been catalogued in previous research as playing a role in host pathogen interactions in mycobacterial pathogens.

Whole-genome sequencing of the 700 M. bovis isolates from studies implemented in NZ, UK, and the USA identified a total of 6,847 SNPs (Supplementary Table 1). Even though the genomes were sequenced on different platforms, in the 692 samples that were sequenced on Illumina HiSeq 2500 (197) and Illumina MiSeq (495), the number of SNPs, depth per base pair, and average depth per genome did not deviate much across platforms (Supplementary Figure 1). Additionally, the number of SNPs that were detected in coding versus non-coding regions of the genome shared similar distributions between all isolates sequenced on the differing platforms (Supplementary Figure 2), suggesting that the choice of sequencing technology did not influence the downstream results. Isolates sequenced on a NextSeq platform (8) were not included in the violin plots constructed for Supplementary Figures S1, S2 due to an insufficient number of isolates characterized with this platform.

When comparing the phylogenetic trees created from the original alignment and the one outputted by Gubbins, we found that there was a Robinson Foulds distance of 592 between the two trees, meaning that one tree would need 592 changes to be converted to the other. The Gubbins phylogeny possessed better nodal support (75.7% of internal nodes with bootstrap value of 75 or above) than the phylogeny produced from the original alignment (70.1% of internal nodes with bootstrap value of 75 or above; Supplementary Table 2), which suggests that it provides a better approximation of the evolutionary history between the M. bovis isolates in this dataset, and therefore, was used both for the rooting method comparison and for simulating the M. bovis alignment. The best rooting strategy for the phylogeny was the Outgroup rooting method [using the reference AF2122/97 (NC002945.4)], since it led to fewer changes in the topology (12 and 14) when compared to the Midpoint and MAD rooting strategies (Supplementary Figure 3, Supplementary Table 3). The substitution model that best fit the data was the TPM2+F+R4. The M. bovis phylogeny produced clustering patterns that were distinguishable based on geographic region (Figure 3). The population clustering analysis identified eight distinct clusters that matched closely with geographical locations. For instance, Cluster 1 had the same isolates as the USA group. Clusters 6 and 7 were subpopulations of the UK group, and NZ isolates contained the remaining population clusters (Clusters 2, 3, 4, 5, and 8). The population clusters did not show a visible pattern based on host species, indicating that the estimated M. bovis population clusters are determined mostly by geographical location.

The simulated M. bovis alignment was produced using the TPM2+F+R4 model, which was also the best supported model. For the SNP Density analysis, we determined that the combination of a window size of 1,000 bp and a step size of 50 bp were the parameters that identified the most significant windows (Figure 4A). In each window, we counted the number of SNPs that were within that window and represented that data as a Poisson distribution. We estimated the Poisson distribution parameter lambda to be 2.84, which suggests that in our simulated alignment, we can expect to have 2.84 SNPs in a window on average with a variance of 2.84 (Figure 4B). For the selective sweep analysis, we identified that a window size of 5,000 bp was needed to maximize the number of significant results (Figure 4C), and similarly we estimated the parameters of the Gamma null distribution to be 0.944 (the shape parameter estimates the typical omega value in a window), and 0.202 (the scale parameter estimate that variance of the omega value), respectively. Based on the inferred parameter values for both null distributions, the genome-wide hypothesis tests computed the probability of finding the amount of SNPs/evidence of a selective sweep within the tested window. After we merged the overlapping significant windows, we determined 14 unique SDRs and 132 unique SSRs (Supplementary Tables 4, 5). A few SDRs and SSRs overlapped regions with regions containging pe or ppe genes, but since the SNPs falling within these genes were removed from the analysis, these genes did not contribute to the identification of SDRs or SSRs.