Single-nucleotide polymorphism-based epidemiological analysis of Korean Mycobacterium bovis isolates.

BACKGROUND: Bovine tuberculosis (TB) is caused by Mycobacterium bovis, a well-known cause of zoonotic tuberculosis in cattle and deer, and has been investigated in many physiological and molecular studies. However, detailed genome-level studies of M. bovis have not been performed in Korea.
OBJECTIVES: To survey whole genome-wide single-nucleotide polymorphism (SNP) variants in Korean M. bovis field isolates and to define M. bovis groups in Korea by comparing SNP typing with spoligotyping and variable number tandem repeat typing.
METHODS: A total of 46 M. bovis field isolates, isolated from laryngopharyngeal lymph nodes and lungs of Korean cattle, wild boar, and Korean water deer, were used to identify SNPs by performing whole-genome sequencing. SNP sites were confirmed via polymerase chain reaction using 87 primer pairs.
RESULTS: We identified 34 SNP sites with different frequencies across M. bovis isolates, and performed SNP typing and epidemiological analysis, which divided the 46 field isolates into 16 subtypes.
CONCLUSIONS: Through SNP analysis, detailed differences in samples with identical spoligotypes could be detected. SNP analysis is, therefore, a useful epidemiological tracing tool that could enable better management of bovine TB, thus preventing further outbreaks and reducing the impact of this disease.

INTRODUCTION

Bovine tuberculosis (TB) is an infectious disease caused mainly by Mycobacterium bovis and has a negative economic impact on the livestock industry worldwide. Members of the Mycobacterium tuberculosis complex (MTC), including M. bovis, have identical 16S rRNA sequences and very similar whole genome sequences, leading to difficulties in the rapid identification and molecular analysis of these strains [1]. To overcome this challenge, methods, such as spoligotyping and mycobacterial interspersed repetitive unit-variable number tandem repeats (MIRU-VNTR), have been developed to analyze the molecular epidemiology of members of the MTC [23]. These methods are used to identify the epidemiological relationships between outbreaks and the transmission of TB. M. bovis research in Korea has traditionally been conducted using spoligotyping and MIRU-VNTR. However, in Korea, only two spoligotypes (SB1040 and SB0140) have been identified to date. Although 24 VNTR types have been recorded in Korea, five VNTR types account for 70.22% of all Korean isolates [4]. There are limitations to the use of spoligotyping and MIRU-VNTR for the epidemiological analysis of M. bovis isolates from Korea, including the small numbers of types identified to date as well as restrictions in dynamic tracking [5]. New methods capable of more fine-grained tracing of bacterial lineages are therefore urgently needed for the molecular epidemiological analysis of M. bovis transmission in Korea.

Recently, a technique for classifying and identifying MTC strains using next-generation sequencing has been reported [67]. Whole-genome sequencing (WGS) was used to evaluate M. tuberculosis strains isolated from humans in an area where TB had been prevalent for several years to investigate the propagation and epidemiology of TB bacteria [89]. The spread of TB bacteria between wild and farm animals was previously tracked by performing WGS of M. bovis strains obtained from four badgers and 26 cows in the UK, with four significant single-nucleotide polymorphisms (SNPs) identified that could successfully map the genetic divergence and spatiality of M. bovis [10].

SNPs occur at a high frequency and relatively uniformly across genomes, and are therefore used as markers in several genotyping methods [1112]. SNPs may differentiate very limited allelic differences, particularly in the 16S rRNA gene and other highly conserved genes. SNP analysis enables the classification of very similar genomes and detailed classification of genomes that are difficult to distinguish using conventional typing techniques, such as those of the MTC, which has > 99% sequence similarity [1314]. Thus, SNPs have also been used in M. bovis for typing and clustering analysis. However, SNP analysis of Korean M. bovis isolates has not been performed [7151617]. In the present study, we performed SNP analysis of Korean M. bovis isolates via WGS to facilitate typing of M. bovis strains. We identified SNPs on a large scale and used them for molecular epidemiological investigation of bacterial isolates from domestic animals with TB in Korea. We compared the efficacy of SNP typing to that of spoligotyping and MIRU-VNTR typing to identify the best method for continued surveillance of M. bovis outbreaks in Korea.

MATERIALS AND METHODS

Field strain isolation of M. bovis

Between 2015 and 2018, 46 lymph node and lung samples from cattle (n = 33), deer (Cervus nippon, n = 10), wild boar (Sus scrofa, n = 1), and wild goral (Naemorhedus caudatus, n = 2) were collected for use in the present study. Cattle and deer samples were collected from TB-like lesions detected during routine abattoir inspection, or from animals clinically suspected of having TB based on a positive reaction in a skin test or gamma interferon test. Bacteria were isolated from the samples showing TB-positive lesions in wild animals. Tissue samples were collected from M. bovis TB lesions in livestock. The tissue was homogenized with an auto-mill machine (Tokken, Japan) and placed in tubes containing 10 mL of phosphate-buffered saline (PBS), followed by decontamination with an equal volume of 10% oxalic acid for 10 min. The homogenate was pelleted by centrifugation for 10 min at 1,000 × g, and the pellet was resuspended in 2 mL of PBS. Next, 200 μL of each sample was inoculated into Middlebrook 7H11 medium (Becton Dickinson, USA) in mycobacterial growth indicator TB tubes (Mediland, Korea), followed by incubation at 37°C for a minimum of 8 weeks until bacterial colonies were visible. We confirmed that there was no requirement to obtain approval for research conducted on animals slaughtered due to infectious diseases from the Institutional Animal Care and Use Committee of South Korea.

M. bovis DNA extraction and polymerase chain reaction (PCR) amplification

Isolated colony samples suspended in PBS (180 mL) were digested with 200 µL (20 ng/mL concentrate) of lysozyme buffer (Sigma, USA) at 37°C for 3 h, after which 25 μL of Proteinase K (Qiagen, Germany) (60°C, 3 h) was added to lyse the cells. Genomic DNA was extracted using a NucleoSpin Soil Kit (Macherey-Nagel, Germany). To confirm the presence of Mycobacteria in the isolated colonies, a multiplex PCR-based DiaPlexC™ MTB/M. bovis Detection kit (Solgent, Korea) was used according to the manufacturer's instructions. PCR products were separated using electrophoresis on a 1.5% agarose gel (140 V, 30 min). The genomic DNA was used for SNP analysis, spoligotyping, and MIRU-VNTR analysis.

Spoligotyping and MIRU-VNTR analysis

Spoligotyping was performed as described by Kamerbeek et al. [18] and Le Fleche et al. [19]. For MIRU-VNTR profiling, 16 target genetic loci (MIRU4, MIRU16, MIRU26, MIRU27, MIRU31, ETR-A, ETR-B, ETR-C, QUB11a, QUB11b, QUB18, QUB26, QUB3232, QUB3336, QUB2401, and QUB3171) were PCR-amplified (Supplementary Table 1). PCRs were performed using AccuPower® HotStart PCR premixture (Bioneer, Korea) in 20 µL reaction mixtures comprising 2 µL of DNA, 1 µL (10 pmol/µL) of each primer, and 14 µL of nuclease-free pure water. Primer sets and annealing temperate conditions are provided in Supplementary Table 1. PCR was performed under the following conditions: initial denaturation at 94°C for 4 min, followed by 35 cycles of 94°C for 30 s; (annealing temp)°C for 1 min; and 72°C for 1 min, with a final extension at 72°C for 7 min, for all loci except QUB3336. For QUB3336, the PCR conditions were as follows: initial denaturation at 94°C for 4 min, followed by 35 cycles of 94°C for 30 s; 57°C for 30 s; and 72°C for 30 s, with a final extension at 72°C for 7 min. To confirm the VNTR fragment sizes, the PCR products were separated on a 1.5%-agarose gel via electrophoresis (145 V, 35 min).

SNP analysis

SNP screening and selection from whole-genome sequences of Korean field isolates

Among M. bovis field isolates from a Korean TB outbreak on a farm identified in a previous study, the dominant strains M. bovis B-1595 (spoligotype SB0140) (GenBank accession no. CP012095; draft sequence, 50 contigs) and B-3222 (GenBank accession no. NZ_LNOF00000000.1; draft sequence, 132 contigs) were used for WGS [20]. WGS data were generated using the Illumina and Pacific Biosciences (PacBio) sequencing platforms. Adapters and low-quality reads were removed using Trimmomatic v0.32 (Illumina, USA). The two dominant strains were compared with the type strain, BCG Pasteur 1173P2 (complete sequence), to identify SNPs/InDels. After mapping to the BCG Pasteur 1173P2 sequence using Bowtie2 2.1.0, variant calling was conducted to confirm the presence of SNPs and indels using GATK 3.2.2. Annotation was performed based on the contents of the existing reference genome (GenBank accession no. NC_008769) using SnpEff 4.0.

Primer design for SNP site confirmation

The whole-genome sequences of B-1595 and B-3332 were used to select SNP sites. The SNP locations were checked via sequence parallel comparison using the UniProt protein sequence database (https://www.uniprot.org). To screen each virulence factor, gene sites were searched for SNPs. Gene sequences of the 89 virulence factors were collected and aligned with M. bovis B-1595 and AF2122/97 whole-genome sequences. To confirm the existence of SNPs using PCR, a set of 87 primer pairs (primers ndh-1,2 and pks12-3,4 can identify two SNPs with one primer set each) to amplify SNP regions of 173–675 base pairs each was created using CLC Main Workbench 7 (Qiagen) (Supplementary Table 2).

PCR-based sequencing for SNP screening

Prior to confirming all SNP candidates via PCR-based sequencing, primary screening was performed to confirm whether there were significant differences in the 89 SNP sites identified between 11 field isolates of M. bovis covering the major VNTR types and the type strains M. bovis AN5 and M. tuberculosis H37Rv. Genomic DNA was extracted from the field isolates, and amplicon libraries of SNP site gene fragments were amplified from the genomic DNA using the SNP search primer sets. The 13 DNA samples (W-1171, B-5648, D-174, B-4559, B-1595, B-3222, D-151, B-5628, D-133, B-5330, D-145, AN5, and H37Rv) were amplified, and purified using the QIAquick Gel extraction Kit (Qiagen). Template DNA (2 µL), primers (forward, 2 µL; reverse, 2 µL), and distilled water (14 µL) were added to the AccuPower® HotStart PCR premixture (Bioneer) for the 87 pairs of primers, and thermal cycles were performed as follows: initial denaturation at 94°C for 5 min, 37 cycles of 94°C for 30 s, 54°C for 1 min, and 72°C for 30 s, and a final extension at 72°C for 5 min and cooling at 4°C. After confirming the amplification pattern and amplicon sizes/patterns via gel electrophoresis, the PCR reaction mix was submitted for sequencing using an ABI 3730xl System DNA sequencer (Macrogen, Korea). The resultant SNP site gene sequences were aligned with the corresponding sequences in the M. bovis B-1595 whole-genome sequence.

SNP analysis of 46 field isolates of M. bovis

Based on the screening of the 89 SNP sites, 34 SNP sites were selected to establish primer sets to identify SNPs within the Korean M. bovis isolates via PCR (Supplementary Table 3). These sites were selected via domestic field strain separation analysis using additional strains (35 field isolates; Supplementary Table 4). Forty-six field isolate strains and the type strains M. bovis AN5 and M. tuberculosis H37Rv were evaluated, and sequencing was performed as described above.

SNP-based grouping of M. bovis field isolates

The sequences of SNP sites with diversity among isolates were extracted and aligned, and the nucleotide sequences were confirmed by evaluating SNP cluster groups (SCGs). To this end, we constructed an unweighted pair-group method with arithmetic mean dendrogram using the Kimura 80 distance measurement to detect differences between the 48 strains according to their SNPs. Phylogenetic neighbors were identified, and SNP pairwise sequence similarity was calculated using CLC Main Workbench 7.

RESULTS

Spoligotype and MIRU-VNTR analyses

Among the 46 Korean M. bovis field isolates, 19 strains were identified as spoligotype SB1040, and 27 strains as spoligotype SB0140 (Table 1, Supplementary Table 5). In addition, nine VNTR types, FF, P, Y, EE, FF, H, C, J, and GG, were identified. Among the 46 field isolates, the H type was the most abundant (17/46, 36.95%), followed by the FF type (7/46, 15.21%), Y type (6/46, 13.04%), P type (5/46, 10.86%), J type (4/46, 8.69%), F type (3/46, 6.52%), CC type (2/46, 4.34%), and the EE and GG types (1/46, 2.17% each) (Table 1).

Table 1

VNTR and spoligotype profiles of 46 Mycobacterium bovis Korean field isolates

Field isolates	Animal	Isolated tissue	M4	M16	M26	M27	M31	EA	EB	EC	Q11a	Q11b	Q18	Q26	Q3232	Q3336	VNTR2410	VNTR3171	Spoligotype	Type
W-1171	Wild boar	LN	4	2	5	2	3	7	5	4	10	4	3	3	10	10	5	2	SB1040	FF
B-5648	Cattle	Lung	3	2	5	3	3	6	5	4	10	4	3	3	10	10	5	2	SB1040	P
D-174	Deer	Lung	3	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	H
B-4559	Cattle	Lung	3	2	5	3	2	7	4	4	10	4	3	4	10	3	4	3	SB0140	C
B-1595	Cattle	Lung	3	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	H
B-3222	Cattle	LN	3	2	5	3	2	7	5	4	10	4	3	4	8	3	4	3	SB0140	F
D-151	Deer	Lung	3	2	5	3	3	7	5	4	10	4	3	3	10	10	5	2	SB1040	Y
B-5628	Cattle	Lung	4	2	5	2	3	7	2	4	10	4	3	3	10	10	5	2	SB1040	EE
D-133	Deer	Lung	3	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	H
B-5330	Cattle	Lung	3	2	5	3	2	7	5	4	10	4	3	4	10	4	4	3	SB0140	J
D-145	Deer	Lung	4	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	GG
B-5351	Cattle	TL	3	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	H
B-5335	Cattle	TL	3	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	H
B-5624	Cattle	Lung	3	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	H
B-5260	Cattle	Lung	3	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	H
D-171	Deer	Lung	3	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	H
B-5281	Cattle	Lung	3	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	H
B-1579	Cattle	Lung	3	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	H
B-1567	Cattle	LN	3	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	H
D-35	Deer	LN	3	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	H
W-1366	Cattle	TL	3	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	H
B-5253	Cattle	LN	3	2	5	3	2	7	5	4	10	4	3	4	10	4	4	3	SB0140	J
B-5315	Cattle	Lung	3	2	5	3	2	7	5	4	10	4	3	4	10	4	4	3	SB0140	J
15B-3	Cattle	Lung	3	2	5	3	2	7	5	4	10	4	3	4	10	4	4	3	SB0140	J
B-3221	Cattle	Lung	3	2	5	3	2	7	5	4	10	4	3	4	8	3	4	3	SB0140	F
B-3561	Cattle	Lung	3	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	H
B-4600	Cattle	Lung	3	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	H
B-3212	Cattle	Lung	3	2	5	3	2	7	5	4	10	4	3	4	8	3	4	3	SB0140	F
B-4216	Cattle	Lung	3	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	H
B-4412	Cattle	Lung	3	2	5	3	2	7	5	4	10	4	3	4	10	3	4	3	SB0140	H
D-155	Deer	Lung	3	2	5	3	3	7	5	4	10	4	3	3	10	10	5	2	SB1040	Y
B-5286	Cattle	Lung	3	2	5	3	3	6	5	4	10	4	3	3	10	10	5	2	SB1040	P
B-5240	Cattle	LN	3	2	5	3	3	6	5	4	10	4	3	3	10	10	5	2	SB1040	P
B-5282	Cattle	Lung	3	2	5	3	3	6	5	4	10	4	3	3	10	10	5	2	SB1040	P
16B-4	Cattle	Lung	3	2	5	3	3	6	5	4	10	4	3	3	10	10	5	2	SB1040	P
B-4452	Cattle	LN	3	2	5	3	3	7	5	4	10	4	3	3	10	10	5	2	SB1040	Y
B-5276	Cattle	Lung	3	2	5	3	3	7	5	4	10	4	3	3	10	10	5	2	SB1040	Y
B-4724	Cattle	Lung	3	2	5	3	3	7	5	4	10	4	3	3	10	10	5	2	SB1040	Y
B-5343	Cattle	Lung	3	2	5	3	3	7	5	4	10	4	3	3	10	10	5	2	SB1040	Y
15W-23	Wild goral	Lung	4	2	5	2	3	7	5	4	10	4	3	3	10	10	5	2	SB1040	FF
15W-442	Wild goral	TL	4	2	5	2	3	7	5	4	10	4	3	3	10	10	5	2	SB1040	FF
B-5278	Cattle	Lung	4	2	5	2	3	7	5	4	10	4	3	3	10	10	5	2	SB1040	FF
D-132	Deer	Lung	4	2	5	2	3	7	5	4	10	4	3	3	10	10	5	2	SB1040	FF
D-166	Deer	Lung	4	2	5	2	3	7	5	4	10	4	3	3	10	10	5	2	SB1040	FF
D-136	Deer	LN	4	2	5	2	3	7	5	4	10	4	3	3	10	10	5	2	SB1040	FF
B-5259	Cattle	LN	3	2	5	3	2	7	4	4	10	4	3	4	10	3	4	3	SB0140	C

Selection of candidate SNPs in virulence factors through draft annotation

For B-1595, 9,589,125 of the 12,359,672 sequences were mapped to the reference genome, and for B-3222, 6,524,564 sequences from a total of 8,000,000 sequences were mapped. SNPs were analyzed for frame-shift mutations, insertions, and deletions (synonymous [silent], non-synonymous [missense], and amino acid mutations). In B-1595, there were 418 synonymous and 759 nonsynonymous mutations in 1177 SNPs in open reading frames (ORFs), while in B-3222, there were 384 synonymous and 668 nonsynonymous mutations in 1052 SNPs in ORFs. The ratio of synonymous to nonsynonymous mutations was 1:2, similar to the ratio in SNP studies in M. tuberculosis [21]. As a result, 1,400 SNP candidates were identified in B-1595 and B-3222 (Table 2). SNPs in major mycobacterial virulence factors were selected based on draft annotation data from the 1,400 SNPs found in both M. bovis B-1595 and B-3222 compared with the AF2122/97 and BCG Pasteur 1173P2 strains (Supplementary Table 6) [22]. A total of 89 SNPs distributed in 73 virulence factors were selected from the SNP annotation draft results.

Table 2

SNPs/InDels in Mycobacterium bovis B-1595 and B-3222 compared to BCG Pasteur 1173P2

SNPs/InDels	B-1595			B-3222
	Total	ORF	Intergenic region	Total	ORF	Intergenic region
SNPs	1,374	1,177	197	1,213	1,052	161
Insertion	73	62	11	61	53	8
Deletion	55	37	18	51	34	17
*Other	1	1	0	0	0	0

SNP, single-nucleotide polymorphism; ORF, open reading frame.

*Other, difficult to determine.

SNP profiles of 46 M. bovis field isolates

Selection of final SNP candidates for epidemiological analysis of Korean field isolates

A total of 89 virulence-related SNP candidates were identified in 11 M. bovis field isolates and the type strains M. bovis AN5 and M. tuberculosis H37Rv. Of these, 56 SNPs showed no differences between the Korean field isolates, whereas 33 SNPs varied between isolates. In addition to the 33 SNPs selected, an SNP present in HisD (termed HisD-1) was detected in the draft SNP annotation, and was included for analysis. Therefore, a total of 34 SNPs were analyzed in 48 strains (including the 46 M. bovis field strains and the M. bovis AN5 and M. tuberculosis H37Rv reference strains) (Supplementary Table 5).

SNP comparison with classical typing methods

Classification of the 34 SNPs according to the spoligotype pattern revealed that all SNP sites except three, HisD-1, glpK-2, and glpK-3, had SNP profiles identical to those of the spoligotypes. HisD-1 and glpK-3 showed SNP diversity only in some isolates, and glpK-2 showed SNP diversity regardless of the spoligotype. Application of the glpK-2 SNP to a representative VNTR pattern divided the P, Y, FF, and CC types into two subtypes, and the H type was divided into three subtypes because it contained both glpK-2 and glpK-3 sites. The J type was not subdivided based on glpK-2, but was divided into two subtypes because of an SNP in HisD-1 (Table 3, Supplementary Table 7).

Table 3

Comparative analysis of SNPs with spoligotypes and representative MIRU-VNTR types

SNP gene symbol	Spoligotype		VNTR representative type
	SB1040	SB0140	P-1	P-2	Y-1	Y-2	EE	FF-1	FF-2	HH	H-1	H-2	H-3	CC-1	CC-2	J-1	J-2	FF	GG
corA	T	C	T	T	T	T	T	T	T	T	C	C	C	C	C	C	C	C	C
glpK	C	G	C	C	C	C	C	C	C	C	G	G	G	G	G	G	G	G	G
fadD9-2	G	A	G	G	G	G	G	G	G	G	A	A	A	A	A	A	A	A	A
fadE20(2)	C	T	C	C	C	C	C	C	C	C	T	T	T	T	T	T	T	T	T
fadD27	G	A	G	G	G	G	G	G	G	G	A	A	A	A	A	A	A	A	A
glnB	G	A	G	G	G	G	G	G	G	G	A	A	A	A	A	A	A	A	A
HisD	A	G	A	A	A	A	A	A	A	A	G	G	G	G	G	G	G	G	G
HisD-1	A	G	G	G	G	G	G	G	G	G	G	G	G	G	G	A	G	G	G
pks12-1	G	A	G	G	G	G	G	G	G	G	A	A	A	A	A	A	A	A	A
pks6b	T	G	T	T	T	T	T	T	T	T	G	G	G	G	G	G	G	G	G
PPE21	A	G	A	A	A	A	A	A	A	A	G	G	G	G	G	G	G	G	G
rhlE	T	C	T	T	T	T	T	T	T	T	C	C	C	C	C	C	C	C	C
speE	G	A	G	G	G	G	G	G	G	G	A	A	A	A	A	A	A	A	A
sseA	G	A	G	G	G	G	G	G	G	G	A	A	A	A	A	A	A	A	A
Mb0353	A	-	A	A	A	A	A	A	A	A	-	-	-	-	-	-	-	-	-
Mb0393	A	C	A	A	A	A	A	A	A	A	C	C	C	C	C	C	C	C	C
mmpL4	C	T	C	C	C	C	C	C	C	C	T	T	T	T	T	T	T	T	T
Mb0849-2	A	G	A	A	A	A	A	A	A	A	G	G	G	G	G	G	G	G	G
Mb0899c	A	G	A	A	A	A	A	A	A	A	G	G	G	G	G	G	G	G	G
Mb1013	G	A	G	G	G	G	G	G	G	G	A	A	A	A	A	A	A	A	A
Mb1365c	C	T	C	C	C	C	C	C	C	C	T	T	T	T	T	T	T	T	T
ndh-2	G	A	G	G	G	G	G	G	G	G	A	A	A	A	A	A	A	A	A
bcp	G	A	G	G	G	G	G	G	G	G	A	A	A	A	A	A	A	A	A
glpK-3	T	T	T	T	T	T	T	T	T	T	A	T	T	T	T	T	T	T	T
Mb2204c	A	C	A	A	A	A	A	A	A	A	C	C	C	C	C	C	C	C	C
plsB2	T	C	T	T	T	T	T	T	T	T	C	C	C	C	C	C	C	C	C
Mb2661	C	G	C	C	C	C	C	C	C	C	G	G	G	G	G	G	G	G	G
Mb3193	T	C	T	T	T	T	T	T	T	T	C	C	C	C	C	C	C	C	C
Mb3478	C	A	C	C	C	C	C	C	C	C	A	A	A	A	A	A	A	A	A
fusA2b-2	CA	A	CA	CA	CA	CA	CA	CA	CA	CA	- A	- A	- A	- A	- A	- A	- A	- A	- A
metK-1	A	G	A	A	A	A	A	A	A	A	G	G	G	G	G	G	G	G	G
Mb1707	C	G	C	C	C	C	C	C	C	C	G	G	G	G	G	G	G	G	G
otsB1-1	T	C	T	T	T	T	T	T	T	T	C	C	C	C	C	C	C	C	C
gkpK-2	C/-	C/-	C	-	C	-	C	C	-	-	-	-	C	-	C	C	C	C	-

MIRU, mycobacterial interspersed repetitive unit; VNTR, variable number tandem repeat; SNP, single-nucleotide polymorphism.

SCG analysis

SNP cluster grouping of the 48 strains resulted in six clusters, which were divided into two groups according to the spoligotype: SB1040 type into SCGs I and II, and SB0140 type into SCGs III–VI (Fig. 1, Tables 4 and 5). The SCGs divided all VNTR types into two or more subtypes, except for the VNTR F type, which was not divided based on the SCG (one VNTR GG type was excluded). In particular, the VNTR H type, which was the dominant type, was divided into three SCGs, demonstrating the potential for more detailed typing based on SNPs.

Fig. 1

UPGMA tree of virulence factor SNPs based on an alignment of only informative, validated SNPs. Mycobacterium bovis AN5, wild-type (D-145, B-5624, B-5335, B- 4600, B-4559, B-4216, W-1366, B-1595, B-1567, 15B-3, D-174, D-171, D-133, D-35, B-5351, B-5281, B-5260, B-5259, B-5253, B-4412, B-3561, B-3222, B-3221, B-3212, B-1579, B-5330, B-5315, B-5286, B-5282, B-5240, B-4724, 16B-4, 15W-442, D-155, D-151, D-136, B-5278, B-5276, B-4452, 15W-23, W-1171, D-166, D-132, B-5648, B-5628, and B-5343), and Mycobacterium tuberculosis H37Rv.

SNP, single-nucleotide polymorphism.

Table 4

Nucleotide information for Korean Mycobacterium bovis field isolates in SCGs in an UPGMA tree

SNP locus	Nucleotide
	SCG-1	SCG-2	SCG-3	SCG-4	SCG-5	SCG-6
corA	T	T	C	C	C	C
glpK	C	C	G	G	G	G
fadD9-2	G	G	A	A	A	A
fadE20(2)	C	C	T	T	T	T
fadE27	G	G	A	A	A	A
glnB	G	G	A	A	A	A
HisD	A	A	G	G	G	G
HisD-1	G	G	G	G	A	G
pks12-1	G	G	A	A	A	A
pks6b	T	T	G	G	G	G
PPE21	A	A	G	G	G	G
rhlE	T	T	C	C	C	C
speE	G	G	A	A	A	A
sseA	G	G	A	A	A	A
Mb0353	A	A	-	-	-	-
Mb0393	A	A	C	C	C	C
mmpL4	C	C	T	T	T	T
Mb0849-2	A	A	G	G	G	G
Mb0899c	A	A	G	G	G	G
Mb1013	G	G	A	A	A	A
Mb1365c	C	C	T	T	T	T
ndh-2	G	G	A	A	A	A
bcp	G	G	A	A	A	A
glpK-3	T	T	T	A	T	T
Mb2204c	A	A	C	C	C	C
plsB2	T	T	C	C	C	C
Mb2661	C	C	G	G	G	G
Mb3193	T	T	C	C	C	C
Mb3478	C	C	A	A	A	A
fusA2b-2	CA	CA	- A	- A	- A	- A
metK-1	A	A	G	G	G	G
Mb1707	C	C	G	G	G	G
otsB1-1	T	T	C	C	C	C
glpk-2	-	C	-	-	C	C

SCG, single-nucleotide polymorphism cluster group.

Table 5

Subtype profiles classified according to spoligotypes, MIRU-VNTR types, and SCGs

Field isolates	Spoligotype	VNTR type	SCG	Sub lineage group
D-145	SB0140	GG	SCG 3	Sub lineage 1
B-5624	SB0140	H	SCG 3	Sub lineage 2
B-5335	SB0140	H	SCG 3	Sub lineage 2
B-4600	SB0140	H	SCG 3	Sub lineage 2
B-4216	SB0140	H	SCG 3	Sub lineage 2
W-1366	SB0140	H	SCG 3	Sub lineage 2
B-4559	SB0140	C	SCG 3	Sub lineage 3
B-1595	SB0140	H	SCG 4	Sub lineage 4
B-1567	SB0140	H	SCG 6	Sub lineage 5
D-174	SB0140	H	SCG 6	Sub lineage 5
D-171	SB0140	H	SCG 6	Sub lineage 5
D-133	SB0140	H	SCG 6	Sub lineage 5
D-35	SB0140	H	SCG 6	Sub lineage 5
B-5351	SB0140	H	SCG 6	Sub lineage 5
B-5281	SB0140	H	SCG 6	Sub lineage 5
B-5260	SB0140	H	SCG 6	Sub lineage 5
B-4412	SB0140	H	SCG 6	Sub lineage 5
B-3561	SB0140	H	SCG 6	Sub lineage 5
B-1579	SB0140	H	SCG 6	Sub lineage 5
15B-3	SB0140	J	SCG 6	Sub lineage 6
B-5253	SB0140	J	SCG 6	Sub lineage 6
B-3222	SB0140	F	SCG 6	Sub lineage 7
B-3221	SB0140	F	SCG 6	Sub lineage 7
B-3212	SB0140	F	SCG 6	Sub lineage 7
B-5259	SB0140	C	SCG 6	Sub lineage 8
B-5330	SB0140	J	SCG 5	Sub lineage 9
B-5315	SB0140	J	SCG 5	Sub lineage 9
B-5286	SB1040	P	SCG 1	Sub lineage 10
B-5240	SB1040	P	SCG 1	Sub lineage 10
B-5282	SB1040	P	SCG 1	Sub lineage 10
16B-4	SB1040	P	SCG 1	Sub lineage 10
B-4724	SB1040	Y	SCG 1	Sub lineage 11
D-155	SB1040	Y	SCG 1	Sub lineage 11
D-151	SB1040	Y	SCG 1	Sub lineage 11
15W-442	SB1040	FF	SCG 1	Sub lineage 12
D-136	SB1040	FF	SCG 1	Sub lineage 12
B-5648	SB1040	P	SCG 2	Sub lineage 13
B-5278	SB1040	FF	SCG 2	Sub lineage 14
15W-23	SB1040	FF	SCG 2	Sub lineage 14
W1171	SB1040	FF	SCG 2	Sub lineage 14
D-132	SB1040	FF	SCG 2	Sub lineage 14
D-166	SB1040	FF	SCG 2	Sub lineage 14
B-4452	SB1040	Y	SCG 2	Sub lineage 15
B-5276	SB1040	Y	SCG 2	Sub lineage 15
B-5343	SB1040	Y	SCG 2	Sub lineage 15
B-5628	SB1040	EE	SCG 2	Sub lineage 16

MIRU, mycobacterial interspersed repetitive unit; VNTR, variable number tandem repeat; SNP, single-nucleotide polymorphism.

Application of SNP analysis for epidemiological tracing

To confirm the applicability of epidemiological analysis using the SNP sites discovered in the present study, we used isolates obtained from five farms with TB outbreaks in Pocheon-si, Gyeonggi-do, Korea to compare the epidemiological analysis capabilities of SNP genotyping, spoligotyping, and VNTR typing. In the Pocheon-si area, there is a strong concentration of livestock farms and, thus, molecular epidemiological analysis in this area can be used to establish methods for controlling TB, based on the understanding of the influx path and TB transmission between farms (Fig. 2).

Fig. 2

Farms where Mycobacterium bovis-infected cattle were found in Pocheon-si, Gyeonggi-do, Korea. (A) TB propagation pathway results when epidemiological tracking was performed based on variable number tandem repeat type. (B) TB propagation pathway results when epidemiological tracking was performed based on single-nucleotide polymorphism subtype.

TB, bovine tuberculosis.

When conventional VNTR typing was used for epidemiological tracking of TB transmission among the five farms, the VNTR type from a deer suspected of having escaped from a farm (D-133) and those of farm E (deer farm) and farms B and C (nearby cattle farms) were identical, indicating that the transmission was caused by a deer, with the outbreak originating from farm E. However, farm A, which is 6 km away from the other farms and has no direct transmission factors, had the same VNTR types as farms B, C, and E, complicating the epidemiological tracking of transmission factors (Fig. 2A). In contrast, SNP typing showed that the VNTR type (H) of farms A, B, C, and E was divided into different sub-lineages (2, 5) with no direct correlation. Farms C and D, both of which had the VNTR P type, were confirmed to have the same sub-lineage upon SNP typing, confirming the possibility of direct transmission between these two farms (Figs. 2B and 3).

Fig. 3

Radial phylogenetic tree of virulence factor SNPs based on an alignment of only informative, validated SNPs, with spoligotype Mycobacterium bovis AN5, field isolates, and Mycobacterium tuberculosis H37Rv.

SNP, single-nucleotide polymorphism.

DISCUSSION

The present study aimed to identify SNPs that enable the demarcation of Korean M. bovis isolates, with the long-term goal to devise a detailed SNP typing method for epidemiological applications. Through SNP screening using WGS, we identified approximately 1,400 SNPs. From these, 34 SNP sites, primarily located in major Mycobacterium virulence genes, were selected for epidemiological analysis of Korean M. bovis field isolates. Some of the SNP sites identified and studied herein have been previously studied via previously sequenced M. bovis strains [722]. However, since the present study aimed to identify markers for epidemiological analysis of Korean M. bovis strains, independent SNPs were selected and only SNPs representing differences between Korean field isolates were analyzed. We confirmed that some of the previously reported SNPs did not show any differences in Korean field isolates.

An SNP-based method for phylogenetic analysis and subclassification of the MTC was established previously [23]. However, SNP typing in M. bovis field isolates has not been performed previously in Korea; thus, the identity of the genetic lineages contributing to the bTB outbreak remain unknown. To this end, we aimed to discover SNPs that can be used for the systematic analysis of M. bovis isolates in Korea. Ducey et al. [24] developed a multilocus genotyping assay to differentiate Listeria monocytogenes strains to below the species level, and a sub-lineage type was identified and applied to outbreak detection and epidemiological investigations. Based on multi-virulence locus sequence typing of 49 epidemiologically unrelated isolates using 28 SNP sets covering six virulence genes, a minimum of 16 SNPs that could distinguish the epidemic clonal lineages from all unrelated lineages were identified. Furthermore, it has been suggested that SNPs in virulence genes can be used as molecular markers for epidemiological investigations [25]. Liu et al. [26] used a set of 12 highly informative SNPs identified from a multilocus sequence typing database to genotype and trace the geographical origins of E. coli isolates of various subtypes. To date, various studies have been conducted for epidemiological or pathological SNP analysis in M. bovis, and several “significant” SNP sites show differences between countries or regions [717]. However, in this study, for detailed epidemiological analysis through SNPs, we separately extracted SNPs showing significant changes between regions or countries, confirmed their variation patterns, and proposed an analysis method using meaningful SNPs between regions.

Globally, spoligotypes have been identified through epidemiological analysis of strains collected from farms with TB outbreaks, enabling tracing of M. bovis infections at the strain level [27]. Spoligotyping has also been used for molecular epidemiology tracking [2829]. In Europe, spoligotyping and VNTR typing have been combined for a more detailed strain classification to determine the epidemiological relationships underlying bovine TB transmission [303132]. In Korea, only spoligotypes SB1040 and SB0140 have been confirmed, while a small number of MIRU-VNTR types are predominant, and are thus generally limited in revealing epidemiological relationships between herds, despite several genotypes being identified. In the present study, SNP typing confirmed the genetic identity of the two spoligotypes that exist in Korea, and they were further divided into different subtypes based on the SNPs identified. Nine previously identified MIRU-VNTRs were subdivided into 16 subtypes, enabling a more detailed epidemiological tracking (Table 6). A more detailed typing and grouping can even be achieved by combining spoligotyping and MIRU-VNTR typing with SNP typing. In addition, although there are geographic differences, the same spoligotypes and VNTR profiles were subdivided into SNPs to clearly distinguish the propagation patterns according to geographic distances, showing that more detailed regional propagation patterns could be confirmed, and the possibility of using this method for epidemiological analysis was shown. When used for epidemiological analysis, this method can overcome the difficulties of epidemiological tracking caused by the limited spoligotypes and MIRU-VNTR types in Korea. In foreign studies, the transmission of a single clone to a specific genotype in wild animals has been reported [16]. In Korea, wild boar, deer, and badger have been reported as the major wild animals and as potential reservoirs for the transmission of M. bovis [33]. For the results of transmission between farms and wild animals in Pocheon, Gyeonggi-do, this study confirmed the possibility of wild animals serving as a medium for transmitting tuberculosis to farms, but host-specific SNPs in wild animals were not observed.

Table 6

Mycobacterium bovis strains suspected to be epidemiologically related

No	Sample ID	Animal species	Herd	VNTR profiles		Spoligotype	SNP type
1	B-4216	Dairy cattle	A*	3 2 5 3 2 7 5 4 10 4 3 4 10 3 4 3	H	SB0140	Sub lineage 2
2	B-4412	Dairy cattle	B*	3 2 5 3 2 7 5 4 10 4 3 4 10 3 4 3	H	SB0140	Sub lineage 5
3	B-5281	Dairy cattle	C*	3 2 5 3 2 7 5 4 10 4 3 4 10 3 4 3	H	SB0140	Sub lineage 5
4	B-5282	Dairy cattle	C	3 2 5 3 3 6 5 4 10 4 3 3 10 10 5 2	P	SB1040	Sub lineage 10
5	B-5286	Korean cattle	D*	3 2 5 3 3 6 5 4 10 4 3 3 10 10 5 2	P	SB1040	Sub lineage 10
6	D-173	Elk	E*	3 2 5 3 2 7 5 4 10 4 3 4 10 3 4 3	H	SB0140	Sub lineage 5
7	D-133	Elk	†	3 2 5 3 2 7 5 4 10 4 3 4 10 3 4 3	H	SB0140	Sub lineage 5

VNTR, variable number tandem repeat; SNP, single-nucleotide polymorphism.

*A–E, farms in Pocheon-si, Gyeonggi-do, Korea; †Suspected to have escaped from a certain farm in the region.

In the present study, SNP sites in M. bovis were identified using standard PCR, which is time-consuming when evaluating large numbers of SNP sites. The methods used in the present study are unsuitable for large-scale, country-wide screening of M. bovis lineages, and techniques that allow for faster SNP typing are urgently needed. High-resolution melt (HRM) analysis is a real-time PCR technique useful for detecting SNPs based on altered melting temperatures of PCR products [34]. Landolt et al. [3536] applied HRM to SNP typing of MTC and successfully identified a rare subtype of M. canettii in clinical samples. Cancino-Muñoz et al. [37] reported accurate distinction between MTC strains and identified a sub-lineage based on SNP typing using HRM in low-concentration DNA samples. Recently, loop-mediated isothermal amplification (LAMP) has attracted increasing attention for SNP analysis. LAMP uses two pairs of primers targeting six regions on a target sequence and a mesophilic DNA polymerase with strand-displacement activity to achieve specific amplification [38]. LAMP-based technologies have successfully achieved accurate and fast SNP analysis. In addition, allele-specific LAMP, one-step strand displacement-coupled LAMP, and LAMP combined with allele-selective oligonucleotide hybridization have been applied to SNP analysis [39]. The combination of these and current SNP typing methods would enable extremely fast and accurate SNP typing of M. bovis and other important pathogens.

In conclusion, we identified several SNPs in virulence genes that can be applied for detailed M. bovis typing for epidemiological analysis. The SNPs identified here may not be sufficient to establish a typing method, and additional screening of significant SNPs should be performed. Moreover, further studies are needed to develop a technique for screening large numbers of SNP candidates. In future, we plan to screen for additional significant SNP sites in genes other than virulence genes. In addition, we will continue our efforts to establish a fast and accurate HRM- or LAMP-based SNP typing method for epidemiological tracking of pathogens in Korea.