Construction and validation of a first-generation Bordetella bronchiseptica long-oligonucleotide microarray by transcriptional profiling the Bvg regulon

Background Bordetella bronchiseptica is a bacterial respiratory pathogen that infects a broad range of mammals, causing chronic and often subclinical infections. Gene expression in Bordetella is regulated by a two-component sensory transduction system, BvgAS, which controls the expression of a spectrum of phenotypic phases transitioning between a virulent (Bvg+) phase and a non-virulent (Bvg-) phase. Results Based on the genomic sequence and using the freely available software ArrayOligoSelector, a long oligonucleotide B. bronchiseptica microarray was designed and assembled. This long-oligonucleotide microarray was subsequently tested and validated by comparing changes in the global expression profiles between B. bronchiseptica RB50 and its Bvg- phase-locked derivative, RB54. Data from this microarray analysis revealed 1,668 Bvg-regulated genes, which greatly expands the BvgAS regulon defined in previous reports. For previously reported Bvg-regulated transcripts, the gene expression data presented here is congruent with prior findings. Additionally, quantitative real-time PCR data provided an independent verification of the microarray expression values. Conclusion The results presented here provide a comprehensive, genome-wide portrait of transcripts encompassing the BvgAS regulon, while also providing data validating the long-oligonucleotide microarray described here for studying gene expression in Bordetella bronchiseptica.


Background
Bordetellae are Gram negative bacterial respiratory pathogens. Bordetella pertussis and Bordetella parapertussis hu , the causative agents of whooping cough, are human-adapted variants of Bordetella bronchiseptica, which naturally infects a broad range of mammals causing chronic and often asymptomatic infections [1]. The majority of gene expression in Bordetella is regulated by a two-component sensory transduction system encoded by the bvg locus. The bvg locus comprises a histidine kinase sensor protein, BvgS, and a DNA-binding response-regulator protein, BvgA. In response to environmental cues, BvgAS controls expression of a spectrum of phenotypic phases transitioning between a virulent (Bvg + ) phase and a non-virulent (Bvg -) phase. During the virulent Bvg + phase, the BvgAS system is fully active and many of the known virulence factors are expressed, such as filamentous hemagglutinin (FHA), pertactin, fimbriae, adenylate cyclase-hemolysin toxin, and dermonecrotic toxin (DNT), as well as a type III secretion system (TTSS) in B. bronchiseptica [2]. Conversely, BvgAS is inactive during the Bvgphase, resulting in the maximal expression of motility loci, virulence-repressed genes (vrg genes), genes required for the production of urease, and in B. bronchiseptica RB50, a siderophore, alcaligin [3][4][5]. Previous studies involving phase-locked and ectopic expression mutants demonstrated that the Bvg + phase promotes respiratory tract colonization by B. pertussis and B. bronchiseptica [6][7][8][9], while the Bvgphase of B. bronchiseptica promotes survival under conditions of nutrient deprivation [6,10].
The signals that activate BvgAS in nature are unknown. However, in the laboratory, BvgAS is active when the bacteria are grown at 37°C and inactive when grown at temperatures below ~26°C or in medium containing MgSO 4 or nicotinic acid at concentrations in the millimolar range. Although originally identified as a positive regulator of virulence gene transcription [11], it is now known that BvgAS controls expression of over a hundred different genes whose products are either proven or predicted to participate in a wide variety of cellular activities including many basic physiological functions [12][13][14]. Additionally, it is now understood that rather than functioning like an ON/OFF switch, BvgAS functions more like a "rheostat" capable of controlling gene expression of many phenotypic phases in response to subtle differences in environmental conditions [10].
The advent of microarray technology has enabled scientists to investigate biological questions, such as those pertaining to bacterial pathogenesis and host-pathogen interactions, in a global fashion. The cDNA microarray represents a popular array type in which double-stranded PCR products are spotted onto glass slides. However, construction of such microarrays presents a number of challenges, largely related to costs associated with amplicon validation, tracking and maintenance. For example, the laborious and problematic tracking of PCR amplicons leads to an estimated 10-30% misidentification [15]. Another limitation of cDNA microarrays is their inability, due to cross-hybridization, to reliably discriminate expression patterns of homologous genes [16]. With oligonucleotide arrays, problems related to clone tracking, homologous gene discrimination, and failed PCR amplicons are avoided, thus making long-oligonucleotide microarrays a more cost-and management-efficient alternative to cDNA microarrays. Here we present the design and assembly of a long-oligonucleotide B. bronchiseptica gene-specific microarray using the currently available genomic sequence generated by the Sanger Institute [17] and the software package ArrayOligoSelector [18]. This long-oligonucleotide microarray was then tested and validated by evaluating changes in the global expression profiles between B. bronchiseptica strain RB50 and its Bvgphase-locked derivative, RB54.

Results and discussion
To construct a B. bronchiseptica-specific whole genome microarray, the freely available software program, Array-OligoSelector (AOS) [18], was used to generate 70-mer oligonucleotide probes for every ORF in the Bordetella bronchiseptica RB50 genome [17]. The rationale behind designing and utilizing oligonucleotide probes versus PCR amplicons as probes, and subsequently the 70-mer length of the oligonucleotide probes, was chosen for several reasons. Long oligonucleotides are a highly sensitive alternative to PCR products and provide a means to readily distinguish between genes with high degrees of sequence similarity, which is an issue for the B. bronchiseptica genome [17]. For example, except for the extreme 5' and 3' termini, the fhaB and fhaS genes are nearly identical [19]. Additionally, previous results involving an anchored set of oligonucleotides revealed a strong relationship between the oligonucleotide length and hybridization performance [18].
For each B. bronchiseptica ORF, the AOS program optimizes the oligonucleotide selection on the basis of uniqueness in the genome, sequence complexity, lack of self-binding and GC content. Candidate oligos closest to the 3' end of the gene are then chosen [18]. There are a number of missing array elements due to gene duplications, prophage duplications, and ORF assignments missing from the completed genome annotation [17]. These missing array elements and, in the case of gene duplications, their corresponding represented array elements are listed in Additional File 1-see Supplementary Table S1. A list of the final 4975 oligonucleotide array elements or probes representing the entire B. bronchiseptica genome is given in Additional File 1-see Supplementary Table S2.
To test the usefulness of the newly constructed long-oligonucleotide Bordetella bronchiseptica microarray, a direct comparison between the transcriptional profile of B. bronchiseptica RB50 and RB54, a B. bronchiseptica Bvgphase locked derivative of RB50, was performed. The rationale behind performing this comparison to validate and test the B. bronchiseptica long-oligonucleotide microarray is that this comparison will globally identify B. bronchiseptica genes regulated by BvgAS. Previous studies, including cDNA microarray studies, have identified 538 genes controlled by BvgAS [2,13,14], thus providing a substantial reportable reference set to validate gene expression data generated from the newly constructed B. bronchiseptica long-oligonucleotide microarray.
Utilizing the B. bronchiseptica long-oligonucleotide microarray, ratio data collected from microarray analysis involving the direct comparison between B. bronchiseptica RB50 and RB54 revealed 1,668 Bvg-regulated genes (identified using SAM with a false discovery rate of <0.1%). This is a substantial increase in the number of Bvg-regulated transcripts compared to the 538 Bvg-regulated genes recently reported by others using cDNA microarray analysis [14]. A complete list of the fold-change expression values from this comparison, along with dye swap experimental data, is provided in Additional File 1- Supplementary Table S3. One possible effect of using the B. bronchiseptica RB50 strain grown at 37°C to represent the Bvg + phase, as opposed to using a Bvg + phase-locked strain, is that some Bvg-activated genes could be missed in this analysis. Possible explanations for the large increase in Bvg-regulated transcripts detected in this study, compared to previous reports, may include (i) differences between strains used, for example wild-type versus phase-locked derivatives, (ii) differences between gene expression threshold cut-off values used in analysis, and (iii) for microarray studies, differences between array platforms, such as cDNA microarrays versus oligonucleotide microarrays. Gene expression profiles reported here of demonstrated Bvgregulated transcripts are consistent with previous results ( Figure 1).
Historically, numerous Bvg-activated genes have been described, such as TTSS genes, fimbrial genes, filamentous hemagglutinin, pertactin, adenylate cyclase-hemolysin toxin, and dermonecrotic toxin [2]. Using the newly constructed oligonucleotide B. bronchiseptica microarray, all of these transcripts were found to be positively regulated by BvgAS ( Figure 1). Along with other virulence-related transcripts classically characterized as Bvg-activated, three putative adhesion genes fhaS (+10.37-fold), fhaL (+8.04fold), and BB0110 (+8.58-fold) were found to be highly Bvg-activated, as well as a putative novel toxin, BB3242 (+50.72-fold), containing an aerolysin and pertussis toxin domain (Table 1 and see Additional File 1-Supplementary Table S3). Seven ATP-binding cassette (ABC) transporters were found to be Bvg-activated, including BB1593 (+83-fold) ( Table 1 and see Additional File 1-Supplementary Table S3). In other bacteria, these transporters serve as important virulence factors based upon their role in nutrient uptake and in secretion of toxins and antimicrobial agents [20]. Five autotransporters were highly Bvg-activated, including two virulence-associated transcripts vag8 and sphB1 (Table 1 and see Additional File 1-Supplementary Table S3). The autotransporter sphB1 has been implicated in the maturation process of FHA [21]. The Bvgactivated expression profile of these transporters implicates their importance in pathogenesis. Five genes involved in protein folding and ushering were highly Bvgactivated. These include dsbG (+35.29-fold) and mucD Expression of demonstrated Bvg-regulated genes Figure 1 Expression of demonstrated Bvg-regulated genes. Hierarchical clustering of known B. bronchiseptica Bvg-regulated genes performed using MeV [32]. Expression profiles of genes are in rows. Data are mean centered for each array element and averaged from biological replicates. Red, indicates increased expression in RB50; green, decreased gene expression in RB50 (and increased expression in RB54); black, no significant change in gene expression. (+2.75-fold), a DegP family serine protease (Table 1 and  see Additional File 1-Supplementary Table 3). Homologues of both dsbG and mucD are known to serve roles in virulence in other organisms [22,23]. Twenty-seven transcripts identified as serving roles in global regulatory functions were found to be Bvg-activated, including the conventionally characterized bvgAS and the TTSS regulator brpL (Table 1 and see Additional File 1-Supplementary  Table S3). Other Bvg-activated transcripts in this group include the other TTSS regulators (BB1645, BB1646, and BB1642), three proposed two-component response regulatory proteins, and numerous proposed transcriptional regulators. Lastly, BB4228 (+15.70-fold), recently designated BopC and characterized as a novel type-III secretion effector [24], was found to be highly Bvg-activated.
The overall global transcriptional program observed during the Bvgphase supports previous data implicating a crucial role for the Bvgphase in survival during environmental stress conditions [2,25]. Genes annotated to serve global regulatory functions were one of the functional categories with the highest number of Bvg-repressed transcripts. Transcripts in this category include forty-nine transcription factors and/or DNA-binding proteins, and ten two-component systems ( Table 2 and see Additional  File 1-Supplementary Table 3). Numerous genes involved in transport were identified as Bvg-repressed, including twenty-seven ABC transporters and TonB-dependent receptors, such as bfrG BB1294 (-2.86-fold) and bfeA BB1942 (-2.80-fold) ( Table 2 and see Additional File 1-Supplementary Table S3). As mentioned earlier, these transporters serve important biological roles in other organisms, including participation in host-pathogen interactions as known virulence factors [20]. The Bvgrepressed expression profile of these genes highlights nutrient scavenging as a critical ability during Bvggrowth and supports prior data demonstrating that the Bvgphase is optimized for growth under nutrient limiting conditions [2,25].
Consistent with previous results, genes required for the expression of the siderophore alcaligin [26] and urease [5] were preferentially activated in the Bvgphase, along with a putative hemolysin BB1186 (-2.83-fold) ( Table 2 and  see Additional File 1-Supplementary Table S3). Also consistent with previous reports, all of the genes known to be involved in chemotaxis, such as cheD (-10.59-fold) and cheW (-15.21-fold), and motility, such as flgB (-44.49fold) and flgC (-32.10-fold), were Bvg-repressed [3] ( Table  2 and see Additional File 1-Supplementary Table S3). Numerous genes involved in electron transport were Bvgrepressed along with numerous genes involved in protein folding and ushering, including a high number of universal stress proteins, congruent with the Bvgphase role in stress survival [6,25] (Table 2 and see Additional File S1-Supplementary Table S3).
Bordetella pertussis only infects humans and, as mentioned above, is responsible for causing an acute upper-respiratory disease known as whooping cough or pertussis. Due to its high degree of similarity to B. pertussis and its broad host range, including animals conveniently used in laboratory studies, B. bronchiseptica is widely used as a model for Bordetella pathogenesis research. Data arising from the recent completion of the comparative sequencing of three different Bordetella strains has revealed a loss or inactivation of a large number of genes in B. pertussis [17]. Thus an intriguing question is how many of the B. bronchiseptica Bvg-regulated genes are intact in B. pertussis. To determine this, the ORF sequence and oligonucleotide probe sequence for every Bvg-regulated gene identified in this study was used to blast the B. pertussis genome sequence. This analysis resulted in identifying 1172 shared genes that are Bvg-regulated in B. bronchiseptica (see Additional File S1-Supplementary Table S4).
Real-time quantitative PCR was performed to provide an independent assessment of microarray expression measurements for selected genes. Genes were chosen to reflect the full spectrum of fold-changes identified by microarray analysis. Specifically, transcripts identified as having little to no change in gene expression by microarray analysis and selected for qRT-PCR include BB0057 rpoA, BB4989 dnaA, BB1037, and BB3703 eno (Table 3). Bvg-activated genes identified by microarray analysis and selected for qRT-PCR include the TTSS regulatory genes, BB1619 bcrH1 and BB1622 bcrH2 (Table 3). On the polar end of the Bvg spectrum, genes identified as Bvg-repressed and selected for qRT-PCR include BB2522 and BB1315, a putative universal stress protein. Additionally, BB4835 rpoH, sigma-32, identified as Bvg-repressed (-3.09-fold) was also selected for qRT-PCR analysis since it had not been identified as Bvg-regulated by previous cDNA microarray studies [13,14]. Real-time quantitative PCR analysis of this gene set provided data consistent with the quantitative measures by microarray analysis, using the newly constructed B. bronchiseptica oligonucleotide microarray (r 2 = 0.94) ( Figure 2). Therefore, this real-time quantitative PCR data provides independent verification of the microarray results.
Given that BvgAS is the major virulence regulator in Bordetella, the information presented in this report should allow researchers to design future experiments targeting some newly identified Bvg-regulated genes and study their role in the pathogenesis of this organism. A recent study demonstrated that Bordetella bronchiseptica, Bordetella parapertussis and Bordetella pertussis all contain a higher number of transcription factors or regulatory elements  DNA microarray analysis was used to measure mRNA levels present in B. bronchiseptica RB50 compared to mRNA levels present in B. bronchiseptica RB54. Differences in mRNA levels are listed as mean fold-changes ± standard error. Fold-changes were calculated by averaging the data from three biological sample sets.
than the genome size of each strain would predict [27]. The microarray analysis presented here identifies, for the first time, an additional 138 predicted global regulator elements as being a part of the BvgAS regulon, which could be described as a modulon. Moreover, the data indicates that a high number of B. bronchiseptica regulators are under control, either directly or indirectly, of BvgAS, thus highlighting an organized hierarchical network, with multiple layers of control, governing transcriptional regulation in B. bronchiseptica. Distinct phenotypic changes occurring at each end of the Bvg regulatory continuum have been described since 1960 [28]. Since then, the molecular basis underlying these morphological phases has been directly linked to the coordinate expression of the BvgAS regulon. The results presented here provide a comprehensive, genome-wide portrait of transcripts encompassing the BvgAS regulon, while also providing data validating the usefulness of utilizing the long-oligonucleotide microarray described here for studying gene expression in B. bronchiseptica.

Conclusion
Long-oligonucleotide microarrays represent a highly sensitive, cost and management efficient, alternative to cDNA microarrays. Using the complete Bordetella bronchiseptica genome sequence and ArrayOligoSelector software, a 70mer oligonucleotide Bordetella bronchiseptica microarray was designed and assembled. This long-oligonucleotide microarray was then tested and validated by comparing changes in the global expression profiles between Bordetella bronchiseptica RB50 and its Bvgphase-locked derivative, RB54. Data from this microarray analysis revealed 1,668 Bvg-regulated genes, which dramatically increases Comparison of gene expression measurements by micro-array hybridization and quantitative real-time PCR the number of Bvg-regulated transcripts identified to date. Additionally, gene expression profiles of previously demonstrated Bvg-regulated transcripts are consistent with previous results and quantitative real-time PCR data provided an independent verification of the microarray expression values. The results presented here provide a comprehensive, genome-wide portrait of transcripts encompassing the BvgAS regulon, and provide data validating the use of the long-oligonucleotide microarray described here for studying gene expression in Bordetella bronchiseptica.
Bordetella Culture and RNA Isolation B. bronchiseptica strain RB50 was isolated from a naturally infected New Zealand White rabbit [6] and RB54 is a Bvgphase-locked derivative of RB50 [10]. B. bronchiseptica strains were maintained on Bordet-Gengou (BG) agar (Difco) containing 10% defibrinated sheep blood for determination of colony morphology and hemolytic activity. Three independent biological replicates of B. bronchiseptica RB50 and RB54, picked from well-isolated colonies on BG plates, were initially grown in Stainer-Scholte (SS) broth [30] supplemented with 40 μg/mL streptomycin. To ensure similar inocula and growth phase among all biological replicates, bacteria were then subcultured at a starting optical density at 600 nm (OD 600 ) of 0.02 into a 250-mL flask containing 50 mL SS broth and grown at 37°C with shaking at 275 rpm for 24 hours. At this time the OD 600 for the RB50 cultures were approximately 1 and the OD 600 for the RB54 cultures were approximately 3.5. Total cellular RNA was extracted using Trizol reagent (Invitrogen, Carlsbad, CA), treated with RNase-free DNase I (Invitrogen, Carlsbad, CA), and purified using RNeasy columns (Qiagen, Valencia, CA) according to manufacturers' recommended protocols.

Preparation of labeled cDNA
A 2-color hybridization format was used for the microarray analysis. For each biological replicate, RNA extracted from B. bronchiseptica RB50 was used to create Cy5-labeled cDNA and RNA extracted from B. bronchiseptica RB54 was used to create Cy3-labeled cDNA. Conversely, dye-swap experiments were performed analogously, in which RNA extracted from B. bronchiseptica RB50 was used to create Cy3-labeled cDNA and B. bronchiseptica RB54 was used to create Cy5-labeled cDNA. Fluorescently-labeled cDNA copies of the total RNA pool were prepared by direct incorporation of fluorescent nucleotide analogs during a first-strand reverse transcription (RT) reaction as follows: 5 μg total RNA and 4.4 μg of random oligonucleotide hexamers were incubated 2 minutes at 98°C, cooled on ice, combined with SuperScript III RTase buffer, 0.5 mM dATP, dGTP, dCTP, 0.2 mM dTTP, 1.5 nmol Cy3-or Cy5-dUTP (Amersham), and 2 μL SuperScript III RTase (in a volume of 26 μL). This mixture was incubated 10 minutes at 25°C and 120 minutes at 50°C. The two differentially labeled reactions to be compared were mixed and buffer exchange, purification, and concentration was accomplished by microcon-10 (Amicon) filtration.

Microarray hybridization and data analysis
Oligonucleotide microarrays were first prehybridized for 1 hour in 5×SSC, 1% BSA and 0.1% SDS at 42°C, followed by washing with H 2 O and then with isopropanol and dried by centrifugation for 5 minutes at 50×g. Following prehybridization, 45 μL hybridization solution (labeled cDNA, 5 μg tRNA, 2×SSC, 25% formamide, 0.1% SDS) was applied to oligonucleotide microarrays and incubated in a humidified chamber overnight at 50°C. Arrays were subsequently removed from humidified chambers and quickly submerged and washed in 1×SSC and 0.05% SDS for approximately 2 minutes, followed by two additional washes for 2 minutes each in fresh 0.06×SSC. Slides were then dried by centrifugation for 5 minutes at 50×g, scanned using a GenePix 4000B microarray scanner, and analyzed with GenePix Pro software (Axon Instruments, Union City, CA). Spots were assessed visually to identify those of low quality and arrays were normalized so that the median of ratio across each array was equal to 1.0. Automatically and manually flagged spots, spots with the sum of medians (635/532) signal intensity less than or equal to 100, and spots with signal intensity below threshold (sum of median intensities plus one standard deviation above the mean background) were filtered out prior to analysis. Ratio data from the three biological replicates were compiled and normalized based on the total Cy3% intensity and Cy5% intensity to eliminate slide to slide variation. Gene expression data were then normalized to 16S rRNA. The statistical significance of the gene expression changes observed was assessed by using the significant analysis of microarrays (SAM) program [31]. A one-class unpaired SAM analysis using a false discovery rate of 0.063% (<0.1%) was preformed. Hierarchical clustering of microarray data using Euclidean Distance metrics and Average Linkage clustering was performed using MeV software from TIGR [32].