Comparison of Francisella tularensis genomes reveals evolutionary events associated with the emergence of human pathogenic strains

The genomic nucleotide sequence is highly conserved between the three strains but different mutation rates are apparent

Although no official genomic criteria exists to classify strains into species, Konstantinidis and coworkers [31] found that almost all 70 strains in their study set that reside in the same species exhibited greater than 94% average nucleotide identity (ANI). They also showed that the classification based on ANI correlates with classifications performed with 16S RNA sequences, DNA-DNA re-association, and mutation rate. In comparison, the few sequences of the other Francisella species available in Genbank, namely Francisella philomiragia, exhibit an ANI of 91.66% with the genome of U112. The ANI corroborates the proposition that novicida arose by diverging from an ancestor common to the subspecies tularensis and holarctica, and that the subspecies tularensis and holarctica subsequently diverged from a common ancestor [31, 32]. Based on the average level of nucleotide identity between the three genomes, it is possible to estimate the rate of substitution in the genomes of holarctica and tularensis after their divergence. The genomes of holarctica strains are estimated to have evolved at an average rate of 0.55 base pairs (bp)/100 bp from the common ancestor, whereas the genome of Schu S4 diverged at the lower rate of 0.25 bp/100 bp.

Genome reorganization occurred in the human pathogenic F tularensisancestral strain during or after differentiation from the nonpathogenic strain

A recent study using paired-end sequencing [24] indicated that the organization of the genomes of holarctica strains and tularensis strains is not conserved. However, the organization was highly similar for the genomes of the 67 holarctica strains analyzed. Similarly, the genome of holarctica strain OSU18 is collinear with the genome of the holarctica strain LVS, but it is organized differently than the genome of Schu S4 [16]. These findings extend the phylogenetic and molecular evidence that the strains are mostly clonal in the subspecies holarctica and that their genome is relatively stable [18, 32–34]. The subspecies tularensis can be divided into two distinct groups (type AI and AII) [18, 35]. According to amplified fragment length polymorphism and restriction fragment length polymorphism analyses, genomes in the subspecies tularensis are organized differently but are similar within groups [33, 34]. Hence, the genome of LVS is representative of all genomes in the subspecies holarctica, whereas the genome of Schu S4 represents genomes in the type AI group.

Sequence alignment of the U112 and Schu S4 genomes reveals 59 chromosomal segments with the same gene content and gene order in both organisms, but arranged differently throughout both genomes (Figure 1). Chromosomal segments with the same gene content and gene order in two bacterial genomes are hereafter termed 'syntenic regions'. The discrepancy in the order of the chromosomal segments between the two genomes suggests that regions have been moved, in one genome or the other. Hence, there are a total of 118 genomic breakpoints when comparing the two genomes. Similarly, 59 syntenic regions are arranged differently when comparing the genomes of U112 and LVS, and 51 are arranged differently between the genomes of Schu S4 and LVS (Figure 1), which is the same amount as found when comparing Schu S4 and OSU18 genomes [16]. Twenty-eight out of the 59 syntenic blocks (47%) are nearly identical in the genomes of Schu S4 and LVS relative to the genome of U112. However, the order in which the blocks are arranged differs greatly. This suggests that these syntenic blocks formed before differentiation between both human pathogenic subspecies, but moved independently later in one or both genomes. The rest of the syntenic blocks in LVS and Schu S4, in comparison with U112, differ both in content and order (Figure 1), which suggests that they formed after differentiation of the two subspecies.

Localization of IS elements at genomic breakpoints suggests that IS elements are involved in most genomic rearrangements in the human pathogenic strains

Six types of IS elements were identified in the three genomes. Five of them are present in the three genomes at least in a remnant form, whereas one, ISFtu5, is only present in the subspecies holarctica and tularensis. As shown in Table 1, the number of each IS element varies greatly in the three strains. The difference in numbers of ISFtu1 and ISFtu2 elements is particularly large. It suggests that ISFtu1 has transposed and proliferated in the genomes of the subspecies tularensis and holarctica, or in the genome of their common ancestor. ISFtu2 exhibits more proliferation in the holarctica genome. ISFtu1 appears to have been replicated essentially in the ancestor of holarctica and tularensis strains becuase 46 out of 53 elements are bordered by the same sequences in both genomes. Nine ISFtu1 elements exhibit the same bordering regions on both sides in the two subspecies genomes. However, 37 other ISFtu1 elements share only one side with an element in the other genome, indicating rearrangements specific to each subspecies. About 13 ISFtu2 elements may have transposed in the ancestral genome of tularensis and holarctica, as indicated by common bordering sequences, but have undergone subsequent rearrangements because ten ISFtu2 elements have only one common side.

These findings strongly support the proposition that genomic rearrangements occurred in the genomes of the tularensis and holarctica strains by homologous recombination at ISFtu1 and ISFtu2 elements [16]. This proposition is also supported by the fact that 82% of breakpoints of LVS-Schu S4 syntenic blocks are bordered by an IS element within 100 bp (Figure 1). Similarly, 60% of the breakpoints in LVS-U112 and Schu S4-U112 syntenic blocks are bordered by IS elements in the genome of the human pathogenic subspecies (Figure 1). This lower incidence may be due to transposition of IS elements subsequent to the initial rearrangement. IS elements appear to play a prominent role in rearrangement events, further corroborating that these events took place in the ancestor of holarctica and tularensis. Indeed, 88% of the Schu S4-U112 syntenic blocks are bordered by an IS element at one extremity or both in the genome of Schu S4. On the other hand, the location of IS elements in the genome of U112 exhibits association with breakpoints for merely four ISFtu2 elements. This suggests that the IS elements did not play a prominent role in the evolution of the strains that are not pathogenic to humans.

In summary, comparative analysis using the genome of U112 revealed that the complex evolutionary scenario of the three F tularensis subspecies involves the transposition of ISFtu1 (tularensis and holarctica) and ISFtu2 (novicida, tularensis, and holarctica), accompanied by replication of these elements and genomic rearrangements at the location of these elements at distinct steps in genome evolution.

The gene content of F t novicidaU112 reveals a species genome backbone

In the genome of U112, 1,731 protein-coding genes, 14 pseudogenes, and seven disrupted genes encoding an IS element transposase were identified. The coding regions (1,751,817 bp) represent 91.72% of the entire genome. Thirty-eight tRNA genes were identified, representing 30 anticodons encoding the 20 amino acids as well as three operons encoding the 5S, 16S, and 23S ribosomal RNAs and tRNAs for alanine and isoleucine. The same RNA genes and operons are found in the genomes of tularensis and holarctica. Overall, 1,813 distinct genes (excluding IS element genes and 33 hypothetical genes that we believe are noncoding) were found in at least one of the three genomes. Out of these 1,813 genes, a total of 1,572 gene sequences (functional or disrupted) are common to the three genomes. Hence, the core gene set may represent about 86.4% of all distinct genes identified in the three genomes (Additional data file 1).

Human pathogenic strains contain genes that are absent from the nonpathogenic strain U112

Human pathogenic strains have undergone substantial loss of function, but not the non-pathogenic strain

Fourteen pseudogenes have been identified in U112 (Additional data file 1). In contrast, the original annotation of Schu S4 listed 201 pseudogenes [25]. Using the genome of U112 as a reference, 53 additional pseudogenes were predicted in the genome of Schu S4 (Additional data file 1) following a procedure described in Materials and methods (see below), most of which were annotated as multiple open reading frames (ORFs) in the published genome. Because the strain LVS was artificially attenuated, it is expected to contain mutations that are not found in any other holarctica genome. Indeed, 11 pseudogene-causing mutations were found to be specific to the LVS genome [39]. We ignored these 11 pseudogenes for the following comparative analysis, because they do not represent a loss of function in the holarctica subspecies as a whole.

When compared with the genome of U112, analysis of the genome of LVS revealed 303 pseudogenes in addition to those contained in IS elements (Additional data file 1). OK The number of protein encoding genes in the genome of LVS and the subspecies holarctica in general may therefore be about 1,400. The higher mutation rate observed in holarctica genomes as compared with tularensis could explain the greater number of pseudogenes. In addition, at least eight genes present in novicida and holarctica were lost by the strain Schu S4, and ten that were present in novicida and tularensis were lost by LVS. A set of 160 genes were inactivated in both LVS and Schu S4. Taking into account gene deletion and inactivation, U112 encodes 164 functions that are no longer active in both holarctica and tularensis strains. Similarly, 18 functions are specific to the strain Schu S4 and potentially to the subspecies tularensis in general (Table 3).

A reduced set of genes was inactivated in the genome of the strain ancestral to human pathogenic strains

A total of 160 genes are inactivated in the genomes of both subspecies holarctica and tularensis. Upon alignment of their sequences, 53% of pseudogenes common to LVS and Schu S4 exhibit at least one common mutation that may have led to their inactivation, whereas 32% of the pseudogenes common to both subspecies share no common variations. The sequence of the remaining 15% is too divergent to determine a potential common inactivating mutation (Additional data file 1). This indicates that at least 53% have arisen in the genome of the human pathogenic ancestor. These 82 pseudogenes bearing common mutations are more likely to be located directly at breakpoints than the pseudogenes not sharing any common mutation (Figure 2b). In addition, the IS insertion is the only inactivating common mutation found in 19 out the 82 pseudogenes from the ancestral strain. This suggests that IS insertions or subsequent sequence rearrangements contributed to at least 22% of the earliest gene inactivations that took place in the emerging human pathogenic strain.

Contribution of IS elements and other early mutations to genome reduction through initiation of genetic drift

When directly compared with the genome of U112, most pseudogenes in the genomes of Schu S4 and LVS appear to result from small indels (1 or 2 bp) or nonsense mutations. In tularensis and holarctica genomes, genes within 1 kb from a genomic breakpoint are twice as likely to be inactivated as were genes in other genomic locations (Figure 2a). The proportion of genes that are within 1 kb from a genomic breakpoint and are inactivated is 28.5% in the genome of Schu S4 (57 out of 200), whereas the global proportion of inactivated genes is 12.6%. Similarly, 24.9% of genes within 1 kb from genomic breakpoints are inactivated in the genome of LVS, whereas the global proportion of inactivated genes is 16.3%. Figure 2a shows that, to a lesser extent, the genes within 3 kb from a breakpoint are also more likely to be inactivated than are the genes in the rest of the genome. In Schu S4, 15.4% of genes between 1 and 2 kb from a breakpoint are inactivated and 17.1% are between 2 and 3 kb. Similarly in LVS 18.8% of the genes between 1 and 2 kb from a breakpoint and 22.1% between 2 and 3 kb are inactivated. It is unlikely that genomic rearrangements could directly have caused mutations as far as 3 kb from the breakpoints. It is more likely that the rearrangements disrupted the transcriptional unit to which these genes belong. If these genes are no longer transcribed, then their sequences are no longer subjected to selection and evolve by neutral genetic drift, eventually causing the disruption of the ORF through mutation.

In agreement with this conjecture, predicted operons located at breakpoints are more likely to contain more than one pseudogene, in Schu S4 by 4-fold and in LVS by 1.4-fold. An additional argument in favor of the inactivation of some genes by genetic drift is the uneven distribution of pseudogenes across functional categories (Figure 2c). Pseudogenes and absent genes of the holarctica and tularensis genomes have been assigned to functional categories based on the annotation of their functional counterpart in the genome of U112. For example, 41.2% of the genes predicted to be involved in amino acid biosynthesis in the genome of novicida are inactivated in the genome of one or both of the other subspecies. Similarly, 43.1% of the genes predicted to encode transporters are inactivated in the genomes of holarctica and tularensis. Remarkably, the distribution in functional categories is the same for genes inactivated in one genome and those inactivated in both. Likewise, it was previously observed in the genomes of Salmonella typhi and S paratyphi that the pseudogenes were different but appeared to belong to the same pathways and operons [11]. The over-representation of pseudogenes in certain functional categories suggests a loss of function associated with specific pathways, resulting in the decay of multiple genes in these categories [40]. Following the disruption of a biologic process by the inactivation of one gene, other genes involved in this process are no longer subjected to selective pressure.

Inactivation of the leucine and valine biosynthesis pathway illustrates the proposed evolutionary scenario

This example illustrates the proposed model of evolution of Francisella human pathogenic strains: initial inactivation of a gene in the ancestor of the subspecies tularensis and holarctica (potentially pathoadaptive) and further gene inactivation in regions no longer subjected to selective pressure before and after subspeciation.

In the genome of U112, the genes involved in leucine and valine biosynthesis are organized in two operons: one contains leuB, leuD, leuC, leuA, and ilvE; and the other one contains ilvD, ilvB, ilvH, and ilvC. All genes are expressed in rich medium (Rohmer and coworkers, unpublished data). In the tularensis and holarctica strains the leucine, isoleucine, and valine biosynthesis pathway is inactivated. Based on the organization of the two regions depicted in Figure 3, we can infer events that took place in leu and ilv loci. Two ISFtu1 elements are associated with the leu operon in both human pathogenic strains and have the same bordering sequences: the same portions of leuA and the upstream sequence of leuB. Hence, the insertion of two ISFtu1 elements has taken place in the leu operon of the ancestor of the two strains and disrupted leuA and the upstream region of leuB. All sequences of the leu operon are still present in the genome of LVS, but they are scattered to three different locations, all associated with ISFtu1 elements. In the genome of Schu S4, leuB, leuD, and leuC have been deleted and one IS element sits in place of the deletion (Figure 3). It seems therefore that the two ISFtu1 elements inserted in the genome of the ancestor underwent different recombination events in each strain. The ilv operon contains distinct mutations in the genome of LVS and Schu S4; in LVS ilvB (FTL_0913-FTL_0914) and ilvD (FTL_0911-FTL_0912) are inactivated by a 100 bp deletion and a 350 bp deletion, respectively, whereas in Schu S4 ilvC (FTT0643) and ilvB (FTT0641) are inactivated because of a nonsense mutation and a single nucleotide deletion, respectively. The distinct origin of the inactivation of the ilv operon indicates that mutations took place after divergence as well.

Potential virulence factors found in the U112 genome and common to all F tularensisstrains

As described in the Introduction (above), virulence strategies overlap in the three subspecies. Here, we provide a list of virulence factors complementary to those previously predicted [16, 25, 41] using the U112 genome as a reference (Additional data file 2). A variety of protein features are potentially indicative of a role in virulence, such as the presence of a protein domain previously associated with a virulence function, the presence of a eukaryotic domain, or homology to eukaryotic proteins sufficiently high to suggest a role in the host cell [42–44]. A total of 129 proteins in U112 revealed one or more of these features. Interestingly, only 80 of them were present and functional in both of the other genomes (Additional data file 2). This suggests that many of these 129 proteins are not involved in virulence or are not essential for the virulence in humans. It is still conceivable that these proteins confer a capacity to infect hosts or to target functions that the subspecies holarctica and tularensis no longer utilize, or they may even be detrimental to the bacterium in the human host.

The ORF FTN_0921 in novicida U112 (FTT1043 in Schu S4) is homologous to a Legionella macrophage infectivity potentiator. FTN_1151 (FTT1170) contains Sel1 eukaryotic tetratrico peptide repeats and is homologous to EnhC and EnhA of Coxiella burnetii, which promote entry of Coxiella into host cells. These two proteins could contribute to entry of the bacteria into the macrophage. FNU1336 (FTT1332) may be a hemolysin. FTN_0403 (FTT0877c) is only homologous to eukaryotic proteins and, in particular, to a family of membrane-bound proteins with which it shares a pair of repeats, each spanning two transmembrane helices connected by a loop. The PQ motif found on loop 2 was shown to be critical for the localization of cystinosin to lysosomes [45]. FTN_0083 (FTT0243) may interact with the cytoskeleton of the host cell because it contains an α-tubulin suppressor or related RCC1 domain. FTN_0171 (FTT0195) has ankyrin repeats, sometimes present in bacterial virulence factors. Larsson and coworkers [25] pointed out that the genome of Francisella tularensis does not encode any of the secretion systems that are usually associated with pathogenicity (type III and type IV). A protein homologous to toxin secretion ABC transporters (FTN_1693) and HlyD-family secretion proteins (FTN_0029, FTN_0718, and FTN_1276) may play a role in the delivery of virulence factors. It has been shown that a secretion system similar to type II and type IV systems is responsible for the secretion of virulence factors in U112 [29]. TolC appears to play a role in virulence in U112 as well in holarctica strains [46]. Secretion through these systems first requires protein translocation through the bacterial inner membrane via an independent export system. A full and functional sec system was identified in the genome of U112 as well as in the genomes of Schu S4 and LVS. This suggests that some of the proteins that are exported outside the cell may contain a signal peptide, promoting their translocation across the inner membrane via the sec system. Hence, we suggest that there may be proteins that interact with host factors that are yet to be identified among the set of proteins with a predicted signal sequence.

Functions specific to the human pathogenic subspecies holarctica and tularensis

We consider functions to be specific to the human pathogenic subspecies if either their DNA sequence is solely found in these strains, or their counterparts in the nonpathogenic novicida are inactivated. We have found 41 genes whose DNA sequence is specific to holarctica and tularensis and five genes common to these subspecies that are pseudogenes in U112. In addition, there are 20 duplicated genes in Schu S4 and 34 in LVS. Included in this set is the duplicated pathogenicity island, of which there is only one copy in U112 [26]. The duplication of the Francisella pathogenicity island may provide a higher level of expression of the virulence genes it carries, as it is the case for the Shiga toxin genes in Shigella dysenteriae 1 [47]. Potentially, greater expression of these pathogenicity genes could play a role in virulence in humans.

Among the 41 genes found solely in the genome of the holarctica and tularensis subspecies, 24 have no predicted function (Table 2). Some of the 41 genes could be linked to the pathogenicity of the human pathogenic strains. Six genes involved in the biosynthesis of the O-antigen of lipopolysaccharide in type A and type B strains have no counterparts in U112. The U112 subspecies carries a different set of genes for this function. This could explain the difference noted in the structure of the O-antigen of U112 as compared with those of tularensis strains [48]. The difference in the O-antigen part of the lipopolysaccharide structure could contribute to the difference in host range observed between the three subspecies. In addition to sequence-specific genes, five U112 pseudogenes are functional in both holarctica and tularensis. It may be that inactivation of these genes impairs the virulence of the strain U112 in humans, but the functions they encode do not suggest this possibility. Two of these genes encode nicotinamide ribonucleoside (NR) uptake permease family proteins (FTT0707 and FTT1090), but four other genes found in the U112 genome encode proteins of this family and some of their counterparts have become pseudogenes in the genome of holarctica and tularensis strains. Hence, these genes may have been inactivated because of functional redundancy. FTT0666c (homologous to some methylpurine-DNA glycosylases), inactivated in U112, may be involved in DNA repair following DNA damage induced by stress. FTT1076 (hipA), a protein that potentially is involved in persistence after exposure to antimicrobial products or other stressful conditions [49], is also inactivated in U112. It is therefore possible that U112 may be less resistant to human responses than the holarctica and tularensis strains. Finally, FTT1450c, wbtM on the O-antigen gene cluster, encodes a dTDP-D-glucose 4,6-dehydratase. Because some components of lipopolysaccharide are missing in U112, it is possible that FTT1450c in U112 has degenerated over time because of lack of selection. It would be interesting to examine the state of these five genes in the novicida strains isolated in humans [20, 21, 50].

Some of the functions specific to F tularensis subspecies tularensisSchu S4 may promote the high virulence of type A strains

Comparison between the three genomes reveals regions encoding nine proteins specific to Schu S4 and potentially to the subspecies tularensis. The RD8 11.1 kb specific region [37] carries six functional genes and two pseudogenes (FTT1066 to FTT1073). Three genes in this region suggest that it could be a phage remnant: a type III restriction-modification system restriction enzyme that is apparently nonfunctional (FTT1067); a DNA helicase, which is also nonfunctional (FTT1070); and a predicted antirestriction protein (FTT1071). The five other proteins have no predicted function. This region is bordered on each side by ISFtu1 elements. Because it is specific to all type A strains and exhibits properties of genomic islands (low G+C content and proteins related to mobile elements), the region may be a pathogenicity island that contributes to the virulence of tularensis. FTT1580c, a hypothetical protein, was detected in the region of difference RD1 [37] as specific to the subspecies tularensis. Two hypothetical proteins, namely FTT1308c and FTT1791, were also determined to be specific to Schu S4 in the three-way comparison. They were not detected in the regions of difference obtained by Broekhuijsen and coworkers [37] and Svensson and colleagues [38], and so it is possible that these genes are not specific to tularensis strains or are not present in all tularensis strains. Alternatively, the differences are not detectable with the techniques used by the authors.

In addition to the sequence-specific functions, some functions (encoded by 20 genes) are specific to Schu S4 because they are pseudogenes or absent in the genomes of U112 and LVS. Table 3 lists these 20 genes. A predicted O-methyltransferase (FTT1766) is only functional in Schu S4, and could influence the composition of the bacterial surface. FTT0939, an adenosine deaminase, is only functional in type A strains. This enzyme is predicted to be involved in purine salvage. This could be important to consider for vaccine design, because inactivation of the purine biosynthesis pathway of a type A strain may not result in the significant reduction of fitness that has been observedin type B [51] and novicida strains (data not shown).

Loss of function specific to holarctica may be responsible for the lower level of virulence of these strains when compared with tularensisstrains

Eight additional genes involved in regulation are inactivated in the genome of holarctica alone (Additional data file 1). Six of these genes belong to the LysR transcriptional regulator family. The regulators of the LysR family have diverse targets, including virulence genes and genes that are involved in response to a specific environment. The genome of holarctica strains also exhibits a higher number of pseudogenes in the functional category 'motility, attachment, and secretion structure'. Although three genes encoding potential pilins are inactivated in both subspecies, the holarctica genome underwent inactivation of four additional genes encoding pilins and two predicted to encode membrane fusion proteins. Attachment and motility are key aspects of pathogenicity, and inactivation of these genes may lower the efficiency of infection of humans by holarctica strains. In addition, six genes that are potentially involved in DNA repair are solely inactivated in holarctica (including one encoding a photolyase that repairs mismatched pyrimidine dimers, and one that encodes the protein mutT, which is involved in removing an oxidatively damaged form of guanine). This could explain the higher rate of mutation in holarctica strains than in tularensis strains, and may indirectly be responsible for the inactivation of genes that are important for the pathogenicity of holarctica strains.

Loss of function common to tularensis and holarcticaprovide clues to possible pathoadptation and to the properties of the environmental niches they occupy during their life cycle

Our data suggest that more than half of the pseudogenes in the human pathogenic strains appeared relatively late in their evolution, after the subspeciation. If pathoadaptive mutations occurred, then it is more likely that they took place before the divergence of the pathogenic strains, rather than twice, independently in each pathogenic subspecies. The 84 pseudogenes in the two human pathogenic strains that have arisen in the genome of their common ancestor are listed in Additional data file 1. Significantly, the gene pepO is part of these early mutants in the human pathogenic strains. This gene is active in U112, but a strain U112 in which pepO (FTN_1186) is inactivated spreads more to systemic sites [29]. Similarly, the system used to secrete pepO and other proteins [29] was also altered in the ancestor of the human pathogenic strains (FTN_0306 and FTN_0389). The distribution of the early pseudogenes across functional categories is similar to the distribution of the entire set of pseudogenes (data not shown). However, although eight independent pathways of amino acid biosynthesis are inactivated in one or both human pathogenic strains (24 genes), only one biosynthesis pathway is inactivated in the ancestral strain: the biosynthesis pathway for leucine, isoleucine, and valine. This suggests that the biosynthesis of most amino acids is not required in the current niche of tularensis and holarctica subspecies, but also that only leucine/isoleucine/valine biosynthesis may have played a role in preventing virulence in the human niche. Three transcriptional regulators are inactivated in both genomes: two regulators of the LysR family, and kdpD and kdpE, which form a two-component regulator. Numerous genes encoding transporters are also inactivated. Hence, it is apparent that the tularensis and holarctica subspecies have lost their ability to adapt to or exploit some conditions, and perhaps have undergone niche restriction.