Genomic Analysis of Isolates From the United Kingdom 2012 Pertussis Outbreak Reveals That Vaccine Antigen Genes Are Unusually Fast Evolving

Abstract

A major outbreak of whooping cough, or pertussis, occurred in 2012 in the United Kingdom (UK), with nearly 10 000 laboratory-confirmed cases and 14 infant deaths attributed to pertussis. A worldwide resurgence of pertussis has been linked to switch to the use of acellular pertussis vaccines and the evolution of Bordetella pertussis away from vaccine-mediated immunity. We have conducted genomic analyses of multiple strains from the UK outbreak. We show that the UK outbreak was polyclonal in nature, caused by multiple distinct but closely related strains. Importantly, we demonstrate that acellular vaccine antigen–encoding genes are evolving at higher rates than other surface protein–encoding genes. This was true even prior to the introduction of pertussis vaccines but has become more pronounced since the introduction of the current acellular vaccines. The fast evolution of vaccine antigen–encoding genes has serious consequences for the ability of current vaccines to continue to control pertussis.

pertussis, genomics, evolution, vaccine

Whooping cough, or pertussis, is caused primarily by the bacterium Bordetella pertussis. In England and Wales, 9711 laboratory-confirmed cases were recorded in 2012, leading to 14 deaths of infants <3 months of age. This was much greater than the previous recent peak year in 2008, in which 902 cases were reported, despite levels of vaccine coverage and diagnostic methods not changing during this period [1, 2]. Similar outbreaks have been reported across the globe [3], contributing to the consensus that pertussis is a resurgent disease that might be no longer effectively controlled by current vaccination programs.

Resurgence has been linked to increased surveillance, better diagnostic techniques, incomplete vaccination of populations, and, primarily, the switch from whole-cell pertussis vaccines (WCVs) to acellular pertussis vaccines (ACVs), which contains 1–5 purified B. pertussis protein antigens: pertussis toxin (Ptx), filamentous hemagglutinin (FHA), pertactin (Prn), and fimbrial type 2 (Fim2) and Fim3. In the United Kingdom (UK), a 5-antigen ACV has been used. ACV-induced immunity appears to be shorter lived than that induced by WCVs, possibly resulting in an expanded pool of carriers, particularly adolescents, and decreased herd immunity [4, 5]. In addition, studies using an infant baboon model revealed that, while ACVs protect the individual from disease symptoms, they are less able to prevent colonization of and transmission from the vaccinee, compared with WCVs. Increased transmission of B. pertussis in populations using ACVs, compared with those using WCVs, is proposed to contribute to the resurgence [6]. Finally, it has been proposed that vaccine-induced escape mutants are arising, because ACV-induced immunity is focused on just a few antigens and because changes in these antigens might result in strains that are less well recognized by this immunity [7].

The frequency of different alleles of vaccine antigen–encoding genes among strains has changed over time [8–11]. The most common allelic profile among currently circulating strains (ptxA1-ptxP3, prn2, fim3-2, and fim2-1) is different from that of strains used for vaccine manufacture [12, 13], and isolates that do not express Prn are increasingly common [14–16]. PtxP refers to alleles of the ptx promoter. PtxP3 is now dominant worldwide [17], and some studies suggest that ptxP3 strains may have increased virulence, compared with ptxP1 strains [18].

The study of genetic changes in B. pertussis over time was hindered by the high levels of homogeneity among B. pertussis and the lack of fine-resolution tools. Recently, the genome sequences of a large panel of B. pertussis strains collected from around the world and across many decades were generated and analyzed [19]. This provided detailed information about the population structure and evolution of B. pertussis, revealing significant genetic changes among strains over the last 50 years. A lack of geographical clustering of strains suggested rapid strain flow between countries. However, this panel of strains did not contain isolates collected more recently than 2008, except for 3 isolates from the Netherlands, which were collected in 2009 and 2010. In addition, the panel did not include a substantial number of isolates from a specific outbreak, leaving the genetic make-up of isolates involved in such events largely unknown. Here, we analyze a large panel of UK strains, with a focus on strains from the recent UK outbreak, to understand the clonal structure of the outbreak and determine whether there is evidence for vaccine-mediated immunity driving the evolution of these strains.

METHODS

Accession Numbers

Genome sequence data has been deposited in the European Nucleotide Archive (http://www.ebi.ac.uk/ena/; Supplementary Table 1).

B. pertussis Strains

One hundred B. pertussis isolates were obtained from the National Reference Laboratory, Respiratory and Vaccine Preventable Bacteria Reference Unit, at Public Health England (Supplementary Table 1). Five strains were collected between 1920 and 1956 (hereafter referred to as the “prevaccine era”), 6 strains collected between 1957 and 2000 (hereafter, the “WCV era”), and 89 strains were collected between 2000 and 2012 (hereafter, the “ACV era”). Serotyping was performed using sera specific for antigens 1, 2, and 3 (89/596, 89/598, and 89/600, respectively; National Institute for Biological Standards and Controls, Potters Bar, UK) as previously described [12]. Tohama I (accession number: BX470248), a strain isolated in Japan in 1954, is the most widely studied strain, provides the reference genome sequence of B. pertussis [20], and is one of the strains used to produce ACVs used in the UK. B. pertussis isolates were grown on charcoal agar for 72 hours at 37°C.

DNA Preparation

Genomic DNA extraction was performed using the Qiagen DNA prep kit according to the manufacturer's instructions.

DNA Sequencing and Single-Nucleotide Polymorphism (SNP) Identification

Twenty four isolates were sequenced previously [19]. For the remainder, multiplex libraries, with fragment sizes between 300 bp and 500 bp, were prepared as previously described [21], with modifications [22]. Reads for each isolate were aligned to the Tohama I reference genome, using SMALT, version 0.7.4 (http://www.sanger.ac.uk/resources/software/smalt/). Base calls were made as previously described [21], using a combination of samtools, mpileup, and bcftools [23], allowing SNPs, and small insertions and deletions relative to Tohama I to be identified. Five strains produced poor-quality sequence and were excluded from the analysis, resulting in 95 strains being taken forward for analysis.

Phylogenetic Analysis

Maximum likelihood phylogenetic analysis was performed on variable sites from across the whole genomes, using RAxML under a general time-reversible evolutionary model and a Γ correction for among site rate heterogeneity [24]. One hundred random bootstrap replicates were run to provide support for relationships identified in the tree.

Analysis of SNP Densities

SNPs were reconstructed on to the phylogenetic tree using parsimony. SNP densities (defined as the number of SNPs per bp) within vaccine antigen–encoding genes (9 genes: fhaB, prn, fim2, fim3, and ptxA–E) or functional cell surface protein–encoding genes (591 genes, as categorized previously [20]) were calculated by counting the number of SNPs per bp of each gene. The difference between the mean per-gene SNP densities of vaccine antigen–encoding genes and cell surface protein–encoding genes was calculated. The significance of this difference was calculated using a nonparametric Monte Carlo simulation. In our randomizations of all the data, preserving relative sample sizes, it was recorded how often a difference as large, or greater than the difference above, was observed, by repeated randomly resampling 2 samples of the same size as above. Under this protocol, if n is the number of observations that have greater than or equal to the observed difference in SNP density and m is the number of simulations (in this case, 10 000), then P = [n + 1]/[m + 1] is the unbiased estimator of the type I error rate.

We performed a similar procedure to compare SNP densities in vaccine antigen–encoding genes between eras. To account for differences in SNP densities between strains from the different eras, the SNP densities of the vaccine antigen–encoding genes were normalized with respect to the SNP densities of all the genes considered (vaccine antigen–encoding and surface protein–encoding genes). A nonparametric Monte Carlo simulation compared the normalized SNP densities in the ACV antigen genes in ACV-era strains with those for pre–ACV-era strains, with P determined as described above.

Allele Typing

The different alleles of prn, ptxA, ptxP, fim3, and fim2 genes have been previously described [11] and were used to identify allele types from DNA sequences.

Analysis of prn From UK50

The prn locus was amplified from UK50 by polymerase chain reaction (PCR) analysis, using primers 5′-CCGCTGATTCGCCACAAG-3′ and 5′-GTGCGGTACTTGCCCTTG-3′. PCR products were cloned using the Gateway system (Invitrogen, Paisley, UK) and sequenced by Eurofins Genomics (Ebersberg, Germany), using standard M13 forward and reverse primers and internal primers 5′-GCGCACGCCTGTCCAAAG-3′ and 5′-TAGCGAGCCAGCACGTAG-3′.

Analysis of Differences in DNA Content Among Strains

To detect gene loss from strains, compared to the DNA content of Tohama I, coverage plots generated using the paired-end reads mapped to this reference genome were used to create a heat map. DNA sequence contigs that did not map to Tohama I were analyzed using Blastn and Blastx (http://blast.ncbi.nlm.nih.gov/Blast.cgi).

RESULTS

Phylogeny of UK Strains From 1920 to 2012

Phylogenetic analysis based on SNPs across the whole-genome sequences was performed to understand the evolutionary relationships between the UK strains analyzed (Figure 1 and Supplementary Table 1). Strains isolated during 1920–1982 form a cluster and are generally separated from strains isolated during 2008–2012. The most distinct clustering separates strains carrying the ptxP1 allele from those carrying the ptxP3 allele, which, as found elsewhere, is the predominant ptxP type among recent strains. This phylogenetic analysis was extended to place the UK strains in the global phylogenetic tree described elsewhere (Figure 2) [19]. This reveals that the UK ptxP3 strains separate into 2 clusters distinguished by the presence of the fim3-2 allele. The UK outbreak strains largely cluster with strains isolated mainly during the early 2000s from a variety of geographical areas, including North America, Europe, and Australia.

Figure 1.

Phylogenetic tree depicting the evolutionary relationships among the United Kingdom (UK) Bordetella pertussis isolates studied here. Maximum likelihood phylogenetic analysis was performed on variable sites from across the whole genomes, using RAxML. Strains are shaded according to their year of isolation and ptxP type.

Figure 2.

Phylogenetic relationships of United Kingdom (UK) strains within a global context. The UK isolates analyzed here are indicated.

Vaccine Antigen Allele Profiles

Previously, ptxP3-ptxA1-prn2-fim3-2 was defined as the dominant allele type circulating in the UK and other countries [17]. Typing of alleles among the outbreak strains reveal no recent change in this profile (Table 1). Numerous isolates deficient for the production of Prn have been reported in other countries, and a number of different mutations in prn responsible for this phenotype identified [14–16]. It has been suggested that loss of Prn expression has been selected by vaccine-mediated immunity pressure. Interestingly, just a single UK strain, UK50, was mutated for prn. This was identified by a lack of sequence reads mapping to a region of prn. The prn locus was amplified by PCR from this strain, and the resulting product was sequenced using Sanger sequencing. This identified that a recombination event between 2 copies of IS1663 resulted in a deletion/insertion mutation in which the 5′ 1326-bp region of the prn coding sequence was deleted. Aberrant mapping was not observed for any other UK strain. In other countries, a common prn mutation arose from insertion of IS481 into prn. We identified paired-end reads in which one read mapped within IS481 but the other did not and, thus, derives from the region flanking IS481. Mapping these reads to the reference genome identified the position of the copies of IS481 within each query strain. No IS481 insertions into prn were identified among UK strains. It is not clear why so few Prn-deficient strains, compared with other countries experiencing pertussis outbreaks, have been identified in the UK.

Table 1.

Frequency of Vaccine Antigen–Encoding Gene Alleles Among United Kingdom Strains, by Vaccine Era

Variable	Prevaccine Era, 1920–1956 (n = 5)	WCV Era, 1957–2000 (n = 6)	ACV Era, 2001–2012 (n = 84)
ptxP
1	100	100	6
3	0	0	94
ptxA
1	20	100	100
2	80	0	0
Prna
1	100	84	5
2	0	0	91
3	0	16	3
4	0	0	1
Fim2-
1	100	100	100
Fim3-
1	100	100	70
2	0	0	29
3	0	0	1
Serotype
1b	20	0	0
1, 2	40	50	37
1, 3	20	17	63
1, 2, 3	20	33	0

Data are percentage of strains tested.

Abbreviations: ACV, acellular pertussis vaccine; WCV, whole-cell pertussis vaccine.

^a Prn allele type was determined for just 76 ACV era strains, owing to poor mapping of reads in this region in 8 strains.

^b Serotype was not determined for one ACV era strain. Thus, frequencies are based on 83, not 84, strains in this era.

SNPs Specific to ptxP3 Strains

ptxP3 strains are the predominate type in current circulation and appear to have infection-associated biologic characteristics that differ from those of ptxP1 strains. The ptxP3 SNP appears to be both a direct cause and a marker for other genetic variations that contribute to this difference [18]. To investigate the genetic traits of UK ptxP3 strains, SNPs specific to this lineage were identified. In total, 22 such SNPs were identified (Table 2). Ten were intergenic, of which 7 were in the direct repeat region of IS elements, which are present in multiple copies in the B. pertussis genome. It is not clear whether these particular IS elements are functional. Twelve SNPs were in coding regions. Of those, 7 were nonsynonymous mutations (NSMs) and 5 were synonymous mutations (SMs). The 7 NSMs were in genes within the functional categories “transport/binding proteins,” “conserved hypothetical,” “pathogenicity,” “unknown,” and “regulation” or as “pseudogenes” as previously described [20]. All of these SNPs were also identified among the global panel of B. pertussis strains [19]. However, of the 22 SNPs identified here as being ptxP3 specific, only 10 were identified as being ptxP3-type specific in the previous study (Table 2); the other 12 SNPs were also identified among non-ptxP3 strains globally.

Table 2.

Single-Nucleotide Polymorphisms (SNPs) Specific to United Kingdom ptxP3 Strains

Locationa	Type	Mutation	Global ptxP3b	Details
36 857	INT	A→G	Yes	93 bp upstream of BP0032 (encoding a putative transport protein), 156 bp upstream of BP0033 (encoding GlyQ-glycyl-tRNA synthetase α chain)
617 083	INT	T→G	No	Within the 5′ repeat region of IS481 (BP0611), 31 bp upstream of transposase start codon.
617 084	INT	C→T	No	Within the 5′ repeat region of IS481 (BP0611), 32 bp upstream of transposase start codon.
1 077 844	INT	C→T	No	Within the 5′ repeat region of IS1663 (BP1035), 139 bp upstream of transposase start codon.
1 170 424	INT	A→G	No	Within the 5′ repeat region of IS481 (BP1114), 31 bp upstream of transposase start codon.
1 222 400	INT	A→C	No	Within the 5′ repeat region of IS481 (BP1157), 31 bp upstream of transposase start codon.
1 635 654	INT	T→G	No	Within the 5′ repeat region of IS481 (BP1557), 31 bp upstream of transposase start codon.
2 259 917	INT	G→C	No	Within the 5′ repeat region of IS481 (BP2135), 98 bp upstream of transposase start codon.
3 263 622	INT	A→C	Yes	193 bp away from BP3062, putative integral membrane transport protein.
3 988 168	INT	G→A	Yes	89 nucleotides away from the start codon of ptxA, ptxP3allele
196 307	NSM	T→C	Yes	BP0194, putative transport protein
299 559	NSM	C→T	Yes	BP0292, pseudogene, conserved hypothetical protein
1 331 840	NSM	G→A	Yes	Pseudogene, BP1261, hypothetical protein
1 547 488	NSM	A→G	No	BP1471, conserved hypothetical protein
2 374 322	NSM	T→C	Yes	BP2249, BscI, type III secretion apparatus protein
2 651 008	NSM	G→A	Yes	BP2502, Hypothetical protein
3 134 458	NSM	G→C	No	BP2946, probable transcriptional regulator
185 405	SM	G→A	No	BP0184, putative periplasmic protein
518 837	SM	T→C	No	BP0507, putative membrane protein
694 521	SM	A→G	Yes	BP0678, Putative peptide chain release factor
3 840 411	SM	G→A	Yes	BP3630, RpsH, 30S ribosomal protein
3 991 376	SM	C→T	No	BP3787, PtxC, Pertussis toxin subunit protein

Abbreviations: INT, SNP is in an intergenic region; NSM, nonsynonymous mutation; SM, synonymous mutation.

^a Tohama I reference genome coordinates (accession no. BX470248).

^b “Yes” denotes that the SNP was also defined as ptxP3 specific in study of the global Bordetella pertussis population [1], whereas “No” denotes that the SNP was not defined as such in the study.

SNP Rates Are High in Vaccine Antigen–Encoding Genes

Previously, it was identified that genes encoding functional cell surface proteins had higher SNP densities than the B. pertussis chromosomal average [19]. However, ACV-mediated immunity is exerting selective pressure primarily on the proteins used in these vaccines and might be driving their evolution. To explore this, the SNP density for the 9 ACV antigen-encoding genes (Ptx comprises 5 different proteins) and for the other 591 genes composing the cell surface protein category was calculated for all strains within each vaccine era and compared. Second, we investigated whether the SNP rate in ACV antigen–encoding genes had increased since the introduction of ACVs.

The difference in mean SNP density across genes within the 2 samples (calculated as the mean SNP density in vaccine antigen–encoding genes minus the mean SNP density in cell surface protein–encoding genes) was calculated. A nonparametric Monte Carlo simulation was used to assess the significance of this difference by determining how often a difference as large or larger than this was derived by randomly resampling, from the pool of vaccine antigen–encoding and cell surface protein–encoding genes, 2 samples of the same size as the 2 samples described above. This revealed that, in each era, vaccine antigen–encoding genes had significantly higher SNP densities than other cell surface protein–encoding genes (P < .05; Table 3), with the difference being greatest among ACV-era strains. This suggests that the vaccine antigen–encoding genes are faster evolving than other surface protein–encoding genes and that they were also faster evolving even before the introduction of widespread vaccination.

Table 3.

Single-Nucleotide Polymorphism (SNP) Rates in Vaccine Antigen–Encoding Genes, Compared With Other Cell Surface Protein–Encoding Genes, Among Strains Isolated in the United Kingdom and Globally During Different Vaccine Eras

Strain Source, Era	Strains, No.	SNPs/bp, No., Mean		Differencea	Normalized Differenceb	P Valuesc
		Vaccine Antigen–Encoding Genes	Cell Surface Protein-Encoding Genes
United Kingdom
Prevaccine era, 1920–1956	5	3 × 10−4	7.8 × 10−5	2.22 × 10−4	2.72	.045
WCV era, 1957–2000	6	4.75 × 10−4	5.9 × 10−5	4.17 × 10−4	6.40	.016
ACV era, 2001–2012	84	1.73 × 10−3	1.55 × 10−4	1.57 × 10−3	8.82	.0004
Global
Prevaccine era, 1920–1956	19	1.45 × 10−3	5.85 × 10−4	8.62 × 10−4	1.44	.012
WCV era, 1957–2000	204	2.62 × 10−3	1.01 × 10−3	1.61 × 10−3	1.56	.002
ACV era, 2001–2012	188	2.91 × 10−3	4.23 × 10−4	2.49 × 10−3	5.41	.0001

Abbreviations: ACV, acellular pertussis vaccine; WCV, whole-cell pertussis vaccine.

^a Calculated as the mean SNP density in vaccine antigen–encoding genes minus the mean SNP density in cell surface protein–encoding genes.

^b Calculated as [the mean SNP density in vaccine antigen–encoding genes minus the mean SNP density in cell surface protein–encoding genes]/mean overall SNP density.

^c For the difference in SNP rate between vaccine antigen–encoding genes and cell surface protein–encoding genes.

To compare SNP densities in vaccine antigen–encoding genes between eras, SNP densities within each era were normalized by dividing by the mean SNP rate across all of the genes concerned (ie, those encoding ACV antigens and cell surface proteins). In comparison to the prior analysis, this has less power, owing to the much smaller sample of ACV antigen–encoding genes, compared with the total number of cell surface protein–encoding genes. Although the normalized SNP density in ACV-era strains was greater than in pre–ACV era strains, the difference was not statistically significant (P = .160). However, the number of pre-ACV strains in this analysis was small. Thus, the same analyses were repeated, using SNP data from the global collection of strains for which the year of isolation was known [19] and incorporating the UK strains sequenced here (Table 3). Again, a significantly greater SNP frequency was found in ACV antigen–encoding genes, compared with other cell surface protein–encoding genes, in all 3 eras. This time, there was also a significantly higher SNP frequency in ACV antigen–encoding genes among ACV era strains, compared with pre–ACV era strains (P = .0177), suggesting that the relative SNP density in ACV antigen–encoding genes has increased since the introduction of ACVs. These results suggest that ACV antigen–encoding genes are intrinsically fast evolving and provide some support for the hypothesis that they have evolved even faster since the introduction of ACVs.

The more rapid evolution in the ACV antigen–encoding genes could be due to either a higher underlying mutation rate or different selection at the protein level. The different selection could be positive selection or weaker purifying selection. To distinguish between these 2 possibilities, SNPs were split into SMs and NSMs. High NSM but not SM rates would suggest altered protein-level selection. A higher rate of synonymous evolution (with possibly a weak nonsynonymous effect) would suggest higher mutation rates. Interpretation here is difficult, owing to well-described but incompletely understood correlation between SM and NSM rates.

Among WCV- and ACV-era global strains but not prevaccine era strains, the SM frequency was significantly higher in ACV antigen–encoding genes, compared with other cell surface protein–encoding genes (Table 4). When comparing ACV-era strains to pre-ACV era strains, the SM frequency in ACV antigen genes was significantly higher (P = .004). NSMs also occurred at a significantly greater frequency in ACV antigen–encoding genes, compared with other cell surface protein–encoding genes (Table 4). The magnitude of this effect is greater than that seen for SMs, suggesting the higher evolutionary rate of ACV antigen–encoding genes, compared with cell surface protein–encoding genes, is largely owing to protein-level selection on the antigens. Evidence for a strong recent increase is less clear-cut. When comparing strains from the ACV-era to pre-ACV era strains, the NSM frequency in ACV antigen genes was on the edge of significance (P = .051). Overall, our results provide support for the hypothesis that the genes encoding antigens chosen for ACVs are intrinsically fast evolving, in part owing to selection pressure on their antigenic products. We cannot discount the possibility that, in the ACV era, there has been an increase in the mutation rate (see “Discussion” section).

Table 4.

Synonymous and Nonsynonymous Mutation Rates in Vaccine Antigen–Encoding Genes, Compared With Other Cell Surface Protein–Encoding Genes, Among Strains Isolated Globally During Different Vaccine Eras

Mutation Type, Era	Strains, No.	SNPs/bp, No., Mean		Differencea	Normalized Differenceb	P Valuesc
		Vaccine Antigen–Encoding Genes	Cell Surface Protein–Encoding Genes
Synonymous
Prevaccine era, 1920–1956	19	1.32 × 10−4	2.4 × 10−4	−1.07 × 10−4	−0.45	.627
WCV era, 1957–2000	204	9.66 × 10−4	4.23 × 10−4	5.43 × 10−4	1.26	.045
ACV era, 2001–2012	188	9.68 × 10−4	1.76 × 10−4	7.92 × 10−4	4.20	.011
Nonsynonymous
Prevaccine era, 1920–1956	19	1.18 × 10−3	3.40 × 10−4	8.38 × 10−4	2.38	.006
WCV era, 1957–2000	204	1.96 × 10−3	5.83 × 10−4	1.37 × 10−3	2.28	.002
ACV era, 2001–2012	188	1.95 × 10−3	2.38 × 10−4	1.71 × 10−3	6.48	.0002

Abbreviations: ACV, acellular pertussis vaccine; SNP, single-nucleotide polymorphism; WCV, whole-cell pertussis vaccine.

^a Calculated as the mean SNP density in vaccine antigen–encoding genes minus the mean SNP density in cell surface protein–encoding genes.

^b Calculated as [the mean SNP density in vaccine antigen–encoding genes minus the mean SNP density in cell surface protein–encoding genes]/mean overall SNP density.

^c For the difference in SNP rate between vaccine antigen–encoding genes and cell surface protein–encoding genes.

Regions of Difference

Deletions have been a major feature of B. pertussis evolution and appear to be ongoing [20, 25]. Compared with the Tohama I reference genome, most of the major deletions observed among the strains analyzed here had been identified previously [25]. Numerous small deletions were found in only a few isolates or just 1 isolate, suggesting that deletion of DNA is common among B. pertussis strains. Interestingly, some deletions appeared to be specific to the UK ptxP3 strains, but no deletions specific to outbreak isolates were detected (Supplementary Figure 1).

Regions from individual strains that were not present in the Tohama I reference genome were investigated by BLAST analyses. These regions were also found within other B. pertussis genomes (BP18323 and CS) or in Bordetella bronchiseptica RB50, similar to that reported in other studies [26]. Thus, there were no novel insertions or gene acquisition among outbreak isolates.

DISCUSSION

The resurgence of pertussis in countries with high levels of vaccination has caused widespread concern. Among other factors, B. pertussis evolution away from efficient control by vaccine-induced immunity has been proposed as a contributor to this. Recently, whole-genome sequencing was used to define global genetic variability among B. pertussis isolates, and it identified genetic changes in the B. pertussis population over time [19].

Here, we have analyzed in detail the genomes of UK B. pertussis isolates, with an emphasis on strains from the 2012 outbreak. We are the first to show that many genetically distinct B. pertussis strains contributed to this outbreak and, importantly, that it was not due to the emergence of a novel, hypervirulent clone or to expansion of an individual lineage. Furthermore, outbreak strains were genetically very similar to those circulating during periods when the incidence of pertussis was low.

The ptxP3 type is the dominant clone worldwide, and UK outbreak strains are also predominantly of this type. Analysis of global isolates identified just 19 SNPs as being ptxP3 specific [19]. Here, 22 SNPs distinguished ptxP3 strains from ptxP1 strains. However, just 10 of these were common to both sets of ptxP3-specific SNPs. If ptxP3 strains have increased fitness or virulence than older isolates, our analysis suggests that very few SNPs are responsible for this, or that particular combinations of SNPs are important, only some of which are ptxP3 specific. Overall, these data argue against large-scale genetic changes being behind the recent resurgence in pertussis.

Changes in alleles of the genes encoding vaccine antigens have been well documented [27], and these findings support the hypothesis that selection pressure from ACV-induced immunity is a driver of B. pertussis evolution. However, definitive studies to demonstrate that allelic variation enhances evasion of vaccine-mediated immunity are lacking and particularly difficult to perform, because of the inability to conduct studies with human hosts, and studies using animal models struggle to detect subtle changes and will not include population-level effects that are certainly important for selection of variants among B. pertussis worldwide. Here, we provide compelling evidence that genes encoding ACV antigens are evolving more rapidly than other cell surface protein–encoding genes (which we consider the most suitable comparator group), containing a significantly higher frequency of SNPs in each of the vaccine eras. Interestingly, this was true even in the prevaccine era. It is likely that, even in the absence of vaccination, the natural immune response to these antigens creates selective pressure, particularly for a pathogen that is restricted to the human respiratory tract. Of particular importance is that we calculated that ACV antigen–encoding gene evolution rates have increased significantly since the introduction of ACVs, the first demonstration of this effect. This might suggest that the use of ACVs has increased selection pressure on ACV antigens, selecting for ACV antigen–encoding gene variants. However, we also calculated that, while the frequency of SMs in ACV antigen–encoding genes was significantly higher in ACV-era strains, compared with older strains, the frequency of NSMs was on the edge of significance (P = .051). In turn, this suggests that selection pressure from vaccine-mediated immunity is not the sole driving force for ACV antigen–encoding gene variation. A different interpretation is that the mutation rate of ACV antigen–encoding genes has increased since the introduction of ACVs. If synonymous sites are under weak purifying selection pressure (ie, not perfectly neutral), then there is a lag between a SNP arising and its elimination by this selection, resulting in an excess of SNPs in the modern era. However, normalizing the ACV antigen–encoding gene SNP rates by the SNP rates for all genes within the era largely eliminates this effect (ie, SMs in cell surface protein–encoding genes should be equally over-represented in the modern era). However, if SMs in ACV antigen–encoding genes and cell surface protein–encoding genes are under different intensities of purifying selection, then our result could be found.

Either way, the more rapid evolution at the protein level (as determined by NSMs) of ACV proteins, compared with other cell surface proteins, across all eras suggests that strains will become increasingly mismatched to those used for vaccine production, and this could lead to decreased vaccine efficacy over time. The ACV antigens were chosen on the basis of their immunogenicity, but it could be that this property has driven the relatively high evolution rates of the genes encoding these antigens. Our results raise fresh concerns over the ability of current ACVs to continue to control disease.

Notes

Financial support. This work was supported by a Public Health England (PhD studentship to N. K. F., A. R. G., A. P.) and the Wellcome Trust (grant 098051 to S. R. H. and J. P.).

Potential conflicts of interest. All authors: No reported conflicts.

All authors have submitted the ICMJE Form for Disclosure of Potential Conflicts of Interest. Conflicts that the editors consider relevant to the content of the manuscript have been disclosed.