Molecular and genomic features of Mycobacterium bovis strain 1595 isolated from Korean cattle.

The aim of this study was to investigate the molecular characteristics and to conduct a comparative genomic analysis of Mycobacterium (M.) bovis strain 1595 isolated from a native Korean cow. Molecular typing showed that M. bovis 1595 has spoligotype SB0140 with mycobacterial interspersed repetitive units-variable number of tandem repeats typing of 4-2-5-3-2-7-5-5-4-3-4-3-4-3, representing the most common type of M. bovis in Korea. The complete genome sequence of strain 1595 was determined by single-molecule real-time technology, which showed a genome of 4351712 bp in size with a 65.64% G + C content and 4358 protein-coding genes. Comparative genomic analysis with the genomes of Mycobacterium tuberculosis complex strains revealed that all genomes are similar in size and G + C content. Phylogenetic analysis revealed all strains were within a 0.1% average nucleotide identity value, and MUMmer analysis illustrated that all genomes showed positive collinearity with strain 1595. A sequence comparison based on BLASTP analysis showed that M. bovis AF2122/97 was the strain with the greatest number of completely matched proteins to M. bovis 1595. This genome sequence analysis will serve as a valuable reference for improving understanding of the virulence and epidemiologic traits among M. bovis isolates in Korea.

Introduction

Mycobacterium (M.) bovis is the causative agent of bovine tuberculosis (bTB), which is responsible for chronic infectious symptoms in livestock, wildlife, and humans [12]. In Korea, field surveillance for bTB relies on the caudal fold skin test using purified protein derivative bovine tuberculin. This is the most widely used method worldwide; however, it is very labor-intensive, because each cow needs to be injected individually, and the investigators are required to visit each farm twice: first for the injection and then again for evaluation of the tuberculin effect [17]. Animals that test positive are compulsorily culled with government compensation, whereas those with a suspected test result are to be re-evaluated after 60 to 90 days from the initial test, and their movement is restricted until the farm achieves bTB-free status [28]. In addition, post-mortem examination at a slaughterhouse is performed by meat inspectors to determine whether the animals have bTB. Although the national eradication program for bTB in Korea has been considered effective, a bTB rate of up to 0.25% has remained stable throughout the country [1417]. Thus, more efficient and effective diagnostic methods and strategies are required for bTB eradication. Accordingly, molecular and genomic studies are required to obtain a better epidemiologic understanding of bTB toward achieving the goal of eradication.

The most common molecular typing methods applied to strains of M. bovis are spoligotyping and mycobacterial interspersed repetitive units-variable number of tandem repeats (MIRU-VNTR) typing. Spoligotyping is performed for genotypic differentiation of Mycobacterium tuberculosis complex (MTC) strains based on polymorphisms of the direct repeat region, which contains a variable copy number of a conserved 36 bp sequence separated by multiple unique spacer sequences [18]. The spoligotyping pattern is specific for strains of a particular evolutionary lineage [33]. MIRU-VNTR typing is adopted to analyze DNA segments containing tandemly repeated sequences with variable copy numbers among strains, and has proven to be a convenient, discriminatory technique for molecular epidemiological studies of M. bovis [2]. VNTR typing has been considered very useful for disease transmission studies or for tracing new infections. The combination of spoligotyping and VNTR is also a very powerful tool to distinguish among epidemiologically related strains [5].

Complete genome sequencing is increasingly used to reveal the functions or virulence factors of a strain. In fact, predicting genomic functions through analysis of genomic DNA or amino acid sequences is very crucial for disease monitoring because experimental methods for elucidating the specific role of genes of interest are currently limited. Therefore, it is first necessary to gain complete understanding of the metabolism and mechanism of pathogenicity of M. bovis strains by conducting sequencing and annotation, which allows for comparison with other strains and can help to explain how genes are maintained in various populations. In particular, previous studies have demonstrated that the mycobacterial cell envelope contains various lipids that function as virulence factors [911]. The thick, waxy mycobacterial cell wall provides protection against harmful environments. M. bovis strains in Korea are of particular concern. In this study, we present the results of the molecular diagnosis and complete genome sequence of M. bovis strain 1595, which was isolated from the laryngopharyngeal lymph node of a bTB-infected cow (Bos taurus coreanae) in Chungcheongbuk-do, Korea. In addition, comparison of this sequence with other related sequences was performed at the genome level.

Materials and Methods

Isolation of bacteria

In 2012, M. bovis strain 1595 was obtained as a field isolate from the laryngopharyngeal lymph node of a cow showing a tuberculosis lesion during routine meat inspection at a slaughterhouse. The cow was raised in an in-house cattle farm in Chungcheongbuk-do province, in the central region of Korea. Handling and sampling of the laryngopharyngeal lymph node were performed by veterinarians in accordance with Korean regulations [30]. The tissue sample was homogenized and decontaminated using 10% (v/v) oxalic acid for 10 min. The supernatant obtained after centrifugation was cultured on Lowenstein-Jensen medium supplemented with 0.1% pyruvate (Sigma-Aldrich, USA) and incubated at 37℃ until bacterial growth was visible. The identity of bacterial colonies was confirmed by performing M. bovis-specific polymerase chain reaction (PCR) using the primers IS1081 and IS6110 (Solgent, Korea). First, the PCR mixture was prepared with 10 pM of each primer pair, 12.5 µL of 2× multiplex premix, and 5 µL of template DNA in a final volume of 25 µL. In each set of reactions, genomic DNA of M. bovis or sterile distilled water was included as the positive or negative control, respectively. The PCRs were performed with a TProfessional TRIO Thermocycler (Biometra, Germany) as follows: Pre-treatment at 50℃ for 3 min and pre-denaturation at 94℃ for 4 min, followed by 35 cycles of denaturation at 95℃ for 30 sec, annealing at 60℃ for 30 sec, extension at 72℃ for 40 sec, and final extension at 72℃ for 5 min. PCR products were analyzed on a 1.5% agarose gel and the 100 bp DNA ladder was used for estimating the size of the PCR products under UV light with RedSafe (Intron, Korea).

DNA extraction

Genomic DNA extraction was performed by applying the enzymatic lysis method [4] with some modifications. In brief, the mycobacterial pellet was harvested and washed once in 5 mL of Tris-ethylenediaminetetraacetic acid (TE) buffer (pH 8.0; Sigma-Aldrich). The cell pellet was frozen at 80℃ overnight and resuspended in 5 mL of TE buffer. Five milliliters of chloroform/methanol (2:1, v/v) was added and rocked at room temperature for 1 h. The suspension was centrifuged at 2,500 × g for 30 min, and then decanted carefully to discard both the organic and aqueous layers. The tube was placed at 55℃ for 15 min, and 5 mL of TE buffer was added for resuspension. The mixture was vortexed vigorously, and 0.5 mL of 1 M Tris-base (pH 9.0; Sigma-Aldrich) and 0.05 mL of lysozyme solution (10 mg/mL; Sigma-Aldrich) were added to the cell suspension and incubated at 37℃ for 12 h. Next, 10 µL of DNase-free RNase (Qiagen, USA) was added and further incubated at 37℃ for 30 min. Subsequently, 0.5 mL of 10% sodium dodecyl sulfate (w/v; Sigma-Aldrich) and 0.05 mL of proteinase K (Qiagen) were added to the cell lysate, and the solution was incubated at 55℃ for 3 h. Then, an equal volume of phenol-chloroform-isoamyl alcohol (25:24:1; Sigma-Aldrich) was added, followed by centrifugation at 12,000 × g for 30 min at room temperature. The aqueous layer was transferred to a sterile tube and centrifuged at 12,000 × g for 30 min at room temperature after adding chloroform-isoamyl alcohol (24:1). Then, 0.5 mL of 3 M sodium acetate (pH 5.2) and 5 mL of isopropyl alcohol were added to the aqueous layer for precipitation of DNA, and the mixture was allowed to stand at 4℃ overnight. The DNA was centrifuged at 12,000 × g for 30 min at 4℃ and washed in 70% ethanol. The pellet of genomic DNA was allowed to air dry and then dissolved in distilled water. To assess the quality of genomic DNA, a Nanodrop spectrophotometer (Thermo Scientific, USA) was used to measure the concentration of DNA. An ideal sample was one in which similar concentration estimates were obtained on both platforms, with A260/280 and A260/230 ratios close to 1.8 and 1.5, respectively, according to the standard expectations of pure DNA.

Spoligotyping and MIRU-VNTR typing

Molecular typing was performed by spoligotyping and MIRU-VNTR typing. Spoligotyping was performed by using a REBA Spoligotyping-Strip (Molecules and Diagnostics, Korea), based on a previously described method [18]. Direct repeat regions were amplified in a total volume of 20 µL containing 10 ng of M. bovis genomic DNA with the PCR mixture included in the kit under the following protocol: Pre-denaturation at 94℃ for 5 min, followed by 35 cycles of denaturation at 94℃ for 30 sec, annealing at 55℃ for 30 sec, extension at 72℃ for 30 sec, and final extension at 7℃ for 7 min. The result was confirmed based on comparison to the international spoligotyping database (Mbovis.org; VISAVET Health Surveillance Centre, Spain). MIRU-VNTR profile identification was performed by using PCR amplification of 14 targeting genetic loci on the basis of the high discriminatory power among M. bovis isolates from Korea (MIRU 4, 16, 26, 27, 31; ETR-A, -B, -C; QUB 11b, 18, 26, 3336; VNTR 2401 and 3171) as previously described [1627]. In brief, the PCR mixtures were prepared with 10 pM of each primer pair and 50 ng of mycobacterial genomic DNA using AccuPower HotStart PCR PreMix (Bioneer, Korea). In each set of reactions, genomic DNA of M. tuberculosis H37Rv was used as the positive control, and distilled water was used as the negative control; the PCR conditions differed according to the annealing temperature of each primer. The size of the amplicons was compared by gel electrophoresis with a 100-bp DNA ladder (Bioneer).

Genome sequencing and comparative genomic analysis

The complete genome of M. bovis strain 1595 was constructed by using sequencing data from Illumina MiSeq (Illumina, USA) and PacBio RS (Pacific Biosciences, USA). The genome sequence of M. bovis strain 1595 was deposited in GenBank under the accession number CP012095, and the genome sequencing, assembly, and annotation were conducted as described in our previous report [21]. In brief, Consed version 19.0 (University of Washington, USA) was used for editing of the assembled sequence, and CLC Genomic Workbench (CLC bio, Denmark) and CodonCode Aligner (CodonCode, USA) were used for mapping and correction. Gene prediction was carried out by using Glimmer version 3 [8], and annotation was performed with the National Center for Biotechnology Information (NCBI) Prokaryotic Genome Annotation Pipeline, using default settings [34]. The transfer RNAs (tRNAs) and ribosomal RNAs (rRNAs) were identified by utilizing tRNAscan-SE version 1.21 [29] and RNAmmer version 1.2 [26], respectively. The annotated open reading frames (ORFs) were edited or corrected by using Artemis (ver. 1.4; Sanger Institute, UK) [6], and clusters of orthologous groups (COGs) of proteins were constructed for functional classifications by using the eggNOG database (ver. 3; EMBL, Germany) [15].

The sequences of mycobacteria closely related to M. bovis strain 1595, including AF2122/97 (GenBank accession No. NC_002945.3), BCG Pasteur 1173P2 (NC_008769.1), W-1171 (JXTK00000000 [22]), and M. tuberculosis H37Rv (NC_000962.3) were obtained from the EzGenome database (ChunLab, Korea) [23] for comparative genomics analysis. The inter-genomic distances among genome sequences were determined from fully or partially sequenced genomes by using the average nucleotide identity (ANI) [10] with the proposed cut-off for species boundary of 95% to 96% [32]. ANI values were calculated by using JSpecies (IMEDEA, Spain), based on the BLAST algorithm (NCBI, USA). From this pairwise distance matrix, an ANI tree was constructed using the unweighted pair group method with arithmetic mean clustering method [7]. In addition, genomic sequences with adjusted start positions were submitted to the NUCmer system (NUCleotide MUMmer version 3.07; GitHub, USA) for alignment [25]. Synteny plots were generated by the MUMmer program using exact matching, clustering, and alignment extension strategies based on the number of identical alignments between two genomes. For results interpretation, a red line in the forward direction indicates the collinearity of two sequences. For amino acid comparison, BLASTP comparison of query ORFs to subject ORFs was performed using annotated genomes. Orthologs and paralogs were differentiated by reciprocal comparison.

Results

Features of M. bovis strain 1595

Rough mycobacterial colonies were observed on the solid medium. Bacterial colonies showed a positive result in the M. bovis-specific PCR with the IS1081 and IS6110 primers. M. bovis strain 1595 was observed to have the spoligotype SB0140 and the MIRU-VNTR type 4-2-5-3-2-7-5-5-4-3-4-3-4-3.

The circularized genome of M. bovis 1595 was observed to contain 4,351,712 bases with a 65.64% G + C content. A total of 4358 protein-coding sequences (CDSs), 3 rRNAs, and 45 tRNAs were predicted. The percentage of protein-coding genes was 93.1% with an average of 929 bp and a 65.86% G + C content. The genomic characteristics of M. bovis 1595 are shown in Fig. 1. Among the 4358 CDSs, 3663 ORFs (84.05%) were classified into COG families comprising 18 functional categories (Table 1). Of these, 1291 ORFs (29.62%) were classified into the “metabolism” (C, G, E, F, H, I, P, and Q) category, almost half of which were related to “lipid transport and metabolism”, “energy production and conversion”, and “amino acid transport and metabolism”. A total of 401 ORFs (9.20%) were classified into cellular processes and signaling-related categories (D, O, M, N, and T), 598 ORFs (13.72%) were classified into information storage and processing categories (J, K, and L), and 1,373 ORFs (31.51%) were attributed to poorly characterized COG groups (R and S); 695 ORFs (15.95%) could not be categorized.

Fig. 1

The entire genomic sequence of Mycobacterium bovis strain 1595. The scale is shown in megabases on the outer black circle. From the outside to the center: RNA features (ribosomal RNAs are shown in blue, and transfer RNAs are shown in red), genes on the forward strand, and genes on the reverse strand (colored according to the clusters of orthologous groups categories). The inner two circles show the GC ratio and GC skew. The GC ratio and GC skew shown in orange and red indicate positive values, respectively, and those shown in blue and green indicate negative values, respectively.

Table 1

Clusters of orthologous groups (COGs) of the genome of Mycobacterium bovis strain 1595

Data are presented as number (%).

Comparative genomics

Comparison of the M. bovis strain 1595 genome sequence with those of four MTC strains showed that its circular DNA molecule of 4351712 bp was similar to that of other related genomes. The genomic characteristics of strains, including strain 1595 and related strains, are summarized in Table 2. Although strain 1595 has a higher number of CDSs, all genomes analyzed in this study showed a similar genome size and G + C content.

Table 2

Genomic characteristics of strains used in this study

Contigs, contiguous; CDSs, coding sequences; rRNA, ribosomal RNA; tRNA, transfer RNA.

To investigate the genomic relatedness among MTC strains, including M. bovis strain 1595 and closely related strains, we analyzed their fully or partially sequenced genomes according to the ANI values (Fig. 2). The resulting phylogenetic tree revealed a strikingly short genetic distance between all of these strains within an ANI value of 0.1%. As shown in Fig. 2, the most closely related strain was M. bovis W-1171. Although the most distant strain was M. tuberculosis H37Rv, the phylogenetic distance was still not greater than 0.1. The MUMmer analysis showed a red line in the forward direction, indicating that the genomes were mostly co-linear (Fig. 3). Except for the genes predicted to be related to general and unknown functions, M. bovis strains such as 1595, W-1171, AF2122/97, and BCG Pasteur 1173P2 showed the same ranking orders in COG analysis, which were “transcription”, followed by “lipid transport and metabolism” and “energy production and conversion”. The top three COG ranking orders of M. tuberculosis H37Rv were “transcription”, “energy production and conversion”, and “lipid transport and metabolism”.

Fig. 2

Genome tree based on the average nucleotide identity (ANI) values showing the relationships among Mycobacterium (M.) tuberculosis complex strains including M. bovis 1595. To convert the ANI value into a genetic distance, its complement to 1 was taken. From this pairwise distance matrix, an ANI tree was constructed using the unweighted pair group method and the arithmetic mean clustering method.

Fig. 3

Nucleotide-based alignments with NUCmer. X-axis: Mycobacterium bovis strain 1595. Y-axis: (A) Mycobacterium bovis W-1171, (B) AF2122/97, (C) BCG Pasteur 1173P2, and (D) Mycobacterium tuberculosis H37Rv. Aligned segments are presented as dots or lines in the NUCmer alignment and were generated by the MUMmer plot script.

In addition, sequence comparison based on BLASTP showed that AF2122/97 shares 83.91% of its complete proteins with strain 1595. The genome of W-1171 ranked second in similarity to strain 1595, with 80.66% of the complete proteins being identical. The most distinct strain was M. tuberculosis H37Rv, which showed the highest percentage of paralogs (3.95%) among all analyzed strains (Table 3).

Table 3

Similarities of open reading frames (ORFs) of various Mycobacterium (M.) tuberculosis complex strains compared to M. bovis 1595 (n = 4,358)

Data present number (%) of ORFs.

Discussion

Bovine tuberculosis is an important contagious disease of livestock and humans that needs to be controlled because it consequently brings a severe economic burden to the agriculture industry and public health worldwide. The objective of this study was to identify the characteristics of M. bovis strain 1595, isolated from the laryngopharyngeal lymph node of Korean cattle, by using bacterial culture, PCR, spoligotyping, and MIRU-VNTR, as well as analysis of the complete genome sequence. Strain 1595 showed the typical morphological characteristics of M. bovis, with rough colonies on solid medium, and was positive to M. bovis-specific IS1081 and IS6110 primers through conventional PCR. Subsequent spoligotype and MIRU-VNTR assays determined that the strain contains the most common and widespread spoligotype found in Korea [20].

Whole-genome sequencing is now commonly adopted, but its application to mycobacteria remains a challenge due to their cellular composition such as the cell wall, which makes DNA extraction complicated [19]. We applied an enzymatic and chemical DNA extraction method and obtained intact DNA, which was used to perform Illumina and PacBio sequencing. The data obtained from an initial version, such as Illumina, can be used alongside data obtained from an updated system, such as PacBio, to correct for the high error rate, in accordance with previous studies [24]. We confirmed that this is the best practice for accomplishing complete bacterial genome assemblies using both short and long reads.

Genome comparison revealed that strain 1595 was very similar to other MTC genomes tested in this study. All bacteria belonging to the MTC group have evolved from a common ancestral strain and diverged as a result of chromosomal deletions without exchanging genes between strains, leading to a high similarity rate [1]. Consistent with this, M. bovis strain 1595 was determined to be closely related with other MTC strains in the genome comparison. The phylogenetic distance among strains evaluated in this study was similar, with all ANI values less than 0.1%, which indicates high genetic relatedness in a pairwise whole-genome comparison between prokaryotic strains. The most closely related strain was M. bovis W-1171, which is likely related to a local geographic association. W-1171 is a Korean isolate from a wild boar, and was previously sequenced in our laboratory [22]. However, all MTC strains showed a genetic distance within 0.1%. Consistent with this, the synteny plot produced by MUMmer analysis showed collinearity between strain 1595 and each MTC genome.

Gene compositions categorized by function were also similar among strains according to the COG results. One particularly notable COG result was that MTC strains tend to contain various kinds of lipid metabolism- and transport-related genes. It is widely known that mycobacterial virulence lipids such as trehalose dimycolate/monomycolate, phthiocerol dimycocerosate, and sulfolipids, and the protein machinery such as mycobacterial membrane protein large (MmpL) and mammalian cell entry protein required for their export are intercalated in the cell wall [13]. Mycobacterial cells have a unique component containing a high proportion of lipids in their envelopes. In particular, an outer layer called the mycomembrane consists of an asymmetric lipid bilayer made of long-chain (C60–C90) mycolic fatty acids in the inner leaflet and free intercalating glycolipids and waxy components on the outer leaflet [3]; this is in contrast to most prokaryotes, which produce fatty acids in the C14 to C18 range [9]. Moreover, mycolic acids are responsible for not only preserving the cell wall structure and function but also for regulating the interaction with the host immune system [9]. Pirson et al. [31] investigated M. bovis AF2122/97, and showed that the production of the cytokines interleukin (IL)-10, IL-12, and IL-6, and tumor necrosis factor-alpha were increased significantly after exposure to virulence lipids. In addition, cell envelope proteins are very important for mycobacterial pathogenesis. For example, fibronectin-binding protein (Fbp), also known as antigen 85 (Ag85), is a complex of FbpA, FbpB, and FbpC2, which are encoded by the fbpA, fbpB, and fbpC2 genes located in M. tuberculosis H37Rv genomic regions Rv3804c, Rv1886c, and Rv0129c, respectively [9]. The Fbp complex can promote the adhesion of mycobacteria to the mucosal surface through fibronectin binding [35]. The proteins of M. bovis 1595, Mb1595_p4218, _p2114, and _p0152 matched with Rv3804c, Rv1886c, and Rv0129c of M. tuberculosis H37Rv, respectively, as identified by BLASTP; these were all found to be well conserved and can be considered orthologous. This can provide evidence that M. bovis 1595 has similar Ag85 function to that of M. tuberculosis H37Rv. Further support for functional similarity relates to MmpL and mycobacterial membrane protein small (MmpS), which mediate the transport of important cell wall lipids across the mycobacterial membrane and are responsible for drug efflux, siderophore export, and heme uptake [3]. According to the BLASTP result, M. tuberculosis H37Rv has 14 MmpL and 5 MmpS proteins, which are all well conserved in M. bovis 1595 (similarity > 95%), although the similarity of Mb1595_p0447 to the Rv0402c (MmpL1) sequence was less than 10%, and Mb1595_p1302 was identified as a paralog of Rv1145 (MmpL13a). Therefore, these proteins may be responsible for lipid transport and lipid biosynthesis, which are essential for mycobacterial pathogenesis. Overall, the nucleotide similarity evaluated by BLASTP showed that MTC strains have more than 70% completely matched proteins; therefore, genes with the same or similar functions are thought to be well conserved. Even though the strain with the greatest similarity (83.91% completely identical proteins) was M. bovis AF2122/97, M. bovis W-1171 also showed a high similarity, with 80.66% completely matching proteins to strain 1595.

The main focus of this study was to provide information about M. bovis strain 1595; however, more data of Korean isolates are needed to fully elucidate the characteristics of M. bovis isolates from Korea. Nevertheless, these findings will be useful for the standardization of M. bovis strains in Korea and should enable gaining a more precise description of the genotypic and phenotypic differences among strains. In further studies, comparison of single nucleotide polymorphisms will be needed to comprehensively investigate the similarity to and genetic evolution of other strains, including isolates from Korea and other countries.