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, a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 inhalation, direct contact, or ingestion of contaminated animal tissue (often rabbits and rodents), or contact with contaminated soil or water [1][2][3][4][5][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. 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.

Strains
Strains of F. tularensis used in this study are listed in S1 Table. 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.

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 MgCl 2 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.

Single-primer genotyping
Single-primer genotyping was performed by PCR with primer Chif1 (5 0 -CTAGGGCTGGTG GG-3 0 ). 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 ( Table 1).

In silico MicrobesGenotyping database
The MLVA data described in this report is included in S1 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]. 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 CO 2 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×10 1 to 6.25×10 2 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.

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]  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 (S1 Table). 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 Table 2. 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. Fig 1 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.
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.  Fig 3). A second subcluster subsequently called M.II comprises the 18 strains  Fig 3 and Fig  4). A single strain isolated in Karakalpakstan (region C Fig 3) forms a separate subcluster M. III. The position of the three subclusters is also shown in Fig 1. Table 3 provides the list of 25 mediasiatica strains investigated here, together with available metadata. Table 4   The behavior of the atypical strain, 60B-57, is reminiscent of the Central Asian strains. Interestingly strain 60B-57appears to define a third group (Fig 2, 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 (Fig 3). 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 MLVA9 Obolensk , 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 Investigation of first Russian isolates of F. tularensis subsp. mediasiatica order to distinguish this set from other selections with an identical number of loci, we call this selection MLVA9 Obolensk . An even simpler MLVA6 assay was proposed by Gürcan et al. [30].  Table 5 shows the composition of the three MLVA assays. The three assays can be compared using the Johansson et al. investigation [25], which constitutes the most complete dataset. Whereas 119 genotypes are resolved among the 192 strains when using all 25 loci, MLVA6, MLVA9 Obolensk and MLVA10 resolve 111, 101 and 116 genotypes, respectively.
Determination of subspecies by single primer amplification. We previously developed a F. tularensis typing assay using a single-primer PCR (Fig 5). 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. Fig 5 shows the results obtained for nine representative strains, including two holarctica, three tularensis, three mediasiatica, and one novicida.
All mediasiatica strains show an identical pattern. Fig 5 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 (S2 Table). 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].  Table). 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][35][36]. The distance from the MRCA to the tips is relatively homogeneous within each subspecies. Within holarctica, this distance goes from   Table 3 and

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) (S3 Table in supplementary materials). S3 Table shows that all mice died between the fifth and eighth day post infection. In Fig 7, the survival curves of mice infected with the minimum dose are presented. The average time to death of mice infected with the minimum dose (about 5-10 CFU per mice, see S3 Table) 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 (Fig 7). 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 S1 Table). 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 S4 Table. The average weight loss (%) through time is shown in Fig 8. 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 (Fig 7). On the fifth day

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. Investigation of first Russian isolates of F. tularensis subsp. mediasiatica 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 (http://7law.info/altajsky/act2i/u298.htm), 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][43][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.