The BvgAS Regulon of Bordetella pertussis

ABSTRACT 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.

ing 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 (virulenceactivated 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 MgSO 4 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 MgSO 4 ) conditions, in which phosphorylated BvgA is undetectable [Bvg(Ϫ) mode], or nonmodulating (no MgSO 4 ) 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 MgSO 4 ; set 2, WT versus ΔbvgAS mutant in the absence of MgSO 4 ; and set 3, ΔbvgAS mutant with or without MgSO 4 (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 MgSO 4 in the absence of bvgAS.
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 MgSO 4 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 MgSO 4 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.
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.
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)(34)(35)(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 MgSO 4 compared to WT with MgSO 4 . While expression of bipA demonstrated only a 2-FD in WT without MgSO 4 versus WT with MgSO 4 , 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 MgSO 4 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). Category 2 contains genes encoding transcriptional regulators, including the ECF sigma factor brpL (also called btrS), the N gene rpoN-1 (BP0133), and several twocomponent 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.
Finally, category 4 (Table S2 and Fig. S1C) includes genes for membrane-associated transporter proteins and hypothetical proteins with unknown functions. The BvgAdependent 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 fim3 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~Pdependent 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 (antisigma 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 brpL transcription was observed only in the presence of BvgA~P, demonstrating that this promoter, like P fhaB , 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. P brpL contains multiple BvgA~P dimer binding sites. To determine whether BvgA~P directly binds to sites within P brpL , we performed iron bromoacetamidobenzyl-EDTA (FeBABE) footprinting using BvgA T194C -FeBABE and BvgA V148C -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 brpL 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 fhaB , served as a positive control. As expected, with either RNAP, we observed the three binding regions for BvgA~P at P fhaB , centered at positions Ϫ44.5, Ϫ66.5, and Ϫ88.5. At P brpL , 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 BvgA T194C -FeBABE, and an additional farther upstream binding site centered at approximately Ϫ114 was also seen. These results indicate that P brpL has a similar BvgA~P binding motif as P fhaB , but the stability and/or oligomerization pattern of BvgA~P on the DNA appears to differ with the two polymerases.

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

Moon et al.
® 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   (67). For instance, the sequence of the perfect inverted repeat is TAGGAAATTTCCTA, whose score is given as "0." 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 brpL . 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 fhaB (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 MgSO 4 (42). After 3 days, cells were scraped from plates and resuspended in PLB medium (42) with or without MgSO 4 to obtain an initial optical density at 600 nm (OD 600 ) of between 0.1 and 0.2. Cells were incubated at 37°C with shaking at 250 rpm until the OD 600 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 fhaB cloned into pTE103, has been described previously (62). 5=-32 P-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 [␥-32 P]ATP. The PCR product was purified on polyacrylamide gels and isolated using the EluteTrap system (Schleicher and Schuell) as described previously (63).
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 MgCl 2 . 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 ϪΔΔCT ) 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 [␣-32 P]UTP at 4.4 ϫ 10 3 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 brpL 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 [␥-32 P]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.