Introduction
============

Bovine tuberculosis is a significant worldwide disease of cattle ([@evx022-B86]). *Mycobacterium bovis*, the causative agent, is also pathogenic for humans and several domestic and wild animals ([@evx022-B87]). *Mycobacterium bovis* belongs to the *Mycobacterium tuberculosis* complex (MTbC), which consists traditionally of *M. tuberculosis*, the main agent of human tuberculosis, plus lineages typical of other vertebrate hosts, *M. africanum*, *M. canettii*, *M. microti*, *M. bovis*, *M. caprae*, *M. pinnipedii*, *M. mungi*, *M. orygis*, and *M. suricattae* ([@evx022-B26]; [@evx022-B79]; [@evx022-B3]; [@evx022-B88]; [@evx022-B70]). All of these species are genomically very similar, with more than 99.9% nucleotide identity (Thoen, LoBue, et al. 2006). It has been hypothesized that MTbC members have evolved from a common ancestor (possibly *M. canettii*) via successive DNA deletions/insertions resulting in the present *Mycobacterium* speciation and associated differences in pathogenicity ([@evx022-B26]; [@evx022-B36]).

Disease control programmes based on regular tuberculin testing, and removal of infected animals ("test-and-slaughter") have been successful in eradicating, or markedly reducing, tuberculosis (TB) from cattle in many countries, but these measures are less effective in countries with wildlife reservoirs of *M. bovis* ([@evx022-B17]; [@evx022-B61]). The gold standard for postmortem diagnosis of bovine TB has been mycobacterial culture, although this method requires a significant amount of time (up to 90 days) ([@evx022-B56]). This fact, together with economic globalization and the pressure from importing markets for a definitive diagnosis of tuberculosis in cattle that exhibit lesions compatible with tuberculosis and the advances in molecular biology, have provided a stimulus for improved TB molecular diagnostic techniques ([@evx022-B85]; [@evx022-B5]).

In order to examine the evolution of *M. bovis* infection in space and time, it is essential to study genetic differentiation of isolates and their phylogenetic relationships, together with cattle movement and wildlife ecology. A clearer understanding of *M. bovis* epidemiology may contribute to better control or to eradication of disease. High-quality epidemiological data increases the efficiency of bovine tuberculosis control programs. With new tools available in molecular epidemiology, there is an increasing possibility of finding answers to issues such as the importance of transmission between individual animals and its risk factors, and the role of wild animals as reservoirs ([@evx022-B41]).

The use of whole genome sequencing (WGS) technology to study infectious bacterial diseases has resulted in unprecedented advances in the ability to resolve epidemiologic data at the global scale ([@evx022-B39]). WGS has provided new insights into within-host replication processes ([@evx022-B30]), and has been used to corroborate the importance of exhaustively identified transmission chains and social drivers of transmission ([@evx022-B33]). WGS of many strains of the same species can provide the necessary statistical power required to find association between genotypes and phenotypes such as virulence, pathogenesis, host range, and antigenic variability.

In order to understand evolutionary processes of bovine tuberculosis, we undertook a large phylogenomic comparative study, addressing aspects of whole genome evolution of *M. bovis*. Our study includes 103 genomes ([supplementary tables S1 and S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, for details), isolated from 14 host species (cattle, pig, sheep, human, bobcat, deer, raccoon, opossum, coyote, jaguar, red deer/elk cross, fallow deer, elk, red deer) in eight countries (Brazil, Argentina, Uruguay, the United States, the United Kingdom, Korea, Mali, and Uganda), and obtained within a 22-year interval (1991--2013). Our goal was to obtain a thorough understanding of genomic characteristics, large-scale gains/losses/duplications, phylogenetic relationships, evolutionary timeframe, positive selection in core-coding genes, and homologous recombination rates. We also highlight important aspects that may be associated with evolution of lineages from the United States, such as patterns of reinfection among host species of *M. bovis*.

We further assess the evolutionary contribution of a specific large gene family typical of MTbC lineages, the PE/PPE family, known to be associated with evolution of these lineages in general, but whose contribution under different evolutionary forces using WGS had still not been contemplated in the particular case of *M. bovis* lineages. The PE/PPE family of genes is large, being distributed across the genome of MTbC, reaching ∼170 copies per genome in *M. tuberculosis* H37Rv ([@evx022-B15]). It is known that many genes within it are associated with antigenic variation and immune evasion, being secreted by the type VII (or ESX) secretion system typical of mycobacteria and some nonmycobacterial actinomycetes, being composed of five copies of the ESX gene cluster, ESX-1 through ESX-5, the former being the most extensively studied ([@evx022-B1]; [@evx022-B8]; [@evx022-B77]; [@evx022-B35]). PE genes have the Pro-Glu motif at positions 8--9, with a conserved stretch of ∼110 amino acids in the N-terminal domain. PPEs embrace a Pro-Pro-Glu sequence at positions 7--9 with conservation of ∼180 AAs in the N-terminal domain. The C-terminal regions of both families are highly variable, also containing repetitive DNA motifs that may vary in copy number among family members ([@evx022-B15]). PE and PPE can be further subdivided by characteristics of their N-terminal regions and phylogenetic relationships ([@evx022-B89]; [@evx022-B59]). PE genes are composed of five subfamilies, the largest being the PE-PGRS (polymorphic GC-rich-repetitive sequence) embracing 65 members in H37Rv, with multiple tandem repeats of Gly-Gly-Ala or Gly-Gly-Asn as typical motifs in the C-terminal region ([@evx022-B59]). PPE genes are also subdivided into five subfamilies, of which the largest sets are of PPE-SVP (for a total of 24 members in H37Rv) and the PPE-MPTR (major polymorphic tandem repeat) subfamilies (with 23 members).

Materials and Methods
=====================

Chromosomal Variations
----------------------

Sibelia ([@evx022-B64]) was used to assess chromosomal structural variations across the five *M. bovis* complete genomes (excluding BCG genomes, as these do not represent natural genomic variation) present in NCBI's RefSeq database ([@evx022-B74]), with comparisons based on raw fasta files to avoid annotation biases. The minimum threshold for detecting large sequence polymorphisms (LSPs) was 2,000 bp segments (this value was chosen to be high enough to uncover several interesting LSPs, but not too low to yield too many less informative polymorphisms).

Data Sets for Phylogenomics
---------------------------

We used two different data sets. Most of the evolutionary analyses were done with a 38-genome data set, which included a total of eight countries ([supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, for details); unless stated otherwise, our analyses refer to this 38-genome data set. We used *M. tuberculosis* H37Rv (a canonical strain from this species) as outgroup. We also present results from a second data set ([supplementary table S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, for details), which had a narrower aim of gaining insights about *M. bovis* evolution within the United States only. This was made possible by the availability of 76 US *M. bovis* genomes from the PATRIC database ([@evx022-B94]).

Regarding the 38-data set, we obtained all non-US genomes from GenBank ([supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, for details). The 11 US genomes within this data set were assembled by the vcf-consensus tool within VCFtools ([@evx022-B19]) using *M. bovis* AF2122/97 (NC_002945.3) as reference, and annotated with the Genome Reverse Compiler ([@evx022-B93]); these genomes can be found at \<<http://pintado.facom.ufms.br/mbovispub/>; last accessed October 3, 2016\>. There is an overlap of those 11 genomes between the 38-genome data set and the 76-genome data set. In the case of the 76-genome data set, we obtained assemblies and annotations for these 11 genomes from PATRIC ([@evx022-B94]), as was the case for all genomes in this data set. We chose not to use our own assemblies and annotations in the analysis of the 76-genome data set in order to keep methods for this data set uniform. In order to check whether differences in assembly/annotation pipelines would introduce systematic biases in the analyses for the 38-data set, we performed *t*-tests of US versus non-US genomes in both the 38 and 76-data set (the latter set with annotations performed by the same tool, RAST, therefore with no annotation bias among genomes) considering different genomic parameters (core-coding size in base pairs, GC-content, and different codon bias measures), with significant differences (*P* ≤ 0.05) between US and non-US genomes being found for all parameters in both data sets (the *t*-tests for the 38-data set specifically are further detailed in Material & Methods, and Results). Therefore, we conclude that the 38-data set was not influenced by heterogeneity in assembly/annotation.

Core- and Pan-Genome Curves
---------------------------

To assess whether the set of genomes used here is representative of the total genetic pool of *M. bovis*, we estimated pan- and core-genome curves ([@evx022-B84]). This was done for the 38-data set using the get_homologues.pl script ([@evx022-B16]) with the "--c" option and 20 resampling steps. At each resampling, genomes were included in an order different from the previous resampling step. Exponential regression lines were then fitted to each curve by the same script.

Inference of Orthologs
----------------------

Inference of the families of homologs was done in the same way for both data sets by Get_Homologues v2.0 ([@evx022-B16]), using the get_homologues.pl script with "-M" flag, which activates OrthoMCL clustering ([@evx022-B55]) with default values. Subsequently, strictly orthologous, unicopy core-genome genes (the *core-coding data set*) were estimated by the same software with the script compare_clusters.pl.

Multiple Alignment
------------------

Individual core-coding DNA gene alignments were obtained in Guidance v2.0 ([@evx022-B78]) anchored by codons using MAFFT v7.0 ([@evx022-B47]) with 100 pseudoreplicate guide trees. Concatenation of the multiple alignment was done by FasConCat ([@evx022-B49]). Manual curation was then performed in Aliview ([@evx022-B53]), and whenever local regions were misaligned we used the software "realign selected block" option and updated the information about each gene start and end. Amount of variable and parsimony-informative sites was obtained in Mega 7.0 ([@evx022-B51]). The fraction of nonsynonymous single nucleotide polymorphisms (nsSNPs), synonymous SNPs (sSNPs), and a reduced matrix containing only 4-fold degenerate sites were also obtained in Mega.

Genomic Characteristics of Core-Coding Data Set
-----------------------------------------------

Dambe v5.5.29 ([@evx022-B95]) was used to look for any patterns of variation across the 38-genome set. For each genome, we obtained values for the effective number of codons (Nc), relative synonymous codon usage (RSCU), codon adaptation index (CAI), GC-content, and total size of core-coding alignment. RSCU (which is an individual measure of the variability of use of each set of codons of an amino acid) per genome was obtained by summing the individual amino acid variances and then extracting the square-root of the result.

Phylogenomic Inference
----------------------

Phylogenetic analysis was done by maximum likelihood (ML) within IQTree v1.3.12 ([@evx022-B69]), with support for each phylogeny calculated from 1,000 UFBoot pseudoreplicates ([@evx022-B63]). An important thing to consider when estimating character-based phylogenies (such as ML-based) is the amount of variation among data partitions, otherwise the resulting topology may be biased ([@evx022-B10]). To test for best partitioning scheme for core-coding genes, we entered the three codon positions as data blocks (hence each block encompasses one-third of the alignment), and then estimated the best partitioning scheme among them and respective model for each partition under IQTree, using BIC as the criterion for data fit. To test for variability of rates along the genome, besides using the typical gamma distribution (+G) and a probability of invariable sites (+I) for each proposed partition, we also tested models with mixture distributions for each site, from two to five such mixtures per model (+MM), which in many cases have better fit ([@evx022-B90]).

Further phylogenetic analyses aimed to test for (and possibly alleviate) systematic biases. As a proxy for phylogenetic reliability of each tree, the retention index (RI; [@evx022-B27]) was used to assess the fit of the character to "sampling locality" (continents used as tip states) across the 38-genome data set, following previous conclusions showing that regionally *M. bovis* samples tend to cluster together (e.g., [@evx022-B38]; [@evx022-B4]; [@evx022-B40]). We ascribed one among five states to each tip according to its respective continental region (South America, North America, Europe, Asia, and Africa). Subsequently, RI was calculated for each of the trees, with higher values indicating better agreement of that tree with geographic location. PAUP v4b10 ([@evx022-B82]) was used to calculate RI across the different topologies, using the command "describetrees".

Maximum parsimony (MP) in Mega v7.0 was used to test for topological differences due to the method employed, with support measured by 1,000 standard bootstrap pseudoreplicates ([@evx022-B28]). Besides testing for phylogenetic method, we aimed to correct for possible data artifacts that could lead to topological biases.

To test if selection, recombination, or the use of nonneutral sites were impacting the RI, reduced data sets from core-coding genes were obtained. To account for positive selection, either all genes under positive selection (see below) were removed, or were maintained but using a codon-based ML analysis to estimate a nonsynonymous to synonymous rate ratio (dN/dS) across a given topology.

In order to account for homologous recombination, genomic blocks shared among all taxa were used when considering the whole core genome (therefore also potentially spanning noncoding genes and nongenic regions) and excluding those blocks with even mild signs of recombination (see below).

To limit the analysis to an enriched set of hypothesized nearly-neutral codons (therefore maximizing true vertical signal), a reduced data set including only codons with 4-fold degenerate positions were obtained (in MEGA).

Aiming at minimizing biases due to conflicting phylogenetic signal among sites (regardless of the source of conflict), Tiger v1.0.2 ([@evx022-B18]) was used, which allocates sites sharing similar relative rates within bins (we used ten bins as default). The largest bin set was used as input for a phylogenetic estimate (termed site-congruence method from here on). Models including ascertainment bias correction (+ASC) were also included for the latter phylogenomic analysis, as it lacked invariable sites ([@evx022-B54]).

Finally, a separate analysis independent of multiple alignment (or alignment-free) was performed by kSNP 3.0 ([@evx022-B32]). Its algorithm works by finding core and noncore SNPs that are present across a subset, or the whole of the taxon set, therefore being very inclusive for assessment of tree topology. This SNP data set was run by ML under IQTree with the same options as described previously.

For the 76-genome data set phylogenomic analysis, we included only the variable characters, performing model-selection, tree search, and assessing confidence as explained above.

Assembly Quality
----------------

After performing selection and recombination analyses, we checked if signals could be the outcome of low quality genome assemblies by remapping paired-end reads of the genome with largest number of contigs (Uganda B2-7505, with 243 contigs) to representative genes indicated to be under these forces in genome AF2122/97. We used Trimmomatic (to filter high-quality reads, dismissing reads that were not paired-ends; [@evx022-B6]), Bowtie 2 (mapping; [@evx022-B52]), Samtools (to order and index .bam files; [@evx022-B37]), and VCFTools (to call and summarize SNVs; [@evx022-B19]).

Detecting Homologous Recombination
----------------------------------

The ten largest locally colinear blocks (LCBs) across the 38 genomes were obtained by the ProgressiveMauve algorithm ([@evx022-B20]). The three tests available in PHIPACK (PHI, NSS, and MaxChi) were used to detect recombination in each of these blocks ([@evx022-B11]). To be conservative, an LCB was considered to be free of detectable recombination if none of three tests were significant under *P *=* *0.1. Blocks with undetected recombination were concatenated and estimated by ML in IQTree as described above. No recombination tests were performed gene-wise due to the scant variability among taxa in each gene alignment (which share \>99% identity on average).

As another way to assess the effect of recombination, ClonalFrameML ([@evx022-B22]) was run by concatenating all LCBs irrespective of recombination signs (the tree with best RI in the phylogenetic estimation phase was fixed). ClonalFrameML estimates relative contributions between homologous recombination and point mutations across the LCBs and for every branch of the tree.

Positive Selection
------------------

The core-coding data set was used to test for positive selection. Due to low sequence variability among *M. bovis* genomes, we had to estimate positive selection from the concatenation of core-coding genes in order to bring more power to the test. The FEL (fixed effects likelihood) method in HyPhy was able to operate on the whole concatenated core-coding data set. FEL estimates site-wise ratios of dN/dS using the entire data for inference of global alignment parameters (such as transition/transversion rates, nucleotide frequencies, etc.). The 38-genome data set used is larger than the minimum of 20 lineages recommended by the method developers ([@evx022-B48]), making the test appropriate for this data set. Three different rate classes for dN and another three for dS were assumed. The best RI phylogeny was constant throughout the analysis.

Clusters of orthologous groups of proteins (COGs; [@evx022-B31]) were inferred using the script CDD2COG.pl v0.1 (available from <https://github.com/aleimba/bac-genomics-scripts/tree/master/cdd2cog>; last accessed October 3, 2016). To detect which COG categories were significantly enriched in positively selected genes, the relative ratio of positively selected genes within a COG to the number of genes in the whole genome within the same COG category was calculated. A Monte Carlo resampling (bootstrap) of the ratios was performed using 1,000 pseudoreplications using an in-house R script ([@evx022-B83]) with a one-tailed *P*-value.

Dating Divergences
------------------

RelTime ([@evx022-B60]) was used to estimate divergence times. It is a likelihood-based algorithm that does not incorporate explicitly rate nor tree priors, whereas providing divergence estimates comparable to the widely used Bayesian software BEAST ([@evx022-B24]; not used here due to lack of convergence among runs having equal priors and data). We used as calibrations (within conservative minima and maxima) two different age estimates for the origin of *M. bovis*: a relatively older interval based on Comas et al. (2013), inspired by the rationale that expansion of *M. bovis* may have been associated with out-of-Africa human expansion (with a conservative upper bound of 88,000 ya obtained by the authors); and a more recent one, based on [@evx022-B7], employing radiocarbon dating of mummies as tip calibrations (with a conservative interval of 800--7,000 ya as outcome; the same minimum was used for Comas' dating prior in our analysis). We estimated divergences using two data sets, the 4-fold degenerate data set and the site-congruence data set, both minimizing biases due to sites under reasonably different rates (due to forces such as selection). The best RI tree was fixed throughout the runs. To assess variance in date estimates, 1,000 bootstraps (generated under RaxML v8.2.4; [@evx022-B81]) for each data set were analyzed separately using the RelTime version within the MEGA-CC suite ([@evx022-B50]). Dating outliers for each of those four run combinations (Bos et al. vs. Comas et al. calibrations; each with either 4-fold or site-congruence data sets) were discarded by calculating the 99% highest density intervals (HDI) on each case using the R library "coda" ([@evx022-B72]). Node date intervals were calculated as \[min(minima), max(maxima)\] across these HDIs, therefore being conservative as it embraced altogether Comas et al. (2013) and [@evx022-B7] calibrations.

Results
=======

General *Mycobacterium bovis* Evolution
---------------------------------------

A total of 6,189 families of homologous genes were obtained. Of these families, 2,550 were inferred as single copy core-coding genes and subsequently multiply aligned (nucleotide sequences) by Guidance (using Mafft), 67 of them being discarded because of alignment issues and/or containing sequences that were not multiple of three (which may be indicative of at least one copy in the alignment being a pseudogene, or the effect of a wrong base call/bad assembly), resulting in 2,483 core-coding groups. The final manually curated multiple alignment of core-coding genes had 2,491,902 bp, with 6,707 variable sites (0.2%), and 1,761 parsimony-informative sites (0.07%). Regarding amino acids, 3,847 variable sites (0.46%), and 1,134 parsimony-informative sites (0.14%) were detected. Pan/core genomic curves are shown in [figure 1](#evx022-F1){ref-type="fig"}.F[ig]{.smallcaps}. 1.---Core- (*left*) and Pan-genome (*right*) curves including only genomes of the 38-data set having \<100 contigs (33 genomes). The core-genome is close to a plateau, indicating the complete core set is included in the 38-data set analysis (below 3,000 genes). The pan-genome curve indicates the pool of genes is quite high even within relatively closely related *M. bovis* genomes. Equations refer to the exponential regression for each curve.

RIs using either core-coding or reduced data sets were lower than the RI for the site-congruence method topology (RI = 0.87), even after changing the phylogenetic method (from ML to MP), or correcting for selection, recombination, nonneutral evolution ([fig. 2](#evx022-F2){ref-type="fig"}; [supplementary fig. S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, for details), or by including both core and non-core SNPs ([supplementary fig. S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, for details). All US genomes clustered together in the site-congruence topology; however, one Uruguayan genome was included within the clade ([fig. 2](#evx022-F2){ref-type="fig"}). Given its better RI, this tree was chosen for downstream analyses.F[ig]{.smallcaps}. 2.---ML tree of the 38-data set matrix obtained using the site-congruence method. *Left*: branch support values (UFBoot); *Right*: same tree with actual branch lengths according to ML. RI for this tree (retention index for the character "sampling locality") = 0.81.

Significant differences (though small in an absolute sense) were observed in genomic characteristics such as codon usage (Nc, RSCU, and CAI), GC-content, and amount of nucleotides in the core-coding data set between US and non-US genomes (all *t*-tests with *P*-value \< 1E-18) ([fig. 3](#evx022-F3){ref-type="fig"}), agreeing with the pattern observed for the best RI tree in identifying a US lineage. Curiously, no synapomorphies were observed for the US clade, either in core-coding, or in the kSNP data sets (which includes core and noncore regions) ([supplementary fig. S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, for details).F[ig]{.smallcaps}. 3.---Variability of genomic characteristics across the 38 core-coding genes obtained in DAMBE. Each index was normalized to a relative average of 1.0 for better visualization. Numbers in *x* axis refer to individual genomes, with the 11 US genomes starting at position \#28. Nc, effective number of codons; RSCU, relative synonymous codon usage; CAI2, corrected version of the codon adaptation index; GC-content, average genomic GC content; Sum(ACGT), size of the core-coding data set in base pairs.

Contrived assessment of read mapping quality (important for testing for spurious outcome in phylogenetic, recombination, and selection results) was done using the 19 PPEs genes in the core-coding, a reasonable choice as PPEs can have more than 100 copies in MTbC genomes and there are also repetitions even within a single PPE paralogue. All single nucleotide variations (SNVs) agreed 100% to SNPs observed in each of the 19 orthologues' multiple alignments (each including the 38 taxa), indicating that phylogenetic inference and genes detected under selection and recombination were not biased to any reasonable extent due to low quality assemblies.

We assessed in ClonalFrameML the amount of recombination in relation to mutations across the branches using the multiple alignment of the ten largest LCBs concatenated (550,634 bp). Rho/theta (ratio of the recombination rate to the mutation rate) and r/m (the ratio of the estimated quantity of polymorphism brought by an event of recombination in relation to an event of mutation) were: Rho/theta = 0.1, and r/m = 0.98, with an average recombination tract length of 67.7 bp. Branch individual values were in general below 1.0 for both ratios except for ten nodes which had Rho/theta ≥ 1.0 and/or r/m ≥ 1.0, seven of them including the US + Uruguay lineage (the branch subtending the clade, the two branches leading to the CO + TX clade, and the other four leading to the MI + Uruguay clade plus all branches within), one subtending the Korean strains, and the other two at more basal positions (the branch leading to the nonAfrican *M. bovis*, and the one leading to all *M. bovis* except the genome from Uganda). The branch leading to the US clade (including the Uruguayan genome) had a homologous recombination import of a PPE family protein segment of 244 bp from Korean strain 1595 (as shown by Blastn (Altschul et al. 1997)). Three different PE-PGRS (a highly polymorphic PE family member; [@evx022-B76]) genes had segments recombined in the UK BAA-935 genome, with best Blastn hits belonging to *M. bovis* genomes. The US genome 12_5092_S3 from Texas had an import of 244 bp of a PPE family protein gene also with best hits belonging to *M. bovis*. A recombinant PPE gene in the UK genome AN5 had a 169 bp import having best hit with *M. canettii*, a species that is probably in a sister group to MTbC ([@evx022-B26]; [@evx022-B36]). The best hit of a segment of the 23S ribosomal RNA from the Uruguayan sample (139 bp) was a distant mycobacterium with high pathogenicity (*M. abscessus*).

A total of 114 genes were inferred to be under positive selection; of these, 102 were categorized into COGs ([table 1](#evx022-T1){ref-type="table"}). Monte-Carlo bootstrapping was conducted to test for COGs significantly enriched under selection. "Cell cycle control, cell division, chromosome partitioning" was the most significant category ([table 1](#evx022-T1){ref-type="table"}), including two copies of an *esx* component (ESX-1 secretion system protein). "Cell motility" was the second most significant; only three genes are within this category, all of them of the PE/PPE family ([table 1](#evx022-T1){ref-type="table"}; [supplementary table S3](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, for details). "Replication, recombination and repair," the third most significant category, had already been described as an enriched category within MTbC ([@evx022-B23]). Another two significant categories were found: "Lipid transport and metabolism" (14 genes), and "Posttranslational modification, protein turnover, chaperones" (the latter being borderline significant). Important genes under selection (other than PE/PPE and ESX-associated products) included cholesterol metabolism-associated products (e.g., all eight copies of acyl-CoA synthetases), CRISPR type III-associated RAMP protein Csm3, and transposition genes, among others. Genes under positive selection in each category can be found in [supplementary table S3](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, for details.Table 1COG Categories Enriched for Genes under Positive SelectionCOG Abbr.COG CategoryAllSel% All% Sel(%Sel)/(%All)ARNA processing and modification100.00040.00000.00BChromatin structure and dynamics000.00000.0000---CEnergy production and conversion14760.06440.05410.84DCell cycle control, cell division, chromosome partitioning3260.01400.05413.85\*\*EAmino acid transport and metabolism16560.07230.05410.75FNucleotide transport and metabolism6940.03020.03601.19GCarbohydrate transport and metabolism10140.04430.03600.81HCoenzyme transport and metabolism12360.05390.05411.00ILipid transport and metabolism188140.08240.12611.53\*\*JTranslation, ribosomal structure and biogenesis11550.05040.04500.89KTranscription16790.07320.08111.11LReplication, recombination and repair95100.04160.09012.16\*\*MCell wall/membrane/envelope biogenesis8940.03900.03600.92NCell motility2230.00960.02702.80\*\*OPosttranslational modification, protein turnover, chaperones10070.04380.06311.44\*PInorganic ion transport and metabolism9850.04290.04501.05QSecondary metabolites biosynthesis, transport and catabolism12950.05650.04500.80RGeneral function prediction only33480.14640.07210.49SFunction unknown18350.08020.04500.56TSignal transduction mechanisms7620.03330.01800.54UIntracellular trafficking, secretion, and vesicular transport2110.00920.00900.98VDefense mechanisms2710.01180.00900.76WExtracellular structures000.00000.0000---YNuclear structure000.00000.0000---ZCytoskeleton000.00000.0000---*Total =*22821111.001.00[^2][^3]

Large-Scale Polymorphisms
-------------------------

We analyzed five complete *M. bovis* genomes occurring in nature available in GenBank (therefore excluding BCG genomes): two from the United Kingdom, one from Brazil, one from the United States, and one from South Korea. We detected 12 large sequence polymorphisms (LSPs) ≥ 2,000 bp, including nine indels, two duplications involving two nonadjacent segments (one of the duplications was probably originally in tandem; see below), and one inversion ([table 2](#evx022-T2){ref-type="table"}). Only LSPs 1 to 4 are newly identified here, the remaining overlapping completely or to a great extent regions previously named as RD3, RDcap_Asia2, RD2, RD5oryx, RDbovis(c)\_fadD18, and RD1 ([@evx022-B9]; [@evx022-B65]). Among the three newly identified indel blocks, two of them include two copies of PE-PGRS family genes (LSP-1 and LSP-4). Further details can be found in [table 2](#evx022-T2){ref-type="table"} and [supplementary table S4](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, for details.Table 2Polymorphic Genomic Blocks ≥2,000 bp Identified across the Five *Mycobacterium bovis* NCBI RefseqsBlockAF212297 (UK)BAA-935 (UK)Strain 1595 (Korea)SP38 (Brazil)Strain 30 (USA)StartEndStartEndStartEndStartEndStartEndLSP-1926157928398---92605292829341763554178596926097928278LSP-21036390105484910298101048269103644410549034286983430544210365261054984LSP-3------1413645141778531647132061114138661418005RD317646451774379b1769178177891367266568240017692031778938RDcap_Asia219885531997632---1993162200224189631590539419930902002169RD221996412210799a220460622157641106529111768722041662215322RD5oryx26050042607181---261005026122271511708151388526094522611629LSP-427623572772053271695427256792767457277715316690191678696cRep1---35971473614668---------Rep2---36146703633236---------Rep1---36332433650816---------Rep2---36508193669385---------RDbovis(c)\_fadD1838860503890304390967539132793891927389617227919082795181cRD142865864296110---429281543023303191235320075942773244286846[^4][^5][^6][^7]

Origin and Divergence of US Strains
-----------------------------------

Individual date intervals were obtained as follows ("US clade" includes the Uruguayan genome). Origin of US clade (4-fold data set): \[167.27, 5209.91\]; Diversification of US clade (4-fold data set): \[172.73, 5209.91\]; Origin of US clade (site-congruence data set): \[149.7, 3492.2\]; Diversification of US clade (site-congruence data set): \[148.7, 3425.1\].

Discussion
==========

Phylogenomics
-------------

The phylogenetic analyses indicated that there may be substantial amounts of conflicting signal in the core-coding set of genes. The conflicting signal remains even after trying a different analytical method (maximum parsimony instead of maximum likelihood), correcting for positive selection (either by removal of genes under selection, or by maintaining them but under a global dN/dS ratio phylogeny), correcting for homologous recombination (by using only the largest collinear blocks without signs of recombination), focusing solely on 4-fold degenerate codons, or else considering the core + noncore regions. Nevertheless, the site-conguence method (an algorithm that is not based on phylogeny) was helpful in pointing at a set of sites evolving at similar rates, and these showed a better RI fit under ML, identifying a clade containing all US samples; moreover, the latter also agreed with US genomic features (codon usage bias, GC-content, and length of core-coding set) being significantly different in comparison to the other genomes. Even though the source of conflicting signal was not identified, a possible hint is the observed length of terminal branches ([fig. 2](#evx022-F2){ref-type="fig"}; [supplementary fig. S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, for details): some of them are relatively long, suggesting that most of the mutations affected external instead of internal branches. A probable cause for the lack of synapomorphies is recent population expansion after a bottleneck, followed by selective sweeps (as suggested by [@evx022-B79]); whenever a selective sweep happens, it tends to hamper variability in the rest of the genome, potentially erasing synapomorphies that led to the evolution of clades, and possibly decreasing the correlation with geographic location.

Impact of Selection and Recombination
-------------------------------------

Homologous recombination is in general assumed to be uncommon within MTbC lineages (except for *M. canettii*) when compared with other organisms ([@evx022-B79]). Mycobacteria lack plasmids in general, and also bear cell wall peculiarities ([@evx022-B73]), which in theory restricts homologous recombination. Our results indicate relatively high amounts of homologous recombination instead. This agrees with a previous study focusing on 24 genomes of different species of MTbC ([@evx022-B67]). The r/m value of 0.98 shows that the magnitude of sites impacted by recombination are on average as much as mutation. Furthermore, it agrees with what has been found recently regarding mycobacterial DNA exchange: individuals may intercross by distributive conjugal transfer, a process in which there can be (multifragment or large-chunk) transfers and internalization by homologous recombination ([@evx022-B92]; [@evx022-B34]). Though deemed rare, when it happens many stretches of the donor DNA may get incorporated in the chromosome at once ([@evx022-B21]).

Positive selection was an important evolutionary force within *M. bovis*, impacting 114 genes according to our analyses. Previously, a study by [@evx022-B23] had shown that many genes under the COG category "Replication, recombination and repair" were under selection. This is of high significance for the *M. bovis* lineages: being able to accumulate sustainable variability may be important if environmental changes start to select among different alleles. In our study, three other categories were shown to be enriched as well ("Cell cycle control, cell division, chromosome partitioning," "Cell motility," and "Lipid transport and metabolism"; and "Posttranslational modification, protein turnover, chaperones" being marginally nonsignificant). We highlight *esx*-associated genes (important constituents for the secretion of virulence factors; [@evx022-B43]); PPE proteins, which are themselves important virulence factors ([@evx022-B76]; [@evx022-B2]; [@evx022-B59]); eight copies of Acyl-CoA dehydrogenases, reinforcing the importance of cholesterol-associated pathways for survival, infection, and persistence ([@evx022-B62]; [@evx022-B91]); esterase/lipase lipF, shown to be important in pathogenesis ([@evx022-B12]; [@evx022-B96]), again reinforcing the importance of virulence-associated genes for *M. bovis* evolution.

Large-Scale Polymorphisms
-------------------------

Although there have been previous reports of polymorphism regarding repetitive markers, SNP polymorphisms, and structural variations within *M. bovis* lineages, this is the first report of LSPs among the five complete genomes occurring in nature. Among the nine indel regions most of the losses occur in one of the UK genomes (BAA-935), except for two that are missing solely from US strain 30 (LSP-4 and RDbovis(c)\_fadD18; [table 2](#evx022-T2){ref-type="table"}). Furthermore, the recombination analysis suggested that there is larger-than-average level of segment imports (five of them between ∼200 and 520 bp) in the BAA-935 lineage (data not shown). These results show that BAA-935 is relatively differentiated, due to both genic losses and homologous recombination. The other available complete genome from the United Kingdom, AF2122/97, is quite different in these respects (low recombination detected, and only one shared loss, LSP-3). Both types of polymorphism (recombination and indels) involve PE/PPE genes, which are known to be important for virulence ([@evx022-B76]; [@evx022-B2]; [@evx022-B59]).

Regarding the two duplicated blocks in UK BAA-935, a parsimonious scenario for their evolution is, first, a duplication of the region encompassing the three regions (Rep1 + 63 kb intervening segment + Rep2) followed by loss of the 63 kbp segment, and then a tandem duplication of Rep1 + Rep2. Regarding the other seven polymorphic LSPs that have been described elsewhere ([@evx022-B9]; [@evx022-B65]), four of them include PE/PPE genes; one of them, RD1, also contains ESX-1 secretion-associated genes.

PE/PPE Genes as Evolutionary Targets within *M. Bovis*
------------------------------------------------------

PE/PPE genes are known to have evolved from duplications and transpositions from MTbC ancestrals ([@evx022-B89]; [@evx022-B68]), whereas being intrinsically associated with *esx* clusters (ESX 1-5, the latter being the more recent cluster), even though other copies exist outside of such clusters. Such association is both in terms of genomic context (i.e., PE/PPE genes can be found within *esx* clusters) and biochemical function (*esx* proteins are largely responsible for PE/PPE secretion; [@evx022-B43]). In fact, *esx*-associated genes (including ESX clusters 1-5) form the TSS7 secretion system typical of Corynebacteriales ([@evx022-B43]).

Although it is accepted that MTbC lineages are largely clonal, it was found that these genes have large rates of recombination ([@evx022-B57]; [@evx022-B45], [@evx022-B46]; [@evx022-B58]; [@evx022-B71]). This same pattern was found in our study, with PE/PPE genes being significantly associated with the evolution of different lineages. A PPE segment of 244 bp present in the US clade (including the Uruguayan genome) was probably a recombination import from a close relative of Korean strain 1595. There were also three other imports in UK BAA-935, one in the United States 12_5092_S3, and another in UK AN5, the latter with best hit with *M. canettii*, which is probably the outgroup of MTbC ([@evx022-B26]; [@evx022-B36]), suggesting genic exchanges with distant lineages may not be uncommon. Our findings once again undermine the belief that *M. bovis* lineages are strictly clonal, and that recombination may be an important force at generating variation that may be useful for immune evasion and/or virulence. This also agrees with the observation that mycobacterial lineages can interact with one another using a particular conjugation system (distributive conjugal transfer; [@evx022-B92]; [@evx022-B34]) in which large portions of an extraneous genome can enter the cell at once, increasing the likelihood of recombination at different genomic regions.

The three genes under positive selection that were found under the "cell motility" COG category were all PE/PPE genes. This observation parallels the pattern of PE/PPE enrichment regarding selection found in closely related *M. tuberculosis* strains ([@evx022-B71]). Finally, two of them were also found within large scale polymorphism indel regions LSP-1 and LSP-4. Overall, these findings evince that PE/PPE genes are targets of different evolutionary processes within *M. bovis*, once again highlighting the importance of the largest mycobacterial family of genes known to date.

It could be argued that, because there are so many PE/PPE copies per MTbC genome (comprising up to ∼7% of the total genome's coding gene set), such genes could have been picked up by chance; in other words, the signals observed here could be simply the effect of sampling bias. This may still be a possibility concerning the identification of PE/PPE genes in the new LSP regions, but we believe this was not the case concerning selection nor recombination, as the signals were identified per site and per gene, respectively, not mattering how many copies the gene family has per genome.

This was the first time that different PE/PPE genes were detected as important targets of both selection and recombination in *M. bovis* under a phylogenomic analysis considering all core-coding genes, providing valuable evolutionary insights, especially considering that many PE/PPE proteins are highly antigenic and associated with immune evasion. These results corroborate and extend previous observations of high levels of synonymous and nonsynonymous substitutions in MTbC lineages that were based on far fewer MTbC genomes (e.g., [@evx022-B29]) or analyses focused specifically on PE/PPE genes without considering the whole set of core-coding genes (e.g., [@evx022-B59]).

Dating and Hypotheses Regarding the Evolution of the US Lineage
---------------------------------------------------------------

Divergence time of origin and evolution of US strains \[148.7, 5209.91\], which included two contrasting sets of calibrations regarding origin of *M. bovis* (Comas et al. 2013, up to 88,000 ya; and [@evx022-B7], at most 7,000 ya), suggested that this lineage evolved at most ∼5,000 ya. We further acknowledge that dating estimates (ours included) in general may suffer from a bias due to change of overall evolutionary rates through time, where short time mutation rate (recent epidemics) changes to a slower longterm fixation rate (e.g., [@evx022-B42]), but we do not know of any available software that corrects for such bias. The US data set was shown to be phylogenetically informative regarding association with sampling locality in certain states, when considering at least three individuals within the same clade ([fig. 4](#evx022-F4){ref-type="fig"}): Michigan (MI), Hawaii (HI), Minnesota (MN), Colorado (CO; three of five samples), and New York (NY; four of five samples). Particularly in the case of MI and HI, it is known that there are wildlife reservoirs of bovine tuberculosis in nature ([@evx022-B61]), but a previous study did not identify monophyly in either population, nor in others ([@evx022-B44]). This is probably because these authors used a subset of SNPs (350, with 306 being genic), whereas the present study used the whole set provided in the core-coding genes (to the exclusion of gene alignments discarded by the program Guidance) amounting to 6,707 SNPs, a set more than 20× larger, which was able to detect the actual clusterings. Therefore, our data add phylogenetic evidence to the hypothesis that MI and HI wildlife populations work as either natural reservoirs or spill-overs (such as deer, racoon, coyote, red fox and others; [@evx022-B61]). MN and NY individuals, which also supposedly include wildlife populations with some *M. bovis* prevalence ([@evx022-B61]), also clustered together by state. The fact that NY forms a monophyletic cluster also agrees with it being a locality prone for human transmission, because it may suggest some prevalence of the disease there. Between 2001 and 2004, 35 human cases were reported in NY (Centers for Disease Control and Prevention, 2005), but apparently these were related to contaminated cheese imported from Mexico; nevertheless, two cases of *M. bovis* were reported in white-tailed deer there ([@evx022-B61]), reinforcing the hypothesis of local sources of reinfection as obtained in our phylogeny. Cattle as host was the pattern observed in terms of more basal nodes within each clade, but support values are not uniformly high and in general the host type composes paraphyletic clades, suggesting reinfections across different species ([fig. 4](#evx022-F4){ref-type="fig"}). It is also possible to observe that there is a large variation in terms of external branches of some US strains, suggesting that some lineages may evolve faster than others ([fig. 4](#evx022-F4){ref-type="fig"}).F[ig]{.smallcaps}. 4.---ML tree of the 76- data set composed of US strains, plus outgroup, based on their core-coding set of genes. *Left:* chronogram with respective ML support. *Right:* same tree with branches in substitutions/site. *Red:* clades of geographically associated genomes (a paraphyletic clade in the case of MN due to inclusion of a sample from TX).

The fact that most clades are paraphyletic in terms of sampling locality indicate that the spread of bovine tuberculosis is relatively fast in the United States, and some lineages evolving much faster than others suggest that different haplotypes may evolve at faster speed. This may be due to selective sweeps, which would tend to erase previous site synapomorphies hence eroding phylogenetic signal in deeper nodes within *M. bovis*, which may have been the case found here regarding the 38-set phylogeny as well ([fig. 2](#evx022-F2){ref-type="fig"}). Overall, these results indicate the possibility of ongoing and localized outbreaks, which can be even more likely given that many of them happened in the United States in the past 20 years. In fact, the country has been considered of moderate risk concerning control of bovine tuberculosis for some time now ([@evx022-B14]). Such fast spreading across different states may be due to commingling of infected cattle with susceptible wildlife reservoirs, as suggested by the paraphyletic clusters in [figure 4](#evx022-F4){ref-type="fig"}, corroborating previous assertions ([@evx022-B61]).

Concluding Remarks
------------------

This was the first evolutionary study focusing on many *M. bovis* genomes and using the largest amount of genomic information available. The gains in resolution of the reconstructed tree and associated evolutionary processes are interesting compared with the use of smaller sets of SNPs, particularly in the case of more recent evolutionary scenarios, as seems to be the case in US strains. Nevertheless, the phylogeny of the 38-data set core-coding alignment was somewhat confounded by conflicting signal among different sites. It may be the case that adding more individuals from each location the synapomorphies start to show up; in any case, it would also be interesting to compare how more customary markers (such as spoligotyping or VNTR) compare in terms of phylogenetic resolution. Going one step further, it is highly likely that the concomitant use of repetitive regions and very large sets of genes will become the most effective way to understand *M. bovis* evolution and population dynamics in detail, both in shallower and deeper levels. The importance of such finer resolution comes from the fact that it may become possible to trace past outbreaks in detail, to further explore hypotheses of wildlife reservoirs, and possibly hint at previously unknown population sources of disease spreading.

Supplementary Material
======================

[Supplementary data](#sup1){ref-type="supplementary-material"} are available at *Genome Biology and Evolution* online.

Supplementary Material
======================

###### 

Click here for additional data file.

This work was partially supported by the Brazilian Agency CAPES (grant number 3385/2013). J.S.L.P. and J.M. were supported by CAPES fellowships. J.C.S. was supported in part by a grant from CNPq (grant number 304881/2015-5). N.F.A. was supported by grants Fundect TO141/2016, TO007/2015, CNPq 305857/2013-4, 473221/2013-6, and CAPES 3377/2013.

[^1]: **Associate editor:** Bill Martin

[^2]: N[ote]{.smallcaps}.---All, all genes; Sel, genes under positive selection; %, relative frequencies; "Total" represents the amount of genes for which a COG was found, and sum of relative frequencies for each COG (=1.0). The last column shows the proportion between (%Sel)/(%All). Test for COG categories significantly enriched for positive selection was done in R by bootstrapping values from this last column and then checking for significance of each category (1,000 pseudoreplicates; α = 0.05; one-tailed critical value = 1.52, i.e., genes with Sel/All \> 152% are significant).

[^3]: Significant under *P*-value = 0.05; \*borderline nonsignificant.

[^4]: N[ote]{.smallcaps}.---UK genome BAA-935 had most of the losses (seven: LSP-1, LSP-3, RD3, RDcap_Asia2, RD2, RD5oryx and RD1), and is the one bearing the two duplications. US strain 30 had two losses (LSP-4 and RDbovis(c)\_fadD18) and the inversion LSP-2, and genome AF2122/97 (UK) had one loss (LSP-3, shared with BAA-935). Gray shading: segment inversion.

[^5]: Only 417 bp.

[^6]: Only 364 bp.

[^7]: Only 500 bp.
