Characterization and Inference of Gene Gain/Loss along Burkholderia Evolutionary History

Collection of genomic sequences

Complete genome sequences from 60 Burkholderia strains (representing 17 species) and five strains of other Burkholderiales species (Supplementary Table S1) were obtained from National Center for Biotechnology Information database (ftp://ftp.ncbi.nih.gov).

Clusters of orthologous genes in Burkholderia

A total of 416,696 open reading frames from the chromosomal regions of all the 65 strains were selected for homology searches. OrthoMCL was used to generate groups of orthologous proteins with default parameters. ¹⁴ This resulted in a set of 52,031 Burkholdria clusters of orthologous genes (BuCOGs; Supplementary Table S2 and Supplementary Table S3), of which 416 were present as single copy genes within the 65 genomes. These single copy genes were selected to reconstruct the organism tree.

Each of the 416 orthologs was concatenated and aligned with MUSCLE using default parameters. ¹⁵ Poorly aligned regions from this alignment were removed by trimAl V1.2. ¹⁶ Phylogenetic trees were inferred by maximum likelihood with PHYML ¹⁷ by 1,000 bootstrap replicates, using JTT substitution model ¹⁸ as the model of amino acid evolution with a gamma distribution with 8 categories. Amino acid usage and the shape parameter were estimated from the data.

Probabilistic model for gene gain and loss

A binary character (0,1) table of gene absence and presence across the studied strains was used to estimate the phyletic patterns. The patterns were analyzed by probabilistic models, using the likelihood framework described by Csürös et al^19,20 to study evolution of introns and gene presence/absence. Compared to the maximum parsimony method used in Prochlorococcus, ²¹ the corresponding probabilistic model uses branch-specific gene gain and loss rate, with different species showing dramatic HGT frequency in the real evolution history. ²²

We used a rates-across-genes Markov model to trait gene evolution, with branch-specific gain and loss rates. Briefly, the procedure was as follows: it was assumed that gene sites evolved independently under a Markov model. ²³ The gene state (encoded by 0 and 1 for absence and presence) changed on each branch b of the phylogeny according to our modified probabilities formulas,

p 0 → 0 ( b ) = ( μ + λ e − ( λ + μ ) ) / ( λ + μ ) p 0 → 1 ( b ) = λ ( 1 + 1 e − ( λ + μ ) ) / ( λ + μ ) p 1 → 0 ( b ) = μ ( 1 − 1 e − ( λ + μ ) ) / ( λ + μ ) p 1 → 1 ( b ) = ( λ − μ e − ( λ + μ ) ) / ( λ + μ )

where λ denotes branch-specific gene gain rate, µ denotes branch-specific loss rate. The gene density at an ancestral node was computed as an expected value conditioned on the observed data, by summing posterior probabilities. ¹⁹ The values of estimated gene gain and loss rate can be affected by several reasons, such as different standard of gene identification in each genome, and different reconstructed methods. In order to overcome the biases of gene numbers of different genome and different gene families, we bootstrapped the gene table, and generated 1,000 random gene tables by selecting genes of the original gene table independently. For each bootstrap table, the likelihood and evolutionary history are calculated. This procedure yields confidence levels for the estimates inferred in this study of about 1%. The mean values of estimate results were considered as final value. We also used the maximum parsimony method to estimate the gene numbers in each branch, and got the same tendency though the gained gene number in recent branches is slightly higher than with our methods.

The individual gain and loss history of gene groups was also analyzed. The posterior probability of gene presence and absence in each node was calculated according to the branch-specific gene gain/loss rate. Only when the gene presence probability was greater than 0.5, it would be considered as being present in the evolutionary node.

Each gene that was predicted to be transferred was compared with the NCBI nr database using BlastP with the expect value cutoff set at 10^–10 to identify homologs in other organisms. For each protein, the first 100 hits were chosen for further potential HGTs analysis. The complete protein sequences of the 100 hits were extracted from the GenBank database, followed by multiple sequences alignment and phylogenetic tree generation by the same criterion described above. Based on tree topology, we identified the donor lineages, which naturally grouped with the Burkholderia genes. The candidate donor genus should have more than three genes and have strong bootstrap support values (>80).

COG classification and gene ontology analysis

The functional annotations of the Burkholderia genes were retrieved from the Burkholderia database ²⁴ except for B. phytofirmans PsJN and Burkholderia glumae BGR1, which were annotated by the NCBI COG database. ²⁵ Gene Ontology annotations are described at Pathema. ²⁶ In order to determine the overrepresented GO terms, we used the GOstats, ²⁷ taking P-values <le^–5 for significant.

Genomic islands (GEIs) in each evolutionary node

In this study, we first used a Markov model to predict the gained genes in 60 Burkholderia strains, which is one of the bacterial genera whose genomes have been most sequenced. Then, the transferred genes from different evolutionary nodes were mapped into the chromosome to investigate the existed GEIs. As we have the information about which genes are transferred, the criterion about GEIs was compared to the supplementary file in B. xenovorans LB400 genome analysis, ⁷ for convenience, we defined minimum continuous four gained genes as a GEI.

Coexpression analysis

A Pearson correlation coefficient was calculated based on 47 arrays of B. pseudomallei K96243 downloaded from NCBI GEO GSE5495 dataset ²⁸ to investigate the function of transferred genes in different evolutionary time (different evolutionary node in phylogenetic tree). Background correction and data normalization was done by the RMA algorithm in Bioconductor ²⁹ and the poor annotated probe-sets were removed. Measurements for unique genes was calculated from means of the probes belonging to the same gene. The gene coexpressed with at least one of either gene in GEIs with the absolute correlation |r| > 0.6, which we expect to be co-regulated in our case, was further selected to analyze the function.

Mapping coexpressed genes onto existing biological pathways

For each gene in GEIs gained in the same node, the top 100 genes in each gained GEIs were selected. These genes were sorted by absolute correlation value |r| referred to as coexpressed genes in descending order. These coexpressed genes with GEIs genes in each node were then mapped onto the pathways from Biocyc. ³⁰ A score of each pathway was calculated by using following formula, which represent the extent of coexpression in each node:

S_pathway = ∑ i G ∑ j N ( 1 − R i j / 100 ) * r i j / G * N , where G is the number of genes in a pathway present in the top 100 of all same node GEIs gene lists, and N is the number of lists in which gene i is present in the top 100 genes. R _ij is the rank of the i^th gene of the specified pathway in a list, r _ij is the absolute correlation value |r|.

Reconstruction of gene gain and loss in the evolution of Burkholderia

Gene gain and loss are ongoing processes in microbial genomes, resulting in diversity in genome sizes, even among closely related strains. We used the probabilistic model (details in Methods) to reconstruct the events that occurred during the evolution of this genus after its divergence from the common ancestor (Fig. 1). By mapping the gene distributions in a phylogenetic context, we could infer the ancestry and dynamics of genes among strains. The reconstruction shows that the gene gain and loss seems to be a prevailing trend in the evolution of Burkholderia. In the Burkholderia ancestor and the ancestral branches of the major clades (the clade definition was presented in Supplementary Table S1) a huge number of gene acquisitions have occurred, whereas gene loss seems to be occurring more frequently in evolution of closely related strains, generating large differences in a short divergence time.

OPEN IN VIEWER

Figure 1

Phylogenese tree of Burkholderiaceae based on concatenated amino acid alignments of one homolog gene and rooted by using proteins from other Burkholderiaceae species as the outgroup. All the evolutionary nodes in this tree were supported by bootstrap values higher than 90%. Reconstruction of gene content evolution based on BuCOGs in Burkholderia is shown. For each species and each internal node of a tree, the inferred number of BuCOGs present (numbers in black), and the numbers of BuCOGs lost (numbers in blue) or gained (numbers in red) along the branch leading to a given node (species) are indicated.

The common ancestor of Burkholderia had an estimated 5,176 genes, losing 322 genes and acquiring 1,335 genes compared with the outgroup of other Burkholderiaceae. Thus, gain of ancestral genes seems to be the prevailing trend in the evolution of Burkholderia (Fig. 1). This result is different from former similar analysis in other bacteria. ³¹ Most of the changes mapped to this ancestral stage of evolution seem to be linked to metabolism. Abundant genes for membrane transport and metabolism were gained via HGT (see Table 1 for the GO over-represented categories). The subsequent evolution of the Burkholderia reveals a considerable number of strain-specific duplications and acquisitions of unique genes as well as ancestral node gene loss (Fig. 1). In order to confirm the ancestral gene aquisitions, we randomly selected 100 ancestral transfer gene families and manually checked the phylogenetic trees. 85 gene families were confirmed nested within other groups, while 8 gene families without enough homologs in other species and 17 other cannot be determined from the trees. For instance, the whole fabG (BuCOG8606) gene family in Burkholderia is clustered with other unrelated bacteria, and the copC (BuCOG142) gene family might have transferred from Enterobacteriaceae (Supplementary Figs. S1 and S2).

OPEN IN VIEWER

Table 1

Over-represented GO categories among BuCOGs gained in the ancestor of the Burkholderia.

GO ID	GO term	Count	Total	P-value	GO level
GO:0051234	Establishment of localization	124	726	1.64e-31	3
GO:0006810	Transport	106	618	3.30e-27	3
GO:0016020	Membrane	84	571	1.38e-15	3
GO:0022891	Substrate-specific transmembrane transporter activity	37	179	2.83e-12	4
GO:0022892	Substrate-specific transporter activity	42	237	8.03e-11	3
GO:0044425	Membrane part	66	497	1.13e-9	3
GO:0022804	Active transmembrane transporter activity	36	200	1.75e-9	4
GO:0044238	Primary metabolic process	145	1521	1.35e-8	3
GO:0005975	Carbohydrate metabolic process	44	290	1.69e-8	4
GO:0044237	Cellular metabolic process	151	1661	1.62e-7	3
GO:0044459	Plasma membrane part	45	327	4.22e-7	4
GO:0005886	Plasma membrane	45	329	4.75e-7	4
GO:0031224	Intrinsic to membrane	47	364	1.92e-6	4
GO:0016746	Transferase activity, transferring acyl groups	23	130	8.15e-6	4

Moore et al ³² have reported that B. thailandensis, a species closely related to of B. pseudomallei, is unable to cause disease. In comparison with B. pseudomallei, the attenuation of B. thailandensis is caused by a functional arabinose biosynthesis operon which is largely deleted in B. pseudomallei. ³² In Burkholderia, the virulence divergence between different strains can be caused by mutation or loss of few genes. ³³ The concerted evolution trends in Burkholderia could be explained by the large-scale ancestral gene acquisition that might facilitate the acquisition of a generally pathogenic lifestyle of Burkholderia, and the strain specific pathogenic lifestyles can be modified by fine-tuning gene loss.

Overall, our results suggested that about 60% of the genes families in Burkholderia genus are inferred to have experienced HGT at least once during their evolution. As a large proportion of genes in Burkholderia are Burkholderia-specific, they could represent new birth of genes or caused by in silico prediction errors. Hence, the HGT percentage might be overestimated and the real HGT genes number in each Burkholderia strains could be less than what we predicted. The candidates for these genes outside the Burkholderia were divided according to taxonomic classification (Fig. 2). With more than 7,000 candidate donors in Proteobacteria, this represents the second largest gene reservoir for Burkholderia. Along with other Betaproteobacteria, there is a clear overrepresentation of Gammaproteobacteria and Alphaproteobacteria (Fig. 2). Most of the Burkholderia species are soil dwellers, and it has shown that the Burkholderia species are ancient symbionts of legumes, ³⁴ which created the opportunity to gain new genes from eukaryotes. ³⁵ Furthermore, more than two third of the gained genes were Burkholderia-only, which means these genes have no homolog outside this genus. This phenomenon may be attributed to the extinction of the donors followed by ancient gene transfer ³⁶ or gene transfer within Burkholderia lineage followed by new gene birth.

OPEN IN VIEWER

Figure 2

Taxonomic distribution of candidate donors for putatively transferred Burkholderia genes. The taxonomy is sorted by the percentages representing numbers of hits relative to the BuCOGs total number (white bars). The percentages relative to the number of phylum-specific proteins in Genbank are shown as black bars. No matches were detected for Dictyoglomi, Elusimicrobia, Synergistetes and Tenericutes (data from 2, 2, 1, and 26 genomes, respectively). Data for Caldiserica, Chrysiogenetes were not available at the time of analysis.

To further understand the evolutionary history of Burkholderia, we selected the representative genome of B. pseudomallei K26943 to elucidate impact of gene gain and loss evolution on different evolutionary nodes. The GEIs of this strain in each evolutionary node is represented in Supplementary Fig. S3. This figure indicates that most of the GEIs in this strain were gained in the ancestor node while only a small fraction of GEIs were gained most recently (Supplementary Fig. S3). We also found that the percentage about transferred chromosome fragments were different between two chromosomes in B. pseudomallei K26943 with 15.58% in the large and 26.88% in the small chromosome, respectively. The result strongly suggests that the small chromosome was much more prone to HGT than the large chromosome in the evolutionary history. It has been reported that the small chromosomes in Burkholderia carry more accessory functions associated with survival and adaptation in different niches, whereas the large chromosomes encode most of the core genes (shared by all strains) involved in central metabolism and cell growth. ⁹ As GEIs play key roles in Burkholderia pathogenicity, our results indicate that the ancient gene transfer played key roles in pathogenicity of Burkholderia.

GEIs are discrete DNA segments which contain a number of genes acquired by HGT that might be beneficial to hosts under certain conditions. ³⁷ Two major methods which can be broadly grouped into sequence-based methods and comparative genomics/phylogeny based methods, have been developed to predict the GEIs computationally.^38–40 The first method can only predict recently acquired GEIs while the prediction accuracy of the latter method is based on the number of related species. In this study, 109 GEIs were found in B. pseudomallei K26943 (Supplementary Table S4). The results indicate that evolutionary history reconstruction is a powerful method to investigate the GEIs, especially if based on a large enough number of closely related species.

Comparison of the number of genes gained or lost on a particular evolutionary scale and the length of the corresponding branch revealed a pattern different from that described for Prochlorococcus (cyanobacteria) and lactic acid bacteria. ²¹ Similar to the Prochlorococcus, the number of gene loss is significantly correlated with the branch length (r = 0.576, P < 8 × 10^–5). ²¹ However, the number of gene gains also has such a correlation (r = 0.500, P < 0.0009). The clock-like behaviour of gene loss and gain could be explained by a large number of small-scale events, which might be randomly distributed along the evolutionary path. However, the gene acquisition and removal model points to a mixed pattern, combining disruption events and complete acquisition/removal. Apparently, multiple type III secretion systems (T3SS) in Burkholderia were acquired from HGT and most of these genes were located as clusters.^8–10,41 The acquisition stage analysis revealed that these genes were gained in the same period, indicating that these T3SS might be integrated through the entire fragment being transferred. Such GEIs probably played a crucial role and could facilitate the evolution by ‘quantum leaps’ as their gain or loss could rapidly and dramatically alter the genome and lifestyle of a bacterium. ³⁷ The transfers of large fragments were frequently associated with particular physiologic adaptation such as virulence, catabolism, or resistance to a toxic compound. ³⁷ Mutated pseudogenes likely accounted for B. mallei being nonmotile or nonflagellated. ³⁴ The presence of numerous insertion sequence elements mediated extensive deletions and rearrangements of the genome. These deletions likely accounted for the large gene variation among closely related strains. For example, comparing the B. mallei SAVP1 and B. mallei ATCC 23344 genomes revealed a large T3SS-encoding fragment lost in B. mallei SAVP1 (Fig. 3). This loss likely accounted for the difference in virulence between B. mallei SAVP1 and 5. mallei ATCC 23344. ⁸ By comparing B. mallei ATCC 23344 with B. pseudomallei K96243, the same event of large T3SS-encoding fragment loss was also observed (Fig. 3).

OPEN IN VIEWER

Figure 3

Profiling of the B. mallei SAVP1 and ATCC 23344 lost genes and regions. (A) Graphs represent the best protein hits scatter plot (black dots) and B. mallei SAVP1 chromosome 2 recent lost genes (blue dots under x axis) in B. mallei ATCC 23344 chromosome 2. (B) Graphs represent the best protein hits scatter plot (black dots) and B. mallei ATCC 23344 chromosome 2 B. mallei lost genes (blue dots under x axis) in B. pseudomallei K96243 chromosome 2. Black scatter plot represents the best protein hits of two genomes chromosome according to the Pathema Scatter Plot Results. The lost gene homologues were plotted onto the another genome's chromosome under the bar with blue dots. The shaded regions represent the lost regions. The black rectangle region represents the T3SS location.

Expression and coexpression analysis

We analyzed the proportion of putative horizontally transferred genes in growth-regulated B. pseudomallei genes by using expression data reported by Rodrigues et al. ²⁸ In three major growth phases (early, log and stationary phase), about 900 genes were differentially expressed, among which 229 were predicted to be horizontally transferred. Nearly 50% of the latter (115 out of 229) displayed differential expression during exponential growth. Rodrigues et al (2006) noted that genes differentially expressed in stationary phase had a close relationship with strain pathogenicity, which confirmed that HGT genes might contribute to their host adaption. A total of 90 acquired genes that were differentially expressed in stationary phase could be assigned to the Burkholderia origin. This result suggests that the ancestor gene acquisition shaped the major characteristics of different Burkholderia species.

As many transferred genes in this strain are functionally unknown, based on co-regulation and coexpression of a set of genes in the same pathway, gene function can be predicted by quantitatively transforming the coexpression correlations as the degree of function similarity. Significant difference was observed after quantitatively mapping the coexpressed genes in different evolutionary nodes onto existing biological pathways (Supplementary Table S5). The transferred genes in each node involved in several important regulation and metabolism pathways, such as sphingolipid metabolism, exhibited a high degree of coexpression based on the correlation coefficients (Supplementary Table S5). By contrast, plotting S_pathway score values to recently transferred and anciently transferred genes, we found that the scores for bacterial secretion system and protein export pathway were several fold higher in case of anciently transferred genes than for recently transferred genes, which means the function of anciently acquired genes are more likely related with the secretion system and export protein pathway (Supplementary Table S5). The result is consistent with the GEIs analysis reported above that the pathogenicity of Burkholderia was mainly determined by ancient gene transfer.

The mapping of coexpressed genes onto biological pathway schemes provides a comprehensive way to identify previously unknown functional patterns in sets of genes with known functions. The results revealed a complex pattern across many biological pathways, indicating that in different evolutional stage, the genes in GEIs are linked to different biological processes. The results presented here indicate that the integration of gene gain history and expression patterns is valuable.

The diversity of the Burkholderia genus is reflected in the diversity of ecological niches occupied by the different species, ranging from soil to aqueous environments, associations with plants, fungi, amoeba, animals, and human, from saprophytes to endosymbionts and from biocontrol agents to pathogens. Our results suggest that HGT events occurred extensively in the adaptive evolution of Burkholderia. As a majority of these acquired genes encode hypothetical proteins or Burkholderia-specific proteins of unknown function, coexpression analysis of these gene products will be instrumental for a better understanding of their role in adaptation and strains divergence. This analysis also suggests that HGT played an important role in adaption and pathogenicity. Most of the transferred genes are located in the small chromosomes and were gained in ancient time, which has resulted in the major differences between different species. However, gene loss and large chromosome fragment rearrangements are also major causes for the diversity and adaption between these closely related strains.