Mycobacteria of the Mycobacterium tuberculosis complex (MTBC) greatly affect humans and animals worldwide. The life cycle of mycobacteria is complex and the mechanisms resulting in pathogen infection and survival in host cells are not fully understood. Recently, comparative genomics analyses have provided new insights into the evolution and adaptation of the MTBC to survive inside the host. However, most of this information has been obtained using M. tuberculosis but not other members of the MTBC such as M. bovis and M. caprae. In this study, the genome of three M. bovis (MB1, MB3, MB4) and one M. caprae (MB2) field isolates with different lesion score, prevalence and host distribution phenotypes were sequenced. Genome sequence information was used for whole-genome and protein-targeted comparative genomics analysis with the aim of finding correlates with phenotypic variation with potential implications for tuberculosis (TB) disease risk assessment and control. At the whole-genome level the results of the first comparative genomics study of field isolates of M. bovis including M. caprae showed that as previously reported for M. tuberculosis, sequential chromosomal nucleotide substitutions were the main driver of the M. bovis genome evolution. The phylogenetic analysis provided a strong support for the M. bovis/M. caprae clade, but supported M. caprae as a separate species. The comparison of the MB1 and MB4 isolates revealed differences in genome sequence, including gene families that are important for bacterial infection and transmission, thus highlighting differences with functional implications between isolates otherwise classified with the same spoligotype. Strategic protein-targeted analysis using the ESX or type VII secretion system, proteins linking stress response with lipid metabolism, host T cell epitopes of mycobacteria, antigens and peptidoglycan assembly protein identified new genetic markers and candidate vaccine antigens that warrant further study to develop tools to evaluate risks for TB disease caused by M. bovis/M. caprae and for TB control in humans and animals.
Mycobacterium tuberculosis has infected more than 2. 5 billion people worldwide with approximately 9 million new tuberculosis (TB) cases reported every year [1]. Animal TB is caused by infection with Mycobacterium bovis and closely related members of the M. tuberculosis complex (MTBC) such as M. caprae. Although cattle are the main concern regarding animal TB in industrialized countries, several other mammals including humans are also infected [2,3]. Eurasian wild boar (Sus scrofa) are a natural reservoir for M. bovis in some regions and thus vaccination strategies are being developed for TB control in this species [4–8]. Several other domestic and wild animals are also infected with M. bovis and may act as reservoir species [9–13]. The life cycle of mycobacteria is complex and the mechanisms resulting in pathogen infection and survival in host cells are not fully understood [14]. Nevertheless, it is generally accepted that after inhalation into the lung or entry to the oropharyngeal cavity, the principal entry routes, mycobacteria of the MTBC are phagocytized by macrophages, which constitute the main host cell. As with other intracellular bacteria, mycobacteria survive inside macrophages by escaping host immune response, which results in the formation of a granuloma that effectively contains infected cells. A change in the host-bacterial equilibrium of granulomas is thought to result in the release of infected cells outside containment and onward transmission of mycobacteria to susceptible hosts [14]. The association between M. bovis spoligotypes and TB lesions in cattle has been used to correlate bacterial genotype with virulence [15]. However, these genotyping methods cover only a small portion of the approximately 4,000 genes contained in the 4. 4 Mb mycobacterial genome [16]. Recently, whole-genome sequencing and comparative genomics analyses have provided new insights into the evolution and adaptation of the MBTC to survive inside the host and explained phenotypic traits related with transmissibility and virulence [16–21]. Although the first M. tuberculosis genome sequence was reported in 1998 [16], it was not until 2003 when the first M. bovis genome was sequenced [19]. Presently, a large number of M. tuberculosis but few M. bovis (except for BCG strains) genome sequences are available [18]. The relative paucity of M. bovis genome sequence information limits the possibility of characterizing mycobacterial evolution and correlation with virulence at the whole-genome level. In this study, the genome of three M. bovis (MB1, MB3, MB4) and one M. caprae (MB2) field isolates with different lesion score, prevalence and host distribution phenotypes were sequenced. Genome sequence information was used for whole-genome and protein-targeted comparative genomics analysis with the aim of finding correlates of phenotypic variation with potential implications for TB disease risk assessment and control.
All animal sampling was post-mortem. Wildlife samples came from hunter-harvested individuals that were shot during the legal hunting season independently and prior to our research while livestock samples were obtained at the slaughterhouse where they were being processed as part of the normal work and submitted to the reference laboratory by the slaughterhouse veterinarian. According to EU and National legislation (2010/63/UE Directive and Spanish Royal Decree 53/2013) and to the University of Castilla–La Mancha guidelines no permission or consent is required for conducting this type of study. Field isolates used come from the EU Reference Laboratory for Animal Tuberculosis (VISAVET). Three M. bovis (MB1, MB3, MB4) and one M. caprae (MB2) field isolates were selected for this study (Table 1). These isolates were originally obtained from wild boar (MB3, MB4), cattle (MB1) and goat (MB2). The study focused on Ciudad Real Province, Spain. This is a high ungulate density area, the west side of the province composed by interspersed game ranges and protected nature areas, with persistent TB infection in extensive livestock farms [22]. Nine hundred MTBC isolates collected from wild ungulates and livestock from 2000 to 2011 were spoligotyped resulting in 62 different spoligotypes ([23] and S1 Fig). The criteria for selection of the MTBC spoligotypes included in the study were based on: Following these three criteria, the M. bovis spoligotypes SB0339 (MB1, MB4) and SB0134 (MB3) were selected. The M. caprae spoligotype SB0157 (MB2) was included in the study as an outgroup but closely related species [29] and for the increased proportion of M. caprae isolated from bovine samples during 2004–2009 [30]. The MB4 isolate was included because although it has the same spoligotype as MB1, it served as a model to characterize possible differences between MTBC isolates otherwise classified with the same spoligotype. The four isolates were grown in 15 ml of Middlebrook 7H9 liquid media supplemented with 0. 36% sodium pyruvate and 10% OADC (Oleic Albumin Dextrose Catalase) for 5 weeks. Chromosomal DNA samples were obtained as described by van Soolingen et al. [31]. Briefly, cultures were centrifuged and pellets were washed twice in 5 ml water. Mycobacteria were heated at 100°C for 15 min to kill the cells. After centrifugation, the cells were resuspended in 5 ml TE buffer (0. 01 M Tris-HCl, 0. 001 M EDTA, pH 8. 0). Lysozyme was added to a final concentration of 1 mg/ml and the tube was incubated over night at 37°C. Eighty hundred and seventy-five microliters of 10% sodium dodecyl sulfate (SDS) with 62. 5 μl of proteinase K (at a 10-mg/ml concentration) were added, and the mixture was incubated for 1 h at 60°C. The extract was transferred to a phase lock gel tube (prime5, Fisher Scientific SL, Madrid, Spain) for a phenol/chloroform DNA extraction. Genomic DNA (2–5 μg) was subjected to mechanical fragmentation using a BioRuptor (Life Technologies, Carlsbad, CA, USA). The number of cycles was adjusted to obtain DNA fragments of a final average size of about 500 pb. Samples were then used to prepare sequencing-amenable TruSeq libraries (NEB-Next, New England Biolabs, Ipswich, MA, USA). Briefly, DNA fragments were made blunt-ended, phosphorylated, adenylated and Illumina-compatible adapters were ligated. After purification, barcoded sequences as well as Illumina-specific sequences were introduced by PCR, followed by quantitation of individual libraries. Libraries were then pooled and quantified again. A quality control of the pooled library made in bioanalyzer is shown in the figure, including an estimation of the percentage of non-overlapping reads that could be obtained using a 2x250 paired-end sequencing protocol. Library was qPCR-quantitated and brought to a final concentration of 10 nM. DNA was then denatured and equilibrated so that a final concentration of 18 pM of library was loaded onto a MiSeq v. 3 flowcell (Illumina, San Diego, CA, USA) and sequenced using a 2x250 paired-end sequencing protocol to obtain more than 400x high quality coverage (1. 9–2. 5 Gb) with 84% of the bases showing a Q30 factor > 30. Reads were finally split according to barcodes and used for bioinformatics analysis. High quality overlapping reads were merged using FLASH (Magoc et al. , 2011) and then assembled using Velvet [32] with k-value = 97 (S1 Table). Contigs were annotated using BG7 [33,34] (S2 Table). For annotation, a set of 191,017 reference proteins was used including (a) all Uniprot proteins from M. bovis and M. tuberculosis, (b) a set of bacterial antibiotic resistance related Uniprot proteins selected using the GO annotation terms “antibiotic resistance” and “response to antibiotic” and a selection of proteins based on similarity to the proteins of ARDB [35], and (c) all Uniprot proteins with Enzyme Code (EC) from MTBC. For whole genome comparative analysis, the 4 genomes were then aligned against the M. bovis reference genome sequence (AF2122/97; http: //www. ncbi. nlm. nih. gov/nuccore/31791177) using Differences program (Turrientes et al. , 2010) that allows comparisons at the whole genome level and particularly the detection of substitutions, insertions or deletions of any length and at any region of the genome. Genome sequence information and annotation was deposited in GenBank under the accession numbers CDHF01000001-CDHF01000049, CDHG01000001-CDHG01000059, CDHH01000001-CDHH01000094 and CDHE01000001-CDHE01000118 for MB1-MB4 isolates, respectively [36]. A phylogenetic tree was constructed based on the SNPs found at the core genome sequence shared with a similarity over a threshold between all genomes included in the analysis using Harvest [37] and visualized with EvolView (v 198. 3) [38]. Harvest defines the SNPs aligning whole assembled genomes and not reads like in other approaches, thus allowing the identification of the gene or the intergenic region where differences are allocated. The SNPs are included in the. vcf file provided by Harvest suite. The genomes used in the SNP phylogenetic analysis using Harvest included MB1-MB4, M. bovis AF2122/97 (http: //www. ncbi. nlm. nih. gov/nuccore/31791177), M. bovis ATCC BAA-935 (http: //www. ncbi. nlm. nih. gov/nuccore/690294709), M. bovis BCG Pasteur 1173P2 (http: //www. ncbi. nlm. nih. gov/nuccore/121635883), and M. tuberculosis H37Rv (http: //www. ncbi. nlm. nih. gov/nuccore/448814763). The multilocus sequence analysis was conducted using 18 genes coding for the proteins linking stress response with lipid metabolism (S3 Table). The nucleotide sequences of the genes were obtained from the genomes of MB1-MB4 isolates used in this study. For comparison, the same sequences were obtained from the reference M. bovis BCG Pasteur 1173P2, M. bovis AF2122/97 and M. tuberculosis H37Rv. M. canettii (NCBI reference sequence NC_019950) was used as outgroup. The nucleotide sequences were concatenated and then aligned with MAFFT (v7), configured for the highest accuracy [39]. After alignment, regions with gaps were removed and 20080 gap-free sites were used in maximum likelihood phylogenetic analysis as implemented in PhyML (v3. 0 aLRT) [40,41]. The reliability for the internal branches was assessed using the approximate likelihood ratio test (aLRT–SH-Like) [41]. Graphical representation and editing of the phylogenetic trees were performed with EvolView (v 198. 3) [38]. The alignments obtained by MAFFT were used to perform codon alignments using HIV database website (www. hiv. lanl. gov; [42]). Non-synonymous (dN) and synonymous (dS) nucleotide substitutions were classified using the SNAP method [43] implemented in HIV database website [42]. SNPs were identified by pairwise comparison of MB1-MB4 and M. bovis AF2122/97 using SNAP. The dN/dS ratio was calculated for the 18 genes coding for the proteins linking stress response with lipid metabolism (S3 Table) and for 81 antigen-coding genes (S4 Table) present on each of the MB1-MB4 isolates included in the study. Under the Datamonkey server (http: //www. datamonkey. org; [44]), the algorithm SLAC [45] was used to detect which nucleotide substitution site were positively or negatively selected. For each dS and non-synonymous dN substitution site, four measurements were made: normalized expected (ES and EN) and observed numbers (NS and NN). The SLAC algorithm then calculated: dN = NN/EN and dS = NS/ES. If dN < dS a codon was negatively selected and if dN > dS a codon was positively selected. A two-tailed extended binomial distribution set at P<0. 05 was used to assess significance of the algorithm. The SLAC algorithm uses a neighbor-joining tree with a maximum likelihood for branch lengths and substitution rates. To identify non-synonymous mutations that may be associated with virulence and/or transmission, the orthologs of peptidoglycan assembly protein locus Rv0050 (H37Rv) in strains MB1-MB4 were compared to the equivalent locus in animal isolates (BCG Pasteur 1173P2, AF2122/97,09–1192) and human isolates (Bz 3115, B2 7505) of M. bovis available in the GenBank. The M. bovis strains were selected among all strains for having distinct Rv0050 locus and/or distinct source of isolation. The presence of a putative signal peptide and cleavage site in the Rv0050 locus was analyzed with a previously validated program for their prediction in M. tuberculosis (SignalP; http: //www. cbs. dtu. dk/services/SignalP-3. 0/) [46]. Results were also confirmed with two other programs, Signal-Blast (http: //sigpep. services. came. sbg. ac. at/signalblast. html) and Phobious (http: //phobius. sbc. su. se/). To confirm selected SNPs identified in the mycobacteria genomes sequenced in this study, sequence-specific oligonucleotide primers were design using reference genomes M. bovis AF2122/97 (BX248333. 1) or M. tuberculosis H37Rv (AL123456. 3) for PCR and sequencing of the amplicons. Selected loci and direct and reverse primers used for analysis included Rv0050 (ponA1; 5´-GACTTTCCCCAAACCGACCGAGG-3’ and 5’-GATCGGTCCCCCGACCACCATT-3’), Rv0589 (MCE2a; 5´-GTGCCAACGCTGGTGACGAG-3’ and 5’-AGAACACGATCAACCCATGA-3’), Rv1198 (ESAT-6; 5´-ATGACCATCAACTATCAATT-3’ and 5’-TCGGCTCCAGCTGGGCCTGA-3’), and Rv1860 (FAP-B; 5´-ATGCATCAGGTGGACCCCAA-3’ and 5’- AGCGGACCTTACCGGCCTGA-3’). The PCR was conducted using 2 μl of DNA with 20 pmol of each primer in a 50 μl reaction PCR Master Mix (Promega, Madison, WI, USA) using a GeneAmp PCR System 2700 thermocycler (Applied Biosystems, Carlsbad, CA, USA). PCR products were electrophoresed on 1. 5% agarose gels to check the size of the amplified fragments by comparison to molecular weight marker GeneRuler 1kb DNA ladder (Thermo Scientific, Waltham, MA, USA). Amplified DNA fragments were purified with a PureLink Quick PCR Purification Kit (Thermo Scientific, Waltham, MA, USA) and sequenced using the reverse primer on each locus. Amplicons from at least two independent PCR reactions were sequenced.
Among the 900 MTBC field isolates analyzed, 62 different spoligotypes were identified suggesting that the study area is one of the regions with the highest diversity of MTBC spoligotypes described in the literature [47]. The high genetic diversity of MTBC in this area supported an important natural scenario where MTBC and particularly M. bovis have diversified, thus offering an interesting epidemiological and evolutionary context where new genotypes can emerge and diversify in terms of adaption to host under a range of environmental and human driven factors. The spoligotypes SB0121, SB0134, SB0339, SB0120, SB1263 and SB0157 were among the most frequent ones found in both livestock and wild ungulates in the study area (S1 Fig). The spoligotypes found at a higher frequency than expected given a predicted mutation rate were identified (Fig 1A) and represented in a hierarchical tree (Fig 1B). The spoligotypes SB0134 (MB3) and SB0339 (MB1 and MB4), which clustered in space and time and showed low mutation rates but high abundance, were selected as emergent under study conditions. The output hierarchical tree suggested a history of mutation events and a relationship between these spoligotypes with different levels of delection for spoligotypes SB0134 (MB3) and SB0339 (MB1 and MB4) (Fig 1B). These spoligotypes were therefore chosen for this study as fulfilling selection criterion (i) and (ii) described above. Additionally, findings in Iberian red deer showed a higher lesion score caused by spoligotype SB0134 when compared to spoligotype SB0339 (Mann-Whitney U test; p = 0. 04) while these spoligotypes and particularly SB0134 and SB0157 suggested an association with high TB severity in Eurasian wild boar thus also fulfilling selection criteria (iii) (Fig 1C). In summary, three M. bovis (MB1, MB3, MB4) and one M. caprae (MB2) field isolates with different prevalence, lesion score and host distribution phenotypes were selected for this study (Table 1). MB3 showed high distribution and lesion score while MB1 and MB4 were highly distributed but with low lesion score. The M. caprae (MB2) isolate had moderate distribution and high lesion score and was selected for comparison with the M. bovis isolates. These phenotypic variations are relevant for pathogen transmission and virulence and could be correlated with genome sequence information with implications for TB disease risk assessment and control. The results of the phylogenetic analysis showed that the MB1 and MB4 isolates with the same spoligotype were the most closely related isolates (Figs 2A and S2). The MB2 isolate (M. caprae) clustered separately from the other isolates, which clustered together with M. bovis sequences. Nevertheless, M. caprae (MB2) was closely related to M. bovis when compared to M. tuberculosis (Fig 2A). The genome of M. bovis BCG Pasteur 1173P2 was the most similar to all sequenced mycobacteria genomes (Figs 2A, 2B and 3). It has been proposed that during evolution, a clone of M. tuberculosis that was originally adapted to cause human TB evolved to infect a non-human mammal and thus began the transition into non-human ecotypes such as M. bovis, which in turn spread to cattle, goats, oryx, seals and pigs [20]. Our results supported a close relationship between M. bovis isolates and suggested that M. caprae is one of the M. bovis-related mycobacteria adapted to infect goats and sheep as well as other hosts such as wild boar, red deer, cattle and humans [29,48,49]. To better understand the relationship between these isolates, a comparative genomics approach was used. The results showed the presence of translocations, deletions of small genomic regions and SNPs between genomic sequences (Figs 2B, 3 and S2). Large-scale polymorphism studies have demonstrated that the MTBC shows a large number of deletions of small genomic regions consistent with the reductive evolution typical of intracellular bacteria [20]. Nevertheless, sequential chromosomal nucleotide substitutions are considered to be the main driver in the M. tuberculosis genome evolution [20]. The results reported here supported these findings for M. bovis isolates and suggested that the isolates with high lesion score, MB2 and MB3, contained the largest number of polymorphisms when compared to the MB1 and MB4 isolates with low lesion score (Figs 2B and 3). However, a clear correlation between phenotype and genome sequence requires a protein-targeted comparative analysis between the different isolates. For protein-targeted analysis, the study was focused on proteins that are known to play an important role in mycobacterial viability or virulence. Protein-targeted comparative analysis was conducted for (a) the ESX or type VII secretion system, (b) proteins linking stress response with lipid metabolism, (c) host T cell epitopes of mycobacteria, (d) antigens and (e) peptidoglycan assembly protein to define possible correlates with bacterial virulence and viability or distribution.
Mycobacteria belonging to the Mycobacterium tuberculosis complex infect humans and animals since pre-history and are a serious health problem worldwide. Whole-genome sequencing and comparative genomics generate information on the evolution and molecular basis of pathogenicity and transmissibility. However, while genomic information is increasingly available for the main human pathogens such as Mycobacterium tuberculosis, little is known about closely related bacteria, Mycobacterium bovis and Mycobacterium caprae. These mycobacteria infect humans causing zoonotic tuberculosis and are the main causative agents of animal tuberculosis. Although human-to-human transmission of zoonotic tuberculosis is limited, the infection often causes extra-pulmonary disease in humans and is still a major public health concern in developing countries, causing not only human disease but also severe effects on livelihoods. In this study, whole-genome sequences and targeted comparative genomics of three Mycobacterium bovis and one Mycobacterium caprae field isolates generated new information on the evolution and phenotypic variation of these mycobacteria. The results identified new genetic markers and candidate vaccine antigens that warrant further study to develop tools to evaluate risks for tuberculosis caused by M. bovis/M. caprae and for disease control in humans and animals.
Comparative Genomics of Field Isolates of Mycobacterium bovis and M. caprae Provides Evidence for Possible Correlates with Bacterial Viability and Virulence