Interrogation of the Burkholderia pseudomallei Genome to Address Differential Virulence among Isolates

Infection by the Gram-negative pathogen Burkholderia pseudomallei results in the disease melioidosis, acquired from the environment in parts of southeast Asia and northern Australia. Clinical symptoms of melioidosis range from acute (fever, pneumonia, septicemia, and localized infection) to chronic (abscesses in various organs and tissues, most commonly occurring in the lungs, liver, spleen, kidney, prostate and skeletal muscle), and persistent infections in humans are difficult to cure. Understanding the basic biology and genomics of B. pseudomallei is imperative for the development of new vaccines and therapeutic interventions. This formidable task is becoming more tractable due to the increasing number of B. pseudomallei genomes that are being sequenced and compared. Here, we compared three B. pseudomallei genomes, from strains MSHR668, K96243 and 1106a, to identify features that might explain why MSHR668 is more virulent than K96243 and 1106a in a mouse model of B. pseudomallei infection. Our analyses focused on metabolic, virulence and regulatory genes that were present in MSHR668 but absent from both K96243 and 1106a. We also noted features present in K96243 and 1106a but absent from MSHR668, and identified genomic differences that may contribute to variations in virulence noted among the three B. pseudomallei isolates. While this work contributes to our understanding of B. pseudomallei genomics, more detailed experiments are necessary to characterize the relevance of specific genomic features to B. pseudomallei metabolism and virulence. Functional analyses of metabolic networks, virulence and regulation shows promise for examining the effects of B. pseudomallei on host cell metabolism and will lay a foundation for future prediction of the virulence of emerging strains. Continued emphasis in this area will be critical for protection against melioidosis, as a better understanding of what constitutes a fully virulent Burkholderia isolate may provide for better diagnostic and medical countermeasure strategies.


Introduction
Melioidosis, the disease caused by Burkholderia pseudomallei, presents with a wide range of non-specific signs and symptoms, including fever, pneumonia, acute septicemia, and chronic localized infection [1][2][3]. Initial infection can also be asymptomatic. Chronic stages of the disease are characterized by abscesses in various organs and tissues, most commonly occurring in the lungs, liver, spleen, kidney, prostate and skeletal muscle [1,3,4]. Melioidosis is community-acquired through bacterial contamination of wounds, inhalation, and ingestion [5]. Research in Thailand and Australia has provided critical information about the clinical epidemiology of the disease; the clinical presentations of melioidosis caused by Thai and Australian strains differ in several ways: 1) parotid abscesses are not prevalent in Australia, but occur in Thailand; 2) prostate abscesses are uncommon in Thailand, but are more commonly seen in Australia [3]; and 3) an encephalomyelitis syndrome is seen in tropical Australia more often than in Thailand [5]. This latter condition was associated with the illnesses caused by B. pseudomallei strains MSHR668 [6] and MSHR305 [7]. However, there is evidence that the same strain can cause different clinical presentations in different individuals, and a number of risk factors, such as diabetes have been identified for melioidosis [8]. Therefore, host factors may be important in determining the severity and duration of disease [9,10].
The high incidences of infection in geographical areas where B. pseudomallei is endemic may be due to its resilience and ability to survive under sometimes harsh environmental conditions. B. pseudomallei can survive nutrient depletion, a wide range of pH differences, salt concentrations, and temperatures [11], detergent solutions [12] and acidic environments [13]. It seems that harsh environmental conditions may confer a selective advantage for the growth of B. pseudomallei [5]. These resilience characteristics may explain why B. pseudomallei can cause persistent infections in the human host that are difficult to cure. Also, B. pseudomallei is naturally resistant to a variety of antimicrobial agents [14,15]. In some cases, there is a latency period before symptoms present that can last for days to years [5]. In other cases, an initial acute infection and extensive antibiotic treatment is followed by a variable period of bacterial persistence, with subsequent recrudescence of the disease months or years after the initial infection [16,17].
Our understanding of B. pseudomallei pathogenesis is further complicated by the natural diversity of its genome. B. pseudomallei is a soil-dwelling bacterium that utilizes lateral gene transfer at a very high rate [18]. As a result, there is substantial variation among B. pseudomallei genomes, which may also contribute to differential virulence. Fortunately, as we now have access to many B. pseudomallei genomes from various geographic locations, it is possible to identify genomic features that the various strains have in common, as well as features that are unique to one or more strains.
Comparative studies of genomes from Australian and Thai B. pseudomallei isolates have revealed genomic differences that contribute to our understanding of this organism. The genomes of B. pseudomallei analyzed so far contain from 16-21 genomic islands (GIs) [7,19,20]. The genome of B. pseudomallei strain K96243 contains 16 GIs [19] that are variably present in other B. pseudomallei genomes [20], and each GI shows micro-evolutionary changes that generate GI diversity [20]. In addition to GIs, the genomes of Thai strains K96243 and 1106a contain a horizontally acquired Yersinia-like fimbrial (YLF) gene cluster, while the comparable region in the Australian strains (MSHR668, MSHR305, DM98, 1655 and 13177) is the B. thailandensis-like flagellum and chemotaxis (BTFC) gene cluster [21]. Previous studies showed that BTFC is dominant in Australian strains, while YLF is dominant in strains from Thailand and elsewhere [21]. In addition, clinical isolates are more likely to belong to group YLF, whereas environmental isolates are more likely to belong to group BTFC [21]. In contrast to these trends, we found that the Australian strain MSHR346 contains the YLF cluster (data not shown), and Tuanyok and colleagues reported that 406e, a clinical isolate from Thailand, has BTFC [21].
Previous studies began to address the question of why different strains show differences in virulence and disease presentation. Many studies have focused on host risk factors such as diabetes and alcoholism; but to date only one study has identified genes associated with different disease presentations [22]. This suggests that virulence factors that are variably present in B. pseudomallei strains may be important for pathogenesis. Taken together with the genomic variation, geographical distribution and differences in environmental habitats [18,21,23], comparative genomic studies suggest that strains associated with human melioidosis may possess an accessory genome that differs from animal and environmental strains [24]. We hypothesize that differences in virulence may be associated with variations in metabolic and regulatory capabilities among B. pseudomallei strains.
In this study we compared three B. pseudomallei genomes, from clinical strains MSHR668, K96243 and 1106a, seeking to identify metabolic characteristics that might explain why MSHR668 is more virulent than K96243 and 1106a in a mouse model of B. pseudomallei infection. Analyses focused on genomic features, including metabolic, virulence and regulatory genes that were present in MSHR668 but absent from both K96243 and 1106a. Features present in K96243 and 1106a but absent from MSHR668 were noted, and we also identified virulence-associated genes that were present in all three genomes. Here we have identified genomic features that may contribute to variations in virulence noted among B. pseudomallei isolates.

Comparative Virulence of B. pseudomallei Isolates
For the purposes of this manuscript, we measured the LD 50 upon intraperitoneal (IP) challenge to assess potential differences in virulence among the three B. pseudomallei strains. While studies evaluating clinical infection are complicated by a range of factors such as host risk factors, exposure routes and dose of exposure, experimental studies using inbred mice were used in an attempt to limit the number of host factors that may contribute to differences. Experiments involving infections of BALB/c and C57BL/6 mice [25] with B. pseudomallei strains K96243, MSHR668, and 1106a revealed differences in LD 50 values among the B. pseudomallei strains. LD 50 values were calculated after 21 and 60 days postchallenge. Differences were more pronounced in the BALB/c model, where the LD 50 values of MSHR668 were 30 to 100-fold lower than those of K96243 and 1106a ( Table 1). The LD 50 values for MSHR668 were also lower in the C57BL/6 model, although the differences were not as great. Since K96243 and 1106a had similar virulence properties in both mouse infection models, we were interested in identifying the genomic features that these strains shared but were not common to MSHR668.

Genome Features
We performed an extensive comparative analysis of the B. pseudomallei genomes to identify genomic features that are common and unique among the various strains, and to begin to address differences in virulence and disease presentation. Because the K96243 genome that we downloaded from NCBI contained nearly 1,500 fewer CDS than the other two genomes, we re-annotated all three genomes using the RAST system [26] to ensure consistent comparisons. Table 2 compares the three complete genomes in terms of their general features. Comparisons of the CDS in each genome identified by RAST annotation compared to the original annotations showed that the number of CDS in the K96243 genome increased by 1,317 (18.7%). The numbers of CDS in the MSHR668 and 1106a genomes were also increased, but by smaller percentages (3.5% and 4.2%, respectively). These analyses provided a common annotation platform from which the ensuing comparisons were made.

Pseudogenes and mobile elements
The number of pseudogenes in each originally annotated genome varied depending on the resource used to identify them. Holden et al. (2004) originally reported that the genome of K96243 contains 26 pseudogenes [19], whereas the IMG system [27] identified 122 pseudogenes in K96243, 5 in MSHR668, and 8 in 1106a. The Pathway Tools [28] identified 136 pseudogenes in K96243, 10 pseudogenes in MSHR668, and 15 pseudogenes in 1106a. Because of this discrepancy, and since we re-annotated the genomes using RAST, which does not include an automatic pseudogene identification step, we identified the potential pseudogenes in each genome using the Psi Phi program [29], which is a comparative method for pseudogene identification. Psi Phi identified no additional pseudogenes in the RAST-predicted CDS of K96243, MSHR668 and 1106a. However, Psi Phi identified a few candidate pseudogenes in the intergenic regions, and there were some CDS with less than full length alignments to known protein sequences in the public databases. Since this report does not focus on pseudogenes, we did not explore these further.
The genomes of K96243 and 1106a contained more genes annotated by RAST [26] as encoding mobile elements (79 and 89, respectively) compared to MSHR668 (S1 Table). The number of genes encoding mobile elements that were identified by RAST annotation of K96243 was greater than the originally reported number of 42 mobile elements in K96243 [19]. This discrepancy is likely due to the higher number of total CDS in the RAST annotation of the K96243 genome.

Chromosome alignments
Individual chromosomes of B. pseudomallei MSHR668, 1106a and K96243 were aligned using Mauve [30], and results showed that they are largely collinear, except for an inversion of the K96243 chromosome 1 and a small gap in between the locally collinear blocks in the inverted region ( Fig. 1). Coding sequence comparisons The protein coding sequences (CDS) in common among the genomes (putative homologs) were identified by a bidirectional best BLASTp hits analysis. This also enabled the identification of unique genes that were only present in each genome or group of genomes. Fig. 2 shows the results of the analyses for each pair of genomes, as well as all three genomes together. A total of 5,808 CDS were shared by all three genomes. The pairwise comparisons showed 5,976 CDS shared between K96243 and MSHR668, 6,192 CDS in common between K96243 and 1106a, and 5,970 CDS shared between MSHR668 and 1106a. The distribution of BLASTp hits to strain K96243 is also displayed in a heatmap in Fig. 3. These comparisons included the two pseudomallei strains plus B. thailandensis, B. mallei and other near neighbors to illustrate overall similarities and differences in percent identities across the genomes. The number of best BLASTp hits in these eight Burkholderia genomes is also summarized at different percent identity cutoffs in Fig. 4.

Gene content comparisons
Although genomic islands (GIs) and their gene content vary greatly among B. pseudomallei strains, a thorough comparison of the GIs in the three genomes was already performed [7]. Therefore to investigate potential virulence and metabolism-related genes, we focused on gene clusters and individual CDS (not found in GIs) that were unique to strain MSHR668 and not present in the genomes of 1106a and K96243 (Tables 3 and 4). Many of the genomic features that were present in strain MSHR668 but absent in the genomes of 1106a and K96243 were also present in the genomes of one or more of the other Australian strains, for example strain MSHR305. This result is particularly interesting because of the similar clinical presentations of disease caused by these Australian strains, involving general septicemic infections and the somewhat rare events of encephalomyelitis caused by strains MSHR668 [6] and MSHR305.
Addressing differences in the gene content of MSHR668 compared to both K96243 and 1106a, Table 3 lists individual genes (CDS) that were present in the MSHR668 genome but not present in the genomes of both K96243 and 1106a. Table 4 compares the gene content of both K96243 and 1106a, listing CDS that were present in both K96243 and 1106a genomes but absent in MSHR668. Most of the individual genes listed in Table 4 have mobile-element related annotated functions.

Metabolic genes and chokepoint reactions
There were some metabolism-related genes in the MSHR668 genome that did not have putative homologs in the K96243 and 1106a genomes ( Table 5). The MSHR668 genome had thirteen genes with annotated functions in metabolism that were not present in the K96243 and 1106a genomes. Only four of the functions listed in Table 5 (cytidine/deoxycytidylate deaminase family protein, beta-glucosidase, putative dienelactone hydrolase, beta-lactamase) had additional copies in the MSHR668 genome. Only one gene (BURPS668_1621, encoding trans-aconitate 2-methyltransferase) was associated with a chokepoint reaction by the Pathway Tools [28]. This enzyme transfers one-carbon groups in the reaction that produces S-adenosyl-L-homocysteine from S-adenosyl-L-methionine [31].
Metabolic genes of interest in the K96243 and 1106a genomes that were not present in MSHR668 (Table 6) included D-glycero-D-manno-heptose 7-phosphate kinase, which is a candidate chokepoint enzyme, a LysM repeat protein and thioredoxin. The thioredoxin function was encoded by additional copies in both K96243 and 1106a genomes. D-glycero-D-manno-heptose 7-phosphate kinase is involved in the biosynthesis of lipopolysaccharide and is a virulence factor and potential protective antigen for B. pseudomallei [32].

Metabolic pathways
Metabolic pathways were identified in the three B. pseudomallei genomes by Pathway Tools [28] and compared using Pathway Tools, MetaCyc [31], KEGG [33], BLAST analysis [34] and IMG [27]. Pathways comprising central carbon metabolism and the main inputs and outputs are listed in S2 Table. All three of the genomes had components of the main pathways of central carbon metabolism and genes encoding transporters systems, anapleurotic reactions, and pathways for amino acid biosynthesis. The genomes of all three strains had complete pathways to make the amino acids and vitamins that humans obtain from diet and that more fastidious host-restricted intracellular pathogens, such as F. tularensis, do not contain. These included histidine, isoleucine, leucine, lysine, methionine, cysteine, phenylalanine, tyrosine, threonine, tryptophan, valine, folate, biotin, lipoic acid, pantothenate, thiamine, riboflavin and vitamin K2 (menaquinone). All of the genomes had genes encoding the cobalamin adenosyltransferase that converts cobalamin to vitamin B12. None of the genomes had genes encoding the enzymes needed to make vitamin K1 (phylloquinone).
Bacterial gene expression is controlled by transcriptional regulators, such as transcription factors and sigma factors. The functions of these proteins in gene expression regulation were first described in Escherichia coli [35] and were previously reviewed in Pseudomonas aeruginosa [36]. There were many transcription and sigma factors, response regulators, and DNA-binding proteins identified in the B. pseudomallei genomes ( Table 2). S3 and S4 Tables list the  differences in regulatory gene numbers, while Tables 5 and 6 compare gene content between MSHR668 and K96243/1106a. Nearly all of regulatory functions listed in these tables were present in additional copies in the genomes, although their exact gene targets are not known. Table 7 lists virulence genes compiled from online databases [37][38][39]and literature [19,40] with annotated metabolic and regulatory functions that were present in all three genomes. At least twenty five of the metabolic genes in Table 7 were identified as potential chokepoints by the Pathway Tools.

Discussion
Experimental infection of mouse models with the three B. pseudomallei strains showed that the K96243 and 1106a strains from Thailand had similar LD 50 values in both BALB/c (more susceptible) and C57BL/6 (more resistant) murine infection models, while the Australian strain MSHR668 was more virulent as measured by LD 50 . Given the incredible amount of genomic diversity among B. pseudomallei strains, we sought to identify candidate genomic differences that may correlate with variations in virulence. We conducted whole genome comparisons focusing on virulence, metabolism and regulation and identified genes in common among all three genomes. We also identified genes that were present in  MSHR668 but absent in K96243 and 1106a (and vice versa). Our findings and the implications on our understanding of melioidosis as a disease are discussed below.
Comparison of the three B. pseudomallei genomes revealed genomic differences that included the previously reported variability in GIs [7,19], which were likely acquired by horizontal transfer [19], as evidenced by their proximity to transposases, integrases, tRNA genes, and the presence of phage-related genes within the GI. This variability in the GI regions may contribute to virulence potential, particularly because these regions can encode a broad array of functions [20]. The intracellular life cycle and adaptation of a pathogen to the host cell environment depends on the expression of virulence factors, which is controlled by regulatory elements, and may be affected by the metabolic state of the pathogen [41]. The genomes of B. pseudomallei MSHR668, K96243, and 1106a contained complete gene sets for the core pathways comprising carbon metabolism. They also contained gene sets encoding transporters and utilization pathways for a wide range of carbon substrates, anapleurotic reactions and fatty acid degradation products (S2 Table), providing many potential targets for metabolic regulation.
An important outcome of the metabolic pathway analysis was identification of chokepoint reactions in the three genomes by the Pathway Tools software [28]. Inhibition of an enzyme that consumes a unique substrate might cause accumulation of the substrate and be potentially toxic to the cell. Conversely, inhibition of an enzyme that produces a unique product might result in starvation for that product, which could cripple essential cell functions. Thus, chokepoint enzymes may be essential to the pathogen and therefore represent potential drug targets. We identified two chokepoint reactions among the lists of genes in Tables 5 and 6, which were differentially present in the three genomes. Among the genes in Table 7, we identified twenty-five candidate chokepoint enzymes in common among the three genomes, involved in a variety of metabolic functions. The complete list of chokepoint reactions, including candidates, totals approximately 1,20021,300 reactions for each genome (data not shown) and requires additional curation and more extensive comparative analysis to determine which ones are the most promising targets. While our findings indicate that there are only a few metabolic differences among the B. pseudomallei genomes, it is becoming increasingly apparent that virulence and metabolism are linked together by complex regulatory interactions occurring between intracellular pathogens and their host cells [41][42][43][44][45]. We did find a few differences in regulatory gene content between MSHR668 and the other two genomes, in particular the K96243 and 1106a genomes contained more predicted sigma factor encoding genes than MSHR668 (S4 Table). Also the MSHR668 genome encoded additional transcriptional regulators, specifically two XRE family, two LysR family and one LuxR family, that were not present in the other genomes ( Table 3). The genomes of both 1106a and K96243 encoded one LysR family and one Cro/CI family regulator that were not present in the MSHR668 genome (Table 4). These results suggest that differences in regulation may contribute to the differences in virulence observed among these strains. Although further work is needed to test this hypothesis, the observed differences in transcriptional regulatory genes may contribute to the differential virulence observed in this study. Increasing evidence indicates that virulence gene expression is regulated by nutrients in the environment surrounding B. pseudomallei [46]. The expression of pathogen genes involved in transport and utilization of nutrients containing carbon and nitrogen is controlled by transcriptional regulators that are activated by the presence of nutrients [47][48][49][50][51]. For example, RpoS is involved in the response of B. pseudomallei to carbon starvation, heat shock, osmotic stress and oxidative stress. The expression of metabolic pathway genes involved in central carbon metabolism is controlled by RpoS, and by RpoS and BpsI co-regulation [52]. Therefore, the inter-regulation of stress response and metabolic genes by RpoS and BpsI may play an important role in B. pseudomallei survival and virulence [52]. RpoS has been reported to play a role in virulence gene expression in Salmonella typhimurium [53], and may also influence host macrophage responses to B. pseudomallei infection [54,55]. The polyamines spermidine and putrescine regulate gene expression at the transcriptional level by affecting regulatory protein binding to DNA. The Fur protein is a positive regulator of peroxidase and iron-containing superoxide dismutase expression, but in response to increased iron concentrations, Fur reduces the transcription of iron-regulated promoters [56].   Several studies have examined transcriptional profiles of B. pseudomallei during infection [46,[57][58][59][60][61]. Results of these efforts support the idea that some virulence functions leading to infection and disease are linked to pathogen metabolism through regulation of gene expression. In some cases, metabolic enzymes may act as virulence factors through their role in providing nutrients to the pathogen during infection. For instance, phosphoserine aminotransferase, encoded by serC, is involved in serine and pyridoxal-5-phosphate synthesis, and may be a virulence factor in B. pseudomallei, as it is co-expressed with other virulence genes and auxotrophic mutation attenuates virulence [59]. Several studies have shown that some genes involved in metabolic processes and virulence are upregulated in B. pseudomallei during infection, while other metabolic genes are downregulated [46,57,58,60,61]. In spite of the increasing evidence linking metabolism and virulence, further work is needed to thoroughly characterize the overlapping roles of virulence factors, regulators and metabolic pathways in B. pseudomallei pathogenicity. Comparative genomic approaches such as those described here can be a key first step in generating hypotheses with respect to the roles of various bacterial factors in virulence. Bacterial pathogens have evolved strategies to alter their lifestyle depending on whether they are in their natural environment or infecting a host, shifting resources from normal cell functions to the production of virulence factors, and altering metabolism to take advantage of the nutrients provided by host cells to facilitate survival and growth [62]. This should be especially true for B. pseudomallei given its presence and survival in a range of soil samples [63][64][65][66][67][68] and ability to cause severe disease in humans. Our comparison of the genomic features of two B. pseudomallei strains from Thailand (K96243 and 1106a) to one strain from Australia (MSHR668) revealed that the genomes are very similar in the repertoires of metabolic and virulence genes that they contain, leading to the conclusion that differential virulence studies on a larger scale are warranted. Detailed experiments will be necessary to characterize the relevance of specific genomic features to B. pseudomallei metabolism and virulence, and particular attention should be focused on the regulatory mechanisms influencing gene expression. Continued emphasis in this area will be critical to protection against melioidosis, as a better understanding of what constitutes a fully virulent Burkholderia isolate may inform better diagnostic and medical countermeasure strategies. The comparative genomic analysis that we present in this report, combined with more detailed functional analyses of metabolic networks, virulence and regulation, shows promise for examining the effects of B. pseudomallei and other intracellular pathogens on host cell metabolism and will lay a foundation for future prediction of the virulence of emerging strains.

Animal Challenges
Mouse challenges and statistical analyses were performed to establish LD 50 values for each strain of B. pseudomallei. The United States Army of Medical Research Institute of Infectious Diseases is compliant with all federal and Department of Defense regulations pertaining to the use of Select Agents. Cultures were initiated by inoculating Glycerol Tryptone broth (GTB-10g/L tryptone, 5 g/L NaCl, 40 ml/ L glycerol) with defrosted freezer stock of B. pseudomallei. The cultures were grown at 37˚C while shaking at 200 RPM for approximately 8-10 hours in order to harvest cells at late logarithmic phase of growth. Challenge doses were prepared according to OD 620 nm values and cultures were plated on sheep blood agar plates to confirm the number of colony forming units per milliliter (CFU/ml). At least 5 dose groups were used and 10 mice were included in each group. BALB/c and C57BL/6 mice were ordered from the National cancer Institute-NCI Frederick and were approximately 7-10 weeks of age at time of challenge. Mice were challenged intraperitoneally with bacterial doses suspended in 200 ml of GTB. Mice were observed at least daily for signs of illness or distress and monitoring frequency increased as indicated by the advancement of clinical signs. Challenged mice were observed at least twice daily for 60 days for clinical signs of illness. Humane endpoints were used during all studies, and mice were humanely euthanized when moribund according to an endpoint score sheet. Animals were scored on a scale of 0-12:0-25no clinical signs; 3-75clinical symptoms; increase monitoring; greater than or equal to 85distress; euthanize. Those animals receiving a score of 8 or greater were humanely euthanized by CO 2 exposure using compressed CO 2 gas followed by cervical dislocation. However, even with multiple checks per day, some animals died as a direct result of the infection. A Bayesian probit analysis was performed for each Burkholderia strain to estimate the lethal dose response curve. Prior distributions for each parameter were assumed to be independent, weakly informative Cauchy distributions with center 0 and scale 10. Using samples from the posterior distributions of the slope and intercept parameters from the probit analysis, the median and 95% credible intervals of the range of dose responses are estimated. Direct comparisons of the posterior samples of the LD50s of each strain permit us to make probabilistic statements about how likely it is that one strain is more or less potent than any other strain, given the observed data.

Genome Analysis
Whole genome sequences were obtained from NCBI (accession numbers NC_006351.1, NC_006350.1, NC_009074.1, NC_009075.1, NC_009076.1, NC_009078.1). To facilitate consistency in genome comparisons, genomes were annotated with RAST [26]. The GenBank format files for the RAST-annotated genomes are included in S1-S3 Files. The numbers of pseudogenes in each genome were obtained through the software package Psi Phi, which was kindly provided by Prof. Lerat [29]. In preparation for running Psi Phi, annotated protein sequences from each query genome were obtained from NCBI and used to query the nucleotide sequences of the other target genomes using tblastn. To identify potential pseudogenes, the Psi Phi software compares protein sequence matches in a query genome to the GenBank file of the target genome. We identified matches having a blast score with E-value ,10 210 and a minimal percentage of protein identity of 80% Matches with 80% to 100% protein sequence identity to the query protein were retained. If a query sequence had two matches in close proximity in the target genome (as might result from frameshifts or insertion), the matches were merged if they were ,300 nt apart [69].
Mobile genetic elements, transcription factors, sigma factors, response regulators, DNA binding proteins and two-component signal transduction systems were identified in each genome by searching the annotated genomes in the SEED [70]. Functional analysis was accomplished through the RegPrecise database [71].
Whole genome alignments were performed with Mauve [30]. To identify putative homologs among the genomes of B. pseudomallei strains K96243, 668 and 1106a, we performed a bidirectional best hits analysis, using BLASTp with an Evalue cutoff of 1e 25 to obtain liberal best hits for the proteins of each genome compared to the others. Genes x and y from genomes 1 and 2 are considered as homologs if y is the best BLASTp hit for x and vice versa. We used the blast2gi program from the Seals package [72] to format the BLAST results in tabular form. Each pair of genomes was subjected to this analysis. To obtain the CDS shared by all three genomes, the sequences in common to each pair of genomes were compared to generate a list of CDS present in all three genomes. Sequences unique to each genome were identified by comparison of the total number of CDS in each genome to the common sequences from each pairwise comparison. We gathered the sequences that were unique to MSHR668 and not found in either K96243 and 1106a, and those that were unique to both K96243 and 1106a but not found in MSHR668. These sets of sequences were compared to the originally annotated genomes from GenBank, to determine whether RAST annotation predicted similar CDS to the previously annotated genomes in GenBank. Predicted CDS were not included in the unique set if there were high identity hits (.95%) in the original annotation. The locus_tags in Tables 3-5 refer to CDS present in both the RAST and original annotations. To create heatmaps comparing CDS from each B. pseudomallei genome to other Burkholderias, we used protein BLAST version 2.2.26+ to compare B. pseudomallei K96243 protein translations against eight other Burkholderia proteomes that we also annotated using RAST. We disabled filtering and set the E-value cutoff to 1e 215 and then saved the best hit to each subject protein. The best hits were binned into groups based on percent identity (100%, 90-99.9%, 80-89.9%, etc) and then displayed as a heatmap (Fig. 3), which was created in R using complete linkage hierarchical clustering with euclidean distances. A matrix showing the numbers of CDS shared in each pairwise comparison and percent identity was created by counting the number of best hits in each bin (Fig. 4).
Virulence gene lists were compiled from [19,[37][38][39][40]73]. Blast analysis was used to compare the virulence gene sequences among the three genomes and between the original and RAST annotations. Metabolic pathways of the original and RASTannotated B. pseudomallei genomes were analyzed using the Pathway Tools version 18.0 [28]. Chokepoint reactions were identified in B. pseudomallei MSHR668, K96243 and 1106a using the chokepoint reaction finder, with human reactions excluded. Supporting Information S1