The BvgAS Regulon of Bordetella pertussis.

Nearly all virulence factors in Bordetella pertussis are activated by a master two-component system, BvgAS, composed of the sensor kinase BvgS and the response regulator BvgA. When BvgS is active, BvgA is phosphorylated (BvgA~P), and virulence-activated genes (vags) are expressed [Bvg(+) mode]. When BvgS is inactive and BvgA is not phosphorylated, virulence-repressed genes (vrgs) are induced [Bvg(-) mode]. Here, we have used transcriptome sequencing (RNA-seq) and reverse transcription-quantitative PCR (RT-qPCR) to define the BvgAS-dependent regulon of B. pertussis Tohama I. Our analyses reveal more than 550 BvgA-regulated genes, of which 353 are newly identified. BvgA-activated genes include those encoding two-component systems (such as kdpED), multiple other transcriptional regulators, and the extracytoplasmic function (ECF) sigma factor brpL, which is needed for type 3 secretion system (T3SS) expression, further establishing the importance of BvgA~P as an apex regulator of transcriptional networks promoting virulence. Using in vitro transcription, we demonstrate that the promoter for brpL is directly activated by BvgA~P. BvgA-FeBABE cleavage reactions identify BvgA~P binding sites centered at positions -41.5 and -63.5 in bprL Most importantly, we show for the first time that genes for multiple and varied metabolic pathways are significantly upregulated in the B. pertussis Bvg(-) mode. These include genes for fatty acid and lipid metabolism, sugar and amino acid transporters, pyruvate dehydrogenase, phenylacetic acid degradation, and the glycolate/glyoxylate utilization pathway. Our results suggest that metabolic changes in the Bvg(-) mode may be participating in bacterial survival, transmission, and/or persistence and identify over 200 new vrgs that can be tested for function.IMPORTANCE Within the past 20 years, outbreaks of whooping cough, caused by Bordetella pertussis, have led to respiratory disease and infant mortalities, despite good vaccination coverage. This is due, at least in part, to the introduction of a less effective acellular vaccine in the 1990s. It is crucial, then, to understand the molecular basis of B. pertussis growth and infection. The two-component system BvgA (response regulator)/BvgS (histidine kinase) is the master regulator of B. pertussis virulence genes. We report here the first RNA-seq analysis of the BvgAS regulon in B. pertussis, revealing that more than 550 genes are regulated by BvgAS. We show that genes for multiple and varied metabolic pathways are highly regulated in the Bvg(-) mode (absence of BvgA phosphorylation). Our results suggest that metabolic changes in the Bvg(-) mode may be participating in bacterial survival, transmission, and/or persistence.

INTRODUCTION

Despite a high vaccination rate, whooping cough, caused by the human pathogen Bordetella pertussis, has resurged for the last 2 decades in many industrialized countries, including the United States and those in Western Europe (1–3). This increasing incidence correlates with a switch from whole-cell pertussis (wP) to acellular component pertussis (aP) vaccines in the mid-1990s (reviewed in references 4 and 5). Some studies have suggested that this may be due, at least in part, to attributes of aP vaccines regarding the length and type of immunity that they engender (6, 7). In baboon challenge studies, aP vaccine, while protecting against disease, was less efficacious than wP vaccine in preventing colonization with or the secondary transmission of B. pertussis from infected transmitters/carriers (8). As efforts turn to developing an improved pertussis vaccine, it is therefore crucial to fully understand how B. pertussis survives within different niches inside the host as well as how it is transmitted from one host to another.

The two-component system BvgAS, comprised of the histidine kinase BvgS and the response regulator BvgA, is the master regulator of B. pertussis virulence (9, 10). This system is coordinated by the intracellular level of phosphorylated BvgA (BvgA~P) (11, 12). It can be modulated by changing growth conditions, which can be divided into three modes [Bvg(+), Bvg(i), and Bvg(−)] (13; reviewed in references 14 and 15). Recent work suggests that the transition between these modes occurs via specific BvgA conformational changes within its periplasmic Venus flytrap domains (VFT1 and VFT2) (16, 17). The nonmodulated mode, Bvg(+), is observed when B. pertussis is grown at 37°C under standard conditions, leading to the full activation of vags (virulence-activated genes) (reviewed in reference 18). Under these conditions, the kinase and phosphotransfer activities of BvgS are high. Therefore, this is considered a default “on-state” since BvgA is phosphorylated without any known chemical stimuli (19). Phosphorylation of BvgA, which is absolutely required for adaptation and colonization in respiratory tracts of experimental animals, results in BvgA~P dimerization and transcriptional activation of key vag promoters (20, 21). In an experimental scenario in which BvgA~P levels are rising, early vags, such as adhesin genes for filamentous hemagglutinin (fha) and fimbriae (fim3 and fim2), are activated first. At this time, BvgA~P also activates the promoter of the early vag bvgR, whose product represses vrgs (virulence-repressed genes) that are induced in the absence of BvgA~P (22). As BvgA~P levels increase, late vags are activated, such as toxin genes cya (adenylate cyclase) and ptx (pertussis toxin) (23). This differential control arises because early vag promoters have fewer and stronger BvgA~P binding sites than late vag promoters. The Bvg(i) mode represents an intermediate state, with an intermediate concentration of BvgA~P where kinase-on and kinase–off BvgS proteins may coexist in equilibrium. It is exemplified by the bipA gene, which is activated by levels of BvgA~P that activate early genes but is repressed by the higher levels required for late gene activation (24). The Bvg(−) mode represents very low levels of BvgA~P and is induced by the addition of modulators such as MgSO4 or nicotinic acid or by growth at low temperatures. It has been proposed that in the presence of modulators, conformational and dynamic changes of the periplasmic domain of BvgS convert BvgS from a kinase-on to a kinase-off state (16, 17). Under these conditions, phosphorylation of BvgA is rapidly lost, BvgR is not produced, and multiple vrgs are expressed (11; reviewed in reference 18).

Despite extensive study, the role of BvgAS regulation in the natural history of B. pertussis is not completely clear. For Bordetella bronchiseptica, a close relative of B. pertussis that has a broad mammalian host range, the Bvg(−) mode has been shown to contribute to survival outside the host environment in laboratory experiments (25, 26). Although B. bronchiseptica has never been cultured from natural environments, this provides at least a plausible rationale for BvgAS-mediated regulation in this species. B. pertussis, on the other hand, does not survive for long periods outside the host respiratory tract, and its natural host range is limited to humans. In addition, the B. pertussis genome has undergone significant contraction and gene inactivation since its evolution from a B. bronchiseptica-like common ancestor, leading one to think of it as a “degraded” species (27). It has been speculated that the Bvg(i) and Bvg(−) modes may represent useless “evolutionary remnants” and that the BvgAS system may represent an evolutionary “fossil” in B. pertussis, necessary only for its role in activating virulence genes but not utilized to adjust virulence gene expression in response to environmental signals (9, 22, 28). Another point of view is that the BvgAS system, together with the genes expressed in the Bvg(−) and Bvg(i) modes, plays a role in B. pertussis different from that in B. bronchiseptica, such as in transmission by the aerosol route (29). In this regard, Bvg(−) mutants have been observed to accumulate among the B. pertussis population in the rhesus monkey nasopharynx, demonstrating that cells in the Bvg(−) mode may not be immediately cleared during infection (30). Given this ambiguity, it is important to definitively identify vags and vrgs within B. pertussis to fully investigate possible functions of these genes.

Global searches in Bordetella species have identified numerous BvgA-regulated genes. Microarray analyses identified ~250 genes regulated in B. pertussis Tohama I BP536 and nearly 500 genes regulated in B. bronchiseptica RB50 (9, 31). Known vags, including those encoding adhesins and toxins, were expressed consistently in both species in the Bvg(+) mode, while more species-specific differences were observed in the Bvg(−) mode. Recently, Ahuja et al. reported the first transcriptome sequencing (RNA-seq) analysis in a Bordetella species, but it was limited to B. bronchiseptica (32). Here, we report the first RNA-seq examination of the BvgAS regulon in B. pertussis. Our work identifies more than 550 genes within the BvgAS regulon with 245 transcripts positively regulated and 326 transcripts negatively regulated by BvgA~P. Our work suggests that B. pertussis has hierarchical regulatory networks under the BvgA phosphorelay. Positively regulated genes include many transcriptional regulators and genes within the type 3 secretion system (T3SS). We demonstrate that expression of the extracytoplasmic function (ECF) sigma factor brpL, which is needed for T3SS expression and is significantly upregulated in the Bvg(+) mode, is directly activated by BvgA~P. In addition, we find that multiple genes involved in metabolic pathways are negatively regulated by BvgAS and thus upregulated in the Bvg(−) mode. Our results are consistent with the idea that metabolic changes in the Bvg(−) mode may be participating in bacterial survival, transmission, and/or persistence.

RESULTS

Analysis of the B. pertussis BvgAS-dependent regulon.

RNA-seq analyses were performed to compare transcriptomes resulting from the presence and absence of BvgA~P by isolating RNA from B. pertussis Tohama I BP536 (wild type [WT]) and a ΔbvgAS derivative. As detailed in Materials and Methods, cells were grown under either modulating (50 mM MgSO4) conditions, in which phosphorylated BvgA is undetectable [Bvg(−) mode], or nonmodulating (no MgSO4) conditions [Bvg(+) mode], in which BvgA is present at higher concentrations and is primarily in the form of BvgA~P (11). Illumina sequencing data were analyzed using CLC Genomic Workbench after normalization with total reads, comparing gene expression of the following experimental sets: set 1, WT with or without MgSO4; set 2, WT versus ΔbvgAS mutant in the absence of MgSO4; and set 3, ΔbvgAS mutant with or without MgSO4 (Table S1A). At least a 1.6-fold change in gene expression with less than an 0.05 false discovery rate (FDR) P value after an empirical digital gene expression (DGE) statistical test in sets 1 and 2 was scored as significant. In addition, set 3 was used to eliminate genes that were affected by MgSO4 in the absence of bvgAS.

(A) Entire set of RNA-seq data. (B) Genes positively regulated by BvgA~P. Genes whose expressions have an FD of ≥1.6 with a P value of <0.05 are considered to be significant. Genes in blue are newly identified in this study; the gene in red is not Bvg(+) as indicated by RT-qPCR. Genes highlighted in yellow were assessed by RT-qPCR. 1Possible function identified by PSI-BLAST in this study. 2Possible function assigned in the work of Ahuja et al. (32). A total of 245 genes are identified as Bvg(+), 237 by RNA-seq and 8 from RT-qPCR. One hundred twenty-one genes are newly identified as Bvg(+) in this study. (C) Genes negatively regulated by BvgA~P. Genes whose expressions have an FD of ≥1.6 with a P value of <0.05 are considered to be significant. Genes in blue are newly identified in this study; the gene in red is not Bvg(−) as indicated by RT-qPCR. Genes highlighted in yellow were assessed by RT-qPCR. 1Possible function identified by PSI-BLAST in this study. A total of 326 genes are identified as Bvg(−), 321 by RNA-seq and 5 from RT-qPCR. Two hundred thirty-two genes are newly identified as Bvg(−) in this study. Download TABLE S1, XLSX file, 6.6 MB.

Tables S1B and C in the supplemental material list the positively and negatively BvgA-regulated genes, respectively, that passed our criteria. Note that (i) in a few instances, genes that did not meet our criteria were included for subsequent reverse transcription-quantitative PCR (RT-qPCR); (ii) most genes in set 3 showed FDR-corrected P values of >0.05, which means that the effect of the addition of MgSO4 was considered to be insignificant; and (iii) tRNAs, rRNAs, and repetitive DNA transcripts were eliminated from our data sets. Overall, mRNA levels of more than 550 genes, or about 15% of open reading frames (ORFs), were affected by the presence or absence of BvgA~P. Of these, 353 are newly identified and were not reported in the previous microarray studies (9, 31). In general, genes known to be directly regulated by BvgA~P, including those encoding virulence factors, had significant changes with more than a 3-fold difference (FD). Although numerous newly identified genes reported in this study were observed to have less significant fold changes (<3 FD), these genes were clearly BvgA dependent since they were affected by the deletion of bvgAS whether MgSO4 was present or absent (set 3).

To verify the expression of many of the genes that we characterized as members of the BvgAS regulon, we performed RT-qPCR analyses (Fig. 1, 2, S1, and S2). RNA was prepared from cultures independently from those used for RNA-seq. Each data set represents three independent cell cultures with at least 2 repeats. Expression levels are indicated by FD with error bars showing the standard deviation.

FIG 1 

Relative expression of selected genes of transcriptional regulators positively regulated by BvgA~P. Values are shown as fold difference using the 2−ΔΔ method (65). Error bars indicate standard deviations among three independent biological replicates for each condition with more than two different trials. Blue bars, FD of WT without MgSO4 versus WT with MgSO4; red bars, FD of WT without MgSO4 versus ΔbvgAS strain without MgSO4. Student’s t test was used for statistical analysis, and genes showing significant differences between WT without MgSO4 and WT with MgSO4 or between WT without MgSO4 and ΔbvgAS strain without MgSO4 are labeled as *** (P < 0.001), ** (P < 0.01), or * (P < 0.05). “ns” represents statistically nonsignificant.

FIG 2 

Relative expression of selected genes involved in metabolic pathways negatively affected by BvgA~P, as measured by RT-qPCR. Values are shown as fold difference using the 2−ΔΔ method (65). Error bars indicate standard deviations among three independent biological replicates for each condition with more than two different trials. Blue bars, FD of WT without MgSO4 versus WT with MgSO4; red bars, FD of WT without MgSO4 versus ΔbvgAS strain without MgSO4. Bracketed genes are located in the same locus of the genome. Student’s t test was used for statistical analysis, and genes showing significant differences between WT without MgSO4 and WT with MgSO4 or between WT without MgSO4 and ΔbvgAS strain without MgSO4 are labeled as *** (P < 0.001), ** (P < 0.01), or * (P < 0.05).

RT-qPCR analyses of genes identified as positively regulated by BvgA~P by the RNA-seq analyses. (A) Well-known vags; (B) T3SS genes; (C) genes of membrane-associated transporter proteins and hypothetical proteins. Values are shown as a fold difference using the 2−ΔΔ method. Error bars indicate standard deviations among three independent biological replicates for each condition with more than two different trials. Blue bar, FD of WT without MgSO4 versus WT with MgSO4; red bar, FD of WT without MgSO4 versus ΔbvgAS without MgSO4. Download FIG S1, EPS file, 4.7 MB.

RT-qPCR analyses of genes identified as negatively regulated by BvgA~P by the RNA-seq analyses. (A) Well-known vrgs; (B) other vrgs, including genes for cold shock membrane proteins and genes with hypothetical functions. (Cold shock genes are bracketed.) Values are shown as a fold difference using the 2−ΔΔ method. Error bars indicate standard deviations among three independent biological replicates for each condition with more than two different trials. Blue bar, FD of WT without MgSO4 versus WT with MgSO4; red bar, FD of WT without MgSO4 versus ΔbvgAS without MgSO4. Download FIG S2, EPS file, 3.1 MB.

As detailed below, many of the BvgAS-regulated genes could be assigned to specific functional classes (positively regulated genes, Tables 1 and S2; negatively regulated genes, Tables S3 and S4). Within these tables, genes in yellow were analyzed by both RNA-seq and RT-qPCR, genes in blue were newly identified in this study, and genes in red (supplemental tables only) were not confirmed by RT-qPCR.

TABLE 1 

Genes of transcriptional regulators positively regulated by BvgA~P

Expression of genes in boldface was checked by RT-qPCR. Underlined genes are newly identified in this study.

Well-known vags, T3SS genes, and genes of membrane-associated transporter proteins and hypothetical proteins. (A) Well-known vags positively regulated by BvgA~P; (B) T3SS genes positively regulated by BvgA~P; (C) genes of membrane-associated transporter proteins and hypothetical proteins positively regulated by BvgA~P. Download TABLE S2, DOCX file, 0.1 MB.

Well-known vrgs and other vrgs, including genes for cold shock and membrane proteins and genes with hypothetical functions. (A) Well-known vrgs; (B) other vrgs, including genes for cold shock and membrane proteins and genes with hypothetical functions. Download TABLE S3, DOCX file, 0.03 MB.

Genes involved in various metabolic pathways negatively regulated by BvgA~P. Download TABLE S4, DOCX file, 0.05 MB.

Genes positively regulated by BvgAS.

Two hundred forty-five genes were identified as positively regulated by BvgA~P (Table S1B). Category 1 (Tables 2 and S2) contains well-known vags, such as virulence determinants and accessory factors, that have been identified by previous genetic and microarray studies (9, 10, 31, 33–36); these include filamentous hemagglutinin/adhesin (fhaB and fhaC), adenylate cyclase toxins (cyaABCDE), fimbriae (fimABCD), pertactin (prn), pertussis toxins (ptxABDEC), and type 4 secretion system proteins (ptlABCDIEFGH). RT-qPCR confirmed the significant increase in RNA levels of several of these vags (Fig. S1A). For example, RNA levels for fhaB, fim2, cyaA, and ptxA increased at least 20-fold in WT without MgSO4 compared to WT with MgSO4. While expression of bipA demonstrated only a 2-FD in WT without MgSO4 versus WT with MgSO4, a dramatic change was observed in the ΔbvgAS mutant in RT-qPCR analyses. We suspect that the low FD for bipA in WT with or without MgSO4 arose from the fact that these samples were taken from each of the extremes of BvgAS activity known to affect its expression, in which low levels of BvgA~P activate the gene while high levels of BvgA~P repress it (12, 37).

TABLE 2 

RNA-seq expression levels of selected well-known vag genes and vrg genes

Expression of genes in bold was checked by RT-qPCR.

FD represents the fold difference in expression of WT without MgSO4 versus ΔbvgAS strain without MgSO4.

Category 2 contains genes encoding transcriptional regulators, including the ECF sigma factor brpL (also called btrS), the σN gene rpoN-1 (BP0133), and several two-component systems, such as kdpDE and BP2934/2935 (Table 1). Expression of brpL was ~100-fold higher in the Bvg(+) mode than in the Bvg(−) mode. RT-qPCR analyses confirmed the BvgA-dependent activation of brpL and many other genes in this set (Fig. 1) as well as several other regulators and putative regulators, which demonstrated less than a 2-FD in our RNA-seq analyses but were previously identified as vags in Bordetella species (Table S5). These included BP0319, BP0498, BP1590, BP2399, BP2878, BP2981, and BP3239 (9). Thus, together with the previous microarray studies, our work reveals that the expression of dozens of transcriptional regulators is upregulated by BvgA~P in B. pertussis. Most of these are newly identified by our analyses.

RT-qPCR analysis for other possible vags and vrgs. Genes were previously identified as vags in the work of Cummings et al. (9) in B. pertussis Tohama I,1 in B. pertussis GMT-1,2 in B. bronchiseptica RB50,3 or in B. bronchiseptica 77.4,5 Genes were previously identified as vrgs in the work of Cummings et al. (9). (A) Possible vag transcriptional regulators, not identified by our RNA-seq; (B) genes reported as vrgs in B. bronchiseptica. Download TABLE S5, DOCX file, 0.02 MB.

Category 3 includes genes within the type 3 secretion system (T3SS) operon (Table S2). As observed before (9), transcripts of the entire T3SS operon are significantly higher in the presence of BvgA~P. RT-qPCR confirmed the significant BvgA-dependent activation of bscF, bapC, bsp22, bopN, and brkA (Fig. S1B).

Finally, category 4 (Table S2 and Fig. S1C) includes genes for membrane-associated transporter proteins and hypothetical proteins with unknown functions. The BvgA-dependent activation of most of these genes is newly identified by this study.

Genes negatively regulated by BvgAS.

A total of 326 genes were negatively regulated by the presence of BvgA~P; 232 of these were newly identified in this study (Table S1C). Category 1 (Tables 2 and S3; Fig. S2A) contained well-known vrgs such as vrg-6, vrg-24, vrg-18, and vrg-73 (35). Their differential expression is consistent with their repression by the product of the bvgR gene, whose expression requires BvgA~P (22, 38). As expected from previous studies, genes within the lipopolysaccharide (LPS) biosynthesis (bpl) and capsular polysaccharide biosynthesis (kpsM-wcbO) loci are in this group (35, 39–41).

The classification of fim3 as a vrg is paradoxical given that the promoters for the fim2 and fim3 genes that encode the Fim2 and Fim3 major fimbrial subunits, respectively, are known to behave as Bvg-activated genes when studied in isolation (20, 42). Given that the fim3 promoter is in the “13C” configuration in the strain used here and thus should not be activated when BvgA~P is present (42), it is reasonable that the level of fim3 RNA is low when BvgA~P is present. However, the results here confirm the surprising finding that the level of fim3 RNA actually increases under growth conditions in which BvgA~P is absent (9). Our recent evidence demonstrates that this arises because of the existence of a strong vrg-like promoter upstream of P and the fim3 gene (unpublished data). Downstream of this promoter is a predicted ORF that we have designated vrgX. Expression of this ORF was not assessed in previous microarray studies as it was not previously annotated as a gene. As seen in Table S3 and Fig. S2, the level of vrgX and fim3 transcripts increased significantly in the Bvg(−) mode, consistent with activation of the vrgX promoter under modulating conditions. However, to date, production of Fim3 protein under these conditions has not been detected.

Category 2 (Fig. 2 and Table S4) includes genes involved in metabolic and catabolic pathways such as fatty acid and lipid metabolism (acyl coenzyme A [acyl-CoA] synthetases/dehydrogenases, enoyl-CoA hydratases/isomerase, and glycerol-3-phosphate dehydrogenase), the tricarboxylic acid (TCA) cycle (sugar and amino acid transporters [BP0120 to BP0122, BP0126, BP0620, and BP0621] and pyruvate dehydrogenase [BP0628/BP0629]), branched amino acid transport systems (BP0301, BP0302, BP0304, and BP0619 to BP0622), the leucine/isoleucine/valine transporter system (BP1273, BP1274, BP1276, and BP1277), the phenylacetic acid catabolic pathway (BP2675 to BP2678, BP2682, and BP2683), peptidoglycan and LPS synthesis genes (glmS [BP0666] and glmU [BP3730]), the glycolate oxidation pathway (glcDEF and glcB), and cyanophycin synthetases (cphA1 and cphA2). The vast majority of genes within category 2 are newly identified by our study. It is not clear why upregulation of these genes was not observed in a previous microarray analysis (9), but differences in the growth medium may be responsible for the dissimilar findings.

Genes for cold shock proteins (BP1770, BP2757, and BP3871) and genes of unknown function, including vrgX, are placed in category 3. Among the three cold shock genes identified in B. pertussis, BP3871 and BP2757 were confirmed to be regulated by BvgA (Table S3 and Fig. S2B). [While the level of the third cold shock gene, BP1770, increased in the Bvg(−) mode in our RNA-seq analysis, this result was not confirmed by RT-qPCR (Fig. S2B).] Taken together with the other metabolic changes observed under the Bvg(−) condition, our observations suggest that B. pertussis undergoes multiple metabolic changes when shifted to modulating conditions. These changes may represent adaptation, even if for a short time, to ensure survival outside the human host, such as in aerosols involved in B. pertussis transmission.

We also used RT-qPCR to quantify the level of transcripts of 5 genes which were reported as upregulated in the BvgA(−) mode in B. bronchiseptica (9) but did not pass our FD cutoff in RNA-seq: BP0966 (sulfate binding protein precursor) and four flagellar synthesis and regulatory genes (fliA, flhD, flhF, and fliO). We found no significant difference in their expression in any of our comparisons (Table S5). The flagellar result is expected as this locus in B. pertussis is disrupted by multiple insertion elements (27).

brpL, which encodes the ECF sigma factor needed for T3SS expression, is directly regulated by BvgA~P.

To investigate whether some of the Bvg(+) genes, whose expression significantly increased in the presence of BvgA~P, are directly regulated by BvgA~P, we performed in vitro transcriptions using B. pertussis RNA polymerase (RNAP) alone or in the presence of BvgA or BvgA~P. The BvgA~P-dependent promoter for fhaB (36, 43) was used as a positive control. Among the newly identified BvgA-regulated genes, we selected BP3441 (phasin-like protein, 7-FD), BP2142 (GntR-like regulator, 6-FD), and BP1005 (hypothetical protein, 5-FD). Among previously identified BvgA-regulated genes, we selected brpL (99-FD), BP2226 (anti-sigma factor antagonist, 7-FD), BP2230 (regulator, 13-FD), and BP2233 (hypothetical protein, 27-FD). DNA regions from approximately −250 to +150 relative to the transcript start from the RNA-seq analysis (details in Materials and Methods) were used as the transcription templates.

As seen in Fig. 3A, P transcription was observed only in the presence of BvgA~P, demonstrating that this promoter, like P, is directly activated by phosphorylated BvgA. Primer extension analysis of the in vitro brpL RNA identified the start site (+1) for in vitro transcription and a reasonable −10 element for a σA-dependent promoter at the expected position of −12 to −7 (Fig. 3B and C). In contrast, transcription of either BP1005 or BP3441 was observed by RNAP alone and was not affected by the presence of BvgA or BvgA~P (Fig. 3A). We assigned promoter sequences for the BP1005 and BP3441 transcripts, based on the size of the detected RNA, and found excellent matches to consensus −35 and −10 elements upstream of the transcription start site (Fig. 3C), consistent with the ability of RNAP alone to transcribe from these promoters without the need of an activator. One possible explanation for the lack of BvgA~P dependence in vitro would invoke the existence of repressors of these genes that are themselves induced under modulating conditions. On the other hand, transcription of BP2142, BP2230, and BP2233 was not detected in vitro (data not shown), suggesting that BvgA~P may indirectly affect expression of these genes and/or that additional transcription factors are required for expression.

FIG 3 

P is directly activated by BvgA~P. (A) Denaturing polyacrylamide gels showing RNA obtained after in vitro transcription of the indicated templates using (as indicated) B. pertussis RNAP alone (−), RNAP plus BvgA (+), or RNAP plus BvgA~P (P). The expected transcripts (nucleotides) are indicated. (B) Primer extension analysis of in vitro RNA to identify the start site of P. A denaturing gel of primer extension products is shown. Left lanes are sequencing ladders; right lanes show primer extension products for RNA obtained in the absence of BvgA (−), the presence of BvgA~P (P), or the presence of nonphosphorylated BvgA (+). (C) Sequence of the brpL promoter, determined from primer extension analyses in panel A, and of the promoters for BP1005 and BP3441, determined from transcript length in in vitro transcription in panel B. Assigned −10 and −35 elements (BP1005 and BP3441) are in blue; +1 is in red.

P contains multiple BvgA~P dimer binding sites.

To determine whether BvgA~P directly binds to sites within P, we performed iron bromoacetamidobenzyl-EDTA (FeBABE) footprinting using BvgAT194C-FeBABE and BvgAV148C-FeBABE (36). As BvgA residue 194 is located near the BvgA~P dimer interface, while residue 148 is located on the outside face of the dimer (36), DNA cleavage by FeBABE conjugated to these residues maps the outside and inside edges of each BvgA relative to the DNA. Our previous work with the BvgA-dependent promoter for fim3 indicated that transcription complexes formed with Escherichia coli RNAP were more stable than those formed with B. pertussis RNAP (20). Consequently, since both RNAPs yielded BvgA~P-dependent transcription from P in vitro (Fig. 3A and S3A), we performed FeBABE footprinting with both RNAPs (B. pertussis, Fig. 4; E. coli, Fig. S3B). In each case, the previously characterized promoter for fhaB, P, served as a positive control. As expected, with either RNAP, we observed the three binding regions for BvgA~P at P, centered at positions −44.5, −66.5, and −88.5. At P, we observed two distinct BvgA~P binding sites centered at −41.5 and −63.5 with E. coli RNAP. Although these sites were again observed with B. pertussis RNAP, overall the cleavage sites were weaker, especially with BvgAT194C-FeBABE, and an additional farther upstream binding site centered at approximately −114 was also seen. These results indicate that P has a similar BvgA~P binding motif as P, but the stability and/or oligomerization pattern of BvgA~P on the DNA appears to differ with the two polymerases.

FIG 4 

FeBABE cleavage sites for P and P. (A) Denaturing gels showing products of FeBABE cleavage using complexes formed with P (left) or P (right), 32P end labeled on either the template or the nontemplate strand, B. pertussis RNAP, and (as indicated) either phosphorylated BvgAT194C-FeBABE or phosphorylated BvgAV148C-FeBABE. The G+A lanes correspond to DNA ladders generated by Maxam-Gilbert sequencing with the positions relative to the transcription start sites indicated (36). Arrowheads denote sites of cleavage. (B) Sequences of P and P with the cut sites generated by BvgAT194C-FeBABE (pink arrowheads) and by BvgAV148C-FeBABE (purple arrowheads). BvgA binding sites are boxed in orange; the −10 regions are shown and boxed in blue. (An alternative −10 element for P is indicated in the dashed purple box.) Transcription start sites (+1) are in red. Numbers above the inverted red arrows indicate the predicted score for BvgA binding half-sites in each promoter according to the algorithm reported by Merkel et al. (67). For instance, the sequence of the perfect inverted repeat is TAGGAAATTTCCTA, whose score is given as “0.”

Cleavage of P by FeBABE-conjugated BvgA~P in transcription complexes made with E. coli RNAP. (A) Portion of denaturing polyacrylamide gels showing P RNA (left) or P RNA (right) obtained after in vitro transcription using, as indicated, E. coli RNAP alone (-), RNAP plus BvgA (+), or RNAP plus BvgA~P (P). (B) Denaturing gels showing products of FeBABE cleavage using complexes formed with P (left) or P (right), 32P end labeled on either the template or the nontemplate strand, E. coli RNA polymerase, and, as indicated, either phosphorylated BvgAT194C-FeBABE or phosphorylated BvgAV148C-FeBABE. The G+A lanes correspond to DNA ladders generated by Maxam-Gilbert sequencing with the positions relative to the transcription start sites indicated. Arrowheads denote sites of cleavage. (C) Sequences of P and P with the cut sites generated by BvgAT194C-FeBABE (pink arrowheads) and by BvgAV148C-FeBABE (purple arrowheads). BvgA binding sites are boxed in orange; the −10 regions are boxed in blue. (An alternative −10 element for P is indicated in the dashed purple box.) Transcription start sites (+1) are in red. Numbers above the inverted red arrows indicate the predicted score for BvgA binding half-sites in each promoter according to the algorithm reported by Merkel et al. (67). Download FIG S3, TIF file, 25.1 MB.

DISCUSSION

The two-component system BvgA (response regulator)/BvgS (histidine kinase) is a master regulator of virulence gene expression in multiple species of Bordetella. In both B. pertussis and B. bronchiseptica, the presence or absence of high levels of BvgA~P activates or represses multiple vags or vrgs, respectively, at multiple locations throughout the genome (9, 10). Among Bordetella spp., B. pertussis is by far the most successful human pathogen and is in fact limited to that host, with no known external reservoirs. In contrast, the highly related veterinary pathogen B. bronchiseptica infects different mammalian hosts and, at least in laboratory experiments, can survive outside the host. Molecular cladistics indicate that B. pertussis evolved from a common ancestor with B. bronchiseptica that was more similar to the latter (44). This evolution, which also resulted in a more acute disease manifestation, was associated with significant genome reduction arising from deletions totaling over 1.2 Mbp as well as inactivation of multiple genes by mutation (27). The BvgAS regulon has been conserved in both pathogens and overlaps extensively. One marked difference is the expression of pertussis toxin in B. pertussis but not in B. bronchiseptica, even though the genes for both the synthesis and the secretion of this toxin are intact in typical B. bronchiseptica lineages. The evolutionary value of such a system of environmentally responsive global regulation of virulence genes is evident for a bacterial pathogen such as B. bronchiseptica with a theoretical external reservoir but is not obvious for an obligate human pathogen such as B. pertussis.

At least two different rationalizations for the presence of the BvgAS system in B. pertussis have been put forward. One holds that the BvgAS regulatory system is an evolutionary fossil, and although it may have played an important role in optimizing gene expression in two different environments in an ancestor, it now simply functions for virulence gene expression (28, 45). Consistent with this view is that B. pertussis does not require the Bvg(−) mode for full virulence in animal models of infection. However, mutants that express virulence genes constitutively can arise easily by single nucleotide changes, and yet this environmental responsiveness has been conserved in evolution. An alternative possibility is that, while the mode of passage between hosts may have become somewhat different and much more abbreviated for B. pertussis, the Bvg(−) or Bvg(i) mode still plays a role in maximizing aerosol transmission. Conceivably, this could occur by allowing bacteria to respond effectively to stresses encountered in the harsh environment of an aerosol created by coughing, by aiding in the preparation for expulsion from the host, or by aiding in the initial establishment of infection after deposition in a new host. While the role of this reciprocal arm of the BvgAS regulon has not yet been rigorously tested in transmission models of pertussis, initial experiments indicate that it contributes to survival of B. pertussis in static aerosols (T. Merkel, personal communication).

Thus, a full characterization of the BvgAS regulon in B. pertussis is needed to further our understanding of virulence gene expression as well as gene expression in the Bvg(−) mode. Here, we report a state-of-the-art analysis of this regulon using RNA-seq analyses coupled with extensive RT-qPCR, which provides a fuller picture of the Bvg(+) and Bvg(−) modes of B. pertussis strain Tohama I.

In the Bvg(+) mode, BvgAS controls dozens of transcriptional regulators.

Our analyses confirm and extend results from previous microarray analyses (9, 31, 32). As expected, we observed the activation of multiple well-known vags, such as genes encoding adhesins and toxins (Table S2), T3SS (Table S2), and transcriptional regulators (Table 1). However, this work doubles the number of transcriptional regulators that were observed to be differentially regulated between the Bvg(+) and Bvg(−) modes, further establishing the role of BvgA~P as a major regulator for networks controlling virulence and metabolism. Among these regulators is brpL, the gene for the ECF sigma factor BrpL, which is responsible for expression of the T3SS operon. Its expression is increased >100-fold by the presence of BvgA~P (Table 1) (32, 46). Our in vitro transcription analyses with B. pertussis RNAP and BvgA~P demonstrate that BvgA~P directly activates P. The locations of the binding sites for BvgA~P at this promoter (centered at −41.5 and −63.5) are similar to those seen at P (36).

The existence of a hierarchical regulatory cascade for a T3SS operon is not novel. Some pathogenic bacteria use such a complex interplay between transcriptional regulators or a regulatory network to control secretion. For instance, PhoP/PhoQ in Salmonella enterica communicates with another two-component system, pmrAB, and also directly activates the expression of ssrB, encoding the response regulator SsrB, in order to express the Salmonella T3SS and other effectors (47; reviewed in reference 48) while CpxRA in Shigella species regulates virF and virB, which activate the T3SS operon via a regulatory cascade (reviewed in reference 49). However, the importance of the T3SS in B. pertussis has not been clear. Although some previous work had failed to detect T3SS proteins in B. pertussis and concluded that the T3SS may be inactive (46), other work has challenged this idea by demonstrating that clinical isolates and laboratory-adapted B. pertussis strains under certain nutrition-limited cultures produce T3SS proteins (50). In addition, Ahuja et al. (32) recently reported that in both B. pertussis and B. bronchiseptica, BtrA (BP2233) functions as an anti-sigma factor for BrpL, the ECF sigma factor responsible for T3SS expression. This observation suggests that T3SS gene expression can be activated by the downregulation of BtrA. Additionally, the same group has demonstrated that the B. pertussis T3SS effector BteA (BP0500) is cytotoxic in cell culture, again implicating the T3SS in B. pertussis pathogenesis. Taken together, these observations reveal a complex B. pertussis regulatory network that is poised to activate T3SS-mediated toxicity under specific conditions.

B. pertussis is likely to undergo significant metabolic changes in the Bvg(−) mode.

The most significant finding from our work is the demonstration that genes for multiple catabolic and metabolic pathways are upregulated in the Bvg(−) mode in B. pertussis. As detailed below, these pathways are varied. However, they all carry the common theme of potentially aiding the bacterium as it undergoes environmental and nutritional stress.

First, multiple genes connected to the TCA cycle are regulated. These include genes encoding sugar and amino acid transporters and pyruvate dehydrogenase that bring different substrates into the TCA cycle, acyl-CoA/enoyl-CoA synthetases and dehydrogenases that participate in lipid and fatty acid metabolism/homeostasis, the Leu/Ile/Val transporter system and other branched amino acid transporter systems, and the phenylacetic acid degradation pathway. The last pathway is a central route by which the catabolisms of many aromatic compounds, such as phenylalanine and 2-phenylethylamine, converge and are directed to the TCA cycle (51; reviewed in reference 52). Interestingly, this pathway is required for full pathogenicity of Burkholderia cenocepacia in Caenorhabditis elegans (53).

The expression of genes in the glycolate/glyoxylate pathway in B. pertussis (glcDEF and glcB) also increases in the absence of BvgA~P. For many bacteria, this pathway provides energetic advantages when under stress. For example, in Escherichia coli the glyoxylate cycle allows cells to utilize simple carbon compounds as a carbon source when complex sources are not available and bypass two oxidative steps of the TCA cycle, preventing the loss of 2 carbon dioxide molecules (54). Interestingly, the glycolate bypass pathway has been shown to be important for the long-term persistence of Mycobacterium tuberculosis infection (55). Even in plants, glycolate metabolism allows seeds to use lipids as an energy source for growth by converting lipids to carbohydrates (56). Taken together, these findings are consistent with the idea that B. pertussis employs glycolate/glyoxylate to conserve energy when under modulation.

Protection from the environment during modulation is suggested by the upregulation of two cold shock genes, of cyanophycin synthase (cphA1), and of glmUS. GlmS converts fructose-6-phosphate to glucosamine-6-phosphate, while GlmU converts N-acetylglucosamine-1-phosphate to uracil-diphosphate-N-acetylglucosamine, a precursor in the syntheses of both peptidoglycan and LPS (57, 58). Thus, in the Bvg(−) mode, B. pertussis may use GlmS/GlmU to modify its cell wall and LPS in order to survive environmental changes. Cyanophycin, a poly(amino acid) generated by CphA using arginine and lysine, is a nitrogen-rich reservoir that functions as storage for nitrogen, carbon, and energy through its deposition into granules in the cytoplasm (reviewed in reference 59). However, poly(amino acids) can also contribute to the formation of capsules that protect bacteria from dehydration or host immune attack. Increased expression of cold shock genes could protect the bacteria from lower temperature outside the host.

Taken together, our results suggest that many genes induced in the Bvg(−) mode could aid survival under Bvg(−) conditions. Although B. pertussis spends almost all of its lifetime within the human respiratory tract, it does need to be transmitted by aerosol droplets from one infected individual to another. At that time, the bacterium will experience a drastic change in its environment, including a significant temperature drop. It is important to note that phosphorylation is rapidly lost, at least in culture, when cells are shifted from nonmodulating to modulating conditions (11), consistent with the idea that the Bvg(−) mode could arise quickly in expelled droplets. Our finding that genes which could aid in this transition are upregulated in the Bvg(−) mode is consistent with the idea that Bvg(−) regulation may be important for transmission. Furthermore, even during infection, Bvg(−) mutants accumulate among the B. pertussis population in the monkey nasopharynx, leading to the speculation that the Bvg(−) mode facilitates bacterial persistence (30). Our identification of new metabolic genes that are induced in the Bvg(−) mode provides a starting point for investigating how these metabolic changes help the pathogen as it proceeds through its complete life cycle.

MATERIALS AND METHODS

Bacterial strains and cell culture.

B. pertussis Tohama I BP536 (60) and its ΔbvgAS derivative (QC3691, this study) were grown on BG agar (42) containing streptomycin (50 μg/ml) at 37°C in the absence (nonmodulating) or presence (modulating) of 50 mM MgSO4 (42). After 3 days, cells were scraped from plates and resuspended in PLB medium (42) with or without MgSO4 to obtain an initial optical density at 600 nm (OD600) of between 0.1 and 0.2. Cells were incubated at 37°C with shaking at 250 rpm until the OD600 reached a value between 0.5 and 0.6. Culture aliquots of 1 ml were then collected by centrifugation for RNA isolation for both RNA-seq and RT-qPCR analyses.

DNAs.

For in vitro transcription templates, PCR fragments harboring promoters of either brpL (−235 to +107) or BP1005 (−182 to +165) were synthesized by PCR using Q5 polymerase (New England Biolabs) and then cloned into the BamHI and SalI sites of the plasmid pTE103 (61) to generate either pGFK3022 or pGFK3023, respectively. Plasmids containing possible promoters for BP2142 (−243 to +57, pGFK3021), BP2226 (−261 to +83, pGFK3026), BP2230 (−262 to +98, pGFK3027), BP2233 (−277 to +63, pGFK3028), and BP3441 (−263 to +160, pGFK3024) were synthesized by GenScript, again using the BamHI and SalI sites of pTE103 (positions are relative to the start of the transcription detected by RNA-seq, except for brpL, whose +1 start site was identified by primer extension). pfha, which contains P cloned into pTE103, has been described previously (62). 5′-32P-labeled DNA fragments harboring either the brpL promoter (−163 to +91) or the fhaB promoter (−138 to +111) were generated by PCR using Q5 polymerase. For production of radiolabeled DNA, the template or nontemplate primer was treated with T4 polynucleotide kinase (Optikinase; Affymetrix) in the presence of [γ-32P]ATP. The PCR product was purified on polyacrylamide gels and isolated using the EluteTrap system (Schleicher and Schuell) as described previously (63).

Proteins.

E. coli RNAP core was purchased from New England Biolabs. B. pertussis RNAP core (20), E. coli σ70 (11), His6-tagged B. pertussis primary σ (42), BvgA (42), BvgAC93A-C103A-V148C (11), and BvgAC93A-C103A-T194C (11) were purified as described previously. (The latter proteins are described in the text simply as BvgAV148C or BvgAT194C.) When indicated, BvgA (or mutants) was phosphorylated by incubation in the presence of 20 mM acetyl phosphate (lithium potassium acetyl phosphate; Sigma-Aldrich) in 20 mM Tris-Cl (pH 8) for 30 min at room temperature. Nonphosphorylated BvgA was incubated for the same period of time in the presence of 20 mM Tris-Cl (pH 8) only. BvgAC93A-C103A-V148C and BvgAC93A-C103A-T194C were conjugated with FeBABE [Fe(III) (S)-1-(p-bromoacetamidobenzyl)EDTA; Dojindo Corp.] as described previously (20).

RNA-seq analyses.

Total RNA, isolated by a modified hot-phenol extraction method as described previously (64), was treated with DNase (Turbo DNA-free kit; Life Technologies, Inc.) for 1 h at 37°C, followed by phenol extraction and ethanol precipitation. After rRNA removal using the RiboZero kit (Illumina), strand-specific DNA libraries were prepared for Illumina sequencing using the ScriptSeq 2.0 kit (Illumina). Libraries were sequenced using a HiSeq 2000 sequencer (Illumina; University of Buffalo Next Generation Sequencing Core Facility). Before mapping, vrgX was annotated into the current B. pertussis genome annotation (NC_002929.2) using GenEdit (Digi Tech). Sequence reads were then mapped to the reference genome and normalized against total reads for further analysis. Differential expression analyses were performed using CLC Biomedical Genomic Workbench version 3.5.4 with default parameters. At least a 1.6-fold change in gene expression with an FDR-corrected P value of ≤0.05 in both set 1 and set 2 was considered significant. tRNAs and rRNAs were eliminated from our analysis. In addition, the multiple insertion sequence (IS) elements and repetitive DNA present in B. pertussis (27) were not analyzed.

RT-qPCR analyses.

Cell pellets were resuspended in 800 μl ice-cold phosphate-buffered saline (PBS), and total RNA was isolated by hot-phenol extraction, treated with DNase, and purified as described previously (65). cDNA was generated using 1 μg of RNA, with a 12 μM final concentration of random hexamers [d(N)6; New England Biolabs] as primers and Moloney murine leukemia virus reverse transcriptase (New England Biolabs) (65).

RT-qPCR analyses were performed using a CFX96 Touch real-time detection system (Bio-Rad). Expression of the primary sigma factor gene (rpoD; BP2184) was used as an internal standard, and SsoAdvanced SYBR green Supermix (Bio-Rad) was used as a signal reporter. Reactions were performed in a 96-well microtiter PCR plate using 1 μl of cDNA at 3.3 ng/μl and sense and antisense primers (0.5 μM each) for amplifying each targeted gene in 1× SsoAdvanced SYBR green Supermix (Bio-Rad) plus 3 μM MgCl2. Cycling conditions were as follows: denaturation (98°C for 30 s), amplification and quantification (95°C for 30 s, followed by 40 cycles of 95°C for 10 s, 53 to 64°C [appropriate annealing temperature] for 30 s, and 72°C for 30 s, with a single fluorescence measurement at both 60 to 53°C and 72°C for 30 s segments), a melting curve program (65 to 95°C with a heating rate of 0.5°C s−1 and continuous fluorescence measurement), and a cooling step to 65°C. Data were analyzed using the CFX manager software (Bio-Rad). The expression level of each sample was obtained by the standard curve for each gene and was normalized by the level of the internal control, primary sigma, rpoD (BP2184). Oligonucleotides used to amplify target genes and the internal control are available upon request. Relative expression values are shown as a fold difference using the threshold cycle (2−ΔΔ) method (65). RT-qPCR experiments were carried out with 3 biological replicates per strain, and each condition and all qPCRs were performed in duplicate. Student’s t test was used for statistical analysis. As with the RNA-seq data, at least a 1.6-fold change in gene expression with an FDR-corrected P value of ≤0.05 was considered significant.

Assignment of genes by PSI-BLAST.

PSI-BLAST (position-specific iterative Basic Local Alignment Search Tool) was used to assign possible functionalities of genes that previously had unknown functions. After a protein sequence was obtained from the DNA sequence using either ExPASY or Artemis, PSI-BLAST was performed at the NCBI Blastp website.

In vitro transcription and FeBABE cleavage reactions.

Transcription reactions were performed as described previously (20) using 0.05 pmol of supercoiled DNA; as indicated, 5 pmol BvgA or BvgA~P; and 0.75 pmol B. pertussis or E. coli RNAP (ratio of core to sigma factor of 1:2.5 to 3) in a total final volume of 10 μl. Reaction mixtures were incubated at 37°C for 15 min before the addition of heparin (final concentration of 200 μg/ml) and nucleoside triphosphates (NTPs) (final concentration of 250 μM [each] ATP, GTP, and CTP and 25 μM [α-32P]UTP at 4.4 × 103 disintegrations per minute [dpm]/pmol). Reactions were stopped, and products were separated on 4% (wt/vol) polyacrylamide-7 M urea denaturing gels as described previously (20). FeBABE reactions were performed as described previously (20).

Primer extension.

Primer extension analysis was performed using the Promega primer extension system as described previously (66), using P RNA (obtained as described above except that the final concentration of each NTP was 250 μM) and the oligonucleotide primer 5′ GATATGCGTTGTCCGGTTTTTC 3′, which is complementary to nucleotides +75 to +96 of the brpL leader region and had been end labeled with [γ-32P]ATP using Optikinase (Affymetrix). Primer extension products were electrophoresed on an 8% polyacrylamide-7 M urea denaturing gel. A DNA sequencing ladder, obtained using the SequiTherm Excel II DNA sequencing kit (Epicentre), was prepared using the in vitro transcription template DNA and the radiolabeled primer.

Accession number(s).

The transcriptomic data are available in the NCBI database (GEO number GSE98155) and in Table S1A.