Genetic diversity of Francisella tularensis subsp. holarctica in Kazakhstan

Abstract

Tularemia is a highly dangerous zoonotic infection due to the bacteria Francisella tularensis. Low genetic diversity promoted the use of polymorphic tandem repeats (MLVA) as first-line assay for genetic description. Whole genome sequencing (WGS) is becoming increasingly accessible, opening the perspective of a time when WGS might become the universal genotyping assay. The main goal of this study was to describe F. tularensis strains circulating in Kazakhstan based on WGS data and develop a MLVA assay compatible with in vitro and in silico analysis. In vitro MLVA genotyping and WGS were performed for the vaccine strain and for 38 strains isolated in Kazakhstan from natural water bodies, ticks, rodents, carnivores, and from one migratory bird, an Isabellina wheatear captured in a rodent burrow. The two genotyping approaches were congruent and allowed to attribute all strains to two F. tularensis holarctica lineages, B.4 and B.12. The seven tandem repeats polymorphic in the investigated strain collection could be typed in a single multiplex PCR assay. Identical MLVA genotypes were produced by in vitro and in silico analysis, demonstrating full compatibility between the two approaches. The strains from Kazakhstan were compared to all publicly available WGS data of worldwide origin by whole genome SNP (wgSNP) analysis. Genotypes differing at a single SNP position were collected within a time interval of more than fifty years, from locations separated from each other by more than one thousand kilometers, supporting a role for migratory birds in the worldwide spread of the bacteria.

Fig 1. Positioning the Kazakhstan strains within the holarctica subspecies by wgSNP and Maximum Parsimony analysis.

Two thousand and twenty-two SNPs were called among 219 strains comprising 39 strains from Kazakhstan and 180 public data sets belonging to holarctica B.4, B.6 or B.12 (B.16 alias japonica strains were used to root the tree and are not shown). The red star shows the position of the MRCA. Branch lengths of more than ten SNPs are indicated. The branch scale is linear. The total tree size is 2034 (homoplasia 0.6%). The circles representing strains from Kazakhstan are colored according to canSNP assignment as indicated. The geographic origin of strains constituting secondary sublineages is indicated using three-letters country codes.

Fig 2. Focus on the 33 Kazakhstan strains belonging to sublineages B.23 and B.42.

Two hundred and fifty-four SNPs were called by mapping on reference SCHU_S4 (assembly accession GCA_000008985). The Maximum Parsimony tree size is 254 (no homoplasy). The red star indicates the position of the MRCA. The blue star indicates the root of the B.66 polytomy shown in a larger context in S2 Fig. Strains are labelled with year of isolation, host, strain Id and CanSNP assignment. The color code reflect the geographic origin within Kazakhstan (region level). Branch lengths longer than one are indicated. The branch scale is linear.

Fig 3. MLVA clustering of the 39 Kazakhstan strains.

The MLVA data for the seven polymorphic loci was clustered using the UPGMA method. The main canSNP sublineages are colored.