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

Sequencing of seven new classical Bordetella genomes

B. bronchiseptica strain 1289 was sequenced and assembled at the Pennsylvania State University using Illumina [16], 454 with GSFLX Titanium chemistry [17, 18], and optical maps [19] from OpGen (Gaithersburg, MD). The other six strains were sequenced at the Sanger Institute. The genomes of Bpp5, 253, MO149, and 18323 were sequenced to approximately 6 to 9-fold coverage with ABI3730 sequencers [20]. The genomes of D445 and Bbr77 were sequenced on Illumina [16] and 454 [17, 18]. Read assembly resulted in finished genomes of MO149, Bpp5, 1289 and 18323, 4 contigs in 253, 11 scaffolds in D445, and 16 scaffolds in Bbr77. The sum of the contig or scaffold lengths ranged from 4.1 Mb to 5.2 Mb, which is comparable to genome sizes of previously sequenced strains.

B. pertussis strains (Tohama I, 18323, and CS) contained no novel genes except transposases or insertion elements, confirming that B. pertussis has evolved from a B. bronchiseptica-like ancestor by genome decay (Figure 1). While we did not observe many genome rearrangements between Tohama I and CS, extensive genome rearrangement distinguishes Tohama I and 18323 (unpublished). More phage-related and membrane protein encoding genes are present in B. bronchiseptica non-human isolates (RB50, 253, and 1289), while more transposases are present in B. bronchiseptica human-isolated strains (D445, Bbr77, and MO149). Notably, a locus novel to the classical Bordetellae was observed in both D445 (D445_1578 – D445_1707) and Bbr77 (Bb77_4113 - Bbr77_4216) strains. This locus is also present in B. petrii and includes a Type IV secretion system-like locus, which is being further investigated (unpublished). B. parapertussis _ov Bpp5 strain has about 150 strain-specific gene families comprised mostly of hypothetical proteins, phage-related proteins, and putative secreted proteins. Comparing the two lineages designated B. parapertussis, there are about 300 gene families that are present in either B. parapertussis _hu or B. parapertussis _ov , but not both. Most of these gene families are shared with B. bronchiseptica strain RB50 and include exported proteins, hypothetical proteins, and membrane proteins. While B. parapertussis _hu 12822 is missing the locus encoding a putative type II secretion system, B. parapertussis _ov strain Bpp5 is missing the locus encoding flagella. Furthermore, Bpp5 has a unique 12 KB plasmid containing genes predicted to encode proteins involved in partition, replication, plasmid conjugal transfer, mobilization and transcriptional regulation. We have not observed this plasmid in any other classical Bordetella strains.

The classical Bordetella pan-genome and core genome

To determine the global gene repertoire (pan-genome) and the universally shared “core” genome, we identified orthologous gene families via OrthoMCL [22, 23] using all of the published Bordetella genomes that have manually curated annotations (RB50, Tohama I, CS, 12822, B. petrii, and B. avium). The pan-genome of these strains has 8,425 gene families, consisting of 1,778 core gene families and 6,647 gene families that are missing in at least one strain (see Additional file 1). While four Bordetella species remain to be sequenced and others are not finished or adequately annotated, these results do provide solid, yet preliminary, insight into the global gene repertoire of Bordetella. Throughout the remainder of this study, all investigations will focus on the evolutionary relationship of the classical Bordetellae. The classical Bordetella pan-genome includes 5,558 gene families, consisting of 2,857 core gene families (Region A in Figure 1) and 2,701 non-core gene families (Regions B, C and D in Figure 1) (see Additional file 2). The core gene families are mainly involved in energy metabolism, central metabolism and information transfer (Figure 2A). The most abundant gene family groups in these genomes encode surface proteins and 40% of these belong to the non-core genome. Variability in these proteins, which are likely surface exposed and antigenic, would be expected to affect host immune recognition [24].

When bacteria specialize to an ecological niche, this adaptation is often accompanied by a loss of genes that are no longer required or beneficial in the new niche [25]. This concept can explain the large-scale genome loss exhibited by B. pertussis and B. parapertussis, which appear to have evolved from a zoonotic generalist into specialized human pathogens [2]. If other lineages had adapted to particular environments (i.e. subsets of hosts), we would expect those lineages to lose genes not required in that environment, but retain a set of shared core genes required for their shared characteristics, such as the ability to infect the respiratory tract of their particular host. B. bronchiseptica complex I non-human isolates (RB50, 253, and 1289; sharing 4,509 gene families) and complex IV human isolates (MO149, Bbr77 and D445; sharing 4,271 gene families) strains have larger core genomes, than B. pertussis strains (Tohama I, 18323, and CS; sharing 3,235 gene families) (Figure 2B) [4]. In fact, the core genome of all available classical Bordetella strains excluding B. pertussis isolates (3,652 gene families) was larger than that of B. pertussis strains, suggesting that this lineage has selectively lost genes/functions that are retained in all others. The non-core genome (2,701 gene families not present in all strains) contains a large variety of gene families and many phage-encoding genes, and likely contributes to the phenotypic variation of this group.

Open pan-genome with limited gene acquisition

To evaluate and compare the total gene pool of all classical Bordetella subspecies (i.e. all eleven sequenced strains), compared to that of all strains excluding B. pertussis strains, and that of B. bronchiseptica strains alone (RB50, 253, 1289, MO149, Bbr77, and D445), we employed the prediction method described by Tettelin et al.[26]. Based on the number of novel gene families per additional genome analyzed, all of the groups (the classical Bordetella subspecies, all the strains excluding B. pertussis strains, and B. bronchiseptica strains only) appeared to have what Tettlin et al. described as open but “slowly closing” pan-genomes [27] (Figure 3A). This indicates that more sequenced genomes have the potential to provide additional genomic diversity but that beyond a couple of hundred sequenced genomes novel genes are less likely to be discovered.

A complementary approach to assess gene acquisition was used to calculate the increase in pan-genome size with each additional genome sequenced (Figure 3B). This expected pan-genome size (n) upon sequencing an Nth genome was modeled by the power law function n = κN ^γ , where an open pan-genome has a γ value greater than zero and less than one, with lower values signifying a more closed genome with fewer acquired genes. The γ value for the classical Bordetella subspecies (0.090) is lower than that for Bacillus cereus (0.43) [28], indicating that the pan-genome of the classical Bordetellae is open but grows more slowly than that of B. cereus with each additional genome sequenced. With 25 additional genomes, the pan-genome of the classical Bordetella subspecies will reach approximately 6,000 gene families, while the pan-genome of the group excluding B. pertussis strains is estimated to contain 6,026 gene families. Furthermore, since a smaller core genome size reflects smaller sets of shared functional genes, the small size of the core genome containing B. pertussis strains is additional evidence that this lineage has lost the functions encoded by approximately 800 genes (Figure 3C). This further supports the previous hypothesis that B. pertussis strains have adapted to their niche by genome reduction [2] and that there is a correlation between the pan/core genome size and the host range.

Variations in virulence factors amongst the classical Bordetellae

To examine the possible basis for differences in virulence phenotypes amongst strains/lineages, we compared the presence/absence of factors thought to be involved in interactions with the host, loosely referred to as “virulence factors”, including filamentous hemagglutinin (FHA), fimbriae (Fims), pertactin (PRN), tracheal colonization factor (TcfA), adenylate cyclase/hemolysin (ACT), dermonecrotic toxin (Dnt), pertussis toxin (Ptx), Bordetella resistance to killing (BrkA) protein [29], O-antigen [30], Type III secretion system (TTSS) [31], and Type VI secretion system (T6SS) [32]. Only genes encoding FHA, PRN, and TTSS locus are conserved in all strains (Figure 4). B. bronchiseptica complex IV strains (D445, Bbr77, and MO149) have the most divergent fha locus (as low as 92% identity) compared to RB50, while B. pertussis strains have the most divergent fim loci (as low as 84% identity) with gene deletions and multiple internal sequence variations. Although tcfA varies in all strains (as low as 84% identity), genes encoding PRN segregate into clear distinctive groups between B. bronchiseptica complex I/B. parapertussis strains and B. bronchiseptica complex IV/B. pertussis strains, as previously described [4]. Together, the presence, absence or divergence of these virulence factors do not always follow the phylogenetic relationships, suggesting some other processes, in addition to sequence divergence by accumulation of point mutations, may contribute to variation of virulence factor genes.

In a prior analysis of the hypovirulent B. bronchiseptica strain 253, we observed that it does not contain the whole ACT locus [13]. Notably, strains 1289, D445, Bbr77, MO149, and 18323 all contain a frameshift mutation within cyaB; however, each strain remains beta-hemolytic and grows within the murine respiratory tracts (unpublished). Therefore, the functional effects of this mutation on ACT secretion and hemolysis appear less dramatic than that of the mutation in strain 253. Strikingly, strains Bbr77 and MO149 do not contain the entire ptx locus and the other human isolate (D445) has a divergent ptx locus (as low as 88% identity) compared to RB50. dnt is divergent (as low as 93% identity) or missing in complex IV strains compared to RB50. One of the complement resistance factors, brkA, is very similar (>97% identity) in all strains except that it is a pseudogene in 12822 and RB50. Combined, these data reveal a discord between the phylogenetic relatedness of strains and the variation in their individual virulence genes.

In addition to small loci (<10 kb) responsible for encoding adhesins and toxins, we compared large loci (>10 kb) containing contiguous genes with highly coordinated functions. The T3SS locus is very similar (>98% identity) among all strains, with the notable exception being a pseudogene present only in B. parapertussis _hu strain 12822 [12, 33]. However, other loci were much more variable. Consistent with previous reports [2, 4, 7], the genetic makeup of the O-antigen locus (0-26 kb) widely fluctuated among strains, with the entire locus missing in all three B. pertussis stains. In strains in which they are present, these genes share as little as 90% identity with their apparent orthologs. In some cases, the locus includes completely different sets of genes. Two pseudogenes are present within the O-antigen locus in Bpp5, suggesting relatively recent loss of function. Similarly, a locus encoding a putative T6SS shows a high degree of variations. For example, B. pertussis strains and D445 are missing parts of the T6SS locus, while both human and ovine B. parapertussis strains either have pseudogenes or are missing subsets of genes within this locus. A predicted pseudogene is also present in this locus of B. bronchiseptica strain 1289. Together, these data reveal the intriguing tendency of virulence loci to be lost or divergent in the human isolates (B. bronchiseptica complex IV strains, B. parapertussis _hu strains, and B. pertussis strains), possibly signifying differential roles in different hosts [34].

Phylogenomic analysis

The observation that there are discrepancies between the phylogenetic relationships of sets of genes suggests that a small set of genes may be limited in its ability to accurately represent phylogenetic relationships among these strains [35]. Thus, to more accurately and robustly evaluate the genetic relationships between Bordetella isolates, we built a phylogenetic tree with all the genomic content, using a previously defined pan-genome matrix in our analysis. This matrix was constructed with genomes and gene families as columns and rows, respectively. Each cell in the matrix presents 1 or 0, depending on the presence or absence of the gene family in each genome. By generating a gene content tree based on this matrix, excluding unique gene families in each strain, we observed many similarities and some notable differences as compared to earlier studies that only used a small set of housekeeping genes. B. pertussis strains form a separate lineage that is isolated from the other two species (Figure 4). This clustering trend can be explained by the drastic genome reduction (absence of genes) in B. pertussis strains. However, B. parapertussis strains are closely related to B. bronchiseptica strains, reflecting a relatively large set of shared core genes. Notably, the B. parapertussis _ov strain (Bpp5) clustered together with B. parapertussis _hu strain 12822, consistent with the species designation but different from the previous minimum spanning tree based on MLST [4], which clustered Bpp5 within B. bronchiseptica complex I strains and separated 12822 from B. bronchiseptica strains.

To include all the additional information contained in the individual SNPs, a phylogenetic tree was generated based on genome-wide SNP candidate sites against the reference genome (RB50) with B. petrii as an outgroup (Figure 5). Similar to the previously published MLST minimum spanning tree, this SNP-based tree suggests that B. pertussis is closely related to B. bronchiseptica complex IV strains [4]. Moreover, building the phylogenetic tree using genome-wide SNPs supports two important genetic relationships that previous phylogenetic analyses did not reveal. First, B. pertussis and B. bronchiseptica complex IV appear to share a more recent common ancestor distinct from that shared by B. parapertussis and B. bronchiseptica complex I strains. Second, B. parapertussis _ov strain Bpp5 is more closely related to B. parapertussis _hu strain 12822 than to B. bronchiseptica complex I strains. The agreement on the closer relationship between the two B. parapertussis clades of both the SNP-based and gene content-based trees suggests that the more distant relationship observed in the previously published MLST-based tree is likely due to the limited amount of information in the seven housekeeping genes that were used.

SNP densities across each genome, compared to the RB50 genome, differ between strains. B. bronchiseptica complex IV strains (MO149, D445, and Bbr77) have the highest overall density of SNP sites (~12 SNPs/1,000 bp), followed by B. pertussis strains (~9 SNPs/1,000 bp), B. bronchiseptica complex I strain 253 (~7 SNPs/1,000 bp), B. parapertussis strains (~4 SNPs/1,000 bp), and complex I strain 1289 (~2 SNPs/1,000 bp), reflecting their overall relatedness to RB50. Of all SNPs, about 61% are synonymous (average dS = 0.0170), while ~26% are non-synonymous (average dN = 0.0023), reflecting the overall pressure of purifying selection. ~12% of SNPs were found within intergenic regions, and approximately 1% are in pseudogenes. Genes encoding phage-related, hypothetical proteins and some virulence factors, including ptx, contained amongst the highest SNP densities (see Additional file 3). One explanation for this observation could be positive selection, for example via the proposed vaccine-driven selection for antigenic variation of genes encoding vaccine components [37, 38]. However, dN/dS ratios are lower than one in most of these genes, indicating overall negative (purifying) rather than positive (diversifying) selection. These results do not support the view that the normal rate of SNP generation, followed by positive selection for variation, resulted in this diversity, and raise the possibility that other mechanisms of diversity generation, such as insertions or deletions, HGT and functional divergence, occurred within these discrete loci.

Horizontal gene transfer and divergent evolution of the classical Bordetellae

To more closely examine the divergence of the ptx/ptl loci, we compared the percent sequence similarity of these genes and their flanking genes in all sequenced strains to that of Tohama I (Figure 6). The ptx/ptl loci of B. pertussis 18323 and CS and B. parapertussis 12822 and Bpp5 are closely related to that of Tohama I (>98% and >95% identity, respectively), while those of B. bronchiseptica non-human strains are much more divergent (as low as 87% identity). Conversely, B. bronchiseptica human isolates, which are much more closely related based on overall SNPs (Figure 5), have either low percentage similarity to Tohama I (D445; as low as 82% identity) or do not contain the entire ptx/ptl loci (Bbr77 and MO149) (Figure 6). There were ten or less SNPs distinguishing the ptx/ptl loci of seven recently published B. pertussis genomes (see Additional file 4). The relatively high SNP densities of the ptx/ptl loci do not appear to be due to positive selection because pair-wise dN/dS ratios are below one, indicating overall purifying selection (see Additional file 5).

An alternative explanation for the high SNP density in the ptx/ptl loci could be the introduction of variation via horizontal gene transfer (HGT). Gerlach et al. previously speculated that the ptx/ptl loci had features of pathogenicity islands (PAIs) based only on the clusters of virulence genes in these loci, the presence of a tRNA (known to be associated with HGT and often a DNA integration target site [39]), and the absence of these loci in some strains [40]. Using the entire sequences of all eleven classical Bordetella strains, we analyzed each genome with Alien_hunter [41], which detects horizontally transferred genome segments. Nine to eighteen percent of each genome was identified as containing potential HGT candidates based on atypical genomic composition (see Additional file 6). The loci encoding phage-related proteins, alcaligin biosynthesis proteins, cytochrome ubiquinol oxidase, NADH-ubiquinone oxidoreductase, ATP synthase and many adhesins were identified as HGT candidates. Among the candidates were a number of genes required for assembly of O-antigen, consistent with our prior prediction of HGT within this locus [7, 42, 43]. Notably, Alien_hunter also predicted that the ptx locus and a part of its associated secretion system (ptl) locus were acquired by HGT.

To further investigate possible HGT of ptx/ptl genes, phylogenetic trees were generated based on individual genes within the ptx/ptl loci; these gene trees (Figure 7 and see Additional file 7) were then compared to the genome-wide SNP tree (Figure 5) to identify incongruence in the tree topology. The majority of the ptx/ptl gene trees are similar to the genome-wide SNP tree in that each strain is clustered based on their species designation, such as B. pertussis, B. parapertussis, and B. bronchiseptica non-human strains. However, the B. bronchiseptica human isolate D445 locus is located on a long branch separate from all the other strains, implying that the ptx/ptl loci in D445 has an evolutionary history different from other Bordetella species. Additionally, B. pertussis clades and B. parapertussis clades clustered together in the individual gene trees unlike in the genome-wide SNP tree. These results suggest HGT being a source of diversity within this locus, although they do not rule out other possibilities.

In addition to the higher overall SNP density, the distribution of SNPs across the ptx/ptl loci appears to be non-random, with areas of low density (consistent with the flanking genes and much of the genome) and areas of much higher density. For example, B. bronchiseptica complex IV strains have accumulated the highest density of SNPs (D445) or have lost the entire locus (MO149 and Bbr77), suggesting a change in the requirement for Ptx in this lineage. However, much of the ptl locus of the non-human B. bronchiseptica strains (RB50, 1289, and 253) have low SNP density, while the ptx locus, most notably the ptxC gene, has much higher SNP density (Figure 6). Interestingly, the large majority of these SNPs across the locus result in silent (synonymous) mutations (dN/dS < 0.6), suggesting the high SNP density is not due to positive selection for variations. It is also interesting to note that sequence similarity differs between the ptx locus (average 96% identity) and the ptl locus (average 98% identity) in all the strains except strain D445, suggesting a different evolutionary history for the ptx and ptl loci. Together, the ptx/ptl loci appear to have been horizontally transferred to the classical Bordetellae and have diverged in different lineages.

Sequencing and assembly of genomes

The genomes of Bpp5 (accession number HE965803, plasmid accession number HE965804), 253 (accession number HE965806), MO149 (accession number HE965807), also known as D444, and 18323 (accession number HE965805) were sequenced to approximately 9-fold coverage (or 6-fold for 18323) by shotgun sequencing from two genomic shotgun libraries, pMAQ1Sac_BstXI (with insert sizes of 6-9 kb) and pOTWI2 (with insert sizes of 4-5 kb; 5-6 kb), that were sequenced using big-dye terminator chemistry on ABI3730 automated sequencers [20]. End sequences from large insert fosmid libraries in pCC1FOS with an insert size of 38-42 kb were used as a scaffold. This generated approximately 0.2- to 0.4-fold coverage. The assembly was generated using phrap2gap. All repeat regions and gaps were bridged by read-pairs or end-sequenced polymerase chain reaction (PCR) products that were again sequenced with big dye terminator chemistry on ABI3730 capillary sequencers. The sequences were manipulated to the ‘Finished’ standard for Bpp5, MO149, and 18323 or to the ‘Improved High Quality Draft’ standard for 253 [54]. The genomes of D445 (accession number HE983627) and Bbr77 (accession number HE983628) were sequenced on both the Illumina GAII platform (54 bp paired end library) and the Roche 454 machine using GSFLX Titanium chemistry (3 kb insert paired end library) [16]. The Illumina data were assembled using Velvet [55] at a k-mer value of 39 or 31, respectively. This assembly was then shredded and combined with 454 data using Newbler [17, 18].

The genome of 1289 (accession number HE983626) was sequenced using the Roche 454 GSFLX sequencing platform [12, 56] and assembled using Newbler [17, 18], which resulted in 32-fold coverage and 133 contigs. The order of these contigs was determined by using RB50 as a reference genome. Gaps were closed by using sequences from PCR amplification of gap regions generated using big dye terminator chemistry on an ABI3730 capillary sequencer, and were assembled using PHRED, PHRAP, and CONSED [57–59], sequences generated from a 76 bp paired end library using a GAIIX [16] (resulting in 188-fold coverage), and sequences generated from a 3 kb insert paired end library using a GSFLX platform [17, 18]. These long-tag paired end reads and optical maps [19] (OpGen, Gaithersburg, MD) were used as scaffolds to identify and correct misassemblies and verify the final assembly. SNP sequencing errors were corrected by mapping the sequences generated by the GAIIX to the final genome assembly using inGAP [60]. The final draft of the 1289 genome consists of seven ordered contigs, covering >99.99% of the genome.

Annotation of genomes

Both RAST [61] and an automated annotation transfer tool at the Sanger Institute (unpublished) were used to annotate the finished or draft genomes of MO149, 1289, 253, D445, Bbr77, Bpp5, and 18323. The annotation tools were also used to annotate the ptx/ptl loci of six B. pertussis genomes (B0558, B1193, B1831, B1834, B1917, and B1920). We used RAST results for the novel regions that were not automatically transferred from the reference genomes, RB50, 12822, and Tohama I, based on their species classification. Each novel gene prediction was also curated with BLAST [62] and FASTA [63] results, and Pfam [64] and Prosite [65] were used to identify protein motifs. Transmembrane domains, signal sequences, and rRNA genes were identified with TMHMM [66], SignalP [67], and BLASTN [62], respectively. ISFinder was also used to detect Insertion Sequence (IS) elements in the genomes [68]. Manual curation was done with these novel regions using Artemis [69] and Artemis Comparison Tool (ACT) [70].

Pan-genome analysis

Coding sequences were extracted from the eleven classical Bordetella genomes as well as B. petrii and B. avium genomes, and orthologous gene families were determined using OrthoMCL [22, 23], which defines putative pairs of orthologs based on reciprocal all-against-all BLASTP [62] with a cutoff E-value of 10^-5, over 70% length coverage, and at least 70% identity. The orthologs are then clustered using Markov cluster algorithm. When gene families were not clustered and had no BLAST hits, they were considered strain specific gene families. The results of OrthoMCL were converted to a pan-genome matrix profile for further analysis. A pan-genome matrix was constructed with each column as a genome and each row as a gene family. Cell (i,j) in the matrix is 1 when gene family “i” is present in genome “j”, or 0 when gene family is absent. A graphical representation of the classical Bordetellae pan-genome (Figure 1) was created using the Circos software [21].

The sequential inclusion of up to eleven strains (eleven columns in a pan-genome matrix) was simulated in all possible combinations (N = 11!/ [(n-1)!*(11-n)!]). The number of new gene families, the core genome, and pan-genome size were estimated using methods adapted from Tettelin et al.[26]. The R [71] function nls was used for non-linear least squares regression on the mean of the new gene families, core gene families, and pan-genome distributions. The number of expected new gene families and core gene families (n) determined by sequencing an Nth genome was modeled by the power law function n = κN^-α, and the pan-genome size (n) by a power law n = κN ^γ . For estimation of novel gene families, the functions were fit to the mean values for all N ≥ 4, when there are more than 6 genomes to avoid the left side bias that could skew the results.

Gene content tree

Hierarchical clustering was performed for the complete pan-genome matrix using pvclust in R [71], modified from Lukjancenko et al.[72]. The dendrogram was generated by hierarchical clustering with complete linkage method and Manhattan distances using R [71]. Bootstrap values were computed by resampling the rows of the matrix 1,000 times.

Virulence factor comparison

Genes and the loci that encode the known virulence factors, filamentous haemagglutinin (FHA), fimbriae (Fims), pertactin (PRN), tracheal colonization factor (TcfA), invasive adenylate cyclase/haemolysin (ACT), dermonecrotic toxin (Dnt), pertussis toxin (Ptx), Bordetella resistance to killing (BrkA) [29], O-antigen [30], Type III secretion system (TTSS) [31], and Type VI secretion system (T6SS) [32], were compared among the eleven genomes via ACT [70]. Percent sequence similarity was calculated based on RB50 sequences with BLASTN [62], and genes that either contain a frame-shift mutation or an in-frame stop codon, or that are absent were highlighted with different colors in the heatmap that was generated by R [71]. The phylogenetic tree (Gene content tree) based on the presence and absence of each gene family in the pan-genome, excluding strain-specific genes, was superimposed on the heatmap.

SNP analysis

All seventeen (B. bronchiseptica strains: 253, 1289, D445, Bbr77, and MO149, B. parapertussis strains: 12822 and Bpp5, B. pertussis strains: Tohama I, 18323, CS, B0558, B1193, B1831, B1834, B1917, and B1920, and B. petrii strain DSM 12804) genomic sequences were randomly shredded into 54 bp long reads and mapped onto the reference genome (RB50), using Ssaha v2.2.1 [73]. High quality candidate SNPs were identified using ssaha_pileup, and 128,752 SNP sites were identified in at least one strain based on RB50. Phylogenetic trees were constructed with RAxML v7.0.4 [74] for all SNP sites in the reference genome, using a General Time Reversible (GTR) model with a gamma correction for among site rate variation and ten random starting trees [75].

dN/dS Analysis

The available Bordetella genomes were aligned onto the reference genome (RB50) using MAUVE aligner with default parameters [76]. Orthologous sequences were extracted from the alignments using gene coordinates defined in RB50 and processed to cover the entire coding regions of each individual genome, excluding pseudogenes. The nucleotide sequence alignments were later refined using their corresponding amino acid sequence alignments. dN and dS values were computed using PAML package [77] with the Nei-Gojobori method [78].

HGT Detection

Putative horizontally acquired regions in the eleven classical Bordetella strains were predicted by Alien_hunter, which is able to detect genomic regions that may originate from foreign sources using Interpolated Variable Order Motifs (IVOMs) [41]. For the phylogenetic tree comparison, multiple alignments of individual genes and the entire ptx/ptl loci were generated by the MEGA5 software [79], and maximum likelihood trees were constructed with a Tamura-Nei model and 1,000 bootstrap replicates.