ABSTRACT

INTRODUCTION

Bordetella pertussis is the causative agent of whooping cough (pertussis), a respiratory disease affecting all age groups, but with the highest disease severity in unvaccinated infants. Whole-cell vaccines against pertussis were introduced in the United States during the 1940s, and greatly reduced disease incidence, but were replaced during the 1990s by acellular vaccines, which produced less-severe side effects (1,–4). In addition to diphtheria and tetanus toxoids, the entire childhood vaccination series of diphtheria–tetanus–acellular-pertussis (DTaP) vaccine contains inactivated pertussis toxin (Ptx) and one or more additional virulence-related bacterial components: filamentous hemagglutinin (Fha), pertactin (Prn), or fimbria (Fim) types 2 and 3. Additionally, a single dose of the adolescent and adult booster Tdap (diphtheria, tetanus, and acellular-pertussis vaccine) was recommended in 2005 to counteract the rise of reported cases in adolescents and adults (3, 5, 6). Despite availability, administration, and high coverage of the acellular vaccine reported in the country as a whole, pertussis cases in the United States over the past decade have increased to record numbers not seen since the 1950s. Many states have experienced epidemic levels of pertussis cases in recent years, specifically California in 2010 and Washington and Vermont in 2012 (7,–11). In 2010, California reported over 9,000 cases (23.4 cases per 100,000 residents), while in 2012 Vermont reported over 600 cases (103 cases per 100,000 residents), over 10 times more than average for that time of year (7, 10). With incidence being highest in unvaccinated infants younger than 6 months and fully vaccinated preadolescents (7 to 10 years) in California, previous studies concluded that unprotected infants were at the highest risk for disease, while waning protection from the childhood acellular vaccine series contributed to disease in preadolescents (7, 12). While there is no single explanation for this resurgence of cases, it has been ascribed to many factors such as pathogen adaptation, improved surveillance and laboratory diagnostics, and waning protection and immune response provided by the acellular vaccine (12,–14).

Since 2010, the U.S. population of circulating B. pertussis strains has increasingly become deficient in Prn, an acellular vaccine component, due to various mutations, primarily the disruption of prn by the mobile genetic element IS481 (11, 14,–18). In 2012, 92% of isolates collected in Vermont did not produce pertactin, suggesting that a selective advantage to prn deficiency may have played a role in that epidemic (14). In separate global locations, two isolates lacking Ptx have been obtained recently, one of which is also Prn deficient (19, 20). The loss of vaccine immunogens identified by current molecular typing methods provides strong evidence that current populations of B. pertussis no longer reflect the genotypic profile of vaccine strains and could be adapting to the immune response elicited by the acellular vaccine (11, 14, 15, 19, 21,–24). This, along with data supporting higher rates of evolution in vaccine antigen genes, has generated concern that other genes encoding vaccine components may harbor mutations not identifiable through current molecular typing techniques (25). Additionally, the recent analysis of a statewide pertussis epidemic in Washington during 2012 revealed that pulsed-field gel electrophoresis (PFGE), a whole-genome restriction digest analysis, was the most powerful indicator of diversity in this population of circulating pertussis strains (11). With evidence suggesting that current strains exhibit genomic diversity that nucleotide variation identified by short-read sequencing fails to explain, it is of great interest to employ long-read technologies to detect mutations and assess the role of genome structure in pertussis resurgence (11, 15, 26).

With hundreds of insertion sequences (ISs) and repeat regions in current circulating strains of B. pertussis, whole-genome sequencing is challenging. There are currently over 400 B. pertussis draft genomes available in public databases, all of which include at least 200 contigs assembled from short-read sequencing. Such fragmented assemblies provide limited information about genome arrangement or gene order and are, therefore, not suitable for evaluating contributions of rearrangement diversity toward pertussis incidence, resurgence, or strain evolution. Complete, de novo genome assembly of B. pertussis is essential in pursuit of this information and now possible with long-read (>1-kb) sequencing and whole-genome mapping (27, 28). The recent release of complete genomes highlights the contrast between modern pertussis strains and current vaccine references, at both the nucleotide and the structural levels (29, 30). In an effort to better understand pertussis resurgence and evolution, here we characterize genome structure and provide direct evidence of large genome rearrangements in annotated, single-contig genome assemblies from 31 isolates collected during two geographically and temporally distinct pertussis epidemics. This depth of resolution facilitates analysis that is not possible through conventional molecular typing or short-read sequencing. The availability of these complete genomes will hopefully aid future transcriptomics and proteomics pursuits toward a deeper understanding of pertussis reemergence and increased incidence.

RESULTS

Large-scale genome rearrangements.

Multiple insertions, deletions, and inversions were seen in all genomes of epidemic isolates relative to E476 and C393 using both optical mapping and whole-genome sequence alignment (Fig. 1). Such genome structural diversity was also observed within epidemic isolates, which included 16 discrete, large-scale architectural profiles. The 13 California genomes comprised 10 rearrangement profiles, and 11 profiles were seen in the 18 Vermont genomes. Further, only five structural profiles were observed in both states, with five and six additional profiles distinct to California and Vermont, respectively (examples of these structural profiles can be found in Fig. S1B in the supplemental material).

Open in a separate window

FIG 1 

Large-scale genome rearrangements within a subset of California (H561 and H563, blue) and Vermont (I496 and I707, red) epidemic isolates compared to vaccine strains E476 and C393 as visualized by whole-genome restriction mapping and alignment with MapSolver (A) and genome sequence alignment with progressiveMauve (B). Connecting lines indicate either conserved restriction fragments (A) or homologous sequence blocks (B). Sequence maps begin at the replication origin, and the approximate replication terminus is indicated with an arrow.

Figure S1 

Genome structural variation among epidemic isolates from California (blue) and Vermont (red). (A) Maximum-likelihood phylogenetic reconstruction based on the order and orientation of conserved, colinear genome fragments identified by whole-genome alignment. (B) Alignment of 14 select genomes with different structural architectures (highlighted with asterisks in panel A). Connected blocks of the same color represent conserved regions, and flags indicate the presence of an IS481 (black), IS1002 (blue), or rRNA operon (red) at the boundary between colinear regions. Download Figure S1, TIF file, 0.4 MB.

Copyright © 2016 Bowden et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

To visualize the dynamic relationship between the vaccine strains and isolates collected during the California and Vermont epidemics resulting from genome structural changes, a Maximum Likelihood for Gene Order (MLGO) tree was constructed from a multiple sequence alignment (Fig. 2). This clustering revealed eight major clades (A to H) that distinguished the epidemic isolates and vaccine strains based on gene rearrangement and gene order (Fig. 2). Isolates with the same or similar PFGE profiles clustered by genome rearrangement (Fig. 2). The vaccine strains formed a distinct clade (H) which is distant from the U.S. epidemic isolates and B1917 (CP009751) and B1920 (CP009752) (29). B1917, the ptxP3 lineage strain corresponding to current circulating strains, clusters with isolates in clade D, while B1920, the ptxP1 lineage strain predominant from 1960 to the 1990s, is no more closely related to I646 and the only ptxP1 epidemic isolate (H375) than it is to the rest of the epidemic and vaccine strain genomes (Fig. 2) (29). Four clades (B, D, E, and F) show clustering of epidemic isolates by state and prn production status. Further, the two isolates collected in Vermont during 2011 (H788 and I669) are diverse and not closely related to each other by genome structure (Fig. 2).

Open in a separate window

FIG 2 

Maximum likelihood for gene order (MLGO) clustering of epidemic strains from the 2010 California (blue) and 2012 Vermont (red) statewide pertussis epidemics. Publicly available genomes (Tohama I, CS, B1917, and B1920) and two resequenced vaccine strains (C393 and E476) were also included (black). Pertactin production and deficiency of each strain are indicated by pink and black boxes, respectively. Clades referenced within the text are indicated by letters. The inset displays banding patterns of the PFGE profiles of all epidemic and vaccine strains sequenced in this study. The dendrogram of PFGE profiles was created using the unweighted pair group method using average linkages (UPGMA) with 1% band tolerance and optimization settings. Genomes for the epidemic isolates I488 and I517 were excluded from the analysis (see Materials and Methods and Data Set S1 in the supplemental material).

Rearrangement boundaries.

To better understand what mediates observed rearrangements, the genome content at each predicted rearrangement boundary was investigated further. When comparing California and Vermont isolates alone, 26 conserved homologous blocks (≥1.5 kb) were identified in these 31 genomes, separated by 25 predicted rearrangement boundaries (Table 2; see also Fig. S1 in the supplemental material). Eighty percent (n = 20) of these boundaries were composed of an IS481 element (Table 2; see also Fig. S1B). Additionally, 12% (n = 3) were composed of an rrn operon, while two boundaries were composed of either an IS1002 element or a combination of IS1002 and IS481 (Table 2; see also Fig. S1B). When two vaccine strains (C393 and E476) were included in the alignment, all additional predicted rearrangement boundaries contained IS481 (Table 2). This larger comparison identified 45 conserved homologous blocks (≥1.5 kb), separated by 44 predicted rearrangement boundaries, 89% (n = 39) of which included copies of IS481 (Table 2).

TABLE 2 

Characteristics of rearrangement boundaries

Characteristic	Value for strain type:
	Epidemic	Epidemic and vaccine
No. of genomes	31a	33a
No. of rearrangement boundaries	25	44
No. of copies/genome
IS481 (≥231)	20	39
rRNA (3)	3	3
IS1002 (4)	1	1
IS1002/IS481	1	1

Open in a separate window

^aI488 and I517 excluded from analysis.

Variant analysis.

The phylogenetic relationships among epidemic isolates were reconstructed from a concatenated alignment of 408 variable sites (Fig. 3). There was greater diversity among the 13 California isolates (Fig. 3), with H375 being most divergent, and no clustering was observed among isolates recovered from fatal cases (Fig. 3). In contrast, Vermont isolates exhibited more homogeneity (Fig. 3). Additionally, four California genomes (H563, H561, H374, and H564) clustered among the earliest collected Vermont genomes, three of which (H788, I475, and I669) were collected from unvaccinated patients (Fig. 3). Pertactin-deficient and pertactin-producing isolates are dispersed throughout the tree (Fig. 3). The two isolates with the prn promoter disruption (I538 and I539) are exclusive to one clade (Fig. 3; Table 1).

Open in a separate window

FIG 3 

Maximum likelihood phylogenetic reconstruction from the concatenated alignment of 408 variable sites (SNVs and MNVs). Vaccine strains resequenced in this study are indicated in black, California isolates are indicated in blue, and Vermont isolates are indicated in red. Yellow circles denote pertactin-deficient strains, and black circles denote strains recovered from fatal cases. The scale bar indicates number of substitutions per site and has been corrected for ascertainment bias based on the nucleotide composition of invariant sites in E476. All epidemic genomes were included in this analysis.

To characterize putative phenotypic effects of observed nucleotide variants, all mutations were classified as noncoding, synonymous, or nonsynonymous, as seen in Fig. 4. Within the 33 epidemic genomes analyzed, a total of 731 variants were discovered and 39.4% (n = 288) were identified as nonsynonymous (Fig. 4A; see also Data Set S3 in the supplemental material). Of these, 34 unique protein-encoding genes were mutated in at least 26 of the 33 epidemic genomes and 17% were found to be involved in inorganic ion transport and metabolism (Fig. 4B; see also Data Set S4). Epidemic genomes were divided by state and characterized in the same manner to identify state-specific mutations. California epidemic genomes contained more variants than Vermont (561 and 413, respectively), with a majority of noncoding mutations (82%) specific to California genomes and a majority of nonsynonymous mutations (53%) specific to Vermont epidemic genomes (Fig. 4A; see also Data Sets S5, S6, S7, and S8). The California isolate H375, the only isolate from the ptxP1 lineage, exhibited the largest number of isolate-specific variants of all epidemic genomes (n = 90) and accounted for seven additional noncoding variants upon exclusion from the California-specific variant analysis (Fig. 4A; see also Data Set S3). The genes affected by nonsynonymous mutations specific to each state were further identified (Table 4). In at least nine of the 12 California epidemic genomes, 63% (n = 5) of the nonsynonymous mutations were seen in five IS1663 elements (Table 4; see also Data Set S9). Nine of the 18 nonsynonymous mutations specific to Vermont affected proteins involved in various metabolic pathways (Table 4; see also Data Set S9). Additionally, no variants were found to be exclusive to prn-deficient isolates outside the prn gene or isolates collected from fatal cases.

Open in a separate window

FIG 4 

Sequence variants within epidemic genomes. (A) Distribution of noncoding (green), synonymous (yellow), and nonsynonymous (blue) variants. (B) Functional classification of predicted protein sequences with observed nonsynonymous mutation in at least 26 epidemic genomes categorized according to EggNOG v.4. Variants were statistically significant by Fisher's exact test.

TABLE 4 

State-specific nonsynonymous variants in California and Vermont epidemic B. pertussis genomes^a

Nonsynonymous variant in state	Functional category	Location in E476 (Tohama I)	Variant within gene	Amino acid change	No. of genomes
California
Transposase for IS1663	Replication, recombination, and repair	116940	947 T→A	316 Leu→Gln	13
Transposase for IS1663	Replication, recombination, and repair	527395	805 T→C	269 Tyr→His	12
Transposase for IS1663	Replication, recombination, and repair	1080959	947 T→A	316 Leu→Gln	9
Transposase for IS1663	Replication, recombination, and repair	1090325	947 T→A	316 Leu→Gln	10
Transposase for IS1663	Replication, recombination, and repair	1447860	947 T→A	316 Leu→Gln	12
Transposase	Replication, recombination, and repair	3457234	947 A→T	316 Gln→Leu	12
Regulator (hypothetical protein)	Signal transduction mechanisms	3614341	11 G→A	4 Gly→Glu	12
Vermont
Zinc transporter ZupT	Inorganic ion transport and metabolism	198405	194 T→C	65 Val→Ala	20
FAD/FMN-containing dehydrogenase	Energy production and conversion	216761	19 T→C	7 Ser→Pro	20
RNase E	Translation, ribosomal structure, and biogenesis	489152	2480 C→T	827 Pro→Leu	20
Endoribonuclease l-PSP (RutC family protein)	Function unknown	528567	186 C→G	62 Asp→Glu	20
Diacylglycerol kinase	Lipid transport and metabolism	868389	97 G→A	33 Asp→Asn	19
Glutamate-ammonia-ligase adenylyltransferase	Posttranslational modification, protein turnover, chaperones	1333208	2519 G→A	840 Gly→Asp	20
Aldo-ketoreductase	Function unknown	1388597	625 G→C	209 Ala→Pro	19
Multidrug efflux pump subunit AcrB	Inorganic ion transport and metabolism	2208591	3221 A→G	1074 Asn→Ser	19
Acyltransferase	Function unknown	2271460	24 G→C	8 Gln→His	17
Hypothetical protein	Function unknown	2276856	158 T→C	53 Val→Ala	20
Transposase	Replication, recombination, and repair	2405180	899 C→G	300 Thr→Ser	19
Potassium-transporting ATPase subunit C	Inorganic ion transport and metabolism	2640489	245 T→C	82 Ile→Thr	19
5-Hydroxyvalerate dehydrogenase (HVD)	Energy production and conversion	3011513	1016 C→T	339 Ala→Val	20
NADPH dehydrogenase	Energy production and conversion	3072916	55 A→G	19 Lys→Glu	20
ADP-dependent (S)-NAD(P)H- hydrate dehydratase	Carbohydrate transport and metabolism	3538308	166 G→A	56 Gly→Arg	11
Putative oligopeptide transporter	Function unknown	1189451-1189452	1765_GC_1766	589 Gly fs	19
ATP-cobalamin adenosyltransferase	Coenzyme transport and metabolism	3304784-3304785	545_CCGG_546	183 Ala fs	17
Hypothetical protein	Function unknown	3490857-3490858	1130_CTA_1131	377 Glu delins Asp*	20

Open in a separate window

^aAbbreviations: fs, frameshift; delins, deletion/insertion; FAD, flavin adenine dinucleotide; FMN, flavin mononucleotide; PSP, perchloric acid-soluble protein; *, produces premature stop codon.

Data Set S3

Seven hundred thirty-one variants in 33 B. pertussis epidemic genomes. Download Data Set S3, XLSX file, 0.1 MB.

Copyright © 2016 Bowden et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

Data Set S4

Thirty-four nonsynonymous variants in at least 26 of the 33 B. pertussis epidemic genomes. Download Data Set S4, XLSX file, 0.1 MB.

Copyright © 2016 Bowden et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

Data Set S5

Five hundred sixty-one variants in California B. pertussis epidemic genomes. Download Data Set S5, XLSX file, 0.1 MB.

Copyright © 2016 Bowden et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

Data Set S6

Four hundred thirteen variants in Vermont B. pertussis epidemic genomes. Download Data Set S6, XLSX file, 0.1 MB.

Copyright © 2016 Bowden et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

Data Set S7

Forty-nine California state-specific B. pertussis epidemic genome variants. Download Data Set S7, XLSX file, 0.1 MB.

Copyright © 2016 Bowden et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

Data Set S8

Thirty-four Vermont state-specific B. pertussis epidemic genome variants. Download Data Set S8, XLSX file, 0.1 MB.

Copyright © 2016 Bowden et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

Data Set S9

State-specific nonsynonymous B. pertussis epidemic genome variants. Download Data Set S9, XLSX file, 0.1 MB.

Copyright © 2016 Bowden et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

To address the issue of vaccine-driven evolution and to determine if certain genes are under selective pressure, the frequencies of nonsynonymous variants and variants within virulence genes across all epidemic genomes were determined. Virulence gene variants had a rate of 0.22 bp/kb compared to 0.24 bp/kb for total epidemic genome variants, with no significant difference between the two. However, nonsynonymous variants had a higher frequency (0.10 bp/kb) than synonymous variants (0.03 bp/kb) (P value = 0.0029).

DISCUSSION

This study illustrates the first effort to evaluate relationships among isolates within and between statewide pertussis epidemics through comparison of complete genome sequences. Prior to the availability of single-contig genomes, PFGE has been the only whole-genome indicator used to fully assess molecular epidemiology and patient transmission of circulating strains of B. pertussis (34). To date, analyses of gene order and genome structure heterogeneity have been limited by the resolution of available methods (35, 36). Additionally, few genomes are available with a complete genome assembly that has been confirmed with a second method such as optical mapping. During the process of complete genome assembly, genome optical maps not only provided structural confirmation of assemblies but proved to be useful as a possible genome typing tool with greater detail than that found with PFGE (with approximate location of restriction cut sites mapped). Depending on future questions, it may be beneficial to develop a hierarchical approach, in which optical mapping is first used to compare genome structures as a typing method, followed by genome sequencing if more detail is required.

Here, we generated high-quality genome assemblies and showed their capacity for detecting discrete genome structural variations among two geographically and temporally independent pertussis epidemics at a resolution never seen before. Furthermore, clustering these genomes according to rearrangement patterns revealed correlations with PFGE profiles and associating prn mutations. However, no strong association was observed between global genome structure and geography, patient fatality, vaccination status, or time of collection (Fig. 2; Table 1). The genome structures of two vaccine strains resequenced as part of this study differed greatly from epidemic isolates, supporting a growing body of evidence that suggests that current circulating strains of B. pertussis have dramatically changed since the isolation of vaccine references (Fig. 2).

Recently, Bart et al. (29) sequenced complete genomes of representative isolates from two predominate global lineages of ptxP3 and ptxP1 (B1917 and B1920, respectively) which exhibited structural architectures similar to those of the U.S. epidemic isolates reported here with comparable molecular typing profiles (Table 1; Fig. 2). The finding that B1920 clustered near the only ptxP1 epidemic isolate may indicate a shift in rearrangement profiles that corresponds to the shift from ptxP1 and ptxP3 lineages. Additionally, B1917 and B1920 strains were isolated in 2000 in the Netherlands, temporally and geographically distant from the epidemic isolates sequenced here, suggesting that rearrangement profiles may be stable over time. Although rearrangements have been identified by others, attempts to measure rearrangement rates in B. pertussis have been limited by the low resolution of PFGE, and reports of gene order changes following multiple laboratory passages, which differ from natural infection cycles, are conflicting (36,–40). Therefore, the historical pattern of genome rearrangement remains unclear. Only through methods like optical mapping and complete genome sequencing, which have been applied to B. pertussis U.S. epidemic isolates here for the first time, can these questions be addressed.

Although a diverse number of PFGE profiles among isolates are observed in this and previous studies, the genomic variation underlying this diversity can now be interrogated (11). To do so, we identified genes located at predicted rearrangement boundaries in all fully assembled genomes in this study. The highly conserved rrn operons and repetitive mobile elements IS481 and IS1002 were found at predicted rearrangement breakpoints in both epidemic and vaccine isolate genomes, the majority of which were composed of IS481 elements (Table 2; see also Fig. S1 in the supplemental material). Almost all predicted rearrangements were in the form of inversions flanked by inverted sequence repeats, often “symmetric” around the origin or terminus of replication, maintaining replichore balance important for genome stability (41, 42). These data support previous evidence that rearrangements in the B. pertussis genome are mediated by insertion sequence (IS) elements through homologous recombination between inverted repeats, presumably during replication, a mechanism long thought to maintain genome plasticity (22, 35, 36, 41, 43, 44).

Aside from rearrangements, mobile elements have been shown to play a role in the creation of pseudogenes, genome reduction, and adaptive evolution toward host specificity (22, 36, 43, 45, 46). Accurate identification of all IS481, IS1002, and IS1663 insertions is necessary to evaluate their contribution to genome evolution, and previous studies using short-read sequencing have been unsuccessful (24, 47). Our pipeline provided more robust sequencing data for proper placement of insertion sequences to construct more accurate assemblies, evident when comparing the corrected number of IS481 elements discovered in the resequenced vaccine strains E476 and C393 and assembly structure errors found when comparing the NCBI CS reference genome with the resequenced C393 (see Data Set S1 in the supplemental material). Alignment of complete genome sequences provides further evidence that IS481 transposition remains active, consistent with previous reports, but also suggests that this expansion of repetitive sequence may provide additional opportunities for rearrangement (48). With the exception of H561, most Vermont and California epidemic genomes contained more IS481 copies than both Tohama I and E476. Additionally, further acquisition of insertion sequence copies in currently circulating strains of B. pertussis may indicate intraspecies host adaptation through gene inactivation that allows B. pertussis to evade protective immunity elicited by the acellular vaccine and cause disease in a vaccinated population. This claim is further supported by studies indicating that the majority of prn-deficient U.S. isolates harbor an IS481 within the prn gene locus (11, 15). However, the results of this study suggest that the role of mobile elements in B. pertussis genome evolution may not be limited to reduction, raising questions about the potential fitness implications of rearrangement, as others have also speculated (40).

In an effort to identify potential correlates of protection not included in the U.S. acellular vaccine, nonsynonymous mutations were investigated in 40 additional virulence-related genes and throughout the genomes of all epidemic isolates. Many virulence-related gene mutations were identified in isolates collected from infants and children with various effects on amino acid sequence within either a single isolate or all epidemic isolates (Tables 1 and ​and3)3) (23, 47, 49,–56). Although statistically insignificant due to the small sample size of this study, rates of SNPs in virulence genes seem to support the idea of genetic drift and are not found to occur more frequently than in other genes. However, the rate of nonsynonymous variants suggests that genetic elements might be under selective pressure. The majority of nonsynonymous mutations (53%) were found to affect proteins either identified as transporters or associated with the transport and metabolism of various cellular molecules (Fig. 4B; see also Data Set S4 in the supplemental material). This overrepresentation of SNPs in transport proteins was also seen in a study of SNP density in the global population of pertussis isolates and may provide insight into the mechanism by which current strains have adapted to the current global vaccination state (24). Although state-specific mutations were discovered in IS elements and metabolic proteins in California and Vermont, respectively, it is difficult to make phenotypic implications based on state specificity, disease presentation, or demographic due to the highly biased selection of the small number of isolates chosen for whole-genome sequencing and the lack of a functional assessment in the study. Transcriptomics, protein expression, and functional assays of proteins affected by nonsynonymous mutations in these genomes are needed to fully understand how these mutations affect the ability of B. pertussis to infect and cause disease, with varying severity, in the host.

It is clear that we do not fully understand the correlates of protection against B. pertussis. However, this study has expanded our understanding of prn deficiency, specifically with genomic data visualizing the true nature of the previously predicted promoter inversion and providing evidence that no conserved single or multiple mutational events may be compensating for prn deficiency in this collection of isolates, calling into question the role that pertactin plays in pertussis disease (11, 15). This analysis in conjunction with proteomics investigation of several of these isolates, which we are actively pursuing, provides a starting point for determining what proteins are required for disease in certain host populations while identifying potential candidates for future vaccine components.

Advances in next-generation sequencing technology, whole-genome mapping, and bioinformatics improve access to complete genomes, which allow comparative analyses at both the nucleotide and the structural levels. In this study, we compared 31 complete genomes from two recent statewide B. pertussis epidemics with vaccine strains and observed the first definitive view of genome structural variation through rearrangement at a nucleotide sequence resolution within this species. These data challenge the previously held view of population clonality, revealing new levels of diversity, even within geographically defined epidemics. As more complete genomic, transcriptomic, and proteomic data are made available for additional strains collected across broader time periods, it will become possible to discern the stability and rate of rearrangements, when certain structural profiles emerged, and whether profiles shift between epidemic and nonepidemic years. Such further analyses are needed to determine how rearrangement events correlate with SNP phylogenetic relationships and aid development of a comprehensive typing methodology that incorporates nucleotide and structural variation to draw meaningful conclusions from associated clinical metadata. It remains unclear if and how rearrangements play a role in adaption and virulence. Genomes of isolates collected from asymptomatic carriers are needed as a baseline, in conjunction with transcriptomic and proteomic analysis, to determine what insertion elements or rearrangements, if any, are specific to increased disease severity in some patients. Although this study cannot link rearrangements to virulence or adaptation, it highlights the need for complete assemblies to fully evaluate the contribution of pathogen evolution, in addition to waning immunity and incomplete protection from current vaccines, toward pertussis reemergence.

MATERIALS AND METHODS

ACKNOWLEDGMENTS

Funding Statement

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

REFERENCES