Russian isolates enlarge the known geographic diversity of Francisella tularensis subsp. mediasiatica.

Francisella tularensis, a small Gram-negative bacterium, is capable of infecting a wide range of animals, including humans, and causes a plague-like disease called tularemia-a highly contagious disease with a high mortality rate. Because of these characteristics, F. tularensis is considered a potential agent of biological terrorism. Currently, F. tularensis is divided into four subspecies, which differ in their virulence and geographic distribution. Two of them, subsp. tularensis (primarily found in North America) and subsp. holarctica (widespread across the Northern Hemisphere), are responsible for tularemia in humans. Subsp. novicida is almost avirulent in humans. The fourth subspecies, subsp. mediasiatica, is the least studied because of its limited distribution and impact in human health. It is found only in sparsely populated regions of Central Asia. In this report, we describe the first focus of naturally circulating F. tularensis subsp. mediasiatica in Russia. We isolated and characterized 18 strains of this subspecies in the Altai region. All strains were highly virulent in mice. The virulence of subsp. mediasiatica in a vaccinated mouse model is intermediate between that of subsp. tularensis and subsp. holarctica. Based on a multiple-locus variable number tandem repeat analysis (MLVA), we show that the Altaic population of F. tularensis subsp. mediasiatica is genetically distinct from the classical Central Asian population, and probably is endemic to Southern Siberia. We propose to subdivide the mediasiatica subspecies into three phylogeographic groups, M.I, M.II and M.III.

Introduction

F. tularensis is a small Gram-negative aerobic coccobacillus. It is a facultative intracellular parasite capable of infecting a wide range of animals and causing a plague-like disease called tularemia. People, particularly those living within endemic foci and working as farmers, hunters, foresters, or meat processing and leather industry workers, are at risk of tularemia infection. This bacterium can be acquired in a variety of ways, including through bites from arthropods, inhalation, direct contact, or ingestion of contaminated animal tissue (often rabbits and rodents), or contact with contaminated soil or water [1-6].

Currently F. tularensis is divided into four subspecies: tularensis (nearctica), holarctica (palaearctica), mediasiatica, and novicida, which differ in their distribution and virulence in humans [2, 7, 8]. Tularensis and holarctica are the most common and well-described subspecies. Subspecies tularensis is mainly found in North America, and shows the greatest virulence in humans [8, 9]. The mortality rate in humans can reach 24% in untreated patients [10]. Subspecies holarctica is commonly found in America and Eurasia, i.e., throughout the Northern Hemisphere. It is less virulent for humans than subsp. tularensis, with a case fatality rate of up to 7% in untreated patients [9].

F. tularensis subsp. novicida is present in North America, but some strains have been isolated in tropical Australia and Southeast Asia (clinical isolates only) [11, 12]. To date, the status of novicida, separate species or subspecies, remains non-standardized and controversial [4, 13–16]. Although this subspecies is almost avirulent in humans, there have been several reports of this subspecies isolated from patients [11, 12, 17–19].

F. tularensis subsp. mediasiatica remains the least studied and understood subspecies. It was found in Central Asia [20, 21], in some regions of Kazakhstan and Turkmenistan. This geographic location is responsible for the subspecies name—in Russian “Middle Asia” is a synonym for “Central Asia”. Nothing is known about the potential for this subspecies to cause disease in humans [22, 23]. The absence of clinical cases is indirect evidence suggesting that virulence in humans is low at best. Virulence in hares is similar to virulence of subsp. holarctica (reviewed in [2]).

In 2013, F. tularensis subsp. mediasiatica was reported to be present in Russia [24]. Three of four Altai-originating strains (designated in this article as A554, A678, and A823) were shown to belong to subsp. mediasiatica. For description of these isolates we used Multiple-locus VNTR (variable number of tandem repeat) analysis (MLVA). This approach has been previously shown to be very efficient as a first line assay able) to distinguish among even very close genetic relatives, including individual strains [2]. MLVA profiles of these strains were similar, and they differed slightly from MLVA-profiles of previously reported Central Asian strains.

In 2013–2014, we received 15 additional strains of F. tularensis from an Altai anti-plague station and Federal Healthcare Service Center for Hygiene and Epidemiology in the Altai region. In the present report, we compare the MLVA and biochemical data from these strains with data from the rest of our collection and from the literature. We show that they all belong to subsp. mediasiatica, but confirm that they constitute a distinct cluster, providing new insight on the intra-subspecies structure. In addition, we investigated the virulence of Altai strains in laboratory mice.

Materials and methods

Strains

Strains of F. tularensis used in this study are listed in . Eighteen strains were isolated within the Republic of Altai and Altai Territory in 2011–2014 by the Federal Healthcare Service Center for Hygiene and Epidemiology in the Altai region and by the Altai anti-plague station. All Altai strains were isolated from ticks and dead rodents, none were clinical isolates.

Bacterial culture

Strains were grown at 37°C on solid (FT-agar) and liquid (FT-broth) nutrient media (SRCAMB, Obolensk, Russia). The composition of the FT-agar was 3.8% erythritol-agar, 1% dried bovine blood, 1% glucose, 0.05% cysteine, and 0.0025% thiamine chloride at pH 7.2. The composition of the FT-broth was 2% casein enzymatic hydrolysate, 1% yeast extract, 1.2% KH2PO4, 1% glucose, 0.001% cysteine, and 0.001% FeCl2 at pH 7.2.

Polymerase chain reaction (PCR)

PCR was performed using a Mastercycler Gradient thermal cycler (Eppendorf, Hamburg, Germany). Each reaction mix contained 5 μl of 5×qPCRmix-HS (Evrogen, Moscow, Russia). This mix corresponds to a final concentration of 200 μM for each dNTP and 3 mM MgCl2 with one unit of Taq DNA polymerase and 10 to 50 ng of DNA per 25 μl reaction. Primers (Syntol, Moscow, Russia) were added at a final concentration of 10pM.

Multiple-locus VNTR analysis (MLVA) genotyping

We used 15 among the 25 VNTR loci described in [25]: FtM3, FtM5, FtM6, FtM7, FtM8, FtM9, FtM10, FtM11, FtM12, FtM13, FtM16, FtM19, FtM20, FtM23, FtM24. In vitro MLVA genotyping was done as previously described [25] except for the use of regular, non-fluorescent primers, agarose gels and monoplex PCR. PCR products size was evaluated using agarose gel-electrophoresis. The PCR products and a 20 bp ladder (Bio-Rad, USA) were electrophoresed at 100 V for 240 min on a 32-cm length 3% agarose gel prepared in 0.5× TBE. The DNA fragments were visualized with ethidium bromide staining and ultraviolet (254 nm) using the Doc-Print gel documenting system and PhotoCaptMw software version 99.04 (Vilber Lourmat, Marne-la-Vallée, France). PCR products larger than 600 bp were reanalyzed on 2% agarose gel for better resolution. Also in these few cases we confirmed the size of amplicon using Experion™ Automated Electrophoresis System (BioRad, Hercules, USA) and by sequencing the fragment, followed by a direct count of the number of repeats.

Cluster analyses were run using BioNumerics version 7.6.2 (Applied-Maths, Belgium).

Single-primer genotyping

Single-primer genotyping was performed by PCR with primer Chif1 (5′-CTAGGGCTGGTGGG-3′). Thermocycling parameters were the following: 94°C for 3for3 min; 30 cycles of 94°C for 30 sec, 54°C for 30 sec, and 72°C for 1 min; and then a final elongation at 72°C for 5 min. The PCR products and a 100 bp ladder (Bio-Rad, USA) were electrophoresed at 100 V for 90 min on a 1% agarose gel prepared in 1× TBE and visualized with ethidium bromide staining and ultraviolet (254 nm) using the Doc-Print gel documenting system (VilberLourmat, France). The particular subspecies was determined by the size of the obtained amplicons ().

In silico MicrobesGenotyping database

The MLVA data described in this report is included in S1 Table and has been imported into the public Francisella database at http://microbesgenotyping.i2bc.paris-saclay.fr/.

Draft whole genome sequencing and whole genome single nucleotide polymorphism (SNP) analysis

Complete genomes and genome assemblies were downloaded via Genbank/NCBI. Sequence Reads Archives (SRA) files were downloaded from the European Nucleotide Archive (ENA). SRA files quality was checked using FastQC. Reads with poor quality towards the ends were trimmed using Fastx_trimmer from the FASTX-Toolkit. Complete genomes and assemblies were converted into artificial 300 bp long reads using wgsim without introducing errors or mutations. Read files were then imported into BioNumerics version 7.6.2 to be mapped on the reference genome. SNPs were called using the BioNumerics SNP analysis pipeline with the “strict SNP clustering (closed SNP dataset)” option.

Biochemical properties determination

Fermentation of glycerol

Acid production during glycerol fermentation was estimated in an FT-broth without glucose or casein hydrolysate, but with the addition of 0.2% glycerol (w/v) and 0.02%(w/v) phenol red as an indicator.

Detection of citrulline ureidase activity

We determined the activity of F. tularensis using a color test with ninhydrin [26] based on ability of citrulline ureidase to degrade citrulline to ornithine, which, on reaction with ninhydrin reagent gives a pink coloration [27].

A 1.0-ml sample of bacterial suspension (1010 bacteria) in 0.1 M phosphate buffered saline (PBS, pH 6.5) was mixed with 1.0 ml of 0.7% (w/v) l-citrulline (Sigma Chemical Co., St. Louis, MO, USA) and incubated for 20 h at 30°C. An aliquot (0.01 ml) of the mixture was removed and added to 0.49 ml of distilled water, 1.0 ml of freshly prepared ninhydrin reagent (625 mg of ninhydrin [Sigma] in 10 ml of 6 M H3PO4 and 15 ml of glacial acetic acid), and 1.5 ml of acetic acid. Samples were boiled for 1 h, and pink colored samples were considered to have citrulline ureidase activity.

Beta-lactamase activity

One microbiological wire loop of each tested strain was suspended into 1 ml of phosphate-buffered saline (PBS). One hundred μl of this suspension were mixed with 10 μl of nitrocefin solution (500 mg/ml).

After incubation at room temperature during one hour, red colored samples were considered to have beta-lactamase activity [28].

Animal experiments

Ethics statement

All protocols for animal experiments were approved by the State Research Center for Applied Microbiology and Biotechnology Bioethics Committee (Permit No: VP-2016/2). They were performed in compliance with the NIH Animal Welfare Insurance #A5476-01 issued on 02/07/2007 and the European Union guidelines and regulations on handling, care, and protection of laboratory animals (http://ec.europa.eu/environment/chemicals/lab_animals/home_en.htm).

Mice

Five-to-eight-weeks-old BALB/C mice (not-specific pathogen free) of both genders, weighing 18–20g (purchased from Laboratory Animals Breeding Center, Shemyakin and Ovchinnikov Institute of Bioorganic Chemistry, Russia), were used in experiments. The mice were housed in polycarbonate cages with space for comfortable movement and easy access to food and water, under constant temperature and humidity conditions (22°C ± 2°C and 50% ± 10%, respectively) and a 12-hour light/12-hour dark cycle. Mice were fed Mouse Mixed Fodder PK-120 (Laboratorkorm, Russia) and provided tap water ad libitum throughout the study.

A minimum number of mice was used for experiments. The mice were randomly divided into experimental groups.

Approved protocols provided scientifically validated humane endpoints, including pre-set criteria for euthanasia of moribund mice by CO2 inhalation. In our animal studies, mice were euthanized when they became lethargic, dehydrated, moribund, unable to rise, or non-responsive to touch.

The health condition of the animals was monitored at least twice a day.

Virulence determination

Mice were inoculated subcutaneously with 0.1 ml F. tularensis cells in PBS (5×101 to 6.25×102 colony forming units (CFU) per animal) in the inner part of the upper thigh.

To determine the actual infectious dose, serial ten-fold dilutions of bacterial cells were plated onto FT-agar, and the colonies were counted after incubation at 37°C for 48 hours. Dead mice were autopsied and subjected to bacteriological studies. The surviving animals were monitored for 14 days.

Vaccination and challenge

BALB/C mice were inoculated subcutaneously with live vaccine strain 15 NIIEG (20 CFU per mouse in 0.1 ml of PBS) in the inner part of the upper thigh. Three weeks later, the mice were challenged with virulent strains of F. tularensis (1000 CFU per mouse). Survival and body weights of immunized mice were monitored every day for 14 days. Mice were weighed in groups of six animals (or less in the event of a death) in a special tray, and the average weight was calculated.

Percent body weight loss was calculated as follows: 100 − final weight/initial weight ×100.

Statistics

Survival curves were generated by plotting the proportion of surviving mice against the time of death (days). For comparison of mouse weights, data was compiled from four independent experiments. The mean ± standard deviation (SD) was calculated for the percentage weight loss and mean time to death. A significant difference in weight loss or time to death between groups of animals was determined using a Student’s t-test with significance set at p<0.05.

Map-containing figures

Map-containing figures were designed using Google-maps service (map data 2017). Maps are the property of Google and partners.

Results

MLVA genotyping of a Russian collection of F. tularensis

Twenty-five polymorphic tandem repeats have been described by Johansson et al. [25]. Different subsets have been subsequently selected to suit specific purposes. We retained the 15 loci with repeat size of 9 bp or more which can be typed confidently, not only on relatively expensive capillary electrophoresis equipment, but also on regular agarose gels (FtM3, FtM5, FtM6, FtM7, FtM8, FtM9, FtM10, FtM11, FtM12, FtM13, FtM16, FtM19, FtM20, FtM23, FtM24). The example of locus size estimating is in S1 Appendix. Low cost and robustness is one major advantage of MLVA as first-line assay and quality check of strain identity prior to running e.g. whole genome sequencing. This is essential for rapid epidemiological investigation and local investigations by regional anti-plague and medical institutions. We first evaluated the validity of this assay by using the MLVA data associated with [25] and accessible via http://microbesgenotyping.i2bc.paris-saclay.fr/. Using all 25 loci, Johansson et al. resolved 119 genotypes among 192 strains. Subspecies were clustered as expected, and additional subdivisions were proposed, tularensis A.I and A.II, holarctica B.I to B.V. In comparison MLVA15 resolves 104 genotypes. The previously defined groups are conserved, except for B.I, B.II and B.III which are not resolved when using this subset (analysis not shown).

Since 2008 SRCAMB acts as the Tularemia Reference Center and has assembled a collection of 173 F. tularensis live strains or DNA samples collected by Russian anti-plague and medical institutions (). The strains were genotyped using MLVA15. One hundred and twelve genotypes are resolved. The Simpson’s diversity index was 0.9894 (95% confidence interval 0.9850–0.9938). The number of detected alleles and Simpson’s diversity index calculated at all used VNTR markers are indicated in In agreement with previous reports, VNTR locus FtM3 is the most variable locus.

The data was merged with MLVA15 data from Johansson et al., to produce a minimum spanning tree. shows the resulting clustering, the color code reflects subspecies and additional subdivisions. The main cluster corresponds to holarctica. Subspecies tularensis subgroups A.I and A.II are clearly identified. All strains of Russian geographic origin are holarctica, with the exception of mediasiatica strains. The mediasiatica group is the most distantly related, with a long branch (six differences) connecting it to the novicida subspecies. Within mediasiatica, two additional long branches are detected.

Minimum spanning tree deduced from MLVA15 data for 365 entries (192 from [25] and 173 from the present investigation).

Branches with length up to four are drawn. Short unlabeled branches correspond to a distance of one locus difference. Branch length of five and more are not drawn.

Although the origin of the Altai strains is not typical for subsp. mediasiatica, clusterisation together with «classical» mediasiatica strains points that these strains do belong to subsp. mediasiatica. Previously, F. tularensis subsp. mediasiatica was believed to be abundant only in Central Asia. Until 2011, there have been no strains isolated in the territory of Russia, other than those belonging to subsp. holarctica. The presence of subsp. mediasiatica strains in Russia can be explained by two alternative hypotheses, recent accidental dissemination beyond its native area, or underestimation of the distribution of this subspecies.

presents a focus on the twenty-nine mediasiatica strains including four from the Johansson et al. investigation. The figure shows available metadata and MLVA genotypes. A first subcluster subsequently called M.I in keeping with previous nomenclature contains the Johansson et al. strains and six additional strains isolated in Central Asia, mostly Kazakhstan (region A ). A second subcluster subsequently called M.II comprises the 18 strains isolated within the Republic of Altai and Altai Territory in 2011–2014 (region B and ). A single strain isolated in Karakalpakstan (region C ) forms a separate subcluster M.III. The position of the three subclusters is also shown in . provides the list of 25 mediasiatica strains investigated here, together with available metadata. illustrates some of the major differences in terms of MLVA genotypes between M.I, M.II and M.III. M.II strains of Altaic origin were characterized by relatively large alleles at loci Ft-M3, Ft-M6, Ft-M7, and Ft-M20.

UPGMA clustering analysis of mediasiatica strains based on MLVA15 data.

The dendrogram was constructed using UPGMA (Unweighted Pair Group Method with Arithmetic Mean). The MLVA data of 12 loci are shown (loci Ft-M16, Ft-M23 and Ft-M24 with an identical allele 1 in all strains are masked). The Table 3# index refers to the location indicated in and numbering in . The cluster cut-off value of 60% defines three groups, labelled M.I (Central Asian origin, red), M.II (Altaic origin, green) and M.III (Karakalpak, blue). Strains from SRCAMB collection are indicated using prefix SRCAMB#.

Table 3

The list of 25 mediasiatica strains investigated in this report.

№	Strain	Region	Geographic Origin	Source	Year
1	A137	Altai	Paspaul village, Choya district, Altai	Tick H. concinna	2014
2	A139	Altai	Choya district, river Achpanak, Altai	Tick H. concinna	2014
3	A187	Altai	Choya district, river Achpanak, Altai	Tick H. concinna	2013
4	A188	Altai	Choya district, river Achpanak, Altai	Tick H. concinna	2013
5	A176	Altai	Karakokcha village, Choya district, Altai	Tick D. Silvarum	2014
6	A174	Altai	Choya district river Pervaya Uba, Altai	Tick H. concinna	2014
7	A82	Altai	Choya district, river Pervaya Uba, Altai	Tick D. silvarum	2014
8	A142 K	Altai	Choya district river Uba, Altai	Tick H. concinna	2014
9	A142 M	Altai	Choya district river Uba, Altai	Tick H. concinna	2014
10	A199	Altai	Choya district, river Uba, Altai	Tick D. Silvarum	2014
11	A116	Altai	Altay district, 25-km, Altai	Tick D. silvarum	2014
12	A84	Altai	Kolovo village, Sovetskoye district, Altai	Tick D. silvarum	2014
13	A678	Altai	Pokrovka village, Pervomayskoye district, Altai	Tick Ix. persulcatus	2011
14	A554	Altai	Martynovo village, Eltsovka district, Altai	Tick H. concinna	2011
15	A373	Altai	Mayma district, river B.Tara, Altai	Tick D. silvarum	2013
16	A99	Altai	Aleksandrovka village, Mayminskiy district, Altai	Tick H. concinni	2013
17	A823	Altai	Molokovo village, Shelabolikha district, Altai	Siberian red vole	2011
18	A343	Altai	Altai, without details	Tick H. concinna	2013
19	117	Central Asia	unknown	unknown	1960
20	120	Central Asia	Bakanas village, Kazakhstan, Central Asia	Hare	1960
21	120 (113)	Central Asia	Bakanas village, Kazakhstan, Central Asia	Hare	1965
22	122	Central Asia	Central Asia, without details	Tick R. pumilio	1973
23	A122 (543)	Central Asia	Almaty region, Kazakhstan, Central Asia	midday gerbil	1965
24	150 (31)	Central Asia	Kumozek, Jambyl Province, Kazakhstan, Central Asia	Hare	1969
25	60B-57	Central Asia, Karakalpak	Chimbay district, Karakalpakstan, Uzbekistan, Central Asia	Tick Hyalomma	1960

Distribution of F. tularensis subsp. mediasiatica.

A, Central Asia region; B, Altai region; C, Karakalpakstan.

Isolation points of Altaic mediasiatica strains.

Strains are designated by numbers: 1—A137; 2—A139; 3—A187; 4—A188; 5—A176; 6—A174; 7—A82; 8—A142 K; 9-A142 M; 10—A199; 11—A116; 12—A84; 13—A678; 14—A554; 15—A373; 16—A99; and 17—A823 as also shown in and .

For instance thirteen of eighteen M.II strains of Altai origin showed more than 30 repeats at locus Ft-M3, whereas M.I Central Asian strains showed less than 27.

The behavior of the atypical strain, 60B-57, is reminiscent of the Central Asian strains. Interestingly strain 60B-57appears to define a third group (, M.III label). This strain is the only strain originating from Karakalpakstan in Western Uzbekistan further illustrating the relation between MLVA genotypes and geographic origin (). This strong anchorage of mediasiatica with geography is suggesting that there is limited spreading between different ecological niches and consequently that the diversity of mediasiatica may be underestimated, given its low impact on human health.

Establishing simple first-line assays for subpopulation assignments

MLVA9Obolensk, as a first line MLVA typing assay

In order to simplify the initial genotyping steps and lower the analysis cost, we have evaluated a minimum panel of VNTRs able to efficiently resolve F. tularensis strains naturally present in Russia and Central Asia. Using the MLVA data from the 144 corresponding strains, the Automated Selection of Typing Target Subsets (AuSeTTS) software [29] provides a list of nine VNTR loci sufficient to achieve the same discrimination as the 15 loci in the tested data set: Ft-M3, Ft-M5, Ft-M6, Ft-M7, Ft-M8, Ft-M12, Ft-M19, Ft-M20, Ft-M24. In keeping with previous similar MLVA assays, and in order to distinguish this set from other selections with an identical number of loci, we call this selection MLVA9Obolensk. An even simpler MLVA6 assay was proposed by Gürcan et al. [30]. Four of the six loci proposed by these authors are included in the MLVA9Obolensk panel (Ft-M3, Ft-M6, Ft-M20, Ft-M24), the last two (Ft-M21 and Ft-M22) with seven and six base-pair repeat units respectively, were not included in the present investigation. Similarly, Vogler et al. [31] proposed an MLVA10 assay, comprising Ft-M2, Ft-M3, Ft-M4, Ft-M5, Ft-M6, Ft-M10, Ft-M20, Ft-M22, Ft-M23, Ft-M24. Table 5 shows the composition of the three MLVA assays. The three assays can be compared using the Johansson et al. investigation [, which constitutes the most complete dataset. Whereas 119 genotypes are resolved among the 192 strains when using all 25 loci, MLVA6, MLVA9Obolensk and MLVA10 resolve 111, 101 and 116 genotypes, respectively.

Table 5

Composition of proposed selections of VNTR loci for simplified MLVA assays and comparison from in silico derived data.

Locus					MLVA assay
	Repeat unit size	Allele number	HGDI1	Standard deviation	MLVA9ObObolensk	MLVA6 [30]	MLVA10 [31]
Ft-M2	6bp	2020	0.45294529	[0.3640,0.5419]			X
Ft-M3	9bp	3030	0.94629462	[0.9358,0.9566]	X	X	X
Ft-M4	5bp	66	0.58375837	[0.5283,0.6391]			X
Ft-M5	16bp	66	0.21352135	[0.1350,0.2920]	X		X
Ft-M6	21bp	66	0.64616461	[0.5985,0.6936]	X	X	X
Ft-M7	16bp	66	0.19481948	[0.1189,0.2707]	X
Ft-M8	16bp	66	0.31633163	[0.2307,0.4019]	X
Ft-M10	16bp	1515	0.44634463	[0.3569,0.5356]			X
Ft-M12	10bp	2	0.03090309	[0.0000,0.0656]	X
Ft-M19	13bp	2	0.40184018	[0.3440,0.4596]	X
Ft-M20	12bp	1515	0.34443444	[0.2558,0.4330]	X	X	X
Ft-M21	7bp	66	0.35323532	[0.2764,0.4300]		X
Ft-M22	6bp	55	0.642642	[0.6115,0.6725]		X	X
Ft-M23	23bp	33	0.17171717	[0.1025,0.2410]			X
Ft-M24	21bp	33	0.49934993	[0.4780,0.5206]	X	X	X

1 the Hunter-Gaston diversity index was calculated using the Johansson et al, [ dataset for 192 strains.

Determination of subspecies by single primer amplification

We previously developed a F. tularensis typing assay using a single-primer PCR (). This method which is essentially a Random Amplified Polymorphic DNA (RAPD) assay [32] is used for the rapid determination of the subspecies of strains arriving in the SCRAMB collection. We observed that this method is able to distinguish all four subspecies by analysing the PCR product using agarose gel electrophoresis. shows the results obtained for nine representative strains, including two holarctica, three tularensis, three mediasiatica, and one novicida.

Single primer—genotyping of representative F. tularensis strains.

The figure shows the distribution of the single-primer PCR amplicons for representatives from all four subspecies. Lines 1, 11– GeneRulerTM 100 bp Plus DNA Ladder; 2– subsp. holarctica 503; 3– subsp. holarctica А1045; 4– subsp. tularensis SCHU S4; 5– subsp. subsp. tularensis BB-399; 6– subsp. tularensis 8859; 7– subsp. mediasiatica 120; 8 –subsp. mediasiatica A678; 9– subsp. mediasiatica A554; 10– subsp. novicida 112.

All mediasiatica strains show an identical pattern. shows the result obtained with M.I strain 120 and M.II strains A678 and A554. Therefore this simple assay appears to be valid also for the new mediasiatica strains.

Phylogeny by whole genome SNP analysis of F. tularensis including one mediasiatica M.II representative

In order to confirm the clustering suggested by MLVA analysis, we sequenced one representative strain from the mediasiatica M.II group, strain A554. Additional public whole genome sequence data including full genomes, contigs, scaffold and reads was downloaded from the European Nucleotide Archive (ENA). Data from 65 strains was retained after removal of datasets of insufficient quality, or representing variants of laboratory strains or duplicates. This final set includes 31 holarctica, 31 tularensis and three mediasiatica strains (). The mediasiatica strains are A554 from M.II (this report) and FSC147 (alias GIEM 543) and FSC148 belonging to M.I, both from Kazakhstan. A total of 9,737 core genome SNPs is retained by the SNPs calling pipeline. The resulting minimum spanning tree has a size of 9,819 indicating a homoplasy of 0.8This observation is in agreement with previous analyses indicating that evolution within Francisella tularensis stricto sensu (not including F. novicida) is clonal [33, 34]. shows the minimum spanning tree of the three clonal subspecies. The Most Recent Common Ancestor (MRCA) of the three subspecies is positioned by using F. novicida and F. hispaniensis to root the tree (). The MRCA is located on the longest branch, which separates holarctica from mediasiatica and tularensis. The topology of the tree is in full agreement with previous reports [6, 31, 34–36]. The distance from the MRCA to the tips is relatively homogeneous within each subspecies. Within holarctica, this distance goes from 2000 in lineage B.V (currently called B.16) up to 2400 in lineage B.IV (currently called B.6). The distances within mediasiatica are approximately 2400 with little variations. The distances within tularensis range from approximately 1130 in A.II up to 1290 in A.I, significantly shorter than in the other two subspecies.

Minimum spanning tree deduced from whole genome SNP analysis.

The 9,737 SNPs constitute a tree with a size of 9,819 SNPs corresponding to an homoplasy level of 0.8%. The red star indicates the position of the MRCA. The tree includes a few laboratory-derived strains: four LVS representatives within holarctica B.I, two ATCC 6223 representatives within tularensis A.II and at least five SCHU S4 variants among A.I. The star-like pattern around the SCHU S4 and short radiating branches may be the result of a very recent outbreak. The correspondence between lineage names defined by MLVA and the more recent names defined by SNPs is indicated, as in B.I = B.12.

Determination of subspecies by biochemical properties

In order to confirm the subspecies assignment, we examined the biochemical properties of the Altaic strains. Subspecies tularensis and mediasiatica may be distinguished from subsp. holarctica by their capacity to ferment glycerol and their possession of the enzyme citrulline ureidase [37, 38]. The main biochemical difference between subspecies tularensis and mediasiatica is beta-lactamase activity: subsp. mediasiatica does not exhibit this activity [39]. We found that all Altaic strains were able to ferment glycerol and showed citrulline ureidase but not beta-lactamase activity. Therefore, we concluded that these strains belonged to subsp. mediasiatica.

Virulence determination

For a determination of the virulence of the Altaic strains, we used a mouse model. BALB/c mice died after subcutaneous injection of Altai-origin F. tularensis (< 10 bacterial cells per mice) ( in supplementary materials). shows that all mice died between the fifth and eighth day post infection. In the survival curves of mice infected with the minimum dose are presented.

Survival curves of naive mice following challenge with F. tularensis subspecies tularensis, holarctica, or mediasiatica.

Step-down survival curves showing the percentage of mice surviving over 10 days post infection are presented. Data is shown for subsp. mediasiatica (strains A678, A823, A554, and 120), subsp. holarctica (strains A1045 and 503), and subsp. tularensis (SCHU S4). The mean time to death is denoted as a bold black cross.

The average time to death of mice infected with the minimum dose (about 5–10 CFU per mice, see ) was 5.8±0.6 days for subsp. mediasiatica, 6.6± 0.8 days for subsp. holarctica, and 5.4±0.5 days for subsp. tularensis (. Student’s t-testing showed no statistically significant differences between the groups (p > 0.05). Therefore, the virulence of the Altaic strains was similar to the virulence of the other subspecies investigated. Because the high pathogenicity of F. tularensis in this animal model did not allow determining differences in the virulence of different subspecies, we decided to use vaccinated mice, which are significantly more resistant to F. tularensis. Mice were immunized subcutaneously with the F. tularensis 15 NIIEG vaccine strain (20 CFU per mouse) commonly used in Russia and some post-Soviet countries (see ). After 21 days, the mice were infected subcutaneously with strains A678 (mediasiatica), SCHU S4 (tularensis), and 503 (holarctica) (1000 CFU per mouse). This experiment was conducted with four independent replicates. Virulence was evaluated by determining the mortality rate and disease severity, which was estimated based on the loss in body weight of infected mice. The weight loss of vaccinated mice after infection with pathogenic strains is shown in . The average weight loss (%) through time is shown in .

Weight loss (%) of immunized BALB/c mice after challenge with highly virulent strains of F. tularensis.

Average weight loss through time in vaccinated mice infected with pathogenic strains of the different subspecies. Red line, weight loss in mice challenged with holarctica strain 503; Blue line, weight loss in mice challenged with tularensis strain SCHU S4; Green line, weight loss in mice challenged with mediasiatica strain A678. Solid line, immune mice; dashed line, non-immune mice (negative controls). All nonimmune mice died by 5 days postinfection. Body weight loss was calculated as the percentage of initial body weight in four experiments.

The mortality rate in mice challenged with mediasiatica strain A678 was 8% (2/24) and in those challenged with tularensis strain SCHU S4 was 12.5% (3/24). All mice challenged with holarctica strain 503 survived. Maximum weight loss was registered between days 5 and 8 after infection. Weight loss caused by challenge with strain A678 was intermediate between the weight loss caused by challenge with strains SCHU S4 and 503 (. On the fifth day postinfection, weight loss was 9.4±1% after infection with strain SCHU S4, 6.5±1.3% with strain A678, and 3.7±1.2% with strain 503 (p<0.05). Thus, we found a statistically significant difference in virulence among the subspecies tested. Virulence of subsp. mediasiatica strain A678 is intermediate between that of subspp. tularensis strain SCHU S4 and holarctica strain 503.

Discussion

The isolation of a number of subsp. mediasiatica strains in the Altai territory suggests the existence of a natural focus of F. tularensis subsp. mediasiatica circulation. This focus is located more than 1500 km eastward from previously identified foci in Kazakhstan and Turkmenistan.

The MLVA data clearly indicates that the Altai population of F. tularensis subsp. mediasiatica is genetically distinct and likely endemic to the Altai region. Thus, currently known strains of subsp. mediasiatica can be divided into three genetic groups, which we propose to call M.I, M.II and M.III. Our data adds to the present knowledge of the intra-specific structure of F. tularensis. To date, intra-subspecific clustering were proposed only for subsp. tularensis and holarctica. Based on MLVA [25], Johansson et al. divided subsp. tularensis into two subgroups (A.I and A.II) and subsp. holarctica into five subgroups (B.I–B.V). We propose dividing subsp. mediasiatica into three groups, M.I, M.II and M.III based on both genotype (MLVA-profile) and geographical location. M.III is defined by a single strain, 60B-57 from Karakalpastan in Western Uzbekistan, no strain from Turkmenistan could be investigated here.

Thus, we have found that Altai has genetically and geographically distinct population of F. tularensis subsp. mediasiatica circulating, and that the distribution of this subspecies is much broader than previously believed. We believe that the lack of earlier reports results mainly from the extremely low human population in Siberia and inaccessibility of some of its regions (mountainous Altai in particular).

These extreme conditions greatly complicate both the identification of pathogens and diagnosis of the disease in this environment. To date, there have been no reported cases of tularemia in humans caused by the subsp. mediasiatica in Altai. Nevertheless, according to the Chief State Sanitary Doctor of the Altai region (), cases of tularemia in humans are identified in this region almost every year. Importantly, when cases of human tularemia are detected, medical and epidemic services usually do not diagnose the subspecies of pathogen—they only identify F. tularensis. Difficulties in transport of the strains to specialized research laboratories make further investigation difficult. Thus, on site research is needed to answer the question of whether tularemia is caused by subsp. mediasiatica or subsp. holarctica. To date, we have conducted focused research on ticks, rodents, soil, and water samples, but not on clinical samples. After the first detection of F. tularensis subsp. mediasiatica in the Altai region in 2013 [24], an intensive search for this microorganism, carried out by the Altai anti-plague station, allowed us to identify 15 strains of this subspecies in 2013–2014 described in this manuscript and six strains (which were not included in this report) in 2015. We hope that the simple genotyping assays proposed here will help develop procedures that can be applied on site at low cost on clinical as well as environmental samples. MLVA is especially interesting as compared to typing of a selection of SNPs because it is an unbiased approach, which will identify a previously unknown clade. The online Francisella tularensis database hosted on the Microbes Genotyping web site provides an easily accessible depository of reference data for comparisons.

Whole genome sequencing will help investigate more precisely the population structure of the M.II group, as illustrated by the investigation of other F. tularensis foci [40] and the evolutionary relationships between subsp. mediasiatica, holarctica and tularensis. Wang et al [41] recently described the presence of all basal holarctica lineages in Western China. The present description of mediasiatica in Altai represents one additional step pointing to similarities between the phylogeography of F. tularensis and Yersinia pestis. The most basal Y. pestis lineages showing no virulence in humans are also found in Altai, at the junction between China, Mongolia, Russia and Kazakhstan [42-44]. The two species share similar ecological behaviors and vectors, it will be interesting to investigate to which extend their phylogeography is coincident.

Conclusion

We showed that the geographic distribution of F. tularensis subsp. mediasiatica is much wider than previously believed, including not only Central Asia but also Southern Siberia and the Altai region of Russia. Because strains of Altaic origin are genetically distinct from strains of Central Asian origin, we propose the division of subsp. mediasiatica into two intra-subspecies phylogeographic groups called M.I (Central Asian origin) and M.II (Altaic origin). In addition, we showed that the virulence of subsp. mediasiatica in a vaccinated mouse model is intermediate between that of subsp. tularensis and subsp. holarctica.

2% agarose- electrophoresis of Ft-M3 PCR-products of 8 F. tularensis strains.

M1-100 bp plus ladder, M2- 20 bp ladder; the 500 bp and 1 kb bands are slightly more intense 1—A99 subsp. mediasiatica (75 repeats, 882 bp), 2- A823 subsp. mediasiatica (63 repeats, 774 bp), 3—A116 subsp. mediasiatica (52 repeats, 675 bp), 4–120 subsp. mediasiatica (19 repeats, 378 bp), 5—A373 subsp. mediasiatica (53 repeats, 684 bp), 6—I346 subsp holarctica (11 repeats, 306 bp), 7—A139 subsp. mediasiatica (26 repeats, 441 bp), 8–319/358 subsp holarctica (24 repeats, 423 bp). Interface of PhotoCaptMw with the calculation of the size of amplicons.

(PDF)

Click here for additional data file.

The strains of F. tularensis used in this study and MLVA15 profiles.

(XLSX)

Click here for additional data file.

(XLSX)

Click here for additional data file.

BALB/c mice death after infection with F. tularensis strains of different subspecies.

(XLSX)

Click here for additional data file.

Body weight (g) of preliminarily vaccinated mice after infection with pathogenic strains A678, SCHUS4, and 503.

(XLSX)

Click here for additional data file.

Table 1

Presence of diagnostically important amplicons in the single-primer typing system for different subspecies of F. tularensis.

approximate size of amplicon, bp	subspecies
	holarctica	tularensis	mediasiatica	novicida
1460*	+	+/-	+/-	+
960	-	+	+	-
560	+	-	-	-
530	-	+	-	-
390	-	+	+	-
290	+	+	+	+

* 1460 bp amplicon is present in all subspecies, but in lower amount in subsp tularensis and mediasiatica strains (designated as+/-)

Table 2

Number of detected alleles and Simpson’s diversity index calculated at 15 VNTR markers using the SRCAMB collection of 173 strains.

Locus	Number of alleles	Simpson’s Index (confidence interval)
FtM3	37	0.95063(0.9407–0.9606)
FtM5	8	0.29074(0.2020–0.3795)
FtM6	7	0.67039(0.6169–0.7239)
FtM7	8	0.28224(0.1934–0.3710)
FtM8	5	0.38855(0.3057–0.4714)
FtM9	4	0.15308(0.0795–0.2266)
FM10	9	0.14384(0.0710–0.2167)
FM11	6	0.27613(0.1915–0.3607)
FM12	3	0.32795(0.2472–0.4087)
FM13	2	0.09920(0.0383–0.1601)
FM16	2	0.05614(0.0084–0.1039)
FtM19	4	0.34576(0.2689–0.4226)
FtM20	14	0.44475(0.3505–0.5390)
FtM23	2	0.02286(0.0000–0.0544)
FtM24	2	0.40695(0.3474–0.4665)

Table 4

Remarkable differences of Tandem repeat copy number between the Central Asian M.I, Altaic M.II and Karakalpak M.III strains.

Locus	Region of origin
	M.I Central Asian = 6+4	M.II Altain = 18	M.III Karakalpakstan (Central Asia)n = 1
Ft-M3	7–27	21–75	23
Ft-M6	2–3	5–6	3
Ft-M7	2	5–9	1
Ft-M20	6	10–14	3