Genomic patterns of pathogen evolution revealed by comparison of Burkholderia pseudomallei, the causative agent of melioidosis, to avirulent Burkholderia thailandensis

Bt is closely related to Bp by 16S phylogeny

The Burkholderia genus is a large bacterial family containing >30 distinct species [16]. We estimated approximate divergence times between Bt and other species in the Burkholderia genus by comparing on a phylogenetic basis 16S rRNA gene sequences from four Burkholderia species and three outgroup species (Pseudomonas aeruginosa LMG1242T, Escherichia coli K12 and Salmonella typhimurium LT2). From this analysis, we found that Bt occupies a branch in the 16S rRNA phylogenetic tree that is highly related, but not identical, to B. pseudomallei (Bp) and B. mallei (supplementary figure one in Additional file 1). The close phylogenetic similarity between Bp and B. mallei is expected given that recent genome analysis suggests that B. mallei is in fact a derivative or clone of Bp [17]. A calibration of the phylogenetic tree against the E. coli and S. typhimurium divergence, estimated at 140 million years (Mya) [18], suggests that Bt diverged from Bp and B. mallei approximately 47 Mya ago. Similar time-frames were also obtained using groEL, another highly conserved housekeeping gene (Y. Yu, data not shown). This analysis suggests that although Bt is avirulent, it is likely to be highly evolutionarily related to virulent Bp, and a good candidate for comparative genomic analaysis. Notably, the Bp/Bt separation times are comparable to that used in other studies comparing pathogenic and non-pathogenic species [19, 20].

Genome sequencing of Bt ATCC700388 and comparison to Bt E264

Genomic DNA from Bt ATCC700388 was sequenced by a shotgun approach to 7x genome coverage, and the computer program ARACHNE was then employed to assemble the shotgun reads into contigs and scaffolds. The median contig length was 40 Kb (86 contigs in total), and 3 scaffolds could be assembled from the contigs, comprising 3.4, 2.9, and 0.35 Mb respectively. A subsequent comparison of these scaffolds to the finished Bt E264 sequence allowed the 1^st and 3^rd scaffold to be collapsed into a single assembly (Chromosome 1). The amount of genome sequence contained in contigs (sum of contig lengths) was 6.66 Mb (Chr 1 – 3.78 Mb, Chr 2 – 2.88 Mb), compared with a total genome length (sum of scaffold lengths, including estimated gaps between contigs) of 6.7 Mb (Chr 1 – 3.80 Mb, Chr 2 – 2.9 Mb). Thus, approximately 99% of the Bt ATCC700388 genome is contained in contigs, with 44 sequence gaps (29 in Chr 1, 15 in Chr 2), of which 15 gaps are larger than 500 bp (largest gap size being 9995 bp).

Although Bt ATCC700388 and Bt E264 represent the same bacterial strain, the specific isolates chosen for sequencing were stored at independent centers (DMERI and USAMRIID), and Bt ATCC700388 exhibits a slightly reduced growth rate compared to Bt E264 (Additional file 2). An analysis of the E264 and ATCC700388 genome sequences revealed that they were near identical except for : I) a chromosomal inversion of 2 Mbp on Chr 1, from position 12442442 (BTH_I1099) to 3328461 (BTH_I2895) based on the E264 genome sequence, II) 4 genes that are absent in ATCC70388, and III) 80 genes with putative protein sequence-altering DNA polymorphisms between the two strains (a full list of these genes is presented in Additional file 3). An initial resequencing analysis of genes in category III) revealed that several of the putative DNA 'polymorphisms' were sequencing artifacts, and thus an extensive analysis of genes in this category was not pursued further (Additional file 3). Of the four ATCC70388 absent genes, two genes (Bth_I1484 and Bth_I1485) were experimentally validated by PCR (Chua HH, data not shown) and encode components of a Type II Oligo-polysaccharide (OPS) synthesis gene cluster [21]. It is possible that the lack of these genes in ATCC700388 may contribute towards the observed growth differences by affecting Bt cell wall and membrane biogenesis. Although this hypothesis is supported by studies of the Bth_I1484 homolog in B. subtilis (TagO, 22), we emphasize that it is still highly speculative and needs to be further investigated. As the genomes of both strains (ATCC700388 and E264) were largely similar, we decided to adopt the finished E264 genome, after incorporating corrections for a few sequencing errors, as a reference Bt genome to compare against Bp.

Genome features of Bt

Similar to Bp, the Bt genome comprises two circular chromosomes of 3.8 and 2.9 Megabases (Figure 1). With a total of 5,645 predicted protein-encoding ORFs (3282 from Chromosome 1 and 2363 from Chromosome 2) (Table 1), the Bt genome ranks in the top twelve largest bacterial genomes amongst the 200-odd bacterial genomes sequenced at the time of writing this report (supplementary figure two in Additional file 1). A striking architectural feature of the Bt genome is the presence of at least fifteen regions exhibiting either atypical GC content or stretches of bacteriophage related genes and phage-like integrases (supplementary table one in Additional file 1). Collectively encompassing 4–5% of the entire genome, we refer to these regions as Bt-genomic islands (Bt-GIs), in keeping with similar regions found in Bp [15]. The Bt-GIs are not found in Bp and B. mallei, consistent with their being acquired by Bt subsequent to the Bp-Bt divergence. Supporting a model of recent acquisition, it has been found that some of these regions (eg. Bt-GI12) may still harbor functionally active bacteriophages [23]. Interestingly, of the ten Bt-GIs on Chromosome 1 and five on Chromosome 2, three Bt-GIs (Bt-GI 8, Bt-GI 9, and Bt-GI 10) occur in the same relative genomic location as GIs in the Bp genome. The presence of different horizontally acquired sequences in the same genomic location in related species suggests that these locations may represent genomic hot spots or landmarks for the acquisition of horizontally acquired sequences.

To experimentally show that the Bt-GIs are likely to represent regions of lateral transfer between different Bt strains, we then performed a series of PCR experiments to confirm the presence or absence of the fifteen Bt-GIs in 29 natural isolates of Bt, obtained from independent soil locations in either Northeast (n = 21) or Central Thailand (n = 8, 9 if E264 is also considered). We found that the distribution of Bt-GIs was indeed highly variable between different Bt strains (Figure 2). Specifically, while Bt-GIs 1, 3, 8 and 10 appear to be unique to Bt E264, GIs 2, 4, 5, 12, 13 were detected in 10–50% of strains and the rest of the GIs were detected in 50–80% of strains. Interestingly, GI 12 was present in 4 out of 9 strains isolated from Central Thailand (including E264), but absent from all Northeast Thailand strains. These results indicate that the Bt genome, like Bp, indeed contains a substantial fraction of genomic material that is highly variable in different strains. It is thus possible that genes on these islands may contribute to the presence of different strain-specific phenotypes.

Genomic synteny and large-scale differences between the Bp and Bt genomes

To assess large-scale genomic similarities between Bp and Bt, we aligned the two genomes and found that both chromosomes were highly syntenic between the two species (Figure 3). Specifically, we identified only four large-scale inversions involving genomic regions greater than 10 Kb; three in Chromosome 1 and one in Chromosome 2. Six out of the eight breakpoints flanking these inversions contain genes involved in DNA recombination such as transposases, phage integrases and recombinases (supplementary table two in Additional file 1). At this large-scale resolution, we did not observe any evidence for exchange of genetic material across the two chromosomes. The high degree of synteny between Bt and Bp is unexpected, given the fundamental plasticity of bacterial genomes, and because a previous analysis has demonstrated a much higher degree of genome rearrangement events between Bp and B. mallei, which is even more closely related [17]. The syntenic conservation of Bp and Bt is also considerably higher than that observed between other virulent and non-virulent bacteria species, despite the latter having smaller genome sizes and similar or shorter divergence times from that of Bt and Bp [5]. This unexpected result raises the possibility that that there may exist certain functional constraints acting as a negative selective forces to preserve the large-scale genomic organization of Bp and Bt.

Metabolism

Although Bp and Bt occupy similar ecological niches, previous biochemical analysis has identified several phenotypic differences between the species, including the ability of Bt, but not Bp, to assimilate the carbon sources arabinose and xylose [24]. Confirming a previous report [25], we found that Bt contains an eight-gene arabinose assimilation operon (BTH_II1626–1633) on Chromosome 2 that is absent in Bp. In Bp, this region has been replaced by a two protein cluster containing one hypothetical protein and one MarR family regulatory protein. A similar situation is likely to explain the differences in xylose metabolism as well. We found that Bt contains a 64 kb region on Chromosome 1 encoding several genes involved in xylose metabolism, including an aldose epimerase, xylose isomerase, xylose ABC transporter and an operon regulator (BTH_I2336–2339). This region, however, is absent in Bp, which instead carries GI 8 in the same location (Figure 4a). This result suggests an evolutionary model where the horizontal acquisition of GI 8 by Bp may have also resulted in deletion of the xylose gene cluster, demonstrating how horizontal transfer events can often result in simultaneous gene acquisition and loss. We also identified a thirteen-gene phosphonate gene cluster (BPSL2848–2860) on Chromosome 1 in Bp that was absent in Bt, consisting of ATP-binding (phnC, phnKL, phnN), phosphonate binding (phnD), metabolism (phnM, phnJ-G), and transport proteins (phnE), and a transcriptional regulator (phnF) (Figure 4b). In other species, gene products of the phn operon are required for the utilization of phosphonate as a phosphorus source [26], and are also involved in the uptake of phosphite, phosphate, and phosphate esters such as phosphoserine [27]. Since phosphorus availability is essential for many aspects of cell physiology, the lack of this operon in Bt suggests that other mechanisms must exist for acquiring phosphate in Bt. Supporting this hypothesis, the Bt genome contains several genes that could function in alternative pathways for phosphorus acquisition, including alkylphosphonate (BTH_I0257), phosphoglycerate (BTH_I3017–3022), and l-phosphonoalanine (BTH_II1908–1910) gene clusters.

Virulence

We attempted to identify genomic regions that might specifically contribute to the differential virulence of Bp and Bt. Variations in surface component proteins have been shown to contribute to virulence in several pathogenic species. For example, fimbriae are short pilus structures that allow bacteria to adhere to environmental surfaces and host cells, and Bp has been shown to be more efficient than Bt in adhering to and invading host cells [28]. We found that Bp contains twice as many fimbrial gene clusters (six) as Bt. Of the three Bp-specific fimbrial gene clusters (BPSL1626–1629, BPSL1799–1801, BPSS0120–0123), one cluster on Chromosome 2 is of particular interest as it contains a fimbrial usher protein, a fimbrial chaperone and two exported proteins with similarity to fimbriae proteins of Yersinia spp, whose family member Y. pestis is the causative agent of plague. This cluster appears to have been horizontally acquired by Bp as its GC content significantly deviates from the rest of the Bp genome, and replaces a flagellar operon (BTH_II0143–0197) in the homologous position in the Bt genome (Figure 5a). In addition to fimbriae, another important class of surface-related proteins that could contribute to virulence are surface polysaccharides. It has been previously shown by other molecular techniques that Bp contains a large gene cluster (BPSL2790–2810) involved in the synthesis and export of capsular polysaccharides, a major determinant of virulence, and that this cluster is absent in Bt [29]. We confirmed this finding and further discovered that the specific location of the Bp capsule gene cluster within the Bp genome is likely to be non-random, as it replaces a pre-existing 10-gene cluster in Bt already dedicated towards the metabolism and processing of polysaccharide structures (Figure 5b). Specifically, the original cluster in Bt contains several genes involved in polysaccharide assembly, such as GDP-6-deoxy-D-lyxo-4-hexulose reductase, a GDP-mannose 4, 6-dehydratase, a methyltransferase and four glycosyltransferases, which are all involved in GDP-D-rhamnose biosynthesis in Pseudomonas aeruginosa [30]. This result expands our depth of understanding surrounding this important acquisition event, and raises the intriguing possibility that a key event in the pathogenic evolution of Bp was to replace a pre-existing or ancestral polysaccharide coat with an alternative pathogenic variant capable of resisting challenges by the immune system of an infected host. In addition to surface molecules, another broad class of molecules that might contribute to virulence is secondary metabolites with toxin activities [31]. There are four polyketide synthase (PKS) and nonribosomal pepeide synthase (NRPS) clusters involved in the production and regulation of secondary metabolites that are present in Bp but absent in Bt, containing genes encoding putative resistance proteins (BPSS1011), efflux proteins (BPSS1268) and lipopeptide antibiotic proteins (BPSS1631–1634), which could be involved in toxin production. Other differences, including the previously reported absence of a type III secretion system (TTS1) in Bt [32], are listed in the supplementary table three and table four (see Additional file 1).

The Bp and Bt genome share similar protein family distributions

The SUPERFAMILY hidden Markov model library [33] is a database representing all proteins of known structure derived using Structural Classification of Proteins (SCOP) defined structural folds. We used SUPERFAMILY analysis to classify and compare evolutionarily related groups of domains between the two organisms. We queried the predicted proteomes of Bp and Bt against the SUPERFAMILY HMM library and obtained structural (and hence implied functional) assignments to protein sequences at the superfamily level for both species. This exercise produced 3706 Bt sequences with 721 unique SCOP superfamilies and 3645 Bp sequences with 705 unique SCOP superfamilies. The top 20 structural families assigned to Bp and Bt protein sequences are provided in supplementary table five (see Additional file 1). We found the distribution of different structural family assignments highly similar between Bp and Bt, supporting our basic conclusion that genomes of Bt and Bp are highly conserved. The most populated superfamily in both organisms were P-loop containing nucleoside triphosphate hydrolases (~7.7%), with a local sequence motif ([AG]-x(4)-G-K-[ST]) for ATP/GTP binding. The second and third most distributed superfamilies are winged helix DNA-binding domain (~6%) and NAD(P)-binding Rossmann-fold domains (~5%), which are involved in DNA-binding and NAD(P)-binding respectively.

Virulence genes in Bp and Bt exhibit increased diversity

Besides the acquisition of large-scale genomic material, alterations in the amino acid composition of protein homologs could also play an important role in the phenotypic differences of Bp and Bt, by altering protein activities and specificities. To perform a systematic analysis of the Bp and Bt proteomes, we identified and compared 4630 orthologous protein pairs (2826 from Chromosome 1 and 1804 from Chromosome 2) between the two genomes, and mapped them to them to their respective metabolic pathways. We found strong conservation between the Bp and Bt proteomes in core metabolic pathways such as amino acid metabolism, cofactor and carrier synthesis, nucleotide and protein biosynthesis, consistent with the ability of Bp and Bt to occupy similar environmental niches [24]. Based upon KEGG annotations, only five out of 1997 genes involved in core metabolic pathways in Bp do not have a clear homolog in Bt, including a putative ATP-binding protein (BPSL2860), a putative acetyltransferase (BPSL1417), a putative RNA 2'-phosphotransferase (BPSL0762) and two ABC transporters (BPSL1824, BPSL2849). This result suggests that the central metabolic machineries utilized by Bp and Bt are likely to be highly similar. Unexpectedly, the Bp and Bt proteomes also appear to share significant similarities in their virulence components as well. Specifically, of 368 known and potential virulence genes in Bp [15], 275 orthologs (71%) are present in Bt at an average similarity of greater than 80% (see supplementary table four in Additional file 1), including two type three secretion systems (TTS2 and TTS3), antibiotic resistance genes, type IV pili-generating proteins, hemolysin-related genes, and several adhesion factors and proteases. Of the remaining 93 Bp-specific virulence-related genes, 20 are located in Bp GIs or regions of atypical GC content, and 73 within the core Bp genome.

We then investigated if specific components of the Bp and Bt genomes might exhibit significantly greater or lesser levels of diversity, by computing synonymous and non-synonymous nucleotide substitution rates (Ka and Ks) for orthologous protein pairs. The average Ka/Ks ratio for Bp and Bt orthologs was 0.069 (0 – 1.14) for Chromosome 1 and 0.099 (0 – 1.20) for Chromosome 2, where low ratios suggest decreased diversity (p < 0.001 for Chr1 vs Chr2). The significantly lower Ka/Ks ratio for Chromosome 1 is consistent with this chromosome encoding many proteins involved in core biological functions, including energy production and conversion genes, genetic information processing, cell cycle control. Supporting this hypothesis, the Ka/Ks ratio of core metabolism genes (0.05) from both chromosomes was significantly lower than the genome average (p < 0.001). In contrast, the average Ka/Ks ratio of ortholgous virulence genes (0.099) was significantly higher than the genome average (p < 0.001, Figure 6a), suggesting that these virulence-associated genes are likely to exhibit increased diversity between Bp and Bt compared to other genes. Intriguingly, genes associated with the Type III secretion cluster TTS3, including BPSS1547 (prgH), BPSS1549 (mxiI/prgJ), BPSS1551 and BPSS1552 (orgA and B), exhibited the highest level of divergence among the virulence related genes between Bp and Bt. In other species, prgH and prgJ proteins are integral members of the Type III secretion needle complex, which establishes a specific physical interface between the pathogen and host [34]. It is tempting to speculate that subtle functional alterations in these specific TTS3 genes could have contributed towards the ability of Bp to successfully infect human and other animal hosts. In addition to the TTS3 members, four genes of unknown cellular function (one from chromosome 1 and three from chromosome 2) also exhibited Ka/Ks ratios greater than one, suggesting that they might be rapidly evolving. However, no predictable function or identified domains could be assigned to these four genes, thus not allowing speculation on their roles in pathogenesis.

Functional biases in Bp and Bt pseudogenes

Besides point differences in amino acid composition, another evolutionary mechanism for generating protein diversity is the introduction of point or frameshift mutations to generate pseudogenes that are subsequently eroded and eliminated by small deletions events [35]. We used the Psi-Phi program to identify pseudogenes in Bp and Bt [36], and found 170 (129 on Chr1 and 41 on Chr2) pseudogenes in Bp and 179 pseudogenes (82 on Chr1 and 97 on Chr2) in Bt. Of the 170 Bp pseudogenes, 44 were previously identified in the original Bp genome annotation [15], which also reported an additional 81 other pseudogenes. These 81 genes lack clear homologs in Bt, and were thus not detected by Psi-Phi, which requires a comparative genomic metric for pseudogene identification. A complete list of these genes is provided in the Additional file 1 (supplementary table six). To better understand the molecular events responsible for these pseudogenes, we grouped the Bp and Bt pseudogenes by their functional inactivation categories, ranging from nonsense mutations, frameshifts, IS insertions, to other insertions or deletions, and truncations (Figure 7). In this analysis, "frameshifts" were classified as very small insertions or deletions (typically only 1 or 2 nt in length), "other insertions" were >10 bp in length, and "truncations" included large deletions at either or both the beginning and end of an ORF. Interestingly, when we considered all recognized pseudogenes including those previously annotated and those newly identified, the Bp and Bt genomes were found to display remarkably different mutational spectra. Specifically, internal stop codons due to nonsense mutations were the most usual disruption (~74%) in the Bt genome, followed by truncations (~21%) and frameshifts (~3%). In contrast, a large proportion of the Bp pseudogenes resulted from truncations (~51%), followed by nonsense mutations (~26%) and frameshifts (~17%). In both species, only a few pseudogenes corresponded to IS element insertions.

Around one third of the pseudogenes in both species corresponded to hypothetical proteins and/or proteins of unknown function. However, among the annotated pseudogenes in the Bt genome, there was a distinct non-random functional trend for proteins exhibiting such inactivation events, as they largely involved membrane-associated and exported proteins (26 Bt pseugodgenes, p < 0.001 by Fisher's Test). This bias may reflect the results of strong selective pressures in the natural environments directly contacted by these bacteria. Intriguingly, in the Bp genome, one frequently involved class of proteins that were subject to inactivation events were genes involved in transcription (21 from Bp; p < 0.001), which could lead to differences in activating upstream stimuli and the selection of target genes. This result suggests that certain aspects of Bp-specific behavior may result from alterations in gene expression, in addition to the other large-scale and protein differences previously discussed.

Comparison of Bp and Bt intergenic regions

Finally, our finding that a number of transcriptional regulators might be targeted for inactivation in Bp made us also investigate the extent to which the transcriptional regulation machinery in Bp and Bt might be similar or different. As a surrogate for gene expression information, we compared intergenic regions of the orthologous protein pairs, which are likely to contain important cis-acting promoter elements and motifs for transcriptional regulator docking and activity. Using a working criterion that an intergenic region should be greater than 30 bp and lie between two distinct ORFs, we identified 1191 (~58%) orthologous intergenic regions in Chromosome 1 and 634 (~49%) in Chromosome 2 for further analysis. When assessed at a global level, the intergenic regions of Bp and Bt appear to be similar, in terms of average sequence identity (93% for integenic pairs in both chromosomes), mean lengths, and size ranges (p > 0.05). The high levels of similarity between the intergenic regions could mean that either the transcriptional regulatory machineries of Bp and Bt are highly similar, or that insufficient time has passed for conserved promoter motifs to emerge against the background mutation rate. To test the genomes for any evidence of increased conservation in the intergenic regions, we asked if there was any relationship between the extent of protein similarities between homologs and the degree of sequence conservation in their 5' intergenic regions. Comparing the orthologous gene pairs, we identified a weak but significant positive correlation between protein similarity and adjacent 5' intergenic sequences for genes associated with core metabolism (r = 0.134, P < 0.01, Spearman's correlation coefficient). This bias became stronger particularly for Chromosome 1 when it was treated as separate from Chromosome 2 (r = 0.147, P < 0.0001 for Chromosome 1; r = 0.084, P < 0.05 for Chromosome 2, Spearman's correlation coefficient). In addition, the percentage similarity for promoters of core genes was significantly higher compared to the genome average (p < 0.001, Figure 6b). The finding that orthologs involved in core metabolism already tend to have more conserved 5' intergenic region suggests that some degree of selection in the integenic regions has already taken place. In contrast, for the virulence genes, such a relationship was not observed (p > 0.05). Thus, for genes that are present in both species, the high levels of similarity between their intergenic regions suggest that the cis-acting promoter elements of Bp and Bt, are likely to be broadly similar.

Phylogenetic analysis

16S ribosomal sequences for Burkholderia (B. pseudomallei K96243, B. thailandensis E264, B. mallei ATCC23343 and B. cepacia J2315), and other spp. (Pseudomonas. aeruginosa LMG1242T, Escherichia coli K12 and Salmonella typhimurium LT2) were aligned and compared using MEGA version 3.0 software [47]. Sequence divergence rates were calculated using a neighbor joining algorithm and a bootstrap value of 2000. P. aeruginosa LMG1242T was used as an outlier subgroup.

Bacterial strains and genome sequencing

Two independent cultures of B. thailandensis E264, derived from the same original isolate but stored at independent centers (DMERI and USAMRIID), were processed for genomic DNA extraction and whole-genome shotgun (WGS) sequencing at separate facilities (The Institute for Genomic Research (TIGR) and the Broad Institute (BI)). At TIGR, sequencing and assembly followed by genome closure was performed as described in Nierman et al. [17]. Identification of B. thailandensis coding sequences (CDs) was performed using GLIMMER (48) with modifications described in Nierman et al. [17]. The BI sequencing approach is highly similar and described in the Main Text and Additional file 4. Our comparative genomic analysis is primarily based on the TIGR E264 genome after incorporating corrections for specific sequencing errors. The TIGR Bt E264 genome sequence has been assigned GenBank accession no. CP000086 and CP0000865. The Whole Genome Shotgun project described in this paper has been deposited at DDBJ/EMBL/GenBank under the project accession AACX00000000. The version described in this paper is the first version, AACX01000000.

Genome alignments and comparative genomics

The genome sequence of virulent B. pseudomallei K96243 was obtained from NCBI [15]. Syntenic regions between the B. thailandensis (Bt) and B. pseudomallei (Bp) genomes were identified and aligned using MUMMER [49], and visualized using the ARGO Genome Browser [50] (Broad Institute, MIT). Protein homologs between the two species were identified using BLASTP at an E-value threshold of < 10^-10 with subsequent confirmation by reciprocal blast. The BLASTP output was also manually curated to generate the final list of orthologous gene pairs between the two species. Pathway comparisons were performed by interrogating the Bp and Bt genomes against the KEGG metabolic database. Fisher's exact test was used to calculate the difference between proportions of functional enrichment. Rates of nucleotide substitution in CDs were calculated by aligning the nucleotide sequences of orthologous gene pairs by Clustalw [51] and determining the number of synonymous nucleotide substitutions per synonymous site (Ks) and the number of nonsynonymous nucleotide substitutions per nonsynonymous site (Ka) by a maximum-likelihood method [52] in PAML [53]. Spearman's rank correlation analysis was performed using GraphPad Prism version 4.0 (GraphPad Inc., Calif.). Intergenic sequences of orthologous gene pairs, defined as the non-coding sequence between the translational start or stop of two successive genes, that were greater than 30 nucleotides between two species were compared by BLASTN. Two-sample t-tests (for normal distributions) and Mann-Whitney tests (for non-normal distributions) were used to determine the difference between means of different populations.

Experimental validation of genomic islands

The presence of absence of 15 GIs was assessed across 29 Bt isolates using PCR. The 29 isolates were collected from independent soil locations in two major areas: Northeast Thailand (21 isolates) and Central Thailand (8 isolates). One target gene and one flanking gene were selected for each genomic island based on Bt E264 sequence data. The target genes chosen were predominantly hypothetical conserved proteins; and genes such as integrases and bacteriophage were avoided. The flanking genes were the genes immediately abutting each putative island. The PCR primers and cycling conditions (supplementary table seven in Additional file 1) for the target and flanking genes were optimized using Bt E264 genomic DNA. PCR amplifications were performed with a PTC-0220 DNA engine (MJ Research, Cambridge, MA) with Platinum Taq polymerase (Invitrogen), and aliquots of reaction mixtures were analyzed by agarose gel electrophoresis.

SUPERFAMILY assignment

SUPERFAMILY (release 1.69) is a library of HMM models based on SCOP and contains 1539 known structural superfamilies [33]. The HMM based assignment tool provides structural assignments to protein sequences at the super-family level using known structural information. SUPERFAMILY was run to assign structural superfamilies and domains for Bt and Bp protein sequences in HMMER mode.

Identifying pseudogenes in Bp and Bt

To identify pseudogenes, we first used protein sequences of Bp (or Bt) genome to query the nucleotide sequence of the Bt (or Bp) genome using TBLASTN. We then applied the Psi-Phi program suite [36] on the BLAST outputs to recover candidate pseudogenes in each genome. Briefly, Psi-Phi is designed to identify pseudogenes by a comparative analysis of related genomes, and retrieves pseudogenes resulting from nonsense mutations, frameshifts generated by small insertions or deletions, large insertions (such as those resulting from transposable elements), and truncations of any specified length as well as any incorrectly annotated spacers resulting from gene degradation [36]. Since Bp and Bt are closely related, a stringent cutoff of E-value < 10^-15 and a minimal percentage of protein identity of 79% was used to identity pseudogenes. The candidate pseudogenes were manually curated, and the disrupting mutations determined by aligning the nucleotide sequences of putative pseudogenes with their functional counterparts using CLUSTALW 1.8 [51]. The IS elements were searched by blasting genome sequences against the IS nucleotide database [54].