High and novel genetic diversity of Francisella tularensis in Germany and indication of environmental persistence

Summary

In Germany tularemia is a re-emerging zoonotic disease. Therefore, we investigated wild animals and environmental water samples for the presence and phylogenetic diversity of Francisella tularensis in the poorly studied Berlin/Brandenburg region. The phylogenomic analysis of three isolates from wild animals revealed three new subclades within the phylogenetic tree of F. tularensis [B.71 from a raccoon dog (Nyctereutes procyonoides); B.74 from a red fox (Vulpes vulpes), and B.75 from a Eurasian beaver (Castor fiber albicus)]. The results from histological, PCR, and genomic investigations on the dead beaver showed that the animal suffered from a systemic infection. Indications were found that the bacteria were released from the beaver carcass into the surrounding environment. We demonstrated unexpectedly high and novel phylogenetic diversity of F. tularensis in Germany and the fact that the bacteria persist in the environment for at least one climatic season. These findings support a broader host species diversity than previously known regarding Germany. Our data further support the assumption derived from previous serological studies of an underestimated frequency of occurrence of the pathogen in the environment and in wild animals. F. tularensis was isolated from animal species not previously reported as natural hosts in Germany.

INTRODUCTION

The highly infectious pathogen Francisella tularensis is the causative agent of tularemia, a zoonotic disease with an unusually wide range of host species [1]. Doses as low as 10 bacteria of F. tularensis may cause disease [2]. Outbreaks in humans may be associated with exposure to different arthropod vectors, infected animals, food, water, fomites, or contaminated aerosols. The clinical manifestations depend mainly on the route of infection. The onset of the disease is typically sudden and influenza-like, with symptoms including high fever, lymph node enlargement, chills, fatigue, headache, and nausea [3]. Two main subspecies of F. tularensis are clinically relevant, namely F. tularensis subsp. tularensis which is found only in North America and F. tularensis subsp. holarctica which occurs in most of Eurasia and North America [1, 3]. Three biovars have been described for F. tularensis subsp. holarctica: an erythromycin-sensitive biovar I (Ery^Sen), an erythromycin-resistant biovar II (Ery^Res), and biovar japonica which is distinguished by its ability to ferment glycerol [1]. More recently, whole genome sequencing and canonical single-nucleotide polymorphism (canSNP) genotyping methods have revealed multiple F. tularensis subsp. holarctica-type clades, including four basal clades B.4, B.6, B.12, and B.16 [4–6]. However, at present it remains unclear how the biovar classification system of F. tularensis subsp. holarctica relates to the phylogeny and current canSNP classification system of the bacterium, except for the B.16 subclade which correlates with biovar japonica and the B.12/B.13 clades correlating with erythromycin susceptibility [4, 7, 8].

In recent years, both multiple tularemia outbreaks and single human cases have been reported from almost the entire area of Germany, and tularemia is now considered a re-emerging zoonotic disease. The knowledge about human clinical cases as well as seroprevalence in humans has increased in recent years [9–14]. Indeed, an unexpectedly high rate (7·5%) of wild animals seropositive for F. tularensis was detected in the Berlin/Brandenburg region [15] and in other regions of Germany [12–14]. The reason for the relatively high (sero)prevalence seen in humans and wildlife is unknown. Therefore, we hypothesized that analysis of the genomic diversity of strains circulating in wild animals in the Berlin/Brandenburg region might indicate whether the studied region is affected with a single genotype of F. tularensis or whether a broader diversity may suggest repeated introductions of the pathogen into this area. Hence, we screened wild animals within the Berlin/Brandenburg region for F. tularensis and investigated the genomic diversity of the isolates.

METHODS

Samples

All samples from wild animals were obtained within the framework of the German rabies surveillance programme [16]. In total, we gathered 30 animal samples from two investigations. First, we collected 29 animals from 11 locations in the Berlin/Brandenburg region during 2012–2014 (Table 1). Our investigation was prompted by the finding of one dead Eurasian beaver (Castor fiber albicus) close to the town of Brandenburg in February 2012. The beaver was found by a voluntary guard of the nature reserve and the local veterinary health authority of the town of Brandenburg (Fig. 1). The region is called the ‘Erdlöcher von Wust’, a small man-made lake close to the river Havel (Fig. 1 c, red circle). This lake is part of a nature reserve. Second, we also included an uncharacterized F. tularensis isolate (isolate A63/63) recovered from the liver of a red fox (Vulpes vulpes; one of 305 animals analysed) found in the same region in 2008.

Fig. 1. Location where the carcass of the tularemic beaver was found. (a) Map of Germany showing the Berlin/Brandenburg region and the town Brandenburg at the Havel (•). (b) Location of sample sites 1–6 in Berlin/Brandenburg. 1, ‘Erdlöcher von Wust/Brandenburg’; 2, Great Plessower lake; 3, Great Wannsee lake; 4, Stölpchensee lake; 5, boar puddle/Wannsee; 6, Teltow sewer (sample sites 7 and 8, not shown, see Methods section). (c) ‘Erdlöcher von Wust/Brandenburg’ – location of the beaver carcass (A, red circle, sample site 1) and the additional water sampling sites (B, C, white circles) of this area are indicated. [(b) and (c), Google Maps/Earth images.]

Table 1. Selected wild animals examined for the presence of tularemia by qPCR

* F. tularensis (Ft) subsp. holarctica was also isolated.

To investigate further the presence of F. tularensis bacteria in the immediate surroundings, we collected water samples from the site where the dead beaver was found (Fig. 1 c, white circles) at ‘Erdlöcher von Wust’ during 2012 and 2013. In 2013, eight sites, but seven additional sites [in addition to Erdlöcher von Wust (site 1)] the Berlin/Brandenburg area of were sampled to screen for the presence of F. tularensis (Fig. 1 b; Erdlöcher von Wust/Brandenburg [1], Great Plessower lake [2], Great Wannsee lake [3], Stölpchensee lake [4], boar puddle/Wannsee [5], Teltow sewer [6], Nordufer sewer/Berlin [7], and Spree river [8]). The five dates of sampling were 1 March 2012, 4 September 2012, 3 March 2013, 6 August 2013, and 28 November 2013.

Organ samples (kidney, liver, spleen, lymph node, urethra, intestine, lung) were stored at −80 °C. We took 200 mg of each organ sample for homogenization using FastPrep tubes containing 1·4 mm ceramic spheres (MP Biomedicals, Germany) immediately after defrosting, and 1 ml sterile PBS was added to each tube. Homogenization was performed three times for 10 s using Precellys 24 (PEQLab Biotechnologie GmbH, Germany). The homogenates were used for cultivation and DNA extraction.

Histology

Tissue samples taken at necropsy from the beaver included the parotid lymph nodes, lung, heart, spleen, liver, kidneys, urethra, and ileum. These were fixed in 10% neutral-buffered formalin and embedded in paraffin wax. Native tissues from the same organs were additionally preserved by freezing at −20 °C. For histopathological examination haematoxylin-eosin (HE) stain was used according to standard protocols.

Culture and erythromycin susceptibility

We cultured 100 µl of each organ homogenate in medium T [17] [1% brain heart infusion broth (Difco Laboratories Inc., USA), 1% bacto tryptone (Difco), 1% technical casamino acids (Difco), 0·005 g MgSO₄ 0·01% FeSO_4, 0·12% sodium citrate, 0·02% KCl, 0·04% K₂HPO₄, 0·06% l-cysteine, 1·5% glucose], on CHAB agar plates [CHA (Difco; Bestbion, Germany), 1% brain heart infusion broth, 1% proteose-peptone, 1% d-glucose, 0·5% NaCl, 0·1% l-cystine, 1·5% agar, 9% sheep blood], or on Neisseria selective medium Plus (Oxoid, Germany) at 37 °C with 5% CO₂ for up to 3 days. Picked clones from water samples were cultivated on CHAB-PACCV (CHABplus 8 × 10⁴ U/l polymyxin B, 2·5 mg/l amphotericin B, 20 mg/l sulfamethoxazol, 100 mg/l cycloheximide, 4 mg/l vancomycin), and on GVPC agar plates (Heipha Dr. Müller GmbH, Germany): ACES-buffered charcoal yeast extract (BCYE) agar plates supplemented with 80 000 IE polymyxin B, 1 mg/l vancomycin, and 80 mg/l cycloheximide. Francisella isolates were grown on CHAB agar plates as described above.

Antimicrobial susceptibility was determined using MIC Evaluator Strips (Oxoid). Isolates were grown on CHAB at 37 °C with 5% CO₂ for 48 h. Four to five single colonies were inoculated in 0·9% NaCl until a turbidity equivalent to a 0·5 McFarland standard occurred. We then transferred 100 µl of the solution onto CHAB agar plates. The test strips were applied onto the inoculated agar surface, and the minimal inhibition concentration of erythromycin was determined after 48 h of incubation.

Molecular analyses

DNA preparation

DNA was extracted from isolates or organ homogenates using the DNeasy Blood and Tissue kit (Qiagen, Germany). Colonies grown on CHAB or 200 µl of the organ homogenates were directly transferred into lysis buffer containing proteinase K. Further steps were performed according to manufacturer's protocol. The purity and concentration of DNA were then determined by using either a Nanodrop spectrophotometer (PEQLab Biotechnologie GmbH) or a Qubit 2·0 fluorometer using the dsDNA BR Assay kit (Life Technologies, USA).

DNA preparation from water samples. A total of 510 ml water was collected at each site and filtered through a 40 µm filter to remove possible debris. Of this sample 500 µl were plated onto two CHAB-PACCV and two GVPC agar plates, respectively. Half of the plates were incubated at 30 °C and the other half at 37 °C for 4–7 days and inspected daily for putative isolates of Francisella sp. Isolates were picked on CHAB_PACCV agar plates and cultivated further for 3–5 days. Next, 500 ml of the sample were centrifuged at 5000 g at room temperature for 30 min. The pellet was resuspended in ~5 ml of sample water and DNA was isolated by using the ‘soil-DNA isolation kit’ from MoBio (Dianova, Germany). For negative control we performed DNA isolation procedure from tap water. We stored 100 µl of total DNA extraction volume at −20 °C until polymerase chain reaction (PCR) analysis.

PCR screening

PCR screening of bacterial colonies from water samples and organ samples for Francisella. A loop of bacteria from picked colonies was suspended in 250 µl water and incubated at 99 °C for 10 min. After centrifugation (10 min, 2500 g ), 10 µl of the supernatant (exhibiting chromosomal DNA) was investigated by PCR using 16s rDNA Francisella-specific primers, F5/F11 [18] to detect putative Francisella isolates. PCR was carried out using a Thermocycler TRIO-Thermoblock (Biometra, Germany) and the HotStar Taq DNA polymerase (Qiagen). Initial denaturation was performed at 95 °C for 5 min with a final extension at 72 °C for 10 min. The cycling conditions (35 cycles) were 94 °C for 1 min, 60 °C for 45 s, and 72 °C for 2 min. PCR amplicons were purified with the QIAquick PCR purification kit according to the manufacturer's instructions (Qiagen). Sequence determination was performed by applying the dye terminator chemistry (Applied Biosystems, Germany), and chromatograms were analysed with Lasergene Software (DNAStar, USA) followed by sequence comparison by BLAST search (http://blast.ncbi.nlm.nih.gov/Blast.cgi).

Quantitative real-time PCR (qPCR) and subtyping

Multiplex real-time PCR (5′ nuclease assay) targeting fopA and tul4 were used to analyse all suspected F. tularensis isolates for confirmation and DNA isolated from water samples. Primer and probe sequences for the detection of fopA were Ft-fopA-F: TTGGGCAAATCTAGCAGGTCA; Ft-fopA-R: ATCTGTAGTCAACACTTGCTTGAACA; and Ft-fopA-TM: FAM-AAGACCACCACCAACATCCCAAGCA-BHQ-1, and for the detection of tul4, Ft-tul4-F: AGATTACAATGGCAGGCTCC; Ft-tul4-R: AGCTGTCCACTT ACCGCTACA; and Ft-tul4-TM: Cy5-TTCTAAGTGCCATGATACAAGCTTCCCAA-BHQ-2. The real-time PCR assays were run in a total volume of 25 µl including 5 µl sample volume or DNA. The reaction mix components were TaqMan Environmental MasterMix 2·0 (6·25 µl), 10 pmol/μl primers (0·75 µl each) and probes (0·25 µl each). Amplification was performed in an ABI7500 (Applied Biosystems, Life Technologies, Germany), each run with an initial denaturation step at 95 °C for 10 min, followed by 35 cycles containing a denaturation step at 95 °C for 15 s, and a combined primer annealing and elongation step at 55 °C for 60 s. An internal amplification control (KoMa) was added [19].

F. tularensis subspecies differentiation was performed by targeting region of difference 1 (RD1) according to protocols described by Broekhuijsen et al. [20]. PCR was carried out using the high fidelity triple master mix system (Eppendorf, Germany) with 15–100 ng template DNA.

Variable number of tandem repeat (VNTR) typing

VNTR markers Ft-M3, Ft-M6, Ft-M20, Ft-M21, Ft-M22, and Ft-M24 were analysed as described by Johansson et al. [21]. PCR was performed using Eppendorf Mastercyclers (Eppendorf) with 50–100 ng bacterial DNA in the presence of 10 pmol/μl primers (1 µl), 2 mm dNTPs (2·5 µl), 25 mm MgCl (1·5 µl), and 1 U/μl Taq DNA polymerase (1 µl) (all by Fermentas, Germany), in a total volume of 25 µl. The cycling conditions were an initial denaturation at 94 °C for 4 min, followed by 35 cycles comprising denaturation at 95 °C for 30 s, annealing at 56 °C or 61 °C for 30 s depending on the primer used, and extension at 72 °C for 30 s with a final extension step at 72 °C for 4 min. Each PCR product was subjected to agarose gel electrophoresis (PEQLab Biotechnologie GmbH).

DNA sequencing

For genome sequencing, indexed sequencing libraries were prepared from 50 ng DNA with the Nextera XT DNA Library Preparation kit (cat. no. FC-131-1024) according to the manufacturer's recommendations (Nextera XT DNA Library Preparation Guide, no. 15031942 rev C, Illumina, USA). A 500 cycle (2 × 250 Paired End) MiSeq Reagent kit v. 2 was used to sequence the libraries on the MiSeq sequencing platform. The genomes were assembled using ABySS v. 1·3·5 [22]. All sequence reads were deposited and are freely downloadable from the Sequence Reads Archive at the National Center for Biotechnology Information (NCBI) under BioProject PRJNA285142 with BioSample accession nos. SAMN03773881: A63/63 (FDC407); SAMN03773882: A271_1 (FDC408); and SAMN03773883: A317 (FDC409).

Phylogenetic analysis

The three draft genomes were aligned with seven representative genomes from previously identified main clades of subspecies holarctica using progressive Mauve with standard parameters. The GenBank assembly accession numbers for the seven published genomes were GCA_000168775·2 (FSC200); GCA_000009245·1 (LVS); GCA_000017785·1 (FTNF0012-00); and GCA_000014605·1 (OSU18), and NCBI biosample SAMN00000315 (OR96-0246); SAMN01047935 (FSC162); and SAMN00092746 (Tul7). The phylogenetic tree was constructed in MEGA6 [23] using the neighbour-joining method, where all positions containing gaps and missing data were eliminated, and the distances were computed using the number-of-differences method [24]. The tree was mid-point rooted. A specific study of the cas-9 gene was made by comparison with the 10 genomes above, and 25 other genomes from Poland, Finland, and Sweden (genomes to be published elsewhere) within the B.39 subclade.

Genotyping of genomes

The three draft genomes were genotyped using canSNPer which is a fast and lightweight tool written in Python for genotype classification in draft genomes of clonal pathogens based on canSNPs [25]. The markers included in the analysis have been described previously: B.1–B.14 by Vogler et al. and Birdsell et al. [6, 26]; B.15–B.25 by Svensson et al. [5]; B.26–B.32 by Chanturia et al. [27]; B.33–B.38 by Gyuranecz et al. [28]; and B.39–B.43 by Karlsson et al. [29]. Newly identified clades were named according to this canSNP classification system for F. tularensis that to date, in addition to the 43 canSNPs used above by canSNPer, also include the subspecies holarctica canSNPs B.67–B.68 [30], and B.69–B.70 [31].

RESULTS

Our investigation in the Berlin/Brandenburg region during 2012–2014 revealed that three (10%) out of a total of 29 examined animals – one found dead (a beaver) and two hunted [a raccoon dog (Nyctereutes procyonoides) and a wild boar (Sus scrofa)] – tested positive by PCR analysis (Table 1). The species distribution of the 29 animals examined for the presence of F. tularensis by qPCR and culture was as follows: six raccoon dogs, five muskrats (Ondatra zibethicus), four red foxes, four beavers, three raccoons (Procyon lotor), four wild boars, one European otter (Lutra lutra), one American mink (Neovison vison), and one nutria (Myocastor coypus) (Table 1). In addition, F. tularensis was isolated from two out of the three PCR-positive animals. One isolate (A317) was obtained from the raccoon dog and three (A271_1-3) from different organs of the beaver. To our knowledge, this is the first report of a F. tularensis subsp. holarctica isolated from a beaver in Germany. We further investigated this case in detail.

The tularemia-positive beaver was an adult male of 17·5 kg weight and 119 cm length found dead in February 2012 (Fig. 1) at the ‘Erdlöcher von Wust’ in a nature reserve close to the river Havel (Fig. 1 c, red circle). The carcass was submitted for necropsy to the Berlin-Brandenburg State Laboratory in Frankfurt (Oder). At necropsy, the animal was found to be emaciated and uraemic. The left parotid lymph node contained a granulomatous and necrotizing lesion of ~8 mm in diameter (Fig. 2 a). Multiple microscopic (pyo-)granulomas and necroses were additionally present in this lymph node (Fig. 2 b, c) and the spleen (Fig. 2 d). Randomly distributed microscopic foci of fresh lytic necrosis were distributed throughout the liver. Cause of post-renal uraemia and death was a rupture of the severely dilated urinary bladder due to multifocal suppurative nephritis and subsequent occlusion of the urethra by a cast composed of a mixture of urinary gravel, cellular debris, and exudate. At necropsy, the abdominal cavity was filled with bloody urine after rupture of the urinary bladder and the urethra was morphologically occluded by a cast. However, a correlation of morphology and function is not always possible in cases of urolithiasis, especially when the urinary bladder is already ruptured at necropsy. Functionally, the occlusion seemed to be only partial at the time of the highest internal pressure shortly before final rupture of the urinary bladder, allowing some urine to leak from the carcass into the surrounding snow. F. tularensis-specific qPCR analysis was positive for the left parotid lymph node, kidneys, lung, and spleen. F. tularensis was cultured from the parotid lymph node, spleen, and lung (isolates A271_1-3).

Fig. 2. Macroscpic and microscopic lesions in the dead Eurasian beaver. (a) Focally extensive necrosis in the left parotid lymph node (arrows). (b) Left parotid lymph node; granulomatous (arrows) and necrotizing (arrowheads) lymphadenitis; H&E stain. (c) Left parotid lymph node; focal pyogranuloma; H&E stain. (d) Spleen; focal pyogranuloma (arrows) adjacent to an intact lymph follicle (arrowheads); H&E stain.

A qPCR specific for F. tularensis was conducted on water samples collected during 2012–2013 in the immediate surroundings of the site where the dead beaver had been found (Fig. 1 c, white circles), and seven more distant sites within the Berlin/Brandenburg region (Fig. 1 b). We detected F. tularensis DNA at sample site A (Fig. 1 c) in March 2012 and at all three samples sites (A–C) in September 2012 (Fig. 1 c, Supplementary Table S1). We could also detect F. tularensis DNA in snow that we collected around the dead beaver which was contaminated with blood-containing urine (data not shown). In contrast, all samples taken 1 year later (2013) were negative for the presence of F. tularensis by qPCR analysis (Supplementary Table S1). Cultures for F. tularensis were negative for all water samples.

For further comparative analysis of the F. tularensis isolates, we also included an uncharacterized F. tularensis isolate (isolate A63/63) recovered from the liver of a red fox (one of 305 animals analysed) from the same region in 2008, hunted in the context of the German rabies surveillance programme [16]. All isolates were assigned to the erythromycin-resistant biovar II type (data not shown) but differed regarding their VNTR profiles (Table 2). In addition, VNTR analysis demonstrated that all three isolates from the beaver (A271_1-3) were identical (data not shown). The genomes of the lymph node isolate A271_1 (beaver), A317 (raccoon dog) and A63/63 (red fox) were sequenced on the Illumina sequencing platform producing draft genome sequences of good quality and sequence depth (Table 3). The canSNPer program [25] using the draft genomes as input, assigned all three isolates to the basal B.12 clade confirming their assignment to biovar II Ery^Res. The beaver isolate A271_1 was assigned to the B.33/34 subclade; the A317 isolate from the raccoon dog to the B.12/13 clade; and the A63/63 isolate from the red fox to B.39. A phylogenetic analysis of the three new genomes and seven previously published reference genomes revealed branch lengths (number of SNPs) for the new branches and enabled the placement of the new genomes in their phylogenetic context (Fig. 3). New subclade names were assigned to the new branches and one representative canSNP per new branch was extracted (Fig. 3). The phylogenetic tree was based on a total of 1 707 622 bp in common between the ten genomes and covered 95% of the total length of the genomes of ~1·89 Mbp. The genomes were highly identical with 1537 base differences (0·09%) in this alignment, thus following a clonal inheriting pattern. In total, we extended the number of canSNPs with five markers that were designated B.71–B.75 according to the canSNP classification system. The branch leading to the beaver isolate A217_1 was designated B.75. The previous B.12/13 clade was divided into two sister clades: the relatively long branch leading to the raccoon dog isolate (A317) was denoted as B.71 and the sister branch was denoted as B.72, leading to all previously known branches within the basal clade B.12, the dominating basal clade in Eastern Europe from Scandinavia to the Black Sea (Fig. 3). The branch leading to the red fox isolate A63/63 was assigned to B.74, while its sister branch was designated as B.73. Examples of canSNP assays that could be developed for typing purposes are B.71 (position 1 330 179, ancestral SNP state G/derived SNP state A), B.72 (337274, G/T), B.73 (198 537, C/T), B.74 (1 304 315 G/A), and B.75 (1 827 247 G/A) in FSC200 (GenBank accession no. CP003862·1) reference genome. Furthermore, genome comparisons revealed that isolate A63/63 contained a B.39 subclade-specific deletion (9 bp, TAGTGATCT) within the 3′ site of the pseudo-cas9 gene [32], near the region Ft-M22 used for VNTR analysis. Preliminary analysis indicated that this was a canonical indel (insertion and deletion) for subclade B.39, based on comparison with 25 draft genomes from three countries that could be targeted in typing assays.

Fig. 3. Neighbour-joining phylogenetic tree showing the relationship between the three new genomes and seven published genomes. The evolutionary distances are in the units of the number of base differences per sequence. The tree is based on 1537 differences identified in an alignment of a total of 1 707 622 bp in common between the ten genomes. CanSNPs B.71–B.75 (circled in green) and strains FDC407, FDC408, and FDC409 are new to this study and are indicated by an asterisk (* this study). B.I, Biovar I; B.II, biovar II.

Table 2. F. tularensis subsp. holarctica isolates from Brandenburg region, Germany, investigated in this study

RKI, Robert Koch Institute; FOI, Swedish Defence Research Agency; VNTR, variable number tandem repeat.

* All three strains showed resistance to erythromycin.

† All three strains showed an identical PCR fragment size for VNTR markers Ft-M21 (396 bp/2 repeats), and Ft-M22 (254 bp/4 repeats).

‡ canSNPer version 1 [25] (http://github.com/adrlar/CanSNPer) (freely available).

§ See Figure 3, identified here.

Table 3. Statistics on the three new genome sequences*

RKI, Robert Koch Institute; FOI, Swedish Defence Research Agency.

* Sequenced on the Illumina MiSeq sequencing platform at FOI with 251 bp read length and paired-end read type.

† Assembled with ABySS v. 1.3.5 [22].

‡ Calculated by (total no. pairs × read length × 2)/1 800 000.

§ The average length of a set of sequences. The N50 length is defined as the length N for which 50% of all bases in the sequences are in a sequence of length L < N.

Altogether, VNTR-typing, erythromycin susceptibility testing, and full genome sequencing of the three F. tularensis isolates from the beaver, red fox and the new branch B.71 (of basal clade B.12) isolate from the raccoon dog were of the erythromycin-resistant biovar II type (Fig. 3, Table 2)

DISCUSSION

Our investigation in the Berlin/Brandenburg region during 2012–2014 revealed that three (10%) out of 29 animals examined were F. tularensis positive by PCR analysis. Despite the limited number of animals, our results are in good agreement with the findings of Kuehn et al. [15], indicating a seroprevalence of 7·5% on average for the Brandenburg area. This is the first report of F. tularensis isolated from a beaver in Germany, and we could corroborate the presence of F. tularensis in raccoon dogs and red foxes previously suggested by serological studies [15]. The unexpectedly diverse animal species found F. tularensis positive indicate that tularemia in wildlife is probably underdiagnosed.

We investigated the case of a dead, tularemia-positive beaver in detail. Culturing from environmental samples from locations close to the site where the dead beaver was found did not identify living F. tularensis in any of the environmental samples during 2012 or 2013. One explanation for this might be that cultivation of F. tularensis from environmental samples is difficult, mainly because the organism grows slowly and fastidiously, and its isolation can be prevented by overgrowth of the background flora despite using selective antibiotics. Another explanation might be the presence of a very low concentration of culturable Francisella bacteria within the environment. In contrast, qPCR analysis indicates that the F. tularensis bacteria were persistent in the aquatic environment during one climatic season, but apparently no further than that (see below). We could not detect significant levels of F. tularensis DNA by qPCR analysis of the distant sites or in the following year, indicating that the entire area investigated is not highly ‘contaminated’ by F. tularensis DNA.

The results of the pathological examination of the beaver (Fig. 2) indicated a chronic systemic infection with final exacerbation leading to massive replication of F. tularensis subsp. holarctica and subsequent contamination of the surrounding environment. We detected F. tularensis DNA in the bloody urine-contaminated snow collected around the dead beaver (data not shown). Contamination of the aquatic environment with F. tularensis by infected voles and beavers (carcasses) has been reported from Finland and the United States [33–36]. In addition, chronic shedding of bacteria by infected voles has been suggested [33]. In the case of F. tularensis contamination of the environment caused by beaver carcasses in the United States, the bacteria were not able to persist in the aquatic environment for a long time [34]. However, the source of infection and the route of transmission in the case of the beaver which died from tularemia remain unclear.

Altogether, this implicates a mechanism where animals with high bacterial load may act as local amplifiers [35] and contaminate the immediate environment, resulting in local hotspots of F. tularensis subsp. holarctica. Thus, mammals would act as local amplifiers to facilitate the spread of the disease, e.g. through contaminated water or by vectors [30, 37, 38].

The three recovered F. tularensis isolates from the beaver, a raccoon dog, and a red fox were of the erythromycin-resistant biovar II type (Fig. 3, Table 2) and were members of the B.33/34 subclade (beaver), B.12/13 subclade (raccoon dog), and B.39 subclade (red fox). Members of the B.33 subclade (i.e. branches B.33–B.38) have previously been found in Germany, Hungary, Austria, and Scandinavia [28, 39], whereas the beaver isolate constitutes a new subclade B.75. The B.39 subclade has not been described previously in Germany but has been found in Sweden, Finland, and Poland [6, 29, 31].

Altogether it is remarkable that our analysis of a small set of animals indicates an unexpectedly high diversity of strains present in the Berlin/Brandenburg area. The finding of a new B.12/13 subclade designated B.71 (Fig. 3) that was distantly related to all other strains previously found in the basal clade B.12 reveals gaps in the knowledge of the phylogeography of F. tularensis subsp. holarctica and shows that the diversity of tularemia is high also in non-Scandinavian countries. Previous assumptions that the high genomic diversity of F. tularensis in Scandinavia was indicative of Scandinavia as the source of the historical spread of tularemia to the rest of Europe [29] could be owing, at least in part, to a sampling bias in Scandinavia and insufficient sampling in other regions of Europe, including Germany. Recent studies of German F. tularensis subsp. holarctica isolates from brown hares (Lepus europaeus) and from an outbreak in non-human primates in animal facilities in Germany also revealed a high genetic and phenotypic diversity of the isolates investigated [13, 40]. In the study of Antwerpen et al., using isolates from animal facilities in Germany, only two strains belong to basal clade B.12, whereas 12 strains belong to the erythromycin-sensitive basal clade B.6 [40]. In the study of Müller et al., 14/31 isolates were assigned to basal clade B.12, 16 to basal clade B.6, and one to basal clade B.4. The identified basal clade B.12 strains all exhibit the same PCR fragment size (480 bp) for VNTR marker Ft-M24 and were localized primarily in Eastern Germany [13]. This is in agreement with our results, showing that all three animal isolates were members of three different B.12 subclades (B.39, B.33/34, and B.12/13). However, the beaver isolate A271_1 in subclade B.33/34 showed a PCR fragment size for VNTR marker Ft-M24 of 459 bp (Table 2) underlining that the exclusive use of VNTR markers for typing may give rise to incongruences in the reconstruction of phylogenetic relationships between strains, leading to potentially wrong interpretation of strain relationships in F. tularensis. However, the VNTR analysis of the three different isolates from the beaver (A271_1-3) demonstrated that all three were identical (data not shown), indicating absence within host diversity. Other recent examples of high Francisella diversity in Germany include an emerging subclade of F. tularensis subsp. holarctica in Ixodes ricinus ticks from south-western Germany [12] and the first German Francisella isolate (W12-1067) not belonging to the species F. tularensis [41].

In conclusion, based on the assumption that high genetic diversity indicates longer presence of a pathogen in a region, one could hypothesize from our data a long presence of tularemia in the study region within East Germany. Moreover, our findings demonstrate that further investigations on Francisella isolates from different animal species, the environment, and from patients are necessary to better understand the occurrence, persistence, and distribution of F. tularensis subsp. holarctica and also the presence of putative further Francisella species in Germany.