Improvement of the MLVA typing scheme for  Burkholderia mallei  strains

Margarita L. Ledenyova; Леденева Маргарита Леонтьевна; Olga S. Bondareva; Бондарева Ольга Сергеевна; Galina A. Tkachenko; Ткаченко Галина Александровна; Dimitriy V. Ustinov; Устинов Димитрий Владимирович; Irina B. Zakharova; Захарова Ирина Борисовна

doi:10.36233/0372-9311-720

Improvement of the MLVA typing scheme for Burkholderia mallei strains

Authors: Ledenyova M.L.¹, Bondareva O.S.¹, Tkachenko G.A.¹, Ustinov D.V.¹, Zakharova I.B.¹
Affiliations:
1. Volgograd Plague Control Research Institute
Issue: Vol 102, No 5 (2025)
Pages: 615-625
Section: ORIGINAL RESEARCHES
URL: https://microbiol.crie.ru/jour/article/view/18873
DOI: https://doi.org/10.36233/0372-9311-720
EDN: https://elibrary.ru/LXJQDU
ID: 18873

Cite item

Full Text

Abstract
Full Text
About the authors
References
Supplementary files
Statistics

Abstract

Introduction. The registration of sporadic cases of glanders in horses in Russia, caused by Burkholderia mallei, highlights the importance of developing genotyping algorithms for this pathogen. The MLVA method (multilocus-variable tandem repeat analysis), based on a comparative analysis of the number of variable tandem repeats (VNTRs), remains a promising genotyping tool. As the number of whole-genome sequences in international databases increases, the informative value of VNTR loci changes, necessitating a revisiting of existing typing schemes.

The aim of this study was to assess the feasibility of including the VNTR locus BPSS1974^#I in the MLVA-6 scheme for intraspecies differentiation of B. mallei.

Materials and methods. The study of 64 strains of B. mallei was conducted in silico and in vitro using MLVA, differentiating region amplification, whole-genome sequencing, and bioinformatic analysis.

Results. Genotyping B. mallei using the MLVA-6 scheme failed to determine in silico VNTR profiles of 13 strains from the GenBank NCBI database for one or more loci due to low read coverage of the corresponding genomic regions or their complete absence (null alleles). The effective number of alleles (ne) and the polymorphic information content (PIC) index for the MLVA-6 scheme loci ranged from 3.842–8.103 and 0.740–0.877, respectively. The potential for including the VNTR locus BPSS1974^#I in this scheme was determined based on the molecular stability of the motif within it and a high values for ne and PIC, which were 4.299 and 0.767, respectively. VNTR profiles of the collection strains at locus BPSS1974^#I were identical to the corresponding strains in the GenBank database. The results of the cluster analysis using a combined MLVA-6 scheme and the BPSS1974^#I locus were consistent with the phylogenetic reconstructions obtained using other molecular genetic methods.

Conclusion. The VNTR locus BPSS1974^#I can be considered a marker, the inclusion of which in the MLVA-6 scheme will improve the accuracy of genotyping and the determination of the regions of origin of newly isolated B. mallei strains.

Keywords

Burkholderia mallei, glanders, genotyping, MLVA, tandem repeats

Full Text

Introduction

The Gram-negative bacterium Burkholderia mallei is the causative agent of glanders, a potentially fatal anthropozoonosis that primarily affects equids: horses, mules and donkeys [1, 2]. Human illness is linked to their professional activities and occurs through close contact with infected animals, for example, among farm and veterinary workers [3]. There are known cases of infection among employees of bacteriological laboratories who work with B. mallei cultures [4, 5]. In Russia B. mallei is classified as pathogenicity group II (hazard) agent, while abroad it is considered a potential group B bioterrorism agent, which is due to the possibility of a public health emergency arising from the use of B. mallei as a biological weapon [6–8].

Currently, outbreaks or sporadic cases of glanders occur in several regions of Asia and the Middle East, North Africa, as well as Central and South America [9, 10]. The registration of imported cases of glanders in countries free of this disease poses a threat of infection to local livestock, humans, as well as the possibility of reintroduction of the infection [11, 12].

Cases of glanders are being reported in countries bordering Russia, which necessitates vigilance from veterinary services regarding this disease, especially in the regions of Siberia and the Far East. The outbreak of glanders in horses in 2023 at the "Chita State Stud Farm and Hippodrome named after Khosaen Khakimov," registered by the World Organization for Animal Health¹, highlights the importance of epidemiological monitoring of glanders in Russia.

Modern methods of intraspecies typing allow obtaining a detailed genetic characterization of the pathogen to address issues such as investigating cases of human or animal disease, and determining the evolutionary and phylogenetic relationships between microbial strains [13–15]. However, in the case of glanders, the search for effective DNA targets is complicated by the high conservation of the B. mallei genome. The high density of microsatellite and minisatellite repeats (Variable Number Tandem Repeats, VNTR) in B. mallei genome determines the prospects of using strain differentiation schemes based on multilocus variable number tandem repeat analysis (MLVA). In foreign sources, a scheme of 23 VNTR loci is most often used for this purpose, it was originally proposed for the causative agent of melioidosis — Burkholderia pseudomallei, a clone of which is supposed to be the ethilogical agent of glanders [16]. However, this scheme is not optimized for the B. mallei. A large number of typing loci increases the overall cost of the analysis.

To reduce the labor intensity and cost of the MLVA method, the development of shortened schemes is relevant, allowing the maintenance of the discriminatory power and reliability of molecular epidemiological investigations using fewer loci. The combination of 6 VNTR loci we previously selected allowed the differentiation of B. mallei strains with high discriminatory power [17]. However, the informative value of VNTR loci changes as the number of new whole-genome sequences in international databases increases, which necessitates a review of existing typing schemes.

In a recent study, we demonstrated the effectiveness of using a scheme of 5 VNTR and 4 SNP loci for establishing the geographical origin of B. pseudomallei strains and determining the clonality of isolates in melioidosis cases [18]. The new VNTR locus BPSS1974^#I, proposed in this study, located within the BPS_RS29560 gene region, was characterized by stable amplification and the presence of amplicons in all B. pseudomallei strains examined. The tandem repeat sequence (motif) consisted of 9 nucleotide pairs and had a degenerate structure that was under stabilizing selection. The presence of an orthologous gene in B. mallei genome determines the possibility of using the BPSS1974^#I locus to improve the 6-locus MLVA typing scheme.

The aim of the study was to assess the feasibility of including the VNTR locus BPSS1974^#I in the MLVA-6 scheme for genetic typing of B. mallei strains.

Materials and methods

The research objects were whole-genome sequences of 56 B. mallei strains from the GenBank NCBI database². Our sample included nucleotide sequences of B. mallei strains isolated in endemic countries (India, Pakistan, Brazil), certain parts of Asia and the Middle East where outbreaks or cases occur sporadically (Turkey, China, Myanmar, Bahrain), as well as in USA and some European countries (Hungary, Yugoslavia).

To assess the stability of DNA loci within genotyping systems, 7 strains of B. mallei from the collection of the Volgograd Plague Control Research Institute were used. Duplicates of the strains P-1, Muksuwar-11, V-120, C-4 and C-5, deposited in the State Collection of Pathogenic Microorganisms and Cell Cultures under the numbers SCPM-O-B-4688, SCPM-O-B-7093, SCPM-O-B-7146, SCPM-O-B-4682 and SCPM-O-B-4683, corresponded to the collection strains. In turn, the collection strains Muksuwar-11, Bogor-37 and Zagreb were matched by strains whose nucleotide sequences have been deposited in GenBank NCBI under accession numbers GCA_033870375.1, GCA_033870395.1 and GCA_033870355.1. The study also included the collection strain B. mallei 16050, isolated from a horse in 2023 during an outbreak of glanders at a state stud farm in Chita. The genome sequence of B. mallei strain 16050 has been deposited in the VGARus database (id vnip002404). A total of 64 strains of B. mallei were analyzed in the study.

For setting up polymerase chain reactions for MLVA typing and amplification of differentiating regions of the genome (Different Region Analysis, DFR), DNA was extracted using the RIBO-prep reagent kit (Central Research Institute of Epidemiology, Moscow). Genomic DNA for sequencing was extracted using the Biolabmix-DU-250 kit (Biolabmix) according to the manufacturer instructions.

Primers and amplification parameters for typing using the MLVA-6 scheme and the BPSS1974^#I locus were described in earlier studies [17, 18]. The VNTR profiles of the strains were determined as the set of allelic variants at each locus and are presented as a numerical pattern in the format "L933k/L3145k/L3652k/L20k/L1217k/S2862k/BPSS1974^#I".

DFR genotypes were determined by real-time hybridization-fluorescence PCR (real-time PCR) for 9 fragments: BmVAT1, BmVAT2, BmVAT3, BmVAT4, BmVAT5, BmVAT6, BmVAT7, BmVAT8, BmVAT9, using the oligonucleotide sequences and genotyping scheme described by O.S. Bondareva et al. [19]. The DFR typing results were converted into a binary matrix, where the presence of an amplicon was denoted by "1" and its absence by "0".

Library preparation for whole-genome sequencing was performed using the Nextera XT library prep kit, and sequencing was carried out on the MiSeq platform using the MiSeq Reagent Kit v3 (all from Illumina Inc.). The raw nucleotide reads were processed using the Trimmomatic utility with standard parameters for Illumina. Processed reads were assembled using SPAdes v3.15.4 with the "-isolate" option and standard command-line parameters [20].

Cluster analysis and dendrogram construction were performed using FAMD v. 1.31 [21] and Mega v. 11.0.11 [22] software, employing the nearest neighbor joining algorithm and the Jaccard genetic distance coefficient. The Hunter–Gaston discriminatory index (HGDI) [23] was used to assess the discriminatory power of the typing schemes. To analyze the informativeness of the microsatellite loci, the number of alleles, the effective number of alleles (ne) and the polymorphic information content (PIC) index were determined.

Results

In silico analysis showed that the entire 692 bp amplicon of the BPSS1974^#I locus, located at positions 1627686–1628377 in the genome of the reference strain B. mallei ATCC 23344 (GCA_033956065.1), corresponded to the 1755 bp BMA_RS16575 gene, located on the second chromosome. Much like in the B. pseudomallei genome, the BMA_RS16575 gene is annotated in the NCBI GenBank database as a collagen-like protein containing a triple helix repeat (WPJ46631.1, 584 bp). The forward primer site occupied positions 1627686–1627707 in the B. mallei ATCC 23344 genome and 28–49 in the BMA_RS16575 gene sequence, while the reverse primer site occupied positions 1628359–1628377 and 701–719, respectively. As in the B. pseudomallei genome, positions 3, 6 and 9 of the motif within the BPSS1974^#I locus were under stabilizing selection.

In silico typing of 56 B. mallei strains from the NCBI GenBank database using primers to the BPSS1974^#I locus revealed 10 allelic variants, the most common of which were repeats with copy numbers 56 and 54 — each accounting for 17 strains. 55 and 53 repeats were found in 9 and 5 strains respectively; the number of repeats was 58 for 2 strains, and the remaining 5 variants were unique. In B. mallei SCPM-O-B-4686, primer binding sites were absent (null allele).

Using the MLVA-6 scheme, null alleles were identified at the L1217k locus in 14 strains from the GenBank NCBI database. Analysis of the molecular nature of the null alleles showed that the absence of amplification products at loci BPSS1974^#I and L1217k is due to the deletion of the corresponding genome fragments. Null alleles at the L1217k locus, identified in 4 complete genomes of strains B. mallei 6, 34, BMQ and 23344, were considered as one of the allelic variants in further work. However, the genomes of B. mallei strains BM 1, BM 5, BM 6, Turkey5, NCTC3709, SCPM-O-B-4682, SCPM-O-B-4683, SCPM-O-B-4684, SCPM-O-B-4688, and SCPM-O-B-4686 are represented by a set of contigs of varying lengths, which does not exclude the absence of a PCR product in in silico analysis due to an incomplete nucleotide sequence. Therefore, these strains were excluded from further analysis.

The MLVA-6 scheme did not allow the in silico determination of the genotypes of B. mallei strains BM-1, BM-3, BM-6 and BM-9 at locus L20k and B. mallei 3076 at locus L3652k, as the target VNTR region was located on different contigs. Due to this, these strains were also removed from the study.

Excluding 13 strains from the sample resulted in the loss of 3 unique allelic variants at locus BPSS1974^#I and one at locus L20k. Table 1 presents the B. mallei strains that were subsequently used for cluster analysis.

Table 1. Results of MLVA and DFR typing of B. mallei strains

B. mallei strain name (GenBank or VGARus accession number)	Location and year of isolation	MLVA-profile*	MLVA-profile/ cluster number**	DFR-type***
6 (GCA_000755845.1)	Turkey, 1950	14/7/8/4,86/–/7/50	1/A	17
Ц-4 (vnip002737)****	Mongolia, 1967	11/7/6/17/–/8/54	2/A	06
Ц-5 (vnip002738)****	Mongolia, 1967	11/7/6/17/–/8/54	2/A	06
BMQ (GCA_000755885.1)	India, 1932	2/4/7/9/–/5/54	1/B	12
23344 (GCA_000755865.1)	Myanmar, 1944	2/4/7/29/–/5/54	2/B	12
NCTC 3708 (GCA_003590195.1)	India, 1932	2/6/10/5/7,9/10/55	1/C	24
16050 (vnip002404)****	Russia, Chita, 2023	5/6/11/17/5/10/54	2/C	06
mongolia_1	Mongolia, 2022	9/5/16/17/5/6/55	3/C	06
В-120 (vnip002736)****	Russia, Ulan-Ude, 1985	10/6/8/17/5/6/55	4/C	06
SCPM-O-B-7146 (GCA_003627695.1)	Russia, Ulan-Ude, 1985	10/6/8/17/5/6/55	4/C	06
BM-2 (GCA_028621615.1)	China, Hebei, no data	6/8/10/19/4,9/6/55	1/D	16
2002721277 (GCA_003590185.1)	USA, 1956	10/8/7/17/4,9/12/55	2/D	16
China5 (GCA_000757315.2)	China, 1956	10/8/8/16/4,9/10,4/55	3/D	16
2000031063 (GCA_000756025.2)	Hungary, no data	10/8/8/16/4,9/12/55	4/D	16
P-1 (vnip002735)****	Yugoslavia, no data	6/5/11/7/-/6/55	1/E	09
SAVP1 (GCA_000015465.1)	India, no data	4/5/9/7/9,9/11/55	2/E	04
2000031066 (GCA_003590125.1)	India, no data	4/5/7/7/13,9/12/54	3/E	21
NCTC 10247 (GCA_000762285.1)	Turkey, Ankara, 1960	5/9/12/7/4,9/7/54	1/F	03
34 (GCA_939576165.1)	no data,1972	10/9/16/7/–/6/54	2/F	20
Bahrain1 (GCA_001729545.1)	Bahrain, 2011	2/10/14/15/5,9/6/54	1/G1	22
BM-4 (GCA_028621665.1)	China, Hebei, no data	6/10/10/18/8/13/54	1/G2	23
11 (GCA_000959405.1)	Turkey, 1949	8/10/10/8/4,9/3/54	2/G2	19
NCTC 10229 (GCA_000015605.1)	Hungary, Budapest, 1961	10/10/14/9/4,9/8/54	1/G3	02
2002734299 (GCA_000959165.1)	Hungary, 1961	10/11/14/9/4,9/8/54	2/G3	02
Ivan (GCA_000986905.1)	Hungary, 1961	10/10/13/9/4,9/8/54	3/G3	02
BM_campo 2.1 (GCA_905359435.1)	Brazil, 2016	10/10/13/9/4,9/8/54	3/G3	02
UFAL2 (GCA_905359425.1)	Brazil, 2017	8/4/9/14/4,9/4/27	1/H	06
PRL-20 (GCA_000169875.1)	Pakistan, Lachor, 2005	5/4/6/6/3,9/7/58	1/I	25
India86-567-2 (GCA_000959465.1)	India, no data	6/4/8/6/8,9/11/58	2/I	01
Turkey1 (GCA_002345985.1)	Turkey, no data	7/4/11/6/6,9/7/56	1/J	15
Turkey2 (GCA_002346025.1)	Turkey, no data	7/4/10/6/6,9/7/56	2/J	15
Turkey3 (GCA_002346065.1)	Turkey, no data	7/4/7/6/6,9/6/56	3/J	15
Turkey4 (GCA_002346085.1)	Turkey, 1960	7/4/11/6/6,9/6/56	4/J	15
Turkey6 (GCA_002346125.1)	Turkey, no data	7/4/11/6/6,9/6/56	4/J	15
Turkey7 (GCA_002346145.1)	Turkey, no data	7/4/11/6/6,9/7/56	1/J	15
Turkey8 (GCA_002346165.1)	Turkey, no data	7/4/11/6/6,9/7/56	1/J	15
Turkey9 (GCA_002346185.1)	Turkey, no data	7/4/12/6/6,9/7/56	5/J	15
Turkey10 (GCA_002346005.1)	Turkey, no data	7/4/11/6/6,9/7/56	1/J	15
Muksuwar-11 (vnip002733)****	India, 1979	5/4/8/17/6,9/12/53	1/K	11
SCPM-O-B-7093 (GCA_003627585.1)	India, 1979	5/4/8/17/6,9/11/53	2/K	07
Mukteswar (GCA_033870375.1)	India, 1996	5/4/8/17/6,9/11/53	2/K	11
Zagreb (vnip002734)****	Yugoslavia, no data	5/4/8/15/6,9/6/53	3/K	07
Bogor-37 (vnip002732)****	Indonesia, 1979	5/4/8/15/6,9/11/53	4/K	07
Zagreb (GCA_033870355.1)	Yugoslavia, 1996	5/4/8/15/6,9/11/53	4/K	07
Bogor (GCA_033870395.1)	Indonesia, 1995	5/4/8/15/6,9/11/53	4/K	07
Kweiyang#4 (GCA_001608335.1)	China, 1942	9/7/4/24,71/6,9/15/56	1/L	22
2002721274 (GCA_002522985.1)	USA, 1956	9/7/4/19/6,9/16/56	2/L	01
ATCC 23344 (GCA_033956065.1)	Myanmar, 1942	8/7/4/31/6,9/15/56	3/L	01
JHU (GCA_002346205.1)	USA, Maryland, 2000	8/7/4/27/6,9/15/56	4/L	01
FMH 23344 (GCA_000755785.1)	Myanmar, 1944	8/7/4/27/6,9/15/56	4/L	01
FMH (GCA_002346045.1)	USA, Maryland, 2000	8/7/4/27/6,9/15/56	4/L	01

Note. *The MLVA profile of the strains is recorded as a numerical pattern of the number of repeats in the scheme "L933k/L3145k/L3652k/L20k/L1217k/S2862k/BPSS1974^#I".

** Each MLVA profile is assigned a serial number within the corresponding cluster formed during the construction of the dendrogram.

***DFR-type numbers have been assigned according to previously published data [19].

****Strains from the collection of the Volgograd Plague Control Research Institute.

As a result of VNTR loci amplification of B. mallei collection strains, PCR products were obtained that appear as single bands on the electropherogram. At the same time, the amplification product for the L1217k locus was absent in 3 strains from the collection (B. mallei P-1, C-4 and C-5). Subsequent amplicon sequencing allowed us to determine the corresponding number of repeats for each locus (Table 1). Subsequently, when calculating the polymorphism indicators of VNTR loci, the MLVA profiles of strains from GenBank NCBI and the collection of the Volgograd Plague Control Research Institute were combined.

When calculating the effective number of alleles (ne) and the polymorphic information content (PIC) index for the BPSS1974^#I locus, the values were found to be 4.299 and 0.767, respectively. In the adjusted sample of strains, the number of alleles per locus in the MLVA-6 scheme was distributed as follows: the minimum number of alleles was 8 and 10 in loci L3145k (ne = 4.257; PIC = 0.765) and L933k (ne = 7.087; PIC = 0.859), respectively; 11 alleles in loci L3652k (ne = 7.368; PIC = 0.864) and L1217k (ne = 3.842; PIC = 0.74); 13 alleles in locus S2862k (ne = 7.628; PIC = 0.869); and the maximum number of alleles was 16 in locus L20k (ne = 8.103; PIC = 0.877). The HGDI did not change when locus BPSS1974^#I was added to the MLVA-6 scheme and remained at 0.981.

For the analysis of genetic polymorphism in strains whose VNTR profiles were reliably determined, the method of amplifying differentiating regions of the genome was additionally used. As a result, 51 strains of B. mallei were distributed into 19 DFR types, of which 6 were identified for the first time (DFR20 — 001110011, DFR21 — 000111110, DFR22 — 111110110, DFR23 — 110110111, DFR24 — 010111010, DFR25 — 111101100).

When comparing dendrograms constructed based on typing results using the MLVA-6 scheme (Figure, a) and when including the BPSS1974^#I locus in this scheme (Figure, b), differences were observed in the arrangement of root clusters and individual strains. Thus, when using the MLVA-6 scheme and the BPSS1974^#I locus in combination, a new cluster was formed, which we designated as G and which consists of subclusters G1, G2 and G3. Subcluster G1 formed the B. mallei strain Bahrain1. Subcluster G2 included strains BM-4 and B. mallei 11, while subcluster G3 included strains B. mallei NCTC 10229, 2002734299, and Ivan, isolated by researchers at the Pecs Institute in Hungary in 1961, and B. mallei BM_campo 2.1. For the G3 subclade strains, a common DFR type (DFR02) was established, and their MLVA profiles differed by 1 repeat at loci L3145k and L3652k, while the VNTR genotypes of strains BM_campo 2.1 and Ivan were completely identical.

Comparative analysis of the typing results of 51 B. mallei strains using the MLVA-6 scheme (a) and the combined use of the MLVA-6 scheme and the BPSS1974#I locus (b).

Identical root clusters are marked in the same color. a: Strains whose positions changed upon inclusion of the BPSS1974^#I locus are marked with a dashed line. The strains from the collection of the Volgograd Plague Control Research Institute are highlighted by a frame

The adjacent branches of the dendrogram contain strains from clusters F and E. For strains from cluster F — B. mallei 34 of unknown origin and NCTC 10247 — the same number of repeats were determined at VNTR loci L3145k, L20k, and BPSS1974^#I, and their DFR profiles were unique and differed at loci BmVAT1 and BmVAT6. Cluster E included strains from India (2000031066 and SAVP1) with identical MLVA profiles at loci L933k, L3145k and L20k, while their DFR types were also unique and differed by the absence of locus BmVAT6 in B. mallei SAVP1. A separate branch within cluster E was formed by the B. mallei P-1 strain from Yugoslavia.

The inclusion of locus BPSS1974^#I in the typing scheme allowed strain B. mallei BM-2 to be placed as a separate branch within cluster D, which consists of strains B. mallei 2002721277, China5, and 2000031063. All strains in this cluster have the same number of repeats at loci L3145k, L1217k and BPSS1974^#I, as well as an identical DFR profile (DFR16).

The genotype at locus BPSS1974^#I, corresponding to 58 repeats, was identified only in the Indian strain B. mallei India86-567-2 and PRL-20, which allowed these strains to be grouped into a separate cluster I. Strain India86-567-2 belonged to type DFR01, while PRL-20 belonged to the newly identified type DFR25, which differed from DFR01 by the absence of the BmVAT5 and BmVAT8 loci.

The position of the strain B. mallei 16050, which we isolated and identified, is particularly interesting. It formed a common cluster with B. mallei mongolia_1 and B. mallei V-120, which are genetically identical in their MLVA profile, and its duplicate SCPM-O-B-7146. All 4 strains were found to have identical DFR profiles (DFR06) and VNTR profiles at loci L20k and L1217k. At the same time, B. mallei mongolia_1 and V-120 were also found to have the same number of repeats at loci S2862k and BPSS1974^#I, and they differed by 1 repeat at loci L933k and L3145k, which allowed these strains to be placed in a separate branch within cluster C. For B. mallei 16050, a closer relationship was established with the B. mallei V-120 strain, as their MLVA profiles contained 6 repeats at the L3145k locus.

A separate cluster group was formed by the UFAL2 strain (Brazil, 2017) with a unique number of repeats at loci L20k (n = 14), S2862k (n = 4), and BPSS1974^#I (n = 27). At the same time, B. mallei UFAL2, along with strains from Russia and Mongolia, belonged to type DFR06.

MLVA profiles of duplicate strains from different collections were identical or differed by only 1 locus, allowing them to remain within the same cluster. Thus, differences of 3 repeats were found at locus L933k between B. mallei C-4 and SCPM-O-B-4682, and of 1 repeat between B. mallei C-5 and SCPM-O-B-4683. The Muksuwar-11 strain differed from the SCPM-O-B-7093 and Mukteswar (GCA_033870375.1) strains by 1 repeat at the S2862 locus. The Zagreb strain from our institute collection differed from Zagreb (GCA_033870355.1) by 5 repeats at the S2862 locus (Table 2).

Table 2. MLVA profiles of B. mallei duplicate strains from various collections

B. mallei strain number (GenBank accession number)	Number of repeats in the loci
B. mallei strain number (GenBank accession number)	L933k	L3145k	L3652k	L20k	L1217k	S2862k	BPSS1974^#I
P-1*	6	5	11	7	–	6	55
SCPM-O-B-4688 (GCA_003627635.1)**	6	5	11	7	–	6	55
Muksuwar-11*	5	4	8	17	6,9	12	53
SCPM-O-B-7093 (GCA_003627585.1)**	5	4	8	17	6,9	11	53
Mukteswar (GCA_033870375.1)***	5	4	8	17	6,9	11	53
В-120*	10	6	8	17	5	6	55
SCPM-O-B-7146 (GCA_003627695.1)**	10	6	8	17	5	6	55
Ц-4*	11	7	6	17	–	8	54
SCPM-O-B-4682 (GCA_003627705.1)**	14	7	6	17	–	8	54
Ц-5*	11	7	6	17	–	8	54
SCPM-O-B-4683 (GCA_003627655.1)**	10	7	6	17	–	8	54
Bogor-37*	5	4	8	15	6,9	11	53
Bogor (GCA_033870395.1)***	5	4	8	15	6,9	11	53
Zagreb*	5	4	8	15	6,9	6	53
Zagreb(GCA_033870355.1)***	5	4	8	15	6,9	11	53

Note. * B. mallei strains from the collection of the Volgograd Plague Control Research Institute. **B. mallei strains deposited in the State Collection of Pathogenic Microorganisms and Cell Cultures.

***B. mallei strains from the Friedrich Loeffler Institute, the German Federal Research Institute for Animal Health

As a result, the addition of locus BPSS1974^#I to the MLVA-6 scheme allowed to identify a closer genetic relationship among a number of strains and the formation of new clusters that correlate with their DFR profiles.

Discussion

The effectiveness of intraspecies typing of infectious agents is largely determined by the availability of their genome nucleotide sequences. Whole-genome sequencing and subsequent assembly of B. mallei genome are complicated by a large number of repetitive sequences and a high GC content. As a result, of the 113 B. mallei genomes deposited in NCBI GenBank as of May 2025, only 33 were complete, while most of the remaining genomes are represented by dozens or hundreds of contigs. At the same time, certain strains were not unique.

For certain strains in the sample, it was impossible to determine the VNTR genotype at individual loci in silico. Low read coverage of the genome regions corresponding to the L20k and L3652k loci prevented the target VNTR region from being localized within a single contig in 5 strains. Insufficient coverage during sequencing can be attributed to the GC content of these regions, which was over 70% in the reference strain B. mallei ATCC 23344.

Another problem when searching for target VNTR loci in incomplete genomes is the absence of corresponding genome regions for these loci (null alleles). In this case, there is no reliable way to determine the cause of the null alleles: it could be the result of both natural evolutionary changes and errors in genome assembly. We detected null alleles at the L1217k locus in 14 strains from the sample, 10 of which were represented by incomplete genomes, leading to their exclusion from further study. When performing PCR with DNA from collection strains, the amplification product for the L1217k locus was absent in 3 strains. However, null alleles in complete genomes increase the risk of false homoplasy, which is why such loci are excluded from analysis in many population studies. However, the probability of registering a null allele at each locus increases with the sample size, which can lead to the exclusion of a large number of loci, and this in turn will negatively impact the quality of the research.

As a result, 13 strains were excluded from the initial sample, and the remaining 51 B. mallei strains were distributed into 39 MLVA types using a 6-locus VNTR typing scheme (HGDI = 0.981). We included the BPSS1974^#I locus in this scheme because of the structure of the motif within it, the degeneracy of which reduces the likelihood of DNA polymerase errors during replication and sequencing. At the same time, the sequence of the motif itself is under the influence of stabilizing selection, which contributes to its preservation in the genome.

In the current study, we found that all VNTR loci had a PIC > 0.5, indicating their high discriminatory power, with the effective number of alleles ranging from 3.842 to 8.103 per locus and averaging 6.381. Based on these indicators, locus BPSS1974^#I occupied an intermediate position between VNTR loci L3145k and L933k, represented by perfect and degenerate minisatellite repeats, respectively.

Successful amplification of the BPSS1974#I locus during PCR with DNA from B. mallei strain 16050, isolated from a sick horse in 2023, indicated the stability of the areas flanking the VNTR region. In turn, the full-length nucleotide sequence of the BPSS1974^#I locus, identified by whole-genome sequencing data analysis of the B. mallei 16050 strain, showed its potential for in silico MLVA typing.

Our cluster analysis results of B. mallei strains based on the determination of repeat numbers in 7 VNTR loci were compared with the core genome multilocus sequence typing (cgMLST) data conducted by S. Appelt et al. [13]. The comparison showed that the composition and mutual arrangement of the clusters had many common features. Thus, according to cgMLST data, strains from China were grouped into two main clusters, corresponding to clusters D and L in the MLVA dendrogram we constructed. These clusters included strains for which the geographical region of origin in GenBank is listed as the USA: cluster L included strains B. mallei 2002721274, JHU and FMH, while cluster D included B. mallei 2002721277. Our clustering was consistent with the cgMLST results and indicated a likely Chinese origin of these strains. Within cluster D, a close relationship was confirmed between strains 2000031063 (Hungary) and China5 (China), whose MLVA profiles differed only at the S2862k locus. In addition to the cgMLST data, the genetic proximity of the strains within MLVA clusters D and L was also confirmed by the results of our DFR typing. Thus, with the exception of B. mallei Kweiyang#4, which lacked an amplification product at the BmVAT6 locus, a common DFR profile was determined for all strains within each cluster.

The strains from Turkey investigated in this study were classified in 3 different clusters. The MLVA-6 scheme, supplemented with locus BPSS1974^#I, allowed for the distribution of strains from these clusters according to their position on a minimum spanning tree constructed based on cgMLST data [13]. Thus, based on typing results using only the MLVA-6 scheme, strain B. mallei PRL-20 was included in the same cluster as strains from Turkey. Determining the number of repeats at the BPSS1974^#I locus allowed B. mallei PRL-20 to be placed on the dendrogram in a separate branch with the B. mallei India86-567-2 strain from India, while the strains from Turkey formed a separate, most numerous cluster J. Two other strains from Turkey, B. mallei 11 and NCTC 10247, were included in sub-cluster G2 and cluster F, respectively.

The inclusion of locus BPSS1974^#I in the MLVA-6 scheme, along with strain B. mallei India86-567-2, allowed us to clarify the position of another Indian strain, B. mallei 3708, which formed a separate branch within cluster C. The clustering of Indian strains established using the developed MLVA-7 scheme was consistent with the data obtained in the study by H. Singha et al. [24]. In the presented study, based on MLVA typing at 23 loci, the strains B. mallei India86 567-2, Mukteswar, SAVP1, BMQ and NCTC 3708 were found to belong to different clusters, which was consistent with the results obtained in this study.

The MLVA and DFR profiles of the Brazilian strains we established confirmed the previously proposed hypothesis of two episodes of introduction of the glanders into Brazil [25]. Therefore, the close genetic relationship we identified within sub-cluster G3 between strains of European origin and BM_campo 2.1 suggested a possible introduction of glanders into Brazil from Europe, likely during colonization or trade.

The inclusion of the BPSS1974^#I locus in the MLVA-6 scheme allowed to localize of B. mallei strain 11 in the neighboring G2 sub-cluster, which was consistent with the data from M.V.D. Falcão et al. [25]. The authors performed typing of strains based on 15 phylogenetically informative single nucleotide polymorphisms using PCR followed by high-resolution melting analysis, which identified the L3B3sB3 branch, including strains from MLVA clusters G2 and G3 formed in this study. Analysis of single nucleotide polymorphisms [25] also showed that strain B. mallei UFAL2 belongs to the L3B2 branch, which includes strains isolated in Brazil. Due to the lack of nucleotide sequences from other Brazilian strains in Genbank, we were unable to determine the MLVA and DFR profiles of the L3B2 lineage strains. At the same time, the unique VNTR profile of strain B. mallei UFAL2 allowed it to be placed in a separate MLVA cluster H, which corresponds to the division of Brazilian strains into different phylogenetic lineages.

Based on the results of MLVA and DFR analysis, a close genetic relationship was established between B. mallei strains V-120 and 16050, isolated 18 years apart from diseased animals in Eastern Siberia, and B. mallei strain mongolia_1, isolated from a diseased horse in Mongolia in 2022. Belonging to the same DFR type (DFR06) and sharing a unique VNTR profile at loci L20k and L1217k, along with the geographical proximity of the territories of these countries where glanders outbreaks were reported, suggests a common source of origin of these strains. At the same time, B. mallei strains C-4 and C-5, isolated in Mongolia in 1967, formed a separate MLVA cluster. The DFR type (DFR06) was also identified in the Brazilian strain B. mallei UFAL2. Common DFR profiles of strains from different MLVA clusters have been identified in several other cases, which may indicate a decrease in the specificity of the selected DFR loci as the number of B. mallei genome sequences in genetic databases increases.

The observed changes in the VNTR profiles of duplicate B. mallei strains from different collections may be related to storage conditions, particularly the methods used for maintaining bacterial cultures. In a study by J.M. U'Ren et al. [26], when studying the stability of 32 VNTR loci proposed for typing melioidosis strains using the serial passage method on plates, changes were found in 12 loci. Mutations were identified in the L933k, L3145k and S2862k loci, which were included in our developed typing scheme for glanders causative agent. The variability of loci L933k and S2862k in the genomes of duplicate strains of B. mallei, as established in this study, indicated a similar pattern of variation in these VNTR loci among pathogenic Burkholderia. The stability of the molecular structure of locus BPSS1974^#I has been determined in strains from different collections, which is important for epidemiological investigations and evolutionary analysis.

Thus, the VNTR locus BPSS1974^#I can be considered a marker, the inclusion of which in the MLVA-6 scheme will improve the accuracy of genotyping and the determination of the regions of origin of newly isolated B. mallei strains. Alongside the improvement of schemes for differentiating B. mallei strains for more effective intraspecies typing, there is a necessity to increase the number of available nucleotide sequences of B. mallei genomes, combined with the addition of metadata about their exact geographical origin.

¹ Russia — Burkholderia mallei (Inf. with) (Glanders) — Follow up report 1 [FINAL].

URL: https://wahis.woah.org/#/in-review/4915

² URL: http://www.ncbi.nlm.nih.gov/genome