Genotyping Brucella canis isolates using a highly discriminatory multilocus variable-number tandem-repeat analysis (MLVA) assay

Differentiation of Brucella canis from other Brucella species are mainly performed through PCR-based methods and multilocus variable-number tandem-repeat (VNTR) analysis (MLVA) procedures. Both PCR-based and MLVA methods are limited in discriminating B. canis strains. A new MLVA-13Bc method for B. canis genotyping was established by combining eight newly-developed VNTRs with five published ones. During 2010 and 2016, 377 B. canis PCR-positives were identified from 6,844 canine blood samples from 22 U.S. states, resulting in 229 B. canis isolates. The MLVA-13Bc method was able to differentiate each of these 229 isolates. The Hunter-Gaston Discriminatory Index of the individual VNTR loci ranged from 0.516 to 0.934 and the combined discriminatory index reached 1.000. Three major clusters (A, B and C) and 10 genotype groups were identified from the 229 B. canis isolates. Cluster A mainly contains genotype groups 1 and 2, and a few group 3 isolates; nearly all Cluster B isolates were from group 6; other genotype groups were classified into Cluster C. Our newly developed MLVA-13Bc assay is a highly discriminatory assay for B. canis genotyping, and can serve as a useful molecular epidemiological tool, especially for tracing the source of contamination in an event of a B. canis outbreak.

for the genotyping of Brucella species and biovars [21][22][23] . Bricker et al. 21 used a MLVA method named HOOF-Prints with eight tandem repeat (TR) loci (MLVA-8) and was able to differentiate Brucella isolates at both the species and biovar levels. A rather comprehensive screening of 107 TRs identified 15 TRs (MLVA-15) that were more informative and were able to differentiate most Brucella species and for some strains even at the biovar level 22 . Al Dahouk et al. 24 added another TR into the MLVA-15 method to form a new assay called MLVA-16 and was used to study genetic diversity of 128 human B. melitensis strains. The study identified 110 genotypes that differentiated most of the 128 strains, yet the MLVA-16 method provided much lower discriminatory power against the eight B. canis and 18 B. ovis strains in the study. Using MLVA-16, B. suis biovar 1, B. suis biovar 2, B. abortus, B. melitensis and B. ovis can be clearly identified. B. suis biovar 5, B. neotomae and the marine mammal strains are closely related strains, and they can be differentiated by this method. However, this method appeared to be insufficient in differentiating strains within the species of B. canis in different studies, especially for those that were from closely related geographic regions 22,24,25 . Whatmore et al. 23 grouped 121 Brucella isolates into 119 genotypes based on 21 VNTR loci. The approach based on 21 VNTR loci provided better strain genotyping information for B. abortus, B. melitensis and B. suis, but was less informative in differentiating strains of B. canis, B. ovis and B. neotomae. Therefore, the goal of this study were: 1) to develop and validate a MLVA genotyping method for the differentiation of B. canis strains; and 2) to study genetic diversity of a collection of 229 B. canis isolates collected from the US in recent years using the newly developed MLVA-13Bc method.

Results
PCR identification and isolation of B. canis strains from canine blood samples. From a total of 6,844 canine blood samples, the duplex diagnostic PCR identified 377 Brucella positives, and all Brucella positives were confirmed to be B. canis strains. The MLVA-13Bc PCR amplifications for 10 genotype groups visualized by QIAxcel were shown in Figure S1. Among the 377 PCR-positive samples, 229 B. canis isolates were obtained. Selected positive samples (n = 20) were further verified by sequencing a region flanking a 976 bp fragment that is deleted only from the B. canis genome, and is intact in all other Brucella species 26,27 . Positive samples were observed in 10 states, mainly from the Midwest region, of the US with an average positive rate of 5.5% (377/6,844). The positive rate for the 6 Figure S2). Based on the average size range and primer annealing temperatures of each VNTR locus, 12 of them were grouped into 6 duplex PCR reactions, namely Reaction 1 (for loci BCTR09 and BCTR06), 2 (BCTR12 and Bruce07), 3 (BCTR03 and Bruce16), 4 (BCTR02 and Bruce09), 5 (BCTR01 and BCTR08), and 6 (BCTR11 and Bruce04), respectively. Reaction 7 is a singular one for locus Bruce18 ( Table 2). All primers were specific to targeted regions and did not amplify other regions of the B. canis genome. The results were compared to individual PCR reactions. No interference was observed while comparing the singular and duplex PCR reactions. There were no overlapping amplicons observed between the 2 loci within each duplex reaction. The result of Sanger sequencing (data not shown) was consistent with the size obtained by the QIAxcel Advanced System.
Genotyping 229 B. canis isolates using MLVA-13Bc method. The newly established method ( Table 2) was then used for genotyping the remaining 229 isolates. Copy numbers of the tandem repeats for each locus were compiled to indicate the genotype for each isolate, and then used to further group the isolate into one of the 10 genotype groups ( Figure S2). Each of the 229 isolates can be differentiated by this newly-developed MLVA-13Bc (MLVA-13 for B. canis) method. To evaluate the discriminatory power of the selected loci, the Hunter and Gaston 28 discrimination index (HGDI) was calculated for the MLVA-13Bc method as a combined method, and for each of the 13 loci used in this study. Among these 229 isolates, the newly developed MLVA assay allowed the classification of all isolates to unequivocal genotypes (HGDI = 1.000). The HGDI value ranged from 0.516 to 0.934 among the 13 individual VNTR loci ( Table 3). Genotyping of the 229 B. canis isolates by the newly developed MLVA-13Bc method identified 10 genotype groups (group 1 -group 10), with the similarity cutoff at 20% ( Figure S2). The US strain 2010009751 was classified into group 1; the Korean strain HSK A52141 was in group 2; US strains ATCC 23365, RM6/66 and 2009004498 were classified into group 3; US strain 2009013648 and the South Africa strain F7/05 A were in group 4; the Argentina strain CNGB 1342 was in group 6; and the Sweden strain SVA13 was in group 9. Group 3 is the largest group with 76 isolates widely spread in 8 states, and it is the predominant group in Colorado (100%) and Kansas (71%). The second largest group, group 6, had strains only from Ohio (71%), Indiana (58%) and Missouri (3%). The 33 isolates in group 2 were mostly obtained from Oklahoma (30); only two from Kansas and one from Missouri. The 22 isolates in group 1 were from Iowa, Missouri, Indiana and Ohio, while the 21 isolates in group 7 were from Idaho, Kansas, Minnesota, Missouri and Indiana, which were widely spread in the US. There were no more than 10 isolates in groups 4, 5, 8, 9 and 10, and each of these groups was observed in only one or two states ( Figure S2). Distribution of the 10 genotype groups were also illustrated in a cluster analysis (Insert in Fig. 1).     (5) as well as most isolates in group 3 (64/76) and 5 (8/9) belonged to Cluster C. The five isolates each in groups 8 and 9 were located in Cluster A and C, respectively (Insert in Fig. 1). The MLVA-16 assay 24 is probably the most used method for genotyping different Brucella spp., and for strains within certain Brucella spp. including B. canis. As the method was not designed specifically to differentiate strains within B. canis, it provided limited discriminatory power against B. canis strains. For example, MLVA-16 differentiated the 29 Chinese B. canis strains into 21 genotypes, as reported by Di et al. 25 . Kang et al. 29 classified 77 Korean B. canis strains into 30 genotypes using the same procedure. As we are using a different procedure, those strains could not be in silico analyzed except for the nine reference strains that have full or near full genome sequences. Our MLVA-13Bc assay differentiated each of the 229 strains used in this study. In comparison, the five most informative loci from the MLVA-16 assay classified the 229 strains into 131 genotypes. To further confirm that the remaining 11 VNTRs in the MLVA-16 assay would not generate useful information for B. canis strains, we have tested the 11 primer pairs on each of the genotype groups, and resulted identical amplicon sizes ( Figure S3). Our data from this study together with other published data suggested that the remaining 11 VNTRs from the MLVA-16 assay do not provide useful discriminatory power for B. canis genotyping.

Distributions of B. canis isolates in different dog breeds and genders.
A survey on canine infection with B. canis was carried out in 22 U.S. states between 2010 and 2016, and the result indicated that the total average positive rate is 5.5% (377/6,844; Table 1). Although B. canis can cause human infections, information of the zoonotic potential of these 229 strains are not available. Due to not all positive samples in years 2010 and 2015, and none in 2014 were subjected to bacterial isolation, the more realistic B. canis isolation rate is reflected by data in years 2011, 2012, 2013 and 2016, which was 79.5% (194/244). This data demonstrated a good correlation between the culture and the molecular detection methods, yet about 20% of PCR positives would not be detected by the traditional culture method, indicating PCR can be more sensitive in B. canis detection.
Dog breed information was available for about half of our strain collection. We identified 24 different dog breeds from 118 out of 229 B. canis isolates. Although we did not see a clear-cut distribution, isolates from certain dog breeds were primarily located in a given cluster (Table S1). Clusters A and B are smaller clusters. Isolates from German Shepherd, French Bulldog, Siberian Husky and Bichon Frise were mostly grouped in Cluster A, and all four isolates from Cocker Spaniel and five out of 20 isolates from Shih Tzu were in Cluster B. We also have dog age information for 60 samples, but there was no correlation between dog age and isolate genotypes observed (Table S1).
There were 133 positive samples (Supplementary Table S1) that were indicated with animal genders. Although we received more samples from females (n = 87) than males (n = 46), it is worth mentioning that we received a significant portion from B. canis-barring male animals that has not been extensively documented, as most of the reports were on abortion cases 9, 30 . Our data indicated that the organism can circulate in blood streams of both Our MLVA genotyping data on 229 B. canis strains may be the largest genotyping study for B. canis isolates in the US. Although the distribution of different genotypes in each state was not clear-cut, some interesting information was revealed. Genotype group 6 was the predominant genotype for both Indiana and Ohio, which are neighbor states. Group 1 isolates were mainly distributed in Missouri, Iowa and Indiana, states that are geographically related. Group 2 isolates were mainly found in Oklahoma and some in Missouri, which are also bordering states. Groups 3 and 7 though were widely distributed in different states.
In a temporal distribution analysis for these 229 B. canis isolates, all 18 samples collected in 2016 were in Cluster C. Most isolates in 2015 (12/13) were in Cluster C, only one was in Cluster A. Most isolates collected from 2011 (28/38) were in Cluster A, some (9/38) were in Cluster C and only one in Cluster B.
In conclusion, most currently-used MLVA assays (MLVA-8, MLVA-10, MLVA-11, MLVA-15 and MLVA-16) were not specifically designed for B. canis at the strain level, although they are good molecular genotyping methods of Brucella at the species, biovar or at strain level for certain Brucella species. Our newly developed MLVA-13Bc method was specifically designed to discriminate strains from within the species of B. canis, and had higher discriminatory power than previous assays 25,31,32 . Although our data is not associated with pathogenesis as such information is not available, the new MLVA-13Bc method genotyped each of the 229 isolates in our collection, indicating that it is a useful molecular epidemiological tool to trace the source of strain introduction in an event of a B. canis outbreak.  (Table 1). All samples were diagnostic samples submitted by clients and no experimental samples were used. All procedures were carried out strictly following institutional biosafety procedures. Samples were freshly collected in sodium citrate (blue-top) or EDTA (purple-top) blood collection tubes, and delivered to KSVDL on ice for next-day delivery. Upon receiving, samples were accessioned and 1 ml blood was transferred into 5 mL Brain Heart Infusion (BHI) broth (HARDY Diagnostics, Santa Maria, CA, USA), mixed, and kept at −80 °C for 3 h (for samples received in the mornings) or at −20 °C for overnight (for those received in the afternoons), then incubated at 37 °C aerobically in a shaker incubator for 40-48 h prior to DNA extraction. A duplex PCR assay was used for B. canis identification. Isolates used in this study were confirmed by culture isolation followed by PCR confirmation on single colonies, and by sequencing on selected isolates (detailed below). Isolation of B. canis from PCR positive samples. Culture isolation of B. canis was achieved by streaking the BHI culture on 5% blood agar plates (HARDY Diagnostics). The inoculated plates were incubated at 37 °C aerobically for 72-120 h. Single colonies morphologically consistent with Brucella were sub-cultured in BHI broth for 24 h, and confirmed by PCR. All isolates confirmed as B. canis were stored with 15% glycerol at −80 °C for long-term use.

Development of the Multiple Locus Variable-Number Tandem Repeat Analysis method for B. canis genotyping.
Published data using MLVA- 15 22 or MLVA-16 24, 29 methods were analyzed with an emphasis on their ability to differentiate B. canis strains. The purpose of analysis on published data was to select those VNTRs that generated polymorphism on B. canis strains and to compile them with new VNTRs (described below) to form a new MLVA method that can be used for genotyping of B canis strains with increased discriminatory power. An effort of identifying additional VNTR loci was performed by searching B. canis whole genome sequences (chromosome I: NC_010103, NC_016778, NZ_CP007758, NZ_CP007629; chromosome II: NC_010104, NC_016796, NZ_CP007759, NZ_CP007630) using Tandem Repeats Finder described by Benson 33 . Additional VNTR loci were selected based on following criteria recommended by Nadon et al. 34 : all loci had no insertions and deletions in the tandem repeats; each tandem repeat was equal or greater than 8 bp; and the loci had conserved flanking sequences from which the PCR primers can be identified. Primers flanking each locus were designed using the online primer design software, Primer 3 (http://bioinfo.ut.ee/primer3-0.4.0/). Twenty B. canis isolates selected from different dog breeds and different colleting states spanning different years were used to evaluate and optimize the method. All individual VNTR were run on the 20 isolates first, and the optimized protocol was used for the remaining isolates. MLVA data analysis. The genotype of each isolate was identified by compiling the genotype generated by each VNTR locus. The phylogenetic tree was generated using an unweighted pair group method with arithmetic mean (UPGMA) method 35 that is embedded in BioNumerics version 7.6 (Applied Maths, Austin, TX). The cluster analysis was performed using the UPGMA with a minimum spanning tree (MST) and distance matrices for categorical data by BioNumerics version 7.6. The Geographic distribution map template was downloaded from www.vectortemplates.com (Graphics Factory CC, Western Cape, Sourth Africa) with permission, and the added pie charts were constructed using BioNumerics version 7.6. Hunter and Gaston discrimination index (HGDI) was calculated using Ridom EpiCompare software version 1.0 (www.ridom.de/epicompare/) to elucidate the discriminatory power of the genotyping methods, which explained the probability of two unrelated and different isolates sampled from the test population grouping as different subtypes by a specific typing method 28 . The 95% confidence intervals (CI) were calculated according to the method previously described 36 . VNTR loci stability of the B. canis genome. In order to investigate the stability of the 13 VNTR loci in the B. canis genome, three B. canis isolates collected from 2011 (D11-106637 #69), 2013 (D13-123417 #42), and 2016 (D16-008400 #3) were subcultured by inoculating 50 μL of stock culture into 5 ml BHI medium and incubated at 37 °C for 2 days. A new 5 ml BHI broth was inoculated with 50 μL of fresh culture, and incubated at 37 °C for 2 days. A total of 20 consecutive subcultures were made for each strain. Portions of each subculture were frozen at −20 °C till use. When 20 subcultures were completed, the original culture, subculture number 10 and number 20 were used for DNA extraction and PCR amplification with the newly developed genotyping method.