Molecular Investigation of Tularemia Outbreaks, Spain, 1997–2008

This disease has reemerged because of persistence of local reservoirs of infection.

studied F. tularensis isolates circulating in Spain during outbreaks with different epidemiologic patterns and investigated whether reemergence of the pathogen after 10 years of no epidemiologic activity was caused by introduction of exotic strains or by establishment of the pathogen in local reservoirs of infection.
All isolates were grown on modified Thayer-Martin agar plates containing 36 g/L GC agar base, 10 g/L soluble hemoglobin powder, and 2 vials/L Vitox supplement (Oxoid, Basingstoke, UK) at 37°C for 2-3 days in aerobic conditions. Biochemical characterization included tests for oxidase and catalase activities, glucose and glycerol fermentation, and urea hydrolysis.

Genetic Characterization
Species and subspecies were identified by real-time and conventional PCRs specific for the fopA gene and the region of difference 1 (RD1) as described (15,16). RD23 was analyzed by PCR for identification of a genetic group of isolates that had been found on the Iberian Peninsula (17).
Pulsed-field gel electrophoresis (PFGE) and multilocus variable number tandem repeat analysis (MLVA) were used to classify isolates into genetic subpopulations. The PFGE protocol described (18), which used restriction enzymes XhoI and BamHI, was optimized to provide major improvements in quality of fingerprint patterns.
Bacterial cells were suspended in SE buffer (25 mmol/L EDTA, 75 mmol/L NaCl, pH 7.5) to an absorbance of 0.5-0.6 at 600 nm. Cells were lysed in agarose plugs and plugs were washed 5 times with Tris-EDTA buffer (10 mmol/L Tris-HCl, 1 mmol/L EDTA, pH 8.0) for 30 min at 50°C. DNA in the plugs was digested with 40 U of XhoI (New England Biolabs, Ipswich, MA, USA) or 40 U of BamHI (New England Biolabs), for 16 and 3 h, respectively, at 37°C following the manufacturer's protocol. DNA fragment sizes were determined by electrophoresis and by comparing bands with a Lambda Ladder PFG Marker (New England Biolabs). MLVA was performed as described for 16 variable number tandem repeat loci (5). To ensure analysis of identical genetic material by PFGE and MLVA, we used DNA from the same culture for both methods. MLVA markers were amplified by using PCR, and sizes of amplification products were determined by electrophoresis on 3.5% high-resolution agarose MS-8 gels (Conda Pronadisa, Madrid, Spain), except for Ft-M3, Ft-M21, Ft-M22, and Ft-M24 MLVA markers, for which the sizes were determined by using capillary electrophoresis. At least 2 alleles were sequenced for each MLVA marker to confirm that size differences observed resulted from the expected variations in numbers of tandem repeats. Forward and reverse sequences were aligned by using MEGA v.4 software (19), and consensus sequences were used to predict the number of tandem repeats in each allele.

Data Analyses
The Simpson index of diversity, which measures the probability that 2 unrelated strains from the test population will be classified into different typing groups (20), was calculated to compare the discriminative power of PFGE typing with that of MLVA for assessing genetic diversity among isolates. The adjusted Wallace coefficient for quantification of agreement between PFGE typing and MLVA results was also calculated. Both analyses were performed by using Comparing Partitions (21). PFGE patterns were analyzed by using Bionumerics v.6.6 (Applied-Maths NV, Sint-Martens-Latem, Belgium) to describe genetic relationships among isolates. Dendrograms were constructed by using the Dice similarity coefficient and the unweighted pair group mathematical average clustering algorithm. MLVA data, expressed as allelic profiles for isolates, were analyzed by using Bionumerics v.6.6. Minimum spanning trees were calculated with priority rules set at first link allelic profiles and maximum numbers of single-locus variants and then maximal numbers of single-locus variants and doublelocus variants. MLVA types were classified as members of a clonal complex if they had the same allele at15 of the 16 MLVA markers. A map of the distribution of isolates showing the geographic origin and number of isolates per province was generated by using Arcgis v.9.2 software (ESRI, Redlands, CA, USA).

Subspecies and Genetic Subclade of F. tularensis Isolates
All isolates were negative for oxidase activity, weakly positive for catalase activity, and positive for acid production from glucose; none of the isolates hydrolyzed urea. Only the reference strain F. tularensis subsp. tularensis Schu (CAPM 5600) produced acid from glycerol. Realtime PCR specific for the fopA gene and size determination at the RD1 region showed that all isolates from Spain and the Czech Republic were F. tularensis subsp. holarctica (online Technical Appendix). All isolates from Spain included in this study had the 1.59-kbp deletion at the RD23 loci, which is characteristic of the F. tularensis subsp. holarctica genetic subclade B.Br:FTNF002-00 (also known as the Iberian clone or the central and western European genetic group).

Characterization by PFGE
All 107 F. tularensis subsp. holarctica isolates showed the same fingerprint pattern by PFGE with the restriction enzyme XhoI, irrespective of their geographic origin, host, or date of isolation (isolate TU41 was not typeable by PFGE analysis). This pattern consisted of ≈20 DNA fragments >70 kbp. The F. tularensis subsp. tularensis strain Schu (CAPM 5600) showed a different banding pattern.
Sixty-nine (63.9%) isolates from Spain clustered into pulsotype A and 26 (24.1%) other isolates from Spain clustered into pulsotype B. All isolates from the Czech Republic, except for isolate CAPM 5538, had the same fingerprint pattern, which was designated pulsotype D (8.3%). Isolate CAPM 5538 showed a pulsotype that clustered with the 2 remaining isolates from Spain (TU8 and TU9) in pulsotype C (1.9%). One isolate from Spain, TU41, could not be genotyped by PFGE despite several attempts (online Technical Appendix).
There were some discrepancies between our findings and those reported for the 37 isolates in Spain from 1997 (18). Improvements in quality of fingerprint patterns enabled us to distinguish between isolates from Spain and those from the Czech Republic by using PFGE and restriction enzyme BamHI. Furthermore, isolates TU3, TU17, TU21, and TU25 were unequivocally assigned to pulsotype B instead of pulsotype A. In instances of discrepancy, analyses were repeated in triplicate with new cultures, and the findings reported were confirmed. Distribution of the 107 F. tularensis subsp. holarctica isolates into 4 pulsotypes resulted in a Simpson index of diversity of 0.522. There was no obvious correlation between the pulsotype of an isolate from Spain and its geographic origin, host, or tularemia outbreak with which it was associated.

Characterization by MLVA
The allele-based analysis of genetic relationships identified 13 MLVA types among the 108 F. tularensis subsp. holarctica isolates and showed that F. tularensis subsp. tularensis strain Schu (CAPM 5600) was more distantly related ( Figure 2). The 10 isolates from the Czech Republic were assigned to 5 MLVA types, which differed from isolates from Spain by ≥2 alleles. Marker Ft-M3 provided the highest number of alleles (6). Six copy numbers were detected among the 109 isolates; for Ft-M6, Ft-M9, and Ft-M20, there were 3 copy numbers. Markers Ft-M5, Ft-M7, Ft-M8, Ft-M10, Ft-M13, FT-M16, Ft-M19, Ft-M21, Ft-M22, Ft-M23, and Ft-M24 each had 2 alleles: these markers with 2 alleles, except for Ft-M24, discriminated only F. tularensis subsp. tularensis strain Schu (CAPM 5600) from all isolates of F. tularensis subsp. holarctica (Table). Ft-M24 had a 464-bp allele that was found in all isolates from Spain analyzed in this study, but was not present in any other isolates. Ft-M24 has been found only in isolates of genetic subclade B.Br:FTNF002-00 (the Iberian clone or the central and western European genetic group). Sequence analysis of the Ft-M24 DNA fragment showed that the unique allele size was caused by deletion of a 16bp sequence adjacent to 2 copies of the Ft-M24 tandem repeat (GenBank accession no. KC696513). For marker Ft-M12, because all 109 isolates had the same copy number, this marker provided no typing resolution. Distribution of F. tularensis subsp. holarctica isolates among 12 MLVA types was uneven; >70% of the isolates in the MLVA types A (49 isolates, 45%) and B (32 isolates, 29.4%) ( Table). Simpson index of diversity, which showed the discriminatory power of MLVA for the 108 F. tularensis subsp. holarctica isolates, was 0.708.
The 98 isolates from Spain were classified into 8 MLVA types, which essentially grouped as 2 closely related clonal complexes that differed at only 1 of the 16 MLVA markers (Figure 2). All MLVA types from Spain were single-locus variants of MLVA type B, which indicated this type was the founder genotype of F. tularensis that caused tularemia in northwestern Spain. No clear relationship was found between genotype and geographic origin (Figure 3), source of infection, or host in Spain. The same genotype was usually isolated from hares, voles, and humans in Spain.
Comparison of isolates from the 2 outbreak periods (37 isolates for 1997-1998 and 61 isolates for 2007-2008) showed that the same F. tularensis genotypes caused tularemia in both outbreaks (Figure 4). Isolates from the second outbreak showed less genetic diversity than those from the first outbreak (Simpson indices 0.62, 95% CI 0.53-0.71 and 0.66, 95% CI 0.57-0.75, respectively; the difference was not significant at the 95% level).
Comparison of allele distribution at the most variable marker (Ft-M3) showed an overall similarity between isolates causing the outbreaks, although the most common copy number was 4 during the first outbreak and 5 during the second outbreak, which might indicate a stepwise increase in copy number over time. Overall, our findings for 37 isolates from the first outbreak were consistent with the data reported by Dempsey et al. (17) although there were 2 discrepancies. First, isolate TU18 had a unique allele with 3 tandem repeats at Ft-M9, which distinguished this isolate from all other F. tularensis subsp. holarctica isolates. Second, isolate TU31 had the same FT-M10 allele as all other isolates. We confirmed our results for these discrepancies in triplicate.

Quantification of Agreement between PFGE Typing and MLVA
Congruence of the 2 methods (PFGE typing and MLVA) for isolate classification was weak for the 107 F. tularensis subsp. holarctica isolates (1 of the 108 isolates was excluded because it was not typeable by PFGE analysis). This finding was true for reverse comparisons of both methods: if 2 isolates were in the same PFGE pulsotype; they had an 18% chance of having the same MLVA type. Conversely, having the same MLVA type was associated with a 40% chance of having the same PFGE pulsotype. The adjusted Wallace coefficient for PFGE versus MLVA was 0.18 (95% CI 0.04-0.32) and that for MLVA versus PFGE was 0.40 (95% CI 0.20-0.58).

Discussion
F. tularensis subsp. holarctica has shown limited genetic diversity worldwide (5,7,22,23). This finding might be the result of a relatively recent bottleneck or clonal expansion event that drastically reduced genetic variation of the bacterial population ( findings, we observed extremely limited genomic diversity among the 98 F. tularensis subsp. holarctica isolates from Spain analyzed by 2 genotyping tools: PFGE with 2 restriction enzymes and MLVA at 16 highly variable tandem repeat loci. PFGE identified 3 genotypes with single band differences (93.3% similarity) (Figure 1). MLVA discriminated these isolates into 8 MLVA types, but in pairwise comparisons they differed at no more than 2 of 16 MLVA markers. Thus, our collection shows extreme genetic homogeneity (Figure 2). An F. tularensis subsp. holarctica lineage in central and western Europe (France, Spain, and Switzerland) has been defined (7,17,24,25). Strains belonging to this lineage have 2 unique genetic traits: a 1.59-kbp genomic deletion at the RD23 locus and a unique 464-bp allele at Ft-M24. All isolates from Spain had these alleles, irrespective of the outbreak, geographic origin, or the host from which they were recovered. Thus, all isolates from Spain analyzed belong to the genetic subclade B.Br:FTNF002-00 (the Iberian clone or central and western European genetic group) (7). Furthermore, all MLVA types for isolates from Spain were single-locus variants of MLVA type B, which suggested that this type might be a founding genotype that has evolved into multiple other genotypes that differ from the founding genotype at a single loci. In this scenario, all strains causing tularemia outbreaks in Spain are linked to this founder (ancestral) genotype.
We found poor congruence between typing results of PFGE and MLVA for 107 F. tularensis subsp. holarctica isolates. However, the 2 methods might indicate different types of genetic variation. PFGE is a suitable approach for detecting rearrangements in a genome, and differences of only 1 band observed among the 98 isolates from Spain are presumably consequences of a single mutation event that might be an inversion, translocation, deletion, or a singlenucleotide polymorphism (26). In contrast, MLVA detects variation in several, fast-evolving, repeated sequences. However, such rapidly evolving sequences are susceptible to homoplasy, and genetic classification of isolates on the basis of a difference at only 1 of 16 genetic loci, as for isolates from Spain, could be biased because of genetic reversion events. Because some of the mutations detected by PFGE or MLVA might not be selectively neutral, we might have observed time-dependent mutations that are transient on evolutionary time scales and will frequently be eliminated by selection pressure acting on them (27). If this hypothesis is true, use of single mutations for genetic discrimination may lead to incorrect phylogenetic inferences. Therefore, poor congruence between PFGE and MLVA for identifying genetic subclade B.Br:FTNF002-00 of F. tularensis subsp. holarctica in Spain might be caused by limited bacterial diversity. Use of more extensive genetic analyses for typing, such as whole genome sequencing, might be useful in subsequent molecular epidemiology studies.
In Spain, tularemia was first reported in late 1997 in association with one of the largest human outbreaks ever described (10). The most common route of infection of humans was by direct contact when hunting and handling hares (L. europaeus). Consistent with infection with F. tularensis through the skin, the most frequent clinical form was ulceroglandular tularemia (55.4%); glandular (15.3%) and typhoid (6.6%) forms of the disease were also observed. A second major human outbreak occurred in the same geographic area in northwest Spain in 2007 and 2008 after 10 years of no epidemiologic activity. The epidemiology of the second outbreak was different from that of the first outbreak. The second outbreak occurred when the population of the common vole (Microtus arvalis) peaked, and >65% of case-patients had typhoidal and pneumonic forms of tularemia (12,13), which is consistent with infection by inhalation.
Few outbreaks of tularemia caused by airborne transmission of the bacteria have been reported. These outbreaks include a notable outbreak of inhalational tularemia in Sweden in 1966-1967 associated with environmental exposure of farmers in which >600 cases were diagnosed (28).There were clusters of outbreaks on Martha's Vineyard (Massachusetts, USA) in 1978 and 2000 (29,30) and cases of tularemia caused by airborne transmission to 53 farmers in northern Finland during 1982. (31). In Germany in 2005, a total of 39 participants in a hare hunt were infected after exposure to contaminated droplets generated by rinsing infected hares (32).
Isolates from Spain obtained during the second tularemia outbreak had the same genotypes as those obtained during the first outbreak (Figures 1, 4). Furthermore, we did not observe any relationships between genotype and geographic origin, or host from which the isolates were recovered, which suggested that in that area, the same clones were circulating in all hosts (Figure 3, panel B). These results are useful because the 2 outbreaks had substantial epidemiologic differences (10,12). Our findings indicate that the outbreak in 2007, after 10 years of no epidemiologic activity, was not caused by introduction of a new strain, but by reemergence of an endemic bacterial population that has been circulating in the region for at least the past 15 years. Furthermore, our findings also suggest that clinical forms of the outbreak are determined by ecologic processes involved in infection (e.g., route of infection, infective dose) rather than by the genotype of the pathogen.
In conclusion, we report genetic characterization of F. tularensis subsp. holarctica isolated in Spain during 2 of the largest tularemia outbreaks worldwide. There were marked epidemiologic differences between the 2 outbreaks, which were separated by 10 years of no epidemic activity. Molecular investigations showed that both outbreaks were caused White sections in circles indicate F. tularensis subsp. holarctica isolates recovered during the first human tularemia outbreak (1997)(1998), and black sections indicate isolates recovered during the second outbreak (2007)(2008). Each circle represents a unique MLVA type and size is proportional to the number of isolates of that type.
by the same group of closely related genotypes in subclade B.Br:FTNF002-00. Therefore, the reemergence of tularemia in 2007 was presumably not caused by introduction of a new strain, but by persistence of local reservoirs of infection. These findings, along with sporadic cases of tularemia in 1998 and 2007, suggest that local foci of tularemia have become established in Spain. Further investigations will help identify these endemic foci and clarify biotic and abiotic factors that have favored establishment of the pathogen in northwestern Spain.