Comparative genomic characterization of Francisella tularensis strains belonging to low and high virulence subspecies

Abstract

Tularemia is a geographically widespread, severely debilitating, and occasionally lethal disease in humans. It is caused by infection by a gram-negative bacterium, Francisella tularensis. In order to better understand its potency as an etiological agent as well as its potential as a biological weapon, we have completed draft assemblies and report the first complete genomic characterization of five strains belonging to the following different Francisella subspecies (subsp.): the F. tularensis subsp. tularensis FSC033, F. tularensis subsp. holarctica FSC257 and FSC022, and F. tularensis subsp. novicida GA99-3548 and GA99-3549 strains. Here, we report the sequencing of these strains and comparative genomic analysis with recently available public Francisella sequences, including the rare F. tularensis subsp. mediasiatica FSC147 strain isolate from the Central Asian Region. We report evidence for the occurrence of large-scale rearrangement events in strains of the holarctica subspecies, supporting previous proposals that further phylogenetic subdivisions of the Type B clade are likely. We also find a significant enrichment of disrupted or absent ORFs proximal to predicted breakpoints in the FSC022 strain, including a genetic component of the Type I restriction-modification defense system. Many of the pseudogenes identified are also disrupted in the closely related rarely human pathogenic F. tularensis subsp. mediasiatica FSC147 strain, including modulator of drug activity B (mdaB) (FTT0961), which encodes a known NADPH quinone reductase involved in oxidative stress resistance. We have also identified genes exhibiting sequence similarity to effectors of the Type III (T3SS) and components of the Type IV secretion systems (T4SS). One of the genes, msrA2 (FTT1797c), is disrupted in F. tularensis subsp. mediasiatica and has recently been shown to mediate bacterial pathogen survival in host organisms. Our findings suggest that in addition to the duplication of the Francisella Pathogenicity Island, and acquisition of individual loci, adaptation by gene loss in the more recently emerged tularensis, holarctica, and mediasiatica subspecies occurred and was distinct from evolutionary events that differentiated these subspecies, and the novicida subspecies, from a common ancestor. Our findings are applicable to future studies focused on variations in Francisella subspecies pathogenesis, and of broader interest to studies of genomic pathoadaptation in bacteria.

Figure 1. Pairwise alignments between five new Francisella genome sequences and a reference genome exhibit >95% sequence conservation.

Genome comparative maps were constructed using CGview software to map pairwise blastn alignments between several Francisella genomes (minimum percent identity = 95 and expected threshold = 1e-5). Specifically, five newly sequenced Francisella genomes (F. tularensis subsp. holarctica FSC257 and FSC022; F. tularensis subsp. tularensis FSC033; and F. tularensis subsp. novicida GA99-3548, and GA99-3549 strains) were aligned to the F. tularensis subsp. holarctica OSU18 reference sequence (outside blue track of genome map). A high degree of similarity between the genomes (>95%) is evident from the continuous blocks of synteny (colored regions).

Figure 2. Dotplot comparison between the Francisella tularensis subsp. holarctica strains showing the occurrence of significant rearrangement events.

Whole genome alignments and dotplot comparisons between Type B strains of Francisella: (A) F. tularensis subsp. holarctica OSU18 and F. tularensis subsp. holarctica FSC022 (reference genome) and (B) F. tularensis subsp. holarctica LVS, F. tularensis subsp. holarctica OSU18 (reference genome), F. tularensis subsp. holarctica FSC257, F. tularensis subsp. holarctica FTA. Alignments were filtered for overlap percentages greater than or equal to 90%. Dotplot (B) shows a nearly linear, overlapping alignment between all of the Type B strains of the main holarctica lineage (not all of the strains of the Type B radiation lineage are shown for clarity), with the notable exception of comparisons between the F. tularensis subsp. holarctica OSU18 and F. tularensis subsp. holarctica FSC257 strains. In contrast, numerous rearrangements are evident from comparisons between the F. tularensis subsp. holarctica FSC022 japonica strain and all of the other Type B strains of the main holarctica lineage (A) (only comparison to the F. tularensis subsp. holarctica OSU18 strain is shown here).

Figure 3. Phylogeny of 20 Francisella strains by analysis of Genome-wide SNP sequences and inferred using Maximum Parsimony method.

The evolutionary history of 20 Francisella strains was inferred from analysis of genome-wide SNP sequences. Genome-wide SNPs at least 20 bp apart were selected for further analysis using MEGA4 software. The tree is drawn to scale, with branch lengths calculated using the average pathway method and are in the units of the number of changes over the whole sequence. Francisella strains sequenced at the Broad Institute are in red. The differentiation of F. tularensis subsp. novicida predated differentiation of the more pathogenic F. tularensis subsp. tularensis and F. tularensis subsp. holarctica subspecies from a common ancestor. Our data indicates differential rates of evolution along the Type A.II and Type A.I lineages, as evident from the branch lengths. The Type B FSC022 strain diverged prior to the radiation of the main holarctica group and also has a greatly reduced branch length, consistent with a much slower rate of evolution. The F. tularensis subsp. mediasiatica FSC147 strain is phylogenetically more closely related to the Type A.II lineage even though its pathogenicity is characteristic of Type B strains.

Figure 4. Multiple alignment of Francisella genomes of the Type B lineage identifies conserved sequence regions with rearrangements; pseudogenes are found proximal to predicted breakpoints.

A comparison of genome rearrangement patterns between the more ancestral F. tularensis subsp. holarctica FSC022 japonica strain and representative strains of the main holarctica group (F. tularensis subsp. holarctica LVS and F. tularensis subsp. holarctica OSU18) was done using MAUVE (A). MAUVE uses an anchored alignment algorithm that permits reordering of the alignment anchors for identification of rearrangements. Colored Local Collinear Blocks (LCBs) are regions of sequence alignment that are free of rearrangements. Each LCB is defined by the anchor regions or predicted sites of rearrangement. A default LCB cutoff of 175, and a filtering for larger blocks containing 10 Kb or larger was done in Mauve. Sequence inversions are denoted by differential positioning of the LCBs relative to a reference axis. A zoomed-in section of the whole-genome alignment is shown so that the annotated ORFs are visible (black outlined boxes). The small red ORFs are rRNA genes. ORFs proximal to predicted rearrangement breakpoints (red circles) have been color coded and labeled (FWD, INV). A summary of the predicted protein products for these genes is provided in (B). These include genes that have been identified as being either disrupted or absent in the F. tularensis subsp. holarctica subspecies.

Figure 5. Comparisons of pairwise sequence alignments provides evidence of pseudogene enrichment within 1 kb of predicted breakpoints in Type B strains.

Pairwise sequence alignments using F. tularensis subsp. holarctica (Type B) strains as a reference were further analyzed for evidence of gene decay proximal to predicted breakpoint sites. Total ORF counts were used to determine the percentage of pseudogenes (y-axis) relative to the distance from the predicted breakpoint (x-axis). Pseudogene enrichment was determined using Fisher's Exact Test to derive p-values. Our findings show the highest percentages of pseudogenes are found within 1 kb of the predicted breakpoints indicating that genome rearrangement events have likely promoted gene decay in strains of the Type B lineage.