Clonality of erythromycin resistance in Francisella tularensis

Abstract

Objectives

We analysed diverse strains of Francisella tularensis subsp. holarctica to assess if its division into biovars I and II is associated with specific mutations previously linked to erythromycin resistance and to determine the distribution of this resistance trait across this subspecies.

Methods

Three-hundred and fourteen F. tularensis subsp. holarctica strains were tested for erythromycin susceptibility and whole-genome sequences for these strains were examined for SNPs in genes previously associated with erythromycin resistance. Each strain was assigned to a global phylogenetic framework using genome-wide canonical SNPs. The contribution of a specific SNP to erythromycin resistance was examined using allelic exchange. The geographical distribution of erythromycin-resistant F. tularensis strains was further investigated by literature search.

Results

There was a perfect correlation between biovar II strains (erythromycin resistance) and the phylogenetic group B.12. Only B.12 strains had an A → C SNP at position 2059 in the three copies of the rrl gene. Introducing 2059C into an rrl gene of an erythromycin-susceptible F. tularensis strain resulted in resistance. An additional 1144 erythromycin-resistant strains were identified from the scientific literature, all of them from Eurasia.

Conclusions

Erythromycin resistance in F. tularensis is caused by an A2059C rrl gene mutation, which exhibits a strictly clonal inheritance pattern found only in phylogenetic group B.12. This group is an extremely successful clone, representing the most common type of F. tularensis throughout Eurasia.

Introduction

The disease tularaemia is caused by the facultative intracellular bacterium Francisella tularensis and occurs naturally as seasonal outbreaks or sporadic cases among humans and other mammals. F. tularensis is classified as a biological threat agent due to its high virulence, its potential for aerosol spread and its history as a biological warfare agent.¹ The F. tularensis species is divided into subspecies by bacterial nomenclature: tularensis (type A), holarctica (type B) and mediasiatica, each with distinct biochemical attributes, geographical distribution and virulence in animal experimental models.² It was recently suggested that Francisella novicida should be included as a fourth F. tularensis subspecies due to its extensive DNA and phenotypical similarity.³ Antimicrobial treatment options for tularemia include quinolones, aminoglycosides and tetracyclines.⁴ Agents of the macrolide group of antimicrobials are generally not recommended because some strains of F. tularensis subsp. holarctica are naturally resistant to these.

Single phenotypical characteristics, including erythromycin susceptibility, have been used to divide F. tularensis strains into biovars. In the 1970–80s F. tularensis subsp. holarctica was subdivided into biovar I, biovar II and biovar japonica; biovar I is susceptible to erythromycin and other antibiotics of the macrolide group but strains of biovar II are resistant to these drugs.^5–8 Biovar japonica strains are susceptible to erythromycin but distinguished from the two other biovars by their ability to ferment glycerol. Characteristics such as antibiotic resistance and other single phenotypes, however, often are not ideal for defining population division within prokaryotes because horizontal transfer of genes (e.g. resistance genes) between different populations of bacteria may occur.⁹ Also, selection pressure (e.g. from an antibiotic) can cause an identical property to evolve independently in unrelated phylogenetic lineages. The latter effect has been clearly demonstrated in Yersinia pestis where the phylogenetic populations do not correspond directly to classical biovars that are based on phenotypic properties.^10,11 Mutations at positions 2058 or 2059 (using Escherichia coli numbering) in domain V of the 23S rRNA gene, rrl, and mutations in the ribosomal genes rplD and rplV encoding the proteins L4 and L22, respectively, are linked to erythromycin resistance in a number of bacterial species.^12,13In silico sequence analysis of a few F. tularensis subsp. holarctica strains has suggested that erythromycin resistance may relate to a nucleotide substitution in rrl position 2059 (A2059C),^14,15 and in vitro experiments have shown that resistance mutations, although not A2059C, are acquired by gradually increasing exposure to erythromycin.¹⁶

There is a modern global phylogenetic framework for F. tularensis based on genome-wide canonical single-nucleotide polymorphisms (canSNPs).^17,18F. tularensis subsp. holarctica displays extremely little genetic variation but can be divided into four major phylogenetic groups by canSNP typing: B.4, B.6, B.12 and B.16.^19–23

We analysed multiple whole-genome sequences and assessed whether the division of F. tularensis subsp. holartica into biovars I and II is associated with one or more specific mutations previously linked to erythromycin resistance. In addition, we investigated the distribution of this resistance trait across different phylogenetic populations of this subspecies.

Materials and methods

Strains and culture conditions

Three-hundred and fourteen F. tularensis subsp. holarctica strains isolated between 1926 and 2011 in 16 countries were included: Austria (n = 1), Bulgaria (n = 1), Czech Republic (n = 4), Finland (n = 1), France (n = 5), Hungary (n = 1), Italy (n = 1), Japan (n = 10), Norway (n = 18), Poland (n = 5), Russia (n = 1), Slovakia (n = 1), Spain (n = 2), Sweden (n = 260), Ukraine (n = 1) and USA (n = 2); see Table S1 for detailed information (available as Supplementary data at JAC Online). In culture experiments described below, F. tularensis was grown on GC II agar with 1% haemoglobin and 1% IsoVitaleX at 37°C in a 5% CO₂ atmosphere under biosafety level 3 laboratory conditions.²⁴

Disc diffusion erythromycin susceptibility analysis

F. tularensis was grown for 24–48 h on 90 mm agar plates and suspended to OD₆₀₀ = 0.13 (∼0.5 McFarland) in 3.0 mL of 0.9% sterile saline solution before being streaked on a new agar plate using a soaked sterile cotton swab. A 15 μg erythromycin disc (Oxoid, Basingstoke, UK) was placed in the centre of the plate and the zone of inhibition was measured after 24 h.

DNA preparation, WGS and canSNP typing

DNA was extracted from bacteria grown on agar for 48 h using the chloroform–phenol protocol described by Larsson et al.²⁵ or a Qiagen EZ1 Advanced instrument with a DNA bacteria card and an EZ1 DNA tissue kit (Qiagen, Hilden, Germany). Three-hundred and two F. tularensis strains were sequenced at SNP&SEQ, Uppsala, Sweden (Illumina Hiseq 2000, 100 bp paired end). Nine strains were sequenced at Source BioScience, LifeSciences, Nottingham, UK (Illumina Genome Analyzer II/IIx, 36 bp single end) and three strains at GATC Biotech, Cologne, Germany (Illumina Genome Analyzer IIx, 76 bp single end). De novo genome assembly was performed using ABySS with a k-value of 79.²⁶ The sequence depth was estimated from the ratio between the number of total bases and the number of bases of the reference genome FSC200 (1 894 157 bp, NC_019551.1/CP003862). The software CanSNPer and assembled draft genomes were used to extract canSNP states.²⁷ The canSNPs used are listed in Table 1.

Sequence analysis of the genes rrl, rplD and rplV

Pairwise alignments between rrl, rplD and rplV in each of the 314 draft F. tularensis subsp. holartica genomes and the reference sequence FSC200 (GenBank CP003862.1; positions 126 515–129 399; 419 484–422 368; 1 123 125–1 126 009) were created using progressiveMauve.²⁸ To evaluate the SNP copy number variation among the three copies of rrl, sequence reads were mapped against the genome of OSU18 (NC_008369.1, position 125 781–129 501) using Bowtie2.²⁹ The SNP positions were then scanned visually in Tablet.³⁰

Allelic exchange of rrl in strain FSC274

The method for allelic exchange in F. tularensis via homologous recombination is described in Golovliov et al.³¹ In brief, a suicide vector containing the 3′ region of rrl including 2059C was constructed using a nucleotide sequence obtained by PCR of the erythromycin-resistant LVS strain (see Supplementary Data). The forward primer was located 40 bp upstream and the reverse primer ∼1200 bp downstream of 2059C and these primers included restriction sites for SalI and SacI, respectively (Eurofins MWG Operon, Ebersberg, Germany). The proofreading Platinum Pfx polymerase was used for PCR amplification according to the manufacturer's protocol (Life Technologies) followed by SalI and SacI digestion (New England Biolabs, Ipswich, MA, USA). Agarose gel electrophoresis and a Qiagen Gel Purification kit (Qiagen) was used to isolate a PCR fragment which was subsequently cloned into a SalI/SacI-cleaved pDM4 vector³² and introduced into E. coli S17-1 λpir by electroporation. The vector was transferred to FSC274 by conjugation.³¹ Clones with plasmid integrated into the chromosome were selected on agar plates containing chloramphenicol (8 mg/L) and polymyxin B (50 mg/L). Legitimate integration was tested using the PCR primers 5′-GTCAGCTAGCACTTATCTATTGACTATATGT-3′ and 5′-GACAACAAGCCAGGGATGTAACGCA-3′.

Erythromycin susceptibility of allelic replacement clones

FSC274 clones were grown on agar supplemented with chloramphenicol (8 mg/L), and WT FSC274 and LVS on agar without chloramphenicol. After 48 h the clones and strains were suspended in 0.9% sterile saline solution to OD_600 = 0.4 (viable counts ∼2.5 × 10⁹ cfu/mL) and streaked on 90 mm agar plates using a soaked sterile cotton swab. An erythromycin 0.016–256 mg/L Etest strip (Biomérieux, Marcy l′Étoile, France) was placed in the centre of the plate. The MIC was monitored at 24, 48 and 72 h.

Sequencing of FSC274 clones and mapping of sequence reads

FSC274 clones and the FSC274 WT were whole-genome sequenced using an MiSeq reagent kit v. 3 (600 cycles, 300 bp paired end) and an Illumina MiSeq (Illumina, San Diego, CA, USA). DNA was extracted using an EZ1 Advanced instrument (Qiagen) as described above, with the bacterial suspensions first incubated at 65°C for 1.5 h followed by RNase A treatment (0.5 μg/μL) (Qiagen) at 37°C for 30 min. DNA libraries were prepared using a Nextera XT kit (Illumina). Sequencing reads were mapped to rrl of OSU18 as described above and the number of reads containing 2059C counted.

In silico search for resistance mutations in all subspecies of F. tularensis and in F. novicida

Any SNPs in the rrl, rplD and rplV genes described in the scientific literature to confer erythromycin resistance in bacteria were examined in multiple genomes of F. tularensis subsp. tularensis (FSC033, FSC043, FSC054, FSC198 and ATCC6223, Bioproject Accession Number PRJNA261819, PRJNA42733, PRJNA73375, PRJNA17375, PRJNA30629) plus in 14 additional unpublished genomes, in subsp. mediasiatica (FSC147 and FSC148, Bioproject Accession Number PRJNA19571, PRJNA73379) plus in one unpublished genome, and in F. novicida (FSC156, FSC159, FSC160 and U112, Bioproject Accession Number PRJNA62751, PRJNA73383, PRJNA73385, PRJNA236529) plus in one unpublished genome.

Literature search for geographical distribution of biovar II and canSNP B.12 strains

PubMed was searched using the search terms ‘Erythromycin resistance AND Francisella’ and ‘Typing AND Francisella’; the database was last accessed on 1 March 2016. Articles in English language with geographical information on erythromycin-resistant or canSNP group B.12 (aka B.Br.012) F. tularensis subsp. holarctica strains were included.

Results

CanSNPer typing of whole-genome sequences

The assembled sequencing reads for the draft genomes of 314 strains provided an average sequence depth of 920× (range 44–4090). Using the CanSNPer software, SNP states for all currently available canSNP markers of F. tularensis subsp. holarctica (n = 43) were retrieved and the genomes were assigned to 17 canSNP groups (Table 1, Figure 1). SNP states per strain are available in Supplementary Data. All previously described major canSNP groups of F. tularensis subsp. holarctica were represented in the data. The majority of the strains were assigned to canSNP groups B.12 (n = 195) and B.6 (n = 100), with smaller numbers of strains assigned to groups B.4 (n = 9) and B.16 (n = 10). Strains belonging to group B.12 were further divided into 11 different subgroups, and those belonging to group B.6 into 4 subgroups (Table 2). Thus, the 314 strains examined in this study covered nearly the entire currently known global diversity of F. tularensis subsp. holarctica, lacking representation only of subgroup B.9 within group B.6 and of subgroups B.27–32, B.37 and B.38 within group B.12. In general, the proportion of strains assigned to each of these major groups is similar to that found in F. tularensis strain collections from the environment and from human cases (i.e. B.12 is the most common group).

Correlation between canSNP group and erythromycin resistance

All 195 strains in group B.12 demonstrated an erythromycin-resistant phenotype with no visible zone of inhibition and hence corresponded to biovar II (Supplementary Data). In contrast, all other 119 strains displayed a zone of inhibition and belonged to groups B.4 (average zone diameter = 34 mm, SD = 4.28), B.6 (31 mm, SD = 3.15) and B.16 (27 mm, SD = 4.54). Overall, there was thus a strict correlation between each of three phylogenetic groups with the three biovars named biovar I (corresponding to phylogenetic group B.4 and B.6), biovar II (group B.12) and biovar japonica (group B.16) (Figure 1).

Sequence analysis of genes previously associated with erythromycin resistance

All the SNPs found in the rrl gene among the 314 F. tularensis subsp. holarctica genomes are shown in Table 2. Only two SNPs, A453G and A2059C (E. coli numbering), were unique to group B.12, corresponding to the biovar II phenotype. A full rrl gene alignment is available in Supplementary Data. The A453G SNP was located within domain I and the A2059C SNP within the peptidyl transferase centre of domain V. An exhaustive search for possible SNP copy number variation among the three copies of the rrl gene among 34 F. tularensis subsp. holarctica genomes revealed that SNPs were universally present in all three copies of the rrl gene (the genomes examined are indicated in Supplementary Data). A search for A453G and A2059C in the rrl genes of additional F. tularensis genomes revealed that there were no such substitutions in multiple genomes of subspecies tularensis (n = 19), mediasiatica (n = 3) or in F. novicida (n = 5).

Alignments of the rplD and rplV genes revealed several SNPs in rplD but none in rplV among the 314 genomes (an rplD gene alignment is available in Supplementary Data). There were two SNPs in the rplD gene, C309 T and A543G (E. coli numbering), but neither of them corresponded directly to the presence or absence of erythromycin resistance.

Erythromycin susceptibility and rrl sequence analysis of allelic replacement clones

Because SNP A2059C in rrl is associated with erythromycin resistance in other bacteria and exhibited a perfect correlation with group B.12, we analysed its contribution to the resistance phenotype by introducing 2059C via allelic exchange in the naturally susceptible F. tularensis subsp. holarctica strain FSC274 of group B.6. Sequence analysis revealed that the allelic exchange was unsuccessful in six of the seven clones as they exhibited no reads with 2059C, similar to the FSC274 WT strain. However, approximately one-third (30.4%; n = 460) of the sequencing reads from the one other clone (clone 2) exhibited 2059C, which is consistent with successful allelic exchange for one out of three copies of rrl (sequencing read frequencies are available in Supplementary Data).

The erythromycin MICs for the six FSC274 clones with no evidence of successful allelic exchange, as well as the FSC274 WT strain, ranged from 1 to 2 mg/L (read at 48 h). In contrast, the naturally resistant control strain LVS, belonging to canSNP group B.12, displayed a MIC ≥256 mg/L and showed no signs of reduced growth rate. For the remaining clone (clone 2), a zone of inhibition was visible at 48 h. The zone was, however, less distinct than the inhibition zone of FSC274 WT (Supplementary Data). After another 24 h the inhibition zone of clone 2 vanished, indicating the capacity of clone 2 to grow at high erythromycin levels (≥256 mg/L), albeit at a reduced growth rate compared with LVS.

Geographical distribution of biovar II alias canSNP group B.12 strains

A PubMed search identified 15 articles published between 1978 and 2016 to be available for analysis of the geographical distribution of biovar II alias canSNP group B.12 strains (article data are compiled in Supplementary Data). A total of 1144 F. tularensis subsp. holarctica strains from 20 countries were described as either biovar II or canSNP group B.12 (when known, data of strain duplicates between articles were excluded). For comparison, 674 strains described in the same articles as either biovar I, or belonging to canSNP group B.4, B.6 or B.16, were also included. The geographical distribution of the strains (n = 1818) is shown in Figure 2.

Discussion

In this work, we provide a genetic basis for the phenotypic erythromycin resistance that has been used to divide strains of F. tularensis subsp. holarctica into biovar I and biovar II. We found that all erythromycin-resistant strains (biovar II strains) belong to a specific major phylogenetic lineage of the subspecies canSNP group B.12, which is a cohesive genetic population. Our findings suggest that erythromycin resistance evolved in an ancestor of the B.12 group and has been maintained in all of the descendants of this group through strict clonal inheritance. We also found that the erythromycin-resistant B.12 group is the most successful phylogenetic group of F. tularensis in Eurasia.

Two SNPs in the ribosomal gene rrl were perfectly associated with erythromycin resistance: A2059C and A453G. These two SNPs were present in all B.12 strains and absent in all other strains. The additional SNPs in rrl and all SNPs found in the gene rplD could be excluded as resistance candidates based on comparison with results of erythromycin susceptibility analyses; i.e. these SNPs were not limited to the B.12 group. We judged the A2059C mutation in rrl to be the most likely candidate to confer resistance because, according to the scientific literature, mutation at position 2058 or 2059 disturbs the antimicrobial effect of erythromycin by blocking its binding to the peptidyl transferase centre of the ribosome in several bacterial species.¹³ In F. tularensis, 2059C was already suspected to confer resistance based on in silico analyses of the genomes of a few strains of F. tularensis subsp. holarctica.^14–16 Our comprehensive genomic evaluation of 314 strains, representing nearly all known genetic diversity of the subspecies, substantially strengthens this suspicion. In addition, we show a gain of resistance phenotype in an erythromycin-susceptible clinical strain by experimentally introducing a nucleotide sequence with 2059C in one gene copy of the rrl gene. Due to the significant challenge of producing allelic replacement mutants in F. tularensis, we did not test the A453G substitution for a possible contribution to resistance. It seems unlikely that this substitution is involved considering our experimental results and its location within domain I of the 23S rRNA that has not previously been described to interact with macrolides.

A gain of resistance phenotype in our susceptible strain was achieved by introduction of the A2059C mutation in one rRNA operon, but this gain appeared to come at a price of reduced growth rate when compared with a WT strain. We found that the allelic replacement clone with only one of three copies carrying the resistance mutation was unlikely to be optimized for growth with regard to rRNA operon gene dosage. Our result is consistent with previous experiments in E. coli that showed marked reduction in growth rate for strains with fewer rRNA operon copy numbers.^33,34 The growth reduction of our experimental clone suggests that one 2059C mutant rRNA operon copy combined with two WT rRNA copies within a single genome result in reduced growth fitness. The WT-resistant genotype carries the signature of 2059C in all three rRNA operon copies and this homogeneity of the rRNA may be required to make a strain fit for normal growth. We can only speculate that a process of gene conversion (wherein one DNA sequence replaces a homologous sequence such that the sequences become identical) has occurred during evolution of the B.12 group. Gene conversion is an rapid evolutionary process conferring, for example, adaptation to changes in selection pressures from antibiotics such as erythromycin by homogenization of favourable rRNA operon mutations among multiple copies.^35–37 With subsequent population bottleneck effects, which are likely to be common in F. tularensis because of its extremely low infectious dose, a resistance mutation may rapidly achieve fixation in all rRNA copies of a genome.

The finding of the A2059C SNP in the rrl gene causing erythromycin resistance and the presence of this SNP restricted to genomes of the phylogenetic group B.12 suggests this mutation has occurred only once during the evolution of F. tularensis. The SNP appears to have been maintained and propagated through strict clonal inheritance within the B.12 group only, consistent with the overall clonal nature of the bacterium.^18,21,38 We found no evidence for horizontal transfer of the resistance trait or recurrent evolution in other genetic lineages. Thus, it appears that as a consequence of the clonality of F. tularensis the erythromycin resistance phenotypic trait is indeed a stable marker for a specific and cohesive genetic subpopulation of F. tularensis subsp. holarctica, supporting the claim by Olsufjev and Meshcheryakova⁷ that erythromycin resistance in F. tularensis is a stable marker that can be used to distinguish bacterial populations below the subspecies level. We also found that the biovar designation japonica, which is based on the ability to ferment glycerol, correlated to the phylogenetic group B.16.

Interestingly, erythromycin resistance coincides with the most successful genetic population of F. tularensis as defined by the number of strains reported in the scientific literature. A reporting bias due to epidemiological factors facilitating isolation of erythromycin-resistant strains in favour of susceptible strains is an intriguing possibility but seems improbable because resistant strains predominate in diverse epidemiological contexts. Resistant strains dominate in tick-borne epizootics among hares in Hungary, in large epidemics among humans infected by contaminated drinking water across Turkey, and in mosquito-borne epidemics among humans in Sweden, Finland and Russia.^6,39–41 We do not think that human activity caused erythromycin resistance in F. tularensis due to the relatively late discovery of erythromycin. The drug was discovered in 1949 in an environmental sample from the Philippines and large-scale production began in the USA during the 1950s.⁴² The erythromycin-resistant type strain of biovar II, strain GIEM 503, was isolated as early as 1949 in Russia.⁸ Erythromycin, however, can be naturally produced by bacteria found in soil, principally actinomycetes.^43–45 For this reason, erythromycin resistance traits are found to be ubiquitous among soil bacteria, including those isolated from pristine soils.⁴⁶ In light of this, there is a possibility that erythromycin resistance confers a selective advantage on F. tularensis during its natural life cycle, including its long-term presence in moist soil and on straw.^47–49

By providing a genetic basis for the biovar I/biovar II division, it was possible to reconcile and compare publications with a whole-genome approach based on canSNPs and publications using the biovar concept. This gave deeper insights into the geographical distribution of erythromycin-resistant F. tularensis subsp. holarctica, i.e. canSNP group B.12 by adding >1800 strains to the analysis. We found the B.12 group to exhibit a wide geographical distribution spanning almost the entire Eurasian continent—from Central Europe all the way to the Primorsky region in the Far East. The B.12 group completely dominates Central and Eastern Europe and the regions around the Black Sea, supporting the idea that this is a very successful genetic subpopulation of F. tularensis.

In conclusion, we have shown that erythromycin-resistant strains perfectly correlate with strains of canSNP group B.12. Our results demonstrate that this strictly clonal erythromycin-resistant population within F. tularensis has been incredibly successful both in terms of its geographical distribution and the numbers of strains reported in the scientific literature.

Funding

This work was supported by the Swedish Civil Contingencies Agency (grant number TA 014-2010-01) and the US Department of Homeland Security's Science and Technology Directorate (award number HSHQDC-10-C-00139) pursuant to the agreement between the US government and the Kingdom of Sweden on Cooperation in Science and Technology for Homeland Security Matters; and also the Swedish Ministry of Defence (FOI project A404014).

Transparency declarations

None to declare.

Acknowledgements

We thank numerous researchers, physicians, veterinarians for over many years providing strains to the Francisella strain collection maintained at the Swedish Defence Research Agency (FOI), Umeå, Sweden. We also acknowledge sequencing support from SNP&SEQ, Uppsala Sweden; a part of the Swedish national infrastructure for large-scale sequencing.

References