Surveillance Data

We analyzed surveillance data from human cases reported in Baden-Wuerttemberg during January 2012–December 2022. We included cases meeting the reference definition for tularemia in Germany (21) and excluded cases in persons with travel history and known risk exposures outside of Baden-Wuerttemberg.

The state veterinary authorities (Veterinary Analysis Agencies) in Baden-Wuerttemberg provided animal data from passive surveillance, primarily consisting of animal carcasses found by hunters, farmers, and hikers (7). We included all animals with available laboratory results for F. tularensis within the observation period.

We described human case characteristics and stratified human cases by meteorological seasons and years for temporal analysis with Poisson regression and calculated median time lag of symptom onset to reporting. For cases with missing onset dates (n = 17 [11%]), we estimated symptom onset by subtracting the median 29-day lag from the notification date. We used sample dates to calculate positivity rates in animal cases and stratified human and animal cases by districts. We used R version 4.1.2 (The R Project for Statistical Computing, https://www.r-project.org) for statistical analysis.

Bacterial Identification, Genome Sequencing, and Data Analysis

Animals identified through passive surveillance were tested for F. tularensis at the Veterinary Analysis Agencies. The isolates were further characterized at the Institute of Bacterial Infections and Zoonoses, Federal Research Institute for Animal Health (Jena, Germany). Cultivation, DNA extraction, and Illumina sequencing of 47 human and animal samples were performed as previously described (4,15,16).

To compare the genomes of strains collected within this study with previously published genomes from Baden-Wuerttemberg, we downloaded raw sequencing data of 7 human (16) and 61 animal (15) strains (Appendix Figure) and analyzed them with the WGSBAC pipeline version 2.2.3 (15). We controlled raw sequencing quality using FastQC version 0.11.7 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc) and calculated raw coverage and performed assembly of raw reads using Shovill version 1.0.4 (https://github.com/tseemann/shovill). We checked the quality of assemblies by using QUAST version 5.0.2 (22). We used Kraken2 version 2.1.1 to test for possible contamination (23). To detect subspecies-specific marker genes, we used in_silico_PCR version 0.5.1 (https://github.com/egonozer/in_silico_pcr) to search for the F. tularensis chromosomal marker RD-1 and checked the amplicon size for the specific length (924 bp) of F. tularensis subsp. holarctica.

We performed genotyping on the basis of predefined canonical SNPs with CanSNPer (24) and CanSNPer2 (https://github.com/FOI-Bioinformatics/CanSNPer2). In addition, mapping-based SNP typing with OSU18 as reference was performed with Snippy version 4.3.6 (https://github.com/tseemann/snippy), which also constructed a core-genome SNP (cgSNP) alignment of all strains. We used Snps-dists version 0.63 (https://github.com/tseemann/snp-dists) to calculate the pairwise SNP distances between strains. Using the cgSNP distance matrix, we performed clustering with the hierClust function version 5.1 of R (https://www.rdocumentation.org/packages/momr/versions/1.1/topics/hierClust). We used a cutoff of 5 cgSNPs to define genomically highly similar strains (15). We constructed phylogenetic trees with RAxML version 8 (25) and visualized with the interactive Tree of Life (iTOL) version 4 web tool (26).

Ethics Statement

Ethical approval was not required because no personal identifying information was processed. Tularemia is a notifiable disease, and the analysis of routinely collected surveillance and sequencing data are part of the health authorities’ responsibilities under the Protection Against Infection Act (German: Infektionsschutzgesetz).

Epidemiology of Notified Human Cases

In Baden-Wuerttemberg, 152 tularemia infections, including 5 deaths, were reported during 2012–2022. Annual cases ranged from 3 to 34 (median 10); annual incidence rates were 0.27–3.1/1 million population. Most cases (n = 34) were reported in 2021, representing a 3-fold increase of cases over the previous year (n = 11). The median lag time from symptom onset to reporting was 29 days; 18 cases were reported after >96 days.

Of the 152 cases, 106 (70%) were in men and 46 (30%) in women (Table 1); the median patient age was 53 years (interquartile range 37–65 years). More than half of patients (66%, n = 85) were hospitalized. The most common laboratory confirmation was serology (67%, n = 102) followed by PCR assay only (22%, n = 33) and PCR plus culture (20%, n = 31). Symptoms were reported in all but 1 case; fever (75%, n = 113) and lymph node swelling (69%, n = 104) were most common, followed by skin ulcer (24%, n = 36), pneumonia (15%, n = 23), and dyspnea (11%, n = 17). The main forms of tularemia were glandular (35%, n = 53) and ulceroglandular (23%, n = 35). Pneumonic tularemia occurred in 17% (n = 26) of the cases, mostly reported in 2019 and 2022 (n = 18); 8 cases occurred in 2021. Most cases with dyspnea (n = 12) were reported during 2019–2022; 9 of those cases were reported in 2021 and 2022. Less than 10% of cases showed symptoms of oropharyngeal tularemia; 12% of cases (n = 18) had unspecific symptoms and unknown exposure risk, possibly indicating typhoidal tularemia.

Of 142 cases for which information on animal exposure was available, 68% (n = 96) had known exposure, primarily to arthropod bites (45%, n = 43), including 28 tick bites. More than half of arthropod-borne infections occurred during 2019–2022 (n = 25). Among those 142 cases, <15% of cases reported contact with wild animals and 35% reported contact with domestic animals. Information on the consumption of contaminated food was available for just 1 case. Four case-patients reported the slaughter of wild animals, but no specific exposure details were available. Of 130 cases with professional exposure information, 31% had potential professional exposure, mainly from hunting and forestry work (43%, n = 17) or farm work (28%, n = 11). Of 84 cases that had other exposure information, 39% (n = 33) reported exposure to potentially contaminated environmental sources, including soil, mud and water sources.

Poisson regression (Figure 1) showed a significant annual increase in tularemia cases by 20% (1.20, 95% CI 1.14–1.27; p<0.001). Periodicity could not be identified in data. Seasonal distribution of cases varied (Figure 2); the lowest case numbers were reported in winter. During 2021 and 2022, when case numbers were higher, cases were almost equally distributed across all seasons except winter.

Epidemiology of Notified Animal Cases

During 2012–2022, the Veterinary Analysis Agencies tested 2,145 animal samples for tularemia. Most samples were from brown hares (86%, n = 1,839), whereas other samples were collected from various animals, including wild boar, rodents, and domestic animals. Laboratory results were available for 2,138 samples, of which 564 (26%) tested positive for F. tularensis. Of those, 551 (98%) tested samples were from brown hares, with the rest from unspecified ground game (n = 7), wild birds (n = 4), 1 zoo animal, and 1 wild boar. Positivity rates fluctuated from 5% in 2012 to 34% in 2016 (Figure 3) but have steadily increased from 16% to 33% since 2017, alongside an increase in the number of samples tested.

Geographic Distribution of Human and Animal Cases

Of the 44 districts in Baden-Wuerttemberg, 37 reported human or animal tularemia cases during 2012–2022 (Figure 4). Most human cases were reported from the districts of Karlsruhe (n = 15) and Ortenaukreis (n = 11). The districts Boeblingen, Enzkreis, Konstanz, and Rhein-Neckar-Kreis reported cases for the first time in 2021, and the city district of Heilbronn reported the first case in 2022.

Most animal cases were reported in the districts Rems-Murr-Kreis (n = 101), Schwaebisch Hall (n = 63), Ostalbkreis (n = 47), and Ludwigsburg (n = 41). In the other districts, the number of laboratory-confirmed animal cases ranged from 1 to 28 cases. Three of the districts reported human cases and no animal cases. The other districts reported animal cases but no human cases.

Genomic Surveillance of F. tularensis subsp. holarctica in Baden-Wuerttemberg

We analyzed genome sequences of 115 F. tularensis subsp. holarctica, including 47 strains sequenced in this study and publicly available genome data of 68 strains (Appendix Table 1). All sequencing data produced in this study were added to the European Nucleotide Archive (Bioproject PRJEB40963). Twenty-four strains were isolated from human cases, whereas animal-derived strains were mainly isolated from brown hares (Lepus europaeus, n = 90); 1 strain was isolated from a Eurasian blue tit (Cyanistes caeruleus).

A minimum of 30-fold coverage was reached for all strains (Appendix Table 2). On average, 99.4% of raw reads belonged to the genus Francisella and 97.6% to F. tularensis subsp. holarctica. The chromosomal marker RD-1 for F. tularensis was detected in all genomes with an amplicon size of 924 bp, corresponding to F. tularensis subsp. holarctica. The mean assembly size was 1.79 Mbp, and the average N50 of the 115 assembled genomes was 42,988 bp (range 13,573–100,618 bp). The mean GC content was 32.2% with little variation.

Canonical SNP-typing uses predefined lists of SNPs to assign major clades and subclades to each strain following a predefined decision tree. Two major clades were detected among the studied strains; most strains belonged to the major clade B.6 (n = 104) for both human isolates and animal-derived strains (Table 2; Appendix Table 1). Using CanSNPer1, we detected 14 canonical subclades. Although 6 subclades (B.33, B.34, B.45, B.49, B.61, and B.63) were shared between human and animal strains, 8 additional subclades were found only in animals. CanSNPer2 further subbranched many clades detected by CanSNPer1into 34 subclades. Within the dataset, 5 subclades (B.182, B.259, B.287, B.305, and B.91) were found only in human-derived strains, 20 were found only in animals, and 9 subclades were detected among strains from humans and animals. The strain from a Eurasian blue tit was isolated from a dead bird in the city of Freiburg im Breisgau. The sample was identified as F. tularensis subsp. holarctica by PCR and matrix-assisted laser desorption/ionization time-of-flight mass spectrometry and was typed as clade B.6, subclade B.49, which was also detected in several strains isolated from hares and humans in Baden-Wuerttemberg.

The variable core-genome alignment contains all genomic positions of the reference genome, which exist in all of the strains and are variable in >1 strain. This alignment consisted of 24,359 bp and enabled the construction of a phylogenetic tree (Appendix Figure). The phylogeny of F. tularensis in Baden-Wuerttemberg constitutes 2 major branches in accordance with the major canonical clades B.6 and B.12, whereas clusters of genetically similar strains mostly follow the canonical subclades. Strains isolated from patients are closely related to strains obtained from animals, because human-derived strains do not form outgroups and they cluster with animal strains. The pairwise number of cgSNPs between strains ranged from 0 to 23,595 cgSNPs with a median of 8 cgSNPs (Appendix Table 1). Hierarchical clustering was performed to group strains into genomically highly similar clusters. Using a cutoff of 5 cgSNPs (1), 83/115 (72%) strains were clustered into 15 groups (Appendix Table 1). Although 11 clusters contained strains that were exclusively isolated from animals, 4 clusters (clusters 1, 3, 8, and 10) contained strains from humans and animals (Figure 5; Appendix Table 1). The largest cluster (cluster 1) contained 26 strains (Figure 5, panel A), of which 5 were isolated from patients during 2017–2022 at different locations in Baden-Wuerttemberg. For example, the F. tularensis strain A1174 isolated from a patient in the city of Karlsruhe was 1 cgSNP different from 4 isolates from hares in southern Baden-Wuerttemberg (Appendix Table 1). Among those 6 isolates, 4 case-patients had been exposed to arthropod bites (3 tick bites and 2 mosquito bites). One of the case-patients had no known risk exposure but reported owning a cat. One case was in a farm worker who was reportedly exposed to sewage through work. Cluster 3 contained 6 strains, including the human isolate A-1667 from 2020 in Heilbronn and the hare-derived strain 22T0158 from the same city (Figure 5, panel B). The strain 16T0054 was isolated from a hare in Stuttgart and 1 cgSNP different from A-1667. No exposure risks were reported in that case, and the mode of transmission remained unclear. Two strains from patients and 1 strain from a hare, all isolated in 2021, grouped within cluster 8 (Figure 5, panel C). Strain A1871 was isolated from a case-patient without exposure risk at work but who was active in a private garden and, therefore, potentially had contact with rodent feces. Its genome was 1 cgSNP different from the hare-derived strain 21T0101 (Appendix Table 1). Strain A1893 was isolated from a farmer who had contact with mice and feces from rodents. Cluster 10 contained 1 human and 1 animal-derived strain. The case-patient reported contact with a dog but no other risk factors (Figure 5, panel D).

In general, canonical (sub)clades of CanSNPer1 are either identical when using CanSNPer2 or split into several subclades (Appendix Table 1). Using CanSNPer1, we found 51/115 strains (44%) belong to subclade B.45. Of those 51 strains, 15 belong to B.45 when using CanSNPer2, whereas 36 belong to 10 different subclades. Subclades are comparable to hierarchical clusters. For example, hierarchical clusters 3, 4, 5, 7, and 11 are completely identical to respective subclades identified by CanSNPer2. Cluster 1 is split into 5 subclades, whereas subclade B.63 is split into hierarchical clusters 6 and 8.