Population Genomics of Francisella tularensis subsp. holarctica and its Implication on the Eco-Epidemiology of Tularemia in Switzerland

Whole genome sequencing (WGS) methods provide new possibilities in the field of molecular epidemiology. This is particularly true for monomorphic organisms where the discriminatory power of traditional methods (e.g., restriction enzyme length polymorphism typing, multi locus sequence typing etc.) is inadequate to elucidate complex disease transmission patterns, as well as resolving the phylogeny at high resolution on a micro-geographic scale. In this study, we present insights into the population structure of Francisella tularensis subsp. holarctica, the causative agent of tularemia in Switzerland. A total of 59 Fth isolates were obtained from castor bean ticks (Ixodes ricinus), animals and humans and a high resolution phylogeny was inferred using WGS methods. The majority of the Fth population in Switzerland belongs to the west European B.11 clade and shows an extraordinary genetic diversity underlining the old evolutionary history of the pathogen in the alpine region. Moreover, a new B.11 subclade was identified which was not described so far. The combined analysis of the epidemiological data of human tularemia cases with the whole genome sequences of the 59 isolates provide evidence that ticks play a pivotal role in transmitting Fth to humans and other vertebrates in Switzerland. This is further underlined by the correlation of disease risk estimates with climatic and ecological factors influencing the survival of ticks.

Introduction

Classification of organisms according to inherent characteristics of their genome has become an indispensable principle in molecular biology. Over the years, progress in sequencing technologies led to an explosion in the amount of available sequencing data and an ever-increasing taxonomic resolution. This development culminated in the advent of massive parallel sequencing methods that allow the classification and characterisation of an organism on a whole genome scale, opening new perspectives in molecular forensics, ecology, and epidemiology. The benefits of whole genome sequencing (WGS)-based classification are most pronounced when applied to genetically monomorphic organisms, for example Mycobacterium tuberculosis (Lee et al., 2015; Stucki et al., 2016), Bacillus anthracis (Girault et al., 2014) and Francisella tularensis (Dwibedi et al., 2016).

Tularaemia is caused by the gram negative, facultative intracellular bacterium F. tularensis. Due to its very low infective dose and high mortality when inhaled as aerosol, the organism is listed as a category A biothreat agent (Rotz et al., 2002). Two subspecies are responsible for most of the tularemia cases. F. tularensis subsp. tularensis (Ftt), and F. tularensis subsp. holarctica (Fth). In Europe, the less virulent Fth occurs and generally causes sporadic cases, of which the first were documented in the 1930s (Jusatz, 1961). Fth, however, occurs throughout the northern hemisphere. Since the late 1990s, case numbers are increasing and outbreaks exceeding 1,000 cases were reported in the Czech Republic, Hungary, Spain (Pérez-Castrillón et al., 2001), Sweden and Finland (Maurin and Gyuranecz, 2016). It is thought that in the United States the majority of F. tularensis infections are due to transmission from tick bites (Petersen et al., 2009) whereas in Europe, the situation is less clear. In Scandinavia, Fth is transmitted predominantly by mosquitoes and it is believed that the pathogen persists in an aquatic life cycle with mosquitoes, mosquito larvae and rodents (Desvars et al., 2015). The terrestrial lifecycle, with arthropods as reservoirs and small terrestrial rodents and lagomorphs, as susceptible hosts, is predominant in most European countries including Switzerland. In these countries, direct contact with infected rodents and lagomorphs seems to be the main route of infection. Concerning the transmission through arthropod vectors, the estimated percentage of tularaemia cases due to tick bites varies between 12% (Slovakia) and 26% (France) (Maurin and Gyuranecz, 2016). Since tularaemia is a notifiable disease in Switzerland, epidemiological data is collected from patient reports sent from the initial point of care to the health authorities. According to this information, tick bites could be associated with 47% of the tularaemia cases reported during the last 10 years. However, the quality of questionnaire-based data is often limited. Additionally, the confirmation of a causal connection between arthropod bites and tularaemia cases using a method providing adequate discriminative power has been lacking. In this study the epidemiology, routes of transmission, and phylo-geographic properties of tularaemia in Switzerland are delineated based on 59 sequenced genomes of Fth isolated from humans, animals and ticks (Table 1).

TABLE 1

Table 1. Overview of the 59 whole genome sequenced Fth strains from Switzerland used in this study.

Materials and Methods

Collection of Ticks

Between 2009 and 2015 a total of 120,000 questing ticks were collected at 165 collection sites throughout Switzerland by flagging low vegetation. Ticks were randomly identified based on morphological characteristics (Keirans, 1997) and immediately stored at −80°C. Subsequently, ticks were washed once in 75% ethanol and twice in deionized water, dried on paper towels, and sorted into pools of 10 nymphs or 5 adult male or female I. ricinus ticks. Pooled samples were stored at −80°C until further processing (Gäumann et al., 2010).

Homogenization of Ticks

Tick homogenates were prepared from frozen tick pools of 10 nymphs or 5 adult female or male ticks as follows: 600 μl of cold PBS was added to each frozen tick pool. Samples were immediately homogenized using the TissueLyser system (Qiagen). Briefly, one 3-mm tungsten carbide bead (Qiagen) was added to each tube (collection microtubes; Qiagen), and tick pools were homogenized for 4 min at 30 Hz. After centrifuging the samples for 5 s at 3,220 g, supernatants were collected and split for further use: One 200 μl aliquot was used for DNA extraction and four 150 μl aliquots were immediately frozen in liquid nitrogen and stored at −80°C. To increase success rates of subsequent cultivation efforts, 40 μl glycerol (100%) was added as a cryoprotectant to two aliquots.

Screening of Tick Homogenates for the Presence of Fth

DNA was automatically extracted using a magnetic bead based protocol (QiaSymphony, QIAGEN). In total, 16,000 tick pools were screened for the presence of Fth using a one-step real-time RT-PCR assay based on the VNTR marker Ft-M19 (Byström et al., 2005). The Ft-M19 PCR positive samples were confirmed by an in-house real-time RT-PCR assay targeting the fopA gene (Table 2).

TABLE 2

Table 2. PCR assays targeting the Francisella Outer Membrane Protein A used for F. tularensis ssp. screening.

Isolation of Fth From Homogenized Tick Samples

Aliquots of the PCR positive samples were cultivated. Briefly, frozen aliquots were quickly thawed and 50 μl each was inoculated in 10 ml Broth Medium T (Becker et al., 2016), as well as streaked on Chocolate Agar PolyViteX (BioMérieux) and Thayer Martin VCNT Neisseria Selective Agar (Oxoid). Plates and liquid cultures were incubated at 37°C using GENBox C02 (BioMérieux) for 2–6 days. Presumptive Fth cultures were confirmed by real-time PCR as described above (fopA; Table 2).

DNA Isolation for WGS

All manipulations with live cultures were performed in a BSL 3 containment laboratory (approval number: A110502/3). Isolates of Fth were grown on Chocolate Agar PolyViteX (BioMérieux) for 2 days and colonies were harvested and suspended in 500 μl AVL Buffer (Qiagen). Lysates were heat inactivated for 15 min at 100°C and DNA was subsequently extracted using the EZ1 or the EZ1 Advanced robot (Qiagen) according to the manufacturer's instructions (EZ1, Tissue Kit, Bacteria Card, Qiagen).

Human Samples

In its function as the Swiss national reference center for tularemia in humans, the Spiez Laboratory maintains a culture collection of Fth isolates derived from routine clinical diagnostics. The patient and epidemiological data are collected by the Federal Office of Public Health according to the Epidemics Law and Ordinance (SR 818.101.126).

Animal Samples

Strains isolated from animals were selected from the strain collection of the Institute of Veterinary Bacteriology, Vetsuisse Bern, Switzerland, which is the national reference center for tularemia in animals.

Whole Genome Sequencing

The 59 Fth isolates were sequenced on the Illumina HiSeq 2500 platform using the TruSeq paired end chemistry for library preparation. The average read length was 126 bp and the average insert size was 350 bp (standard deviation = 87.593, lower quantile = 272, upper quantile = 488). Prior to the assembly the ~50,000,000 sequencing reads were quality trimmed with Trimmomatic V0.32 (Bolger et al., 2014).

The quality trimmed reads were assembled using SPAdes-3.10.1 (Bankevich et al., 2012) resulting in ~100 contigs with an average N50 of 27,000 bp. The coverage of the Fth genomes was >1,000 × in all samples. In order to use the pan genome pipeline Roary (Page et al., 2015), gff3 annotation files were generated with Prokka (Seemann, 2014) using the Franco-Iberian Fth FTNF002-00 (NC_009749) strain as reference. Based on these gff3 files, the core- and pan-genomes were inferred using Roary by applying a 95% blastp identity threshold. A given gene had to be present in 100% of the samples to be included into the core-genome. Mafft (Katoh et al., 2002) was used for core genome alignment. Ambiguously aligned regions and Indels were removed from the Mafft alignment by Gblocks (Castresana, 2000) and SEQfire (Ajawatanawong et al., 2012). Maximum likelihood phylogenies were inferred with PhyML (Initial tree: BioNJ; Tree topology search: NNIs; Model of nucleotide substitution: TN93; Bootstrap replicates: 1000) (Guindon et al., 2010). Mutations in the sequenced strains relative to the reference genomes were identified with breseq (Deatherage and Barrick, 2014).

To designate the sequenced strains to the established canSNP node nomenclature the software canSNPer (Lärkeryd et al., 2014) was used. Naming of the novel subtypes found in this study is in agreement with the canSNP nomenclature scheme assigning an increasing ordinal number to each newly identified SNP. To implement the novel subtypes in the canSNPer tool, the text files containing the SNP positions and the tree topology were adjusted accordingly.

Calculation of average nucleotide identity using standard MUMmer algorithm (ANIm) from whole genome was done using a python program pyani (v0.2.4). Hierarchical cluster analysis of ANI correlation matrix was done using package hclust (R v3.4.3).

Spatial Statistics

Phylogeographic interpolation was performed with the R package Phylin (Tarroso et al., 2015). Calculation of spatial disease clusters was performed (carried out) with the R package DCluster (Hornik et al., 2003).

Results

Prevalence of Fth in Swiss Ticks

A total of 120,000 questing ticks were analyzed and only 25 tick homogenates were positive for Fth by PCR. This corresponds to a prevalence of ~0.02%.

Fourteen Fth isolates were successfully recovered from the 25 positive tick homogenates. The recovery from glycerol preserved samples was slightly improved and the isolation was often easier from Neisseria Agar than from Chocolate Agar due to reduced contaminating flora.

Phylogenetic Context of Swiss Fth Isolates

ANIm is a standard in silico method for bacterial species delimitation with a threshold value of 94–96%. Here, ANIm was applied to validate that 59 isolates belong to F. tularensis subsp. holarctica. Whole genomes were compared to type strains: F. tularensis subsp. tularensis NIH B-38, F. tularensis subsp. holarctica FTNF002-00, F. tularensis subsp. mediasiatica FSC147, and F. tularensis subsp. novicida U112. In this small dataset, we found that an ANIm threshold of 99.5% could distinguish subspecies of F. tularensis. All 59 isolates showed >99.9% ANI identity to F. tularensis subsp. holarctica FTNF002-00 and F. tularensis subsp. holarctica FSC200.

To put the 59 Swiss Fth whole genome sequences in a global phylogenetic context we used a canonical SNP (canSNP) approach (Lärkeryd et al., 2014). According to this typing method, 54 strains (91.5%) were assigned to the B.11 clade, which is predominant in Western Europe including Switzerland (Pilo et al., 2009; Vogler et al., 2009; Dwibedi et al., 2016). Three animal and two human isolates (11%) belong to the Northeast European B.12 clade (Figure 1).

FIGURE 1

Figure 1. Overview of the current canSNP typing scheme. Green blobs indicate the canSNP types found in the 59 sequenced Swiss Fth isolates.

To increase the phylogenetic resolution a core genome SNP analysis based on a 1477202bp alignment was applied and revealed a similar branching pattern of the B.11 clade as described in the European survey by Dwibedi et al. (2016). Clade B.45 was represented by 44 (18 human, 14 tick, 12 animal) and B.46 by four isolates (three human, 1 animal). In addition, six isolates (2 human, 4 animal) were assigned to a clade, which was previously represented by one Swiss isolate only. This novel B.11 subclade is distinguished by 5 SNPs (B.82/B.83/B.84/B.85/B.86; Figure 2).

FIGURE 2

Figure 2. Phylogenetic relationships based on the core genome alignment of 59 sequenced Fth strains. The clades and subclades are labeled according to the canSNP nomenclature. The ending letters of the tip labels indicate the source of the isolates (H, Human; A, Animal; T, Tick).

Except of the clades B.48 and B.52 exclusively found in Spain, as well as the French B.54 clade, all subclades derived from B.45 were represented among the 59 sequenced isolates. In addition, five novel subclades deriving from B.45 could be identified (B.87/B.88, B.89, B.90, B.91, B.92; Figure 2).

The spatial dependence of the genetic variation among the isolates is represented by the semi-variogram shown in Figure 3. The good fit (R² = 0.96) of the applied spherical model clearly indicates a statistically significant dependence between genetic variation and geographic distance. The Kriging interpolation (Krige, 1951), based on the model derived from the semi-variogram, shows that clade B.45 is spread over the majority of the Swiss territory below 1,500 m of altitude. The subclades of B.45 that show a higher spatial coherence are visible as focal hotspots (Figure 4). The clades B.46 and B.86 show a geographically more focused distribution (Figure 5). Three animal and one human isolate assigned to the Northeast European B.12 clade originated from areas close to the northern and eastern Swiss border. One human isolate was isolated in central Switzerland.

FIGURE 3

Figure 3. Semivariogram showing the spatial dependence of the genetic similarity between two strains as a function of distance. The size of the circles is proportional to the number of strains used for variance calculation. The red line indicates the fit of the spherical model used to describe the spatial structure.

FIGURE 4

Figure 4. Spatial interpolation of phylogenetic clusters by krigin. Due to the high genetic diversity of Fth, the geographic distribution of the strains is shown in two separate maps. This Figure shows the interpolated distribution of the B.45 subclades. Figure 5 comprises the data for subclades B.33, B.86 and B.46.

FIGURE 5

Figure 5. Same as Figure 4.

Genetic Variation in the Swiss Fth Population

Breseq analysis revealed 259 polymorphisms within the 59 sequenced Fth genomes relative to the genome FTNF002-00. We identified 36 deletions, 28 insertions and 195 SNPs. Seventy-nine SNPs were non-synonymous. Of the 116 synonymous SNPs, 46 were located in coding regions of the genome. Regarding the mutations that were exclusively found in the newly defined B.86 subclade, 9 of the in total 17 SNPs were non-synonymous. Two non-synonymous mutations were found in genes associated with the heme (hemA) and cytochrome c (FTA_0224) biosynthesis. One mutation alters a member (FTA_0090) of the major facilitator superfamily (MFS) which is involved in the transport of small molecule substrates across all classes of organisms and plays an important role in sustaining the osmotic balance between the cytosol and the environment. The remaining 5 non-synonymous SNPs are found in genes coding for proteins involved in the carbohydrate (FTA_0510, FTA_1026, FTA_1126) and fatty-acid (FTA_0617) metabolism of Fth (Table 3).

TABLE 3

Table 3. List of the non-synonymous SNPs identified with the breseq pipeline.

All five strains assigned to the B.12 clade showed the SNPs in the rrl (23S rRNA) gene (A2059C and A453G) which were found to be responsible for the erythromycin resistance, which is characteristic for this clade (Karlsson et al., 2016). These findings are in agreement with the antibiotic susceptibility testing and MLVA typing results of Origgi et al. performed on the same isolates (Origgi et al., 2014).

All genetic variations found in the core genome alignment with the strain FTNF00-002 as reference are included in a variant call file (vcf) as Supplementary Material (Supplementary Table 1).

Epidemiology of Tularemia in Switzerland

Since 2004, the incidence of tularemia in Switzerland has risen ~15-fold from 0.04 to 0.69 per 100,000 per year (Figure 6). With the exception of patients >75 years of age, males were generally more affected by tularemia. We observe a bimodal distribution of case numbers per age group with a narrow peak between 10 and 12 years of age and a broader peak spanning the ages from 27 to 75 years (Figure 7).

FIGURE 6

Figure 6. Increase in the absolute case numbers and incidence of tularemia in Switzerland from 2004 to 2016.

FIGURE 7

Figure 7. Bimodal age distribution of tularemia patients. Men are more likely to be affected than women.

In agreement with a terrestrial lifecycle, Fth isolates from humans, animals and ticks are represented in the same terminal subclades. In subclade B.51, the tick derived isolates FT40/FT47 show no differences in the 1,474,944 bases core genome alignment with the human FT75 isolate. The same applies for the isolates FT05 (human) and FT58 (animal) in clade B.90.

The majority of the 276 patients reported an insect or arthropod bite prior to the occurrence of tularemia symptoms (33% tick-bite, 24% insects—total 57%). Contact with wild animals was assumed as the cause of 22% of the infections. In 23% of the cases, the source was unknown. The proportion of the cases with a suspected vector-borne source of tularemia infection correlated with the clinical manifestation of the disease. The glandular/ulceroglandular type typically associated with the bite of a hematophagous arthropod was observed in 60% of the cases. The pulmonary type was represented in 24% and the abdominal/oropharyngeal type in 14% of the patients (Figure 8). Hence, 38% of the 276 reported infections can be attributed to direct or indirect contact with infected wild life.

FIGURE 8

Figure 8. Classification of tularemia in Switzerland between 2004 and 2016 according to the clinical manifestation.

Population corrected risk estimation on the spatial distribution of tularemia cases shows an uneven distribution. Areas with elevated relative risk are predominately located in the northeastern area of Switzerland (Figure 9). Our data suggest that ticks are the predominant source of infection in Switzerland. Therefore, we correlated the risk estimations with climatic as well as ecological factors known to influence the survival and vector competence of ticks. Relative humidity (RH) and saturation deficit (SD) are known to be of pivotal importance for tick abundance and survival. Consequently, we compared the surface normalized averages of waterlogging levels in areas with elevated risk for acquiring tularemia vs. low risk areas and found a highly significant effect of soil moisture (Figures 10, 11).

FIGURE 9

Figure 9. Spatial distribution of cases by postal code area (first two digits). The estimated relative risk is shown, which is defined by the ratio of observed cases to expected ones. Most areas have less cases than expected, e.g., the relative risk is below 1 (blank zone). The remaining areas are shaded from orange to red according to their relative risk value.

FIGURE 10

Figure 10. Geographical overview of the water saturation of the soil (waterlogging). In general, the soil in the north-eastern part of Switzerland shows a higher degree of water saturation compared to other regions.

FIGURE 11

Figure 11. Regions with an elevated relative risk of tularemia show significantly higher degrees of waterlogging compared to low risk regions (p = 1.2 * 10⁻¹⁵).

The reduction of biodiversity due to urbanization influences the competence of vectors by shifting the ratio of competent to non-competent hosts. Most of the regions with an elevated risk for acquiring tularemia are highly fragmented and show below average mesh sizes (Figure 12).

FIGURE 12

Figure 12. Geographical overview of the degree of landscape fragmentation of the Swiss cantons. Most of the regions with a high risk for acquiring tularemia are also highly fragmented.

Discussion

Prevalence, Genetic Diversity, and Population Genomics of Swiss Fth Isolates

The determined prevalence of 0.2‰ is low compared to other pathogens found in ticks (Oechslin et al., 2017) and confirms principally a previous smaller survey of the Swiss Army (Wicki et al., 2000).

The successful cultivation of Fth from ticks substantiate the role of I. ricinus as important vector of the pathogen. Furthermore, the availability of isolates was a prerequisite for the WGS inferred high-resolution phylogeny.

As described elsewhere (Origgi et al., 2014) the B.13 and B.11 strain are co-circulating in Switzerland with a predominance of the latter. Strikingly, the geographic origin of four of the five B.13 isolates is close to the northern and eastern Swiss border. Since the B.11 and B.13 strains are equally established in the German hare population (Müller et al., 2013) and B.13 prevails in Austria (Tomaso et al., 2005), our data suggests a sporadic intrusion of the B.13 strain from neighboring countries.

In accordance with recent publications regarding the dispersal of Fth in Europe (Dwibedi et al., 2016), Switzerland seems to harbor an exceptionally high genetic diversity of the pathogen compared to other European countries. The phylogenetic structure and proportion of the B.11 subclades derived from Swiss isolates are almost identical to the one reported in the European survey. In this aspect, Switzerland can be considered as a small-scale model for Fth population structure. Moreover, the identification of six additional B.11 subclades underlines the notion that the alpine region harbors an evolutionary older and highly diverse founder population (Dwibedi et al., 2016).

A major factor that is positively associated with the biodiversity of an ecosystem is the topographic complexity of a habitat (Zhou et al., 2015). Together with climate-related factors, topography is known to influence a wide array of environmental parameters including hydrology, nutrient dispersion, soil structure, and the microclimate. Furthermore, mountain ranges on the scale of the Swiss Alps form a relevant barrier for animal migration. In this light, the observed diversity of the B.11 clade in Switzerland and the statistically significant correlation between genetic- and geographic distance may reflect small-scale topographic fragmentation of the habitat of the susceptible hosts and pathogen vectors.

Analysis of a Newly Identified B.11 Subclade

In agreement with the European survey, the B.45 subclade of B.11 was found to be the most successful subclade in Switzerland in terms of geographical distribution and prevalence. Noticeably, all of the 14 tick isolates from five different geographic regions were assigned to the B.45 clade, whereas B.46 and B.86 were solely comprised by human and animal isolates. One possible explanation for this observation could be that the strains of the B.45 subclade are better adapted to the arthropod vector and thus elevate the vector competence of the ticks, compared to ticks carrying B.46/B.86 strains. A hint for this line of argumentation is provided by our breseq data. Seven non-synonymous SNPs were found exclusively in the B.86 subclade of B.11, which might modulate the function of proteins involved in intracellular survival. Members of the MFS are involved in the transport of small molecule substrates across all classes of organisms and play an important role in sustaining the osmotic balance between the cytosol and the environment. In F. tularensis, phagosomal transport proteins that are a sub-family of MFS were identified as factors essential for lethality to adult fruit flies (Akimana and Kwaik, 2011). Furthermore, members of the MFS class are discussed as a target for attenuation and vaccine development (Marohn et al., 2012). The tolerance to oxidative stress and the acquisition of iron are fundamental aspects for the survival of a pathogen in a host environment. Besides its pivotal role in cellular respiration, the iron containing porphyrin-ring heme is involved in the function of a variety of enzymes like catalases and nitric oxide synthase (Choby and Skaar, 2016). In this light, the non-synonymous mutation of hemA that is involved in the first steps of heme biosynthesis may affect the tolerance of clade B.86 strains to oxidative stress (Ezraty et al., 2017). Additionally, we identified a mutation in another component associated with the electron transport chain: cytochrome c-type biogenesis proteins are membrane-bound proteins that may play a role in the guidance of apocytochromes and heme groups for their covalent linkage by the cytochrome-c-heme lyase. In summary, the mutations of proteins involved in the adaption to environmental conditions may modify the fitness of the B.86 subclade and may prevent persistence in ticks. This would lead to a diminished host range and thus restricted geographical dispersal. However, this interpretation should be treated with caution and needs to be substantiated further by functional studies.

Investigation of Eco-Epidemiological Aspects

According to the reported tularemia cases, ticks are the predominant vector for disease transmission in Switzerland over the last 10 years. The importance of arthropod vectors is also reflected in the clinical manifestation of the disease where the glandular / ulcero glandular form prevails. These findings together with the sporadic occurrence and the low number of cases are in agreement with a terrestrial life cycle of Fth involving rodents, lagomorphs and ticks as main source for human infections. Despite the lacking evidence of a transovarial transmission of Fth in Ixodes ricinus (Genchi et al., 2015), it is known that the pathogen is propagated transracial from larvae to nymphs and adults in the 3-year life cycle of the vector. Therefore it was suggested that ticks can be regarded as a de facto reservoir of tularemia (Petersen et al., 2009; Foley and Nieto, 2010).

Our spatial data from tick, human, and vertebrate isolates, as well as from reported clinical cases suggest that the Fth population is not randomly distributed in Switzerland. The highest risk for acquiring tularemia coincides with regions where Fth was found in the ticks, a hint to the importance of their role as a biological niche for Fth. Since the initial tick survey was conducted with the aim to assess the prevalence of TBEV in all regions of Switzerland below 1,500 m altitude (Gäumann et al., 2010), the identification of Fth foci within the TBEV sample pool was not biased by the number of reported tularemia cases in a given region.

Analysis of Habitat Characteristics and Abiotic Factors Influencing the Tick Abundance

An important factor influencing the distribution of the pathogen in Switzerland is the complex topography of the country where 50% of the territory lays above an altitude of 1,000 m above sea level. Only one case of tularemia was reported in a patient living at 1,700 m altitude and contaminated water or food was suspected to be the source of infection (Ernst et al., 2015). In principle, the incidence of tularemia shows the same altitude dependence as was shown for borreliosis, where the prevalence of the pathogen in tick nymphs correlated negatively with altitude (Gern et al., 2008). There are multitudes of additional environmental parameters that are known to have an influence on the transmission dynamics of tick-borne pathogens. Abiotic factors like temperature, wind, RH and SD are recognized to be of pivotal importance for tick abundance and survival (Vial, 2009; Alonso-Carné et al., 2015). It was shown that a RH of 85% is required for the survival of Ixodes ricinus in its non-parasitic stage (Keymer et al., 2009). Studies, which assessed the usability of satellite-based remote sensing data as indices for geospatial disease monitoring, identified soil moisture as an important factor for the establishment of tick habitats (Beck et al., 2000; Barrios et al., 2012).

Comparing surface normalized averages of waterlogging levels in areas with elevated risk for acquiring tularemia vs. low-/expected risk areas, we find a highly significant effect of soil moisture (Figures 10, 11). Thus, the elevated prevalence of Fth in ticks and the higher incidence of tularemia in the northeastern part of Switzerland may at least be partially explained by the above average soil moisture in this region, favoring tick survival during dry periods. Another observation that fits in this line of argumentation is the congruence of high-risk tularemia regions with areas of elevated risk of Tick Born Encephalitis (TBE) (Figure 13).

FIGURE 13

Figure 13. Depiction of localized clusters of Tick borne encephalitis (TBE). It is based on information from the mandatory reporting system between 2007 and 2016.

Habitat Fragmentation and Urbanization in Switzerland

Besides micro-climatic factors, anthropogenic changes of the ecosystem, and social factors are known to have a fundamental impact on the host-vector-pathogen dynamics (Bayles et al., 2013; Lou et al., 2014). As seen in other European countries, the fragmentation and urbanization of the Swiss landscape is on the rise which is reflected in a linear decrease of the “effective mesh size” during the last 70 years. The mesh size is a metric that describes the probability that two random points in the landscape can be connected without the interference of artificial structures, for example, transportation routes, buildings or developed land (Jaeger, 2000, 2006). The more barriers fragmenting the landscape, the lower the probability that two points are connected, and the lower the effective mesh size. In an ecological view, the mesh size can be interpreted as the probability that two animals of the same specie find each other in the landscape. Therefore, the reduction of the mesh size due to progressing landscape fragmentation may reach a limit where a given species is no longer able to sustain a stable population, which reduces the biodiversity of the area (Di Giulio et al., 2009). Most of the regions with an elevated risk for acquiring tularemia are highly fragmented and show a below average mesh sizes (Figure 12). The reduction of biodiversity may alter the ratio of hosts vs. low or non-competent hosts, thereby reducing the pathogen dilution effect. A diminished pathogen dilution effect may in turn elevate the vector competence of the ticks and thus the probability of disease transmission (Schmidt and Ostfeld, 2001; Pfäffle et al., 2013; Dantas-Torres, 2015).

Risk Factors for Human Disease

The bimodal distribution of the tularemia cases per age group and the predominance of male patients is an indicator of social factors influencing the transmission probability of the disease and is in accordance with the incidence in other countries (Desvars et al., 2015). According to clinicians treating patients in the elevated risk areas, the peak in the age group of the 8–12 year old children is due to the joining of youth organizations like the scouts that are combing through the forests away from fixed paths seeking for adventure. The second peak with a maximum at 50 years reflects the demographic age structure of the Swiss population.

Conclusion

The combination of high-resolution whole genome phylogenies with epidemiological and ecological parameters allows a profound insight into the transmission modalities of tularemia. The exceptionally high genetic diversity of the Swiss Fth population allows us to study the population and disease dynamics of the pathogen in a small geographical context. In our work, we found evidence for the spatial correlation of disease incidence and ecological/anthropogenic factors that influence the survival and vector competence of the local tick population.

Based on these findings, we conclude that I. ricinus plays a pivotal role in the establishment of the disease in a given ecological niche and in sustaining the transmission cycle.

Author Contributions

MW analyzed the genomic data, performed the spatial statistics, and drafted the manuscript. EA compiled the epidemiological data and calculated the spatial relative disease risk. PP provided the animal isolates and expertise in the field of tularemia. SG performed the functional SNP analysis and provided the expertise in whole genome phylogeny. CB, RA-G designed and organized the tick sampling and provided expertise in tick borne diseases. UK provided the expertise in the clinical aspects of tularemia. DJ, RG gave theoretical and practical advice for bacterial isolation and provided expertise in the eco-epidemiology of tularemia. NS coordinated and supervised the presented work. All authors contributed to manuscript revision, read, and approved the submitted version.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We like to thank Kerstin Myrtennäs from the Swedish Defense Research Agency (FOI) for fruitful discussions regarding the phylogeny of Fth and canSNP nomenclature. Furthermore, we thank the NBC Defense Laboratory 1 of the Swiss Armed Forces for the extensive tick sampling campaigns, ADMED (Analyses et Diagnostics Medicaux) for supplying DNA and last but certainly not least Sandra Paniga, Susanne Thomann, and Fritz Wüthrich for excellent technical assistance. This work was partially supported by the Health program 2014–2020, through the Consumers, Health, Agriculture and Food Executive Agency (CHAFEA, European Commission); EMERGE Joint Action grant number:677066.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcimb.2018.00089/full#supplementary-material

Supplementary Table 1. Genetic variations found in the core genome alignment with the strain FTNF00–002.

References

Ajawatanawong, P., Atkinson, G. C., Watson-Haigh, N. S., Mackenzie, B., and Baldauf, S. L. (2012). SeqFIRE: a web application for automated extraction of indel regions and conserved blocks from protein multiple sequence alignments. Nucleic Acids Res. 40, W340–W347. doi: 10.1093/nar/gks561

CrossRef Full Text | Google Scholar

Akimana, C., and Kwaik, Y. A. (2011). Francisella-arthropod vector interaction and its role in patho-adaptation to infect mammals. Front. Microbiol. 2:34. doi: 10.3389/fmicb.2011.00034

PubMed Abstract | CrossRef Full Text | Google Scholar

Alonso-Carné, J., García-Martín, A., and Estrada-Peña, A. (2015). Assessing the statistical relationships among water-derived climate variables, rainfall, and remotely sensed features of vegetation: implications for evaluating the habitat of ticks. Exp. Appl. Acarol. 65, 107–124. doi: 10.1007/s10493-014-9849-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Bankevich, A., Nurk, S., Antipov, D., Gurevich, A. A., Dvorkin, M., Kulikov, A. S., et al. (2012). SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J. Comput. Biol. 19, 455–477. doi: 10.1089/cmb.2012.0021

PubMed Abstract | CrossRef Full Text | Google Scholar

Barrios, J. M., Verstraeten, W. W., Maes, P., Clement, J., Aerts, J. M., Farifteh, J., et al. (2012). Remotely sensed vegetation moisture as explanatory variable of Lyme borreliosis incidence. Int. J. Appl. Earth Obs. Geoinf. 18, 1–12. doi: 10.1016/j.jag.2012.01.023

CrossRef Full Text | Google Scholar

Bayles, B. R., Evans, G., and Allan, B. F. (2013). Knowledge and prevention of tick-borne diseases vary across an urban-to-rural human land-use gradient. Ticks Tick Borne Dis. 4, 352–358. doi: 10.1016/j.ttbdis.2013.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Beck, L. R., Lobitz, B. M., and Wood, B. L. (2000). Remote sensing and human health: new sensors and new opportunities. Emerg. Infect. Dis. 6, 217–227. doi: 10.3201/eid0603.000301

PubMed Abstract | CrossRef Full Text | Google Scholar

Becker, S., Lochau, P., Jacob, D., Heuner, K., and Grunow, R. (2016). Successful re-evaluation of broth medium T for growth of Francisella tularensis ssp. and other highly pathogenic bacteria. J. Microbiol. Methods 121, 5–7. doi: 10.1016/j.mimet.2015.11.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30:2114. doi: 10.1093/bioinformatics/btu170

PubMed Abstract | CrossRef Full Text | Google Scholar

Byström, M., Böcher, S., Magnusson, A., Prag, J., and Johansson, A. (2005). Tularemia in denmark: identification of a Francisella tularensis subsp. holarctica strain by real-time PCR and high-resolution typing by multiple-locus variable-number tandem repeat analysis. J. Clin. Microbiol. 43, 5355–5358. doi: 10.1128/JCM.43.10.5355-5358.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Castresana, J. (2000). Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol. Biol. Evol. 17, 540–552. doi: 10.1093/oxfordjournals.molbev.a026334

PubMed Abstract | CrossRef Full Text | Google Scholar

Choby, J. E., and Skaar, E. P. (2016). Heme synthesis and acquisition in bacterial pathogens. J. Mol. Biol. 428, 3408–3428. doi: 10.1016/j.jmb.2016.03.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Dantas-Torres, F. (2015). Climate change, biodiversity, ticks and tick-borne diseases: the butterfly effect. Int. J. Parasitol. Parasites Wildl. 4, 452–461. doi: 10.1016/j.ijppaw.2015.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Deatherage, D. E., and Barrick, J. E. (2014). Identification of mutations in laboratory-evolved microbes from next-generation sequencing data using breseq. Methods Mol. Biol. 1151, 165–188. doi: 10.1007/978-1-4939-0554-6_12

PubMed Abstract | CrossRef Full Text | Google Scholar

Desvars, A., Furberg, M., Hjertqvist, M., Vidman, L., Sjöstedt, A., Rydén, P., et al. (2015). Epidemiology and ecology of Tularemia in Sweden, 1984-2012. Emerg. Infect. Dis. 21, 32–39. doi: 10.3201/eid2101.140916

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Giulio, M., Holderegger, R., and Tobias, S. (2009). Effects of habitat and landscape fragmentation on humans and biodiversity in densely populated landscapes. J. Environ. Manage. 90, 2959–2968. doi: 10.1016/j.jenvman.2009.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Dwibedi, C., Birdsell, D., Lärkeryd, A., Myrtennäs, K., Öhrman, C., Nilsson, E., et al. (2016). Long-range dispersal moved Francisella tularensis into Western Europe from the East. Microb. Genomics 2:e000100. doi: 10.1099/mgen.0.000100

PubMed Abstract | CrossRef Full Text | Google Scholar

Ernst, M., Pilo, P., Fleisch, F., and Glisenti, P. (2015). Tularemia in the Southeastern Swiss Alps at 1,700 m above sea level. Infection 43, 111–115. doi: 10.1007/s15010-014-0676-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Ezraty, B., Gennaris, A., Barras, F., and Collet, J.-F. (2017). Oxidative stress, protein damage and repair in bacteria. Nat. Rev. Microbiol. 15, 385–396. doi: 10.1038/nrmicro.2017.26

PubMed Abstract | CrossRef Full Text | Google Scholar

Foley, J. E., and Nieto, N. C. (2010). Tularemia. Vet. Microbiol. 140, 332–338. doi: 10.1016/j.vetmic.2009.07.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Gäumann, R., Mühlemann, K., Strasser, M., and Beuret, C. M. (2010). High-throughput procedure for tick surveys of tick-borne encephalitis virus and its application in a national surveillance study in Switzerland. Appl. Environ. Microbiol. 76, 4241–4249. doi: 10.1128/AEM.00391-10

PubMed Abstract | CrossRef Full Text | Google Scholar

Genchi, M., Prati, P., Vicari, N., Manfredini, A., Sacchi, L., Clementi, E., et al. (2015). Francisella tularensis: no evidence for transovarial transmission in the tularemia tick vectors Dermacentor reticulatus and Ixodes ricinus. PLoS ONE 10:e0133593. doi: 10.1371/journal.pone.0133593

PubMed Abstract | CrossRef Full Text | Google Scholar

Gern, L., Morán Cadenas, F., and Burri, C. (2008). Influence of some climatic factors on Ixodes ricinus ticks studied along altitudinal gradients in two geographic regions in Switzerland. Int. J. Med. Microbiol. 298, 55–59. doi: 10.1016/j.ijmm.2008.01.005

CrossRef Full Text | Google Scholar

Girault, G., Thierry, S., Cherchame, E., and Derzelle, S. (2014). Application of high-throughput sequencing : discovery of informative SNPs to Subtype Bacillus anthracis. Adv. Biosci. Biotechnol. 5, 669–677. doi: 10.4236/abb.2014.57079

CrossRef Full Text | Google Scholar

Guindon, S., Dufayard, J.-F., Lefort, V., Anisimova, M., Hordijk, W., and Gascuel, O. (2010). New algorithms and mehtods to estimate maximum-likelihood phylogenies: asessing the performance of PhyML 2.0. Syst. Biol. 59, 307–321. doi: 10.1093/sysbio/syq010

CrossRef Full Text | Google Scholar

Hornik, K., Leisch, F., Zeileis, A., Gómez-Rubio, V., Ferrándiz, J., and López, A. (2003). Detecting Clusters of Diseases with R. Available online at: http://www.ci.tuwien.ac.at/Conferences/DSC-2003/ (Accessed June 19, 2017).

Jaeger, J. (2000). Landscape division, splitting index, and effective mesh size: new measures of landscape fragmentation. Landsc. Ecol. 15, 115–130. doi: 10.1023/A:1008129329289

CrossRef Full Text | Google Scholar

Jaeger, J. (2006). Degree of Landscape Fragmentation in Switzerland. Federal Office for the Environment FOEN. DIV-8409–E

Jusatz, H. J. (1961). “The geographical distribution of tularemia throughout the world, 1911–1959,” in Welt-Suchen Atlas, ed E. Rodenwaldt (Hamburg: Falk-Verlag), 7–12.

Karlsson, E., Golovliov, I., Lärkeryd, A., Granberg, M., Larsson, E., Öhrman, C., et al. (2016). Clonality of erythromycin resistance in Francisella tularensis. J. Antimicrob. Chemother. 71, 2815–2823. doi: 10.1093/jac/dkw235

PubMed Abstract | CrossRef Full Text | Google Scholar

Katoh, K., Misawa, K., Kuma, K.-I., and Miyata, T. (2002). MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 30, 3059–3066. doi: 10.1093/nar/gkf436

PubMed Abstract | CrossRef Full Text | Google Scholar

Keirans, J. E. (1997). Ticks of north-west Europe. By Paul D. Hillyard. Med. Vet. Entomol. 11, 102–103. doi: 10.1111/j.1365-2915.1997.tb00296.x

CrossRef Full Text | Google Scholar

Keymer, D. P., Lam, L. H., and Boehm, A. B. (2009). Biogeographic patterns in genomic diversity among a large collection of Vibrio cholerae isolates. Appl. Environ. Microbiol. 75, 1658–1666. doi: 10.1128/AEM.01304-08

PubMed Abstract | CrossRef Full Text | Google Scholar

Krige, D. (1951). A statistical approach to some basic mine valuation problems on the Witwatersrand. J. Chem. Metall. Min. Soc. South Afr. 52, 201–203.

Google Scholar

Lärkeryd, A., Myrtennäs, K., Karlsson, E., Dwibedi, C. K., Forsman, M., Larsson, P. R., et al. (2014). CanSNPer: a hierarchical genotype classifier of clonal pathogens. Bioinformatics 30, 1762–1764. doi: 10.1093/bioinformatics/btu113

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, R. S., Radomski, N., Proulx, J.-F., Levade, I., Shapiro, B. J., McIntosh, F., et al. (2015). Population genomics of Mycobacterium tuberculosis in the Inuit. Proc. Natl. Acad. Sci. U.S.A. 112, 13609–13614. doi: 10.1073/pnas.1507071112

PubMed Abstract | CrossRef Full Text | Google Scholar

Lou, Y., Wu, J., and Wu, X. (2014). Impact of biodiversity and seasonality on Lyme-pathogen transmission. Theor. Biol. Med. Model. 11:50. doi: 10.1186/1742-4682-11-50

PubMed Abstract | CrossRef Full Text | Google Scholar

Marohn, M. E., Santiago, A. E., Shirey, K. A., Lipsky, M., Vogel, S. N., and Barry, E. M. (2012). Members of the Francisella tularensis phagosomal transporter subfamily of major facilitator superfamily transporters are critical for pathogenesis. Infect. Immun. 80, 2390–2401. doi: 10.1128/IAI.00144-12

PubMed Abstract | CrossRef Full Text | Google Scholar

Maurin, M., and Gyuranecz, M. (2016). Tularaemia: clinical aspects in Europe. Lancet Infect. Dis. 16, 113–124. doi: 10.1016/S1473-3099(15)00355-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Müller, W., Hotzel, H., Otto, P., Karger, A., Bettin, B., Bocklisch, H., et al. (2013). German Francisella tularensis isolates from European brown hares (Lepus europaeus) reveal genetic and phenotypic diversity. BMC Microbiol. 13:61. doi: 10.1186/1471-2180-13-61

PubMed Abstract | CrossRef Full Text | Google Scholar

Oechslin, C. P., Heutschi, D., Lenz, N., Tischhauser, W., Péter, O., Rais, O., et al. (2017). Prevalence of tick-borne pathogens in questing Ixodes ricinus ticks in urban and suburban areas of Switzerland. Parasit. Vectors 10:558. doi: 10.1186/s13071-017-2500-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Origgi, F. C., Frey, J., and Pilo, P. (2014). Characterisation of a new group of Francisella tularensis subsp. holarctica in Switzerland with altered antimicrobial susceptibilities, 1996 to 2013. Eurosurveillance 19:20858. doi: 10.2807/1560-7917.ES2014.19.29.20858

PubMed Abstract | CrossRef Full Text | Google Scholar

Page, A. J., Cummins, C. A., Hunt, M., Wong, V. K., Reuter, S., Holden, M. T. G., et al. (2015). Roary: rapid large-scale prokaryote pan genome analysis. Bioinformatics 30, 3691–3693. doi: 10.1093/bioinformatics/btv421

CrossRef Full Text | Google Scholar

Pérez-Castrillón, J. L., Bachiller-Luque, P., Martín-Luquero, M., Mena-Martín, F. J., and Herreros, V. (2001). Tularemia epidemic in Northwestern Spain: clinical description and therapeutic response. Clin. Infect. Dis. 33, 573–576. doi: 10.1086/322601

PubMed Abstract | CrossRef Full Text | Google Scholar

Petersen, J. M., Mead, P. S., and Schriefer, M. E. (2009). Francisella tularensis: an arthropod-borne pathogen. Vet. Res. 40:7. doi: 10.1051/vetres:2008045

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfäffle, M., Littwin, N., Muders, S. V., and Petney, T. N. (2013). The ecology of tick-borne diseases. Int. J. Parasitol. 43, 1059–1077. doi: 10.1016/j.ijpara.2013.06.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Pilo, P., Johansson, A., and Frey, J. (2009). Identification of Francisella tularensis cluster in central and western Europe. Emerg. Infect. Dis. 15, 2049–2051. doi: 10.3201/eid1512.080805

PubMed Abstract | CrossRef Full Text | Google Scholar

Rotz, L. D., Khan, A. S., Lillibridge, S. R., Ostroff, S. M., and Hughes, J. M. (2002). Public health assessment of potential biological terrorism agents. Emerg. Infect. Dis. 8, 225–230. doi: 10.3201/eid0802.010164

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmidt, K. A., and Ostfeld, R. S. (2001). Biodiversity and the dilution effect in disease ecology. Ecology 82, 609–619. doi: 10.1890/0012-9658(2001)082[0609:BATDEI]2.0.CO;2

CrossRef Full Text | Google Scholar

Seemann, T. (2014). Prokka: rapid prokaryotic genome Annotation 30, 2068–2069. doi: 10.1093/bioinformatics/btu153

PubMed Abstract | CrossRef Full Text | Google Scholar

Stucki, D., Brites, D., Jeljeli, L., Coscolla, M., Liu, Q., Trauner, A., et al. (2016). Mycobacterium tuberculosis lineage 4 comprises globally distributed and geographically restricted sublineages. Nat. Genet. 48, 1535–1543. doi: 10.1038/ng.3704

PubMed Abstract | CrossRef Full Text | Google Scholar

Tarroso, P., Velo-Anton, G., and Carvalho, S. B. (2015). Phylin: Spatial Interpolation of Genetic Data. Available online at: https://cran.r-project.org/package=phylin

Tomaso, H., Al Dahouk, S., Hofer, E., Splettstoesser, W. D., Treu, T. M., Dierich, M. P., et al. (2005). Antimicrobial susceptibilities of Austrian Francisella tularensis holarctica biovar II strains. Int. J. Antimicrob. Agents 26, 279–284. doi: 10.1016/j.ijantimicag.2005.07.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Vial, L. (2009). Mise au point biological and ecological characteristics of soft ticks (ixodida: argasidae) and their impact for predicting tick and associated disease distribution. Parasite 16, 191–202. doi: 10.1051/parasite/2009163191

CrossRef Full Text | Google Scholar

Vogler, A. J., Birdsell, D., Price, L. B., Bowers, J. R., Beckstrom-Sternberg, S. M., Auerbach, R. K., et al. (2009). Phylogeography of Francisella tularensis: global expansion of a highly fit clone. J. Bacteriol. 191, 2474–2484. doi: 10.1128/JB.01786-08

PubMed Abstract | CrossRef Full Text | Google Scholar

Wicki, R., Sauter, P., Mettler, C., Natsch, A., Enzler, T., Pusterla, N., et al. (2000). Swiss Army Survey in Switzerland to determine the prevalence of Francisella tularensis, members of the Ehrlichia phagocytophila genogroup, Borrelia burgdorferi sensu lato, and tick-borne encephalitis virus in ticks. Eur. J. Clin. Microbiol. Infect. Dis. 19, 427–432. doi: 10.1007/s100960000283

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, T., Chen, B.-M., Liu, G., Huang, F.-F., Liu, J.-G., Liao, W.-B., et al. (2015). Biodiversity of Jinggangshan Mountain: the importance of topography and geographical location in supporting higher biodiversity. PLoS ONE 10:e0120208. doi: 10.1371/journal.pone.0120208

PubMed Abstract | CrossRef Full Text | Google Scholar