Continuing evolution of Burkholderia mallei through genome reduction and large-scale rearrangements

Abstract

Burkholderia mallei (Bm), the causative agent of the predominately equine disease glanders, is a genetically uniform species that is very closely related to the much more diverse species Burkholderia pseudomallei (Bp), an opportunistic human pathogen and the primary cause of melioidosis. To gain insight into the relative lack of genetic diversity within Bm, we performed whole-genome comparative analysis of seven Bm strains and contrasted these with eight Bp strains. The Bm core genome (shared by all seven strains) is smaller in size than that of Bp, but the inverse is true for the variable gene sets that are distributed across strains. Interestingly, the biological roles of the Bm variable gene sets are much more homogeneous than those of Bp. The Bm variable genes are found mostly in contiguous regions flanked by insertion sequence (IS) elements, which appear to mediate excision and subsequent elimination of groups of genes that are under reduced selection in the mammalian host. The analysis suggests that the Bm genome continues to evolve through random IS-mediated recombination events, and differences in gene content may contribute to differences in virulence observed among Bm strains. The results are consistent with the view that Bm recently evolved from a single strain of Bp upon introduction into an animal host followed by expansion of IS elements, prophage elimination, and genome rearrangements and reduction mediated by homologous recombination across IS elements.

Keywords: bacterial evolution; bacterial virulence; comparative genomics; genome erosion.

FIG. 1.—

Multigenome alignment of eight Bp and seven Bm strains. Each circle represents a genome as presented in Materials and Methods. All genomes are aligned with Bp K96234 genome as a reference, which appears as the outermost multicolored circle. The Bp genomes are the eight outermost circles, and Bm genomes are internal. Areas in each color represent homologies between the subject genome and the reference. Areas in black in the reference chromosome (outermost circle) are regions present in K96243 but absent in query genome. Areas in black in each of the concentric circles are regions present in the query genome but absent from K96243. Representative Bp GIs are shown with red arrows. Representative clusters of Bp-specific genes absent from all Bm genomes (black on the K96243 ring) are highlighted with a yellow arrow.

FIG. 2.—

Pan-genome analysis of seven Bm and eight Bp strains. The CDSs in all Bm genomes (blue line) and Bp genomes (red line) were compared, and the number of new genes was plotted against the number of genomes used. The blue dashed line represents the extrapolated number of Bm strain-specific genes. The red dashed line represents the extrapolated minimum number of new genes discovered with each Bp genome.

FIG. 3.—

IS407A rearrangement of whole genomes. (A) Relative occurrence of the nucleotides in the 4-bp direct repeat of IS407A element insertion is shown as bar graphs for each position in the box below. (B) Four fully sequenced Bm genomes were aligned using WebACT. Red lines denote homology between chromosomes organized in the same orientation. Blue lines show homology but inverse orientation in each chromosome. Yellow lines show the presence of IS407A elements. Regions with no homology are shown by the absence of red or blue lines. (C) Four fully sequence Bp genomes were aligned, as described for Bm.

FIG. 4.—

IS407A mediated rearrangements of rrn and replichores among Bm strains. (A) rrn rearrangements due to IS407A recombinations. The outermost ring corresponds to Bp K96243 but is a representative of all Bp genomes. Green, ATCC23344; orange, NCTC10229; purple, NCTC10247; brown, SAVP1. The brown rrn cluster represents the locus rearranged into chromosome II in Bm. Red bars represent degenerate rrn loci. (B) guanine/cytosine-skew representation of the NCTC10247 genome generated in DNAplotter (Carver et al. 2009). Green represents a negative guanine/cytosine-skew suggesting ORF are oriented in the negative strand and purple represents a positive guanine/cytosine-skew suggesting ORF oriented in the positive strand. The origin of replication for NCTC10247 chromosome I is predicted at around 2.3 Mb and the termination around 1.0 Mb. (C) Alignment of chromosome II of ATCC23344 with chromosome I of BpK96243 as Bp representative. Regions of homology are represented by blue color. For the sake of clarity, only the genomic regions of interest are depicted.

FIG. 5.—

Genomic organization of the fliP locus in Bp and Bm. The wild-type fliP locus is present in all Bp. The fliP CDS is represented by dark purple rectangles. The NCTC10247 allele is interrupted by an IS407A (aquamarine) element. In ATCC23344, an ISBma1 (gray) is located upstream of the IS407A element and an additional 65 kb was inserted at this location. PRL20 had additional IS407A mediated insertions into fliP. Figures are not to scale, and IS407A elements in PRL-20 were made smaller.

FIG. 6.—

Distribution of dN/dS in variable and core genes of Bm genomes aligned with corresponding regions of the reference strain ATCC 23344. dN and dS rates were calculated as described in Materials and Methods. Cumulative data for the seven Bm strains is shown.

FIG. 7.—

Evolutionary tree of Bm showing the number of genes deleted and the evolutionary point of change. In total, 5,686 gene changes can be mapped onto this tree in a manner that assumes only single evolutionary deletion events. Conversely, 997 gene changes require 2 or 3 independent deletions of the same gene. Because we did not compare these genes with Bp, we do not know the ancestral state for 45 of these genes. These 45 genes could be additions or deletions with equal parsimony with mutations occurring along the basal branches of this tree. Letters in red represent the variable regions lost in each branch.