Comparative genomics of the classical Bordetella subspecies: the evolution and exchange of virulence-associated diversity amongst closely related pathogens

Abstract

Background: The classical Bordetella subspecies are phylogenetically closely related, yet differ in some of the most interesting and important characteristics of pathogens, such as host range, virulence and persistence. The compelling picture from previous comparisons of the three sequenced genomes was of genome degradation, with substantial loss of genome content (up to 24%) associated with adaptation to humans.

Results: For a more comprehensive picture of lineage evolution, we employed comparative genomic and phylogenomic analyses using seven additional diverse, newly sequenced Bordetella isolates. Genome-wide single nucleotide polymorphism (SNP) analysis supports a reevaluation of the phylogenetic relationships between the classical Bordetella subspecies, and suggests a closer link between ovine and human B. parapertussis lineages than has been previously proposed. Comparative analyses of genome content revealed that only 50% of the pan-genome is conserved in all strains, reflecting substantial diversity of genome content in these closely related pathogens that may relate to their different host ranges, virulence and persistence characteristics. Strikingly, these analyses suggest possible horizontal gene transfer (HGT) events in multiple loci encoding virulence factors, including O-antigen and pertussis toxin (Ptx). Segments of the pertussis toxin locus (ptx) and its secretion system locus (ptl) appear to have been acquired by the classical Bordetella subspecies and are divergent in different lineages, suggesting functional divergence in the classical Bordetellae.

Conclusions: Together, these observations, especially in key virulence factors, reveal that multiple mechanisms, such as point mutations, gain or loss of genes, as well as HGTs, contribute to the substantial phenotypic diversity of these versatile subspecies in various hosts.

Figure 1

Comparative genome content of eleven classical Bordetellae strains. The outermost circle indicates the predicted functional categories of each gene families in the pan-genome. Species-specific gene families are individually blocked with black line (B. bronchiseptica, B. parapertussis, and B. pertussis). Internal circles indicate the presence (solid color) or absence (unfilled) of each gene family in each strain examined. Circles from outer to inner are started with B. bronchiseptica strains followed by B. parapertussis strains and B. pertussis strains. Gene families that are shared in all strains are in the region A, while gene families that are not conserved in all strains are in the region B, C, and D. Region C indicates the phage encoding gene families that are only present in RB50. This figure was created using the Circos software [21].

Figure 2

Core and non-core genome of eleven classical Bordetellae pan-genome. The number of core and non-core gene families in each predicted functional category of the classical Bordetellae pan-genome (A) and each species pan-genome (B) are summarized (B. bronchiseptica non-human isolates: RB50, 253 and 1289, B. bronchiseptica human isolates: MO149, D445, and Bbr77, B. parapertussis strains: 12822 and Bpp5, and B. pertussis strains: Tohama I, CS and 18323). Functional categories are SF: surface proteins, CH: conserved hypothetical proteins, Misc: miscellaneous information, R: regulators, CI: central/intermediary metabolism, IS: phage/insertion sequence (IS) elements, IT: information transfer proteins, U: unknown proteins, PAC: pathogenicity/adaptation/chaperones, DS: degradation of small molecules, EM: energy metabolism, Pseudo: pseudogenes and DL: degradation of large molecules. * indicates the gene families that more than 80% of gene families are conserved in all strains.

Figure 3

Mathematical estimation of size of novel gene families, pan-genome, and core gene families. The number of novel gene families (A), pan-genome size (B), and the number of core gene families (C) were estimated for the classical Bordetellae (blue: all eleven genomes), B. bronchiseptica strains only (green: RB50, 253, 1289, MO149, D445, and Bbr77), and all the strains except B. pertussis strains (red: RB50, 253, 1289, MO149, D445, Bbr77, 12822, and Bpp5). If n genomes are selected from 11, there are 11!/ [(n-1)!*(11-n)!] possible combinations. Each possible combination is plotted as a point, and the line is fitted to the power law model adapted from the methods of Tetellin et al. [26]. γ in the power law model for the pan-genome size estimation was reported for each group. The numbers shown on the right side of each graph are the number of expected novel gene families, pan-genome size, and core gene families with 25 genomes.

Figure 4

Diversity in virulence factor genes/loci among the classical Bordetellae compared to RB50. The heatmap was generated based on nucleotide percentage identity compared to RB50 for each gene/loci. Absence of a certain gene and presence of a pseudogene are highlighted with white and sky blue color with #, respectively. * indicates the missing nucleotides due to the draft status of the genome. The gene content tree, the dendrogram of hierarchical clustering of the complete pan-genome matrix, was superimposed on top of the heatmap. Manhattan distances, linkage method, and 1,000 bootstrap replicates were used for the clustering.

Figure 5

Maximum likelihood phylogenetic tree of eleven classical Bordetellae with genome-wide SNP sites. The phylogenetic tree was reconstructed with genome-wide SNP sites based on B. bronchiseptica RB50.

Figure 6

Percent sequence similarity of ptx/ptl locus with flanking genes against Tohama I. Percent sequence similarity of the ptx/ptl locus and flanking genes based on Tohama I was plotted between 80% and 100% using zPicture [36]. Intergenic regions, coding regions, and a tRNA were highlighted with red, blue, and green, respectively.

Figure 7

The phylogenetic tree of the ptx/ptl locus. The tree was reconstructed by maximum likelihood methods with 1,000 bootstrap replicates based on the ptx/ptl locus sequences. B1920 represents B1834 and Tohama I represents CS because they have the identical sequences.