The transcriptional response of Pasteurella multocida to three classes of antibiotics

Background Pasteurella multocida is a gram-negative bacterial pathogen that has a broad host range. One of the diseases it causes is fowl cholera in poultry. The availability of the genome sequence of avian P. multocida isolate Pm70 enables the application of functional genomics for observing global gene expression in response to a given stimulus. We studied the effects of three classes of antibiotics on the P. multocida transcriptome using custom oligonucleotide microarrays from NimbleGen Systems. Hybridizations were conducted with RNA isolated from three independent cultures of Pm70 grown in the presence or absence of sub-minimum inhibitory concentration (sub-MIC) of antibiotics. Differentially expressed (DE) genes were identified by ANOVA and Dunnett's test. Biological modeling of the differentially expressed genes (DE) was conducted based on Clusters of Orthologous (COG) groups and network analysis in Pathway Studio. Results The three antibiotics used in this study, amoxicillin, chlortetracycline, and enrofloxacin, collectively influenced transcription of 25% of the P. multocida Pm70 genome. Some DE genes identified were common to more than one antibiotic. The overall transcription signatures of the three antibiotics differed at the COG level of the analysis. Network analysis identified differences in the SOS response of P. multocida in response to the antibiotics. Conclusion This is the first report of the transcriptional response of an avian strain of P. multocida to sub-lethal concentrations of three different classes of antibiotics. We identified common adaptive responses of P. multocida to antibiotic stress. The observed changes in gene expression of known and putative P. multocida virulence factors establish the molecular basis for the therapeutic efficacy of sub-MICs. Our network analysis demonstrates the feasibility and limitations of applying systems modeling to high throughput datasets in 'non-model' bacteria.


Background
Pasteurella multocida is a gram-negative bacterial pathogen that has a unique history of serving as a model species for new discoveries. In his seminal experiments in the late eighteenth century, Louis Pasteur demonstrated both attenuation as well as protective immune response utilizing P. multocida in birds [1,2]. P. multocida is also an excellent model species for studying the effects of antibiotics on gram-negative bacteria because its short lipopolysaccharide O side chains make it more permeable, allowing investigations on the effects of antibiotics on bacterial metabolism [3].
Of the different P. multocida subspecies, multocida has a broad host range and causes diseases in poultry, cattle, pigs, and rabbits [4]. Zoonosis with P. multocida is caused by cat and dog bites and scratches and respiratory infection [5]. Fowl cholera (avian pasteurellosis) caused by P. multocida can be either chronic or acute and is a septicemic disease that causes significant economic loss to the poultry industry. Due to the importance of this pathogen to the avian community, the genome of an isolate recovered from chicken fowl cholera (strain Pm70) was sequenced in 2001 [6]. The genome sequence has enabled subsequent functional genomics research with this species [7][8][9]. It has also enabled investigations on the proteomic and transcriptomic response of P. multocida Pm70 to sub-MICs of antibiotics [10][11][12].
Much attention has been focused on characterization of specific antibiotic resistance mechanisms; however, recent studies are shedding light on non-target effects of antibiotics. In the current study, we used P. multocida as a model species to investigate the effects of three antibiotics with disparate mechanisms of action (amoxicillin, chlortetracycline, and enrofloxacin) on bacterial metabolism. Because biological systems utilize highly complex, interrelated metabolic and signaling networks, a thorough understanding of bacterial physiology requires studying the global response to a given stimulus in the context of interacting gene networks. Computational systems biology facilitates this level of analysis. Systems level analysis can identify key regulatory elements of a molecular interaction network; dynamic changes in these elements govern a response [13].
Application of systems approaches to P. multocida is still in the initial stages, but we have started using these methods to characterize the P. multocida response to antibiotics at the proteome level [11]. Here we describe the use of systems analysis at the transcriptome level to investigate the P. multocida response to sub-MICs of antibiotics.

Results and discussion
Differentially expressed genes Statistical analysis of microarray data revealed that 1/4 MIC of AMX, CTC, and ENR resulted in significant changes in expression of approximately 25% of the genome (525 genes, Additional file 1, supplemental table 1). The differences in gene expression that were determined to be statistically significant ranged from as small as a 5% decrease (relative to control) to as high as an 11.5fold increase in expression (recN, Table 1). Earlier microarray studies with P. multocida arbitrarily chose a 1.5-fold significant change in gene expression as differentially regulated genes. Instead of choosing a pre-determined threshold for determining the biological relevance of changes in gene expression, we considered every significant change in expression to be a valid change based on our rigorous statistical testing. We validated changes as small as a 15% decrease or increase (for example, expression of murA in response to AMX and ENR, respectively) (Additional file 1, supplemental table 1) by qPCR. Figure  1 shows the correlation between the trends observed in microarray and qPCR for the subset of genes with qPCR validations.
Compared to the untreated control, expression of 413, 392, and 473 genes had significantly altered expression in response to AMX, CTC, and ENR, respectively, and the overlap between treatments is shown in the Venn diagram ( Figure 2). The extensive overlap of regulated genes with the three different antibiotics is consistent with the recent transcriptome analysis of ampicillin and ofloxacin effects in E. coli, which showed that these two unrelated antibiotics had a significant overlap of regulated genes [14].

Functional analysis of DE genes: COGs
The three antibiotics used in this study have distinct modes of action; AMX is a cell wall biosynthesis inhibitor, CTC is a protein synthesis inhibitor, and ENR inhibits DNA gyrase and DNA topoisomerase IV activities [15]. At the COG category level, the three antibiotics had different effects on gene expression ( Figure 3). The AMX and CTC treatments resulted in an overall decrease in gene expression in all COG categories. The overall suppression of gene expression could indicate that the antibiotics have a marked detrimental effect on the fitness of P. multocida at doses below MIC, or it is also possible that the overall transcriptional shutdown is a compensatory response by slowing metabolism. ENR had varying effects on gene expression in different COG categories. In particular, ENR caused a pronounced increase in expression of genes in categories I (lipid metabolism), J (translation, ribosomal structure and biogenesis), and L (DNA replication, recombination and repair). All three antibiotics had reduced expression of genes in COGs V (defense mechanisms), G (carbohydrate transport and metabolism), P (inorganic ion transport and metabolism) and Q (secondary metabolites biosynthesis, transport and catabolism).

Effects of antibiotics on known target genes
We observed altered expression of some of the known targets (direct or inferred based on mechanism of action) with each antibiotic. CTC reduced the expression of genes encoding ribosomal proteins S4 and S7, which interact with tetracyclines [16]. One quarter-MIC of ENR resulted in a 1.5 fold increase in the expression of gyrB, which encodes a known ENR target, DNA gyrase B (Table 2). We had previously reported that the expression of RecA protein increased in response to sub-MIC of ENR [10], and in the current study we detected increased expression of recA.
With AMX, except for the expression of ponC, there were no detectable differences in the expression of genes encoding penicillin binding proteins (pbp). There was a 50% decrease in ponC (pbp) gene expression. However, this was a not an AMX specific effect, as all three antibiotics caused a decrease in ponC gene expression.
While regulation of expression of known targets is an expected response to antibiotics, it is not a consistent finding; expression of the target is usually not affected. For example, microarray analysis showed that enrofloxacin, trimethoprim, brodimoprim, and cefquinome had no effects on expression of their targets in P. multocida [12]. However, even though expression of the known antibiotic target is often not affected, bacteria typically have a "signature" gene expression response to specific antibiotics. Indeed, bacterial transcription profiles are now sometimes used to obtain initial indications of the mechanism of action for new compounds [17,18].

Common trends in response to AMX, CTC and ENR
Despite differences in their mechanism of action, the majority of genes with significantly altered expression (282) were common to all three antibiotics ( Figure 2). Beyond the effects on specific target gene expression, antibiotics are known to cause secondary effects on genes involved in general physiology as part of the adaptive response to antibiotic stress. A hallmark of recently described antibiotic mediated bacterial cell death common to bactericidal antibiotics in E. coli is the generation of hydroxyl radicals with subsequent induction of SOS response [19,20]. These effects were noted for MIC as well as sub-MIC concentrations of the antibiotics. There is increasing appreciation for the fact that fine tuning of these responses is unique to each organism and needs to be evaluated as such to facilitate identification of targets that potentiate the bactericidal activity of antibiotics. The P. multocida genome is known to contain genes necessary for antibiotic mediated SOS response [6] (for example, sulA, which is necessary for beta-lactam mediated SOS response).
Antibiotic stress in P. multocida resulted in increased ATP synthesis (Additional file 2, supplemental table 2). While this effect has been described for gyrase inhibitors [21], in P. multocida this response was common to all three antibiotics. In particular, all three antibiotics influenced gene expression involved in de novo nucleotide biosynthesis. The overall effects on expression of purM (increase), purK (decrease), and purF (decrease) could reduce the intracellular concentrations of IMP, while increased expression of pyrG and pyrF with ENR could result in increased concentrations of CTP and UDP.
Both AMX and CTC impaired the expression of hsp90, thus affecting proper protein folding under stress conditions in P. multocida. All three classes of antibiotics adversely affected putrescine (potE), ribose (rbscA and rbscC), and molybdate (modC) transport systems. Alterations to modC expression could have wide ranging effects on bacterial physiology due to effects on the synthesis of molybdoenzymes (Additional file 2, supplemental table 2). Reduced levels of polyamines like putrescine could have profound negative effects on protein translation and impair the ability to cope with oxidative stress [22].
Expression of 121 genes annotated as coding for hypothetical proteins were altered in response to the antibiotics (Additional file 1, supplemental table 1). Lack of functional information for these genes hampers our understanding of their role in adaptation of P. multocida to antibiotic stress and demonstrates the need for continuous functional annotation of genes beyond the initial annotation that is described with the genome sequence.
One of the demonstrated outcomes of antibiotic stress in E. coli is alterations in iron homeostasis; namely, expression of genes involved in siderophore mediated iron transport are affected (for example, fecC) [19]. In P. multo-

AMX ENR
Validation of microarray differential gene expression data by quantitative real time RT-PCR (qPCR) Figure 1 Validation of microarray differential gene expression data by quantitative real time RT-PCR (qPCR). Ratio of treated vs. control were calculated for microarray data (ordinate) and qPCR data (abscissa). A total of six genes that were up or down regulated with amoxicillin, chlortetracycline or enrofloxacin were compared. The correlation coefficient r 2 was 0.84. Changes in the ratio of NAD+/NADH are known to influence the levels of toxic superoxide formation, which ultimately leads to hydroxyl radical formation. NADH I levels are directly linked to the activities of TCA cycle enzymes that generate NADH from NAD+ [19]. ENR treatment caused decreased expression of sucB, and all three antibiotics caused decreased lpdA expression, which could reduce hydroxyl radical formation. Expression of the gene encoding anaerobic respiration electron transfer enzyme molybdoenzyme dimethylsulfoxide reductase was reduced in response to AMX and ENR. AMX treatment reduced the expression of acetyl-CoA carboxylase enzyme (accB, accD), which catalyzes the first step of type II fatty acid biosynthesis as well as the step catalyzed by acyl carrier protein S-malonyltransferase (fabD); by contrast, ENR increased the expression of these genes, which could result in increased fatty acid biosynthesis.
DE genes: network analysis COG analysis of DE genes was useful for identifying overall trends at a broad functional category level. Gene ontology (GO) analysis of P. multocida, which is available through Uniprot, was useful to some extent in adding specific functions to the genes within each of the broad functional categories represented by COGs. However, neither of these two approaches could actually show exactly how the genes within a functional category were interacting with each other to bring about a specific function. To enable network analysis, we used Pathway Studio to build interaction networks for all DE P. multocida genes with AMX, CTC and ENR. Careful analysis of these networks revealed expression cascades for some of the DE genes common to CTC and ENR.
Expression of fnr increased in response to sub-MIC of CTC and ENR (Figure 4). Fnr regulates the transcription of many genes that are involved in either aerobic or anaerobic respiration, including Nap nitrate reductase, which facilitates nitrate reduction in the periplasmic space. Expression of napD and napF increased in response to ENR and decreased with CTC (Additional file 3, supplemental figure 1). NapD is a chaperone protein that is specific for NapA and is required for the stability of NapA protein. NapF stimulates the Nap reductase. NapF promoters regulate the transcription of cmABCDEFGH operon in E. coli, which constitutes the cytochrome biogenesis system. Expression of ccmH increased with CTC and decreased with ENR ( Table 3). The identification of this cascade of gene regulation was possible by combining the network analysis with the functional information available at EcoCyc. The implications of up regulation of enzymes that facilitate anaerobic respiration in response to antibiotics remain unclear at present.
Network analysis also facilitated the identification of differences in the SOS response with different antibiotics in P. multocida ( Figure 5 (ENR) and Additional file 3, supplemental figures 2 (AMX) and 3 (CTC)). In particular, the response to ENR showed the signature gene expression pattern typical of SOS response ( Figure 5). Increased  (Table 1). Expression of recN was increased, which also relieves LexA repression [23][24][25]. Genes involved in fine-tuning the actions of RecA, namely recX and uvrD, also had increased expression. Expression of primosome genes (dnaJ and dnaK) that encode enzymes required for re-starting the stalled replication fork and other replicative enzymes like DNA polymerase III increased with ENR.

Summary of significant changes in P. multocida gene expres-sion grouped by COGs
It is reported that beta lactam mediated cell killing involves SOS response through SulA and DpiBA two component system [26]. Based on the genome annotation of Pm70 [6], all the components of this regulon are not present. There was no induction of the expression of recA or lexA with AMX (Additional file 3, supplemental figure  2). The expression profile of SOS response genes remained either unchanged or showed decreased expression, which indicates that sub-MIC AMX does not activate SOS response in P. multocida. Excluding the increased expression of lexA and gyrB, CTC (Additional file 3, supplemental figure 3) also did not activate SOS response.

DE virulence genes
Productive infection by bacterial pathogens relies on the expression of virulence factors that have wide ranging functions like competence, adherence, capsule synthesis and export, evading host immune responses etc. Transcription profiling of the response of a bovine P. multocida isolate (L386) to MIC of eight different antibiotics identified mostly reduced virulence gene expression [12]. In the current study, transcription profiling of the P. multocida response to AMX, CTC and ENR identified significant changes in the expression of known and putative virulence genes. AMX and CTC reduced the expression of ompA, a known virulence factor of P. multocida involved in binding to host cells [27] (Table 4). AMX and CTC decreased expression of the gene encoding detoxifying enzyme superoxide dismutase, which is one of the virulence genes identified in the P. multocida L386 study [12]. Additional virulence factors (based on L386 study) which had reduced expression in Pm70 are: clpB protease with AMX, capsule transport protein hexA with all three antibiotics, PM1714 regulator with AMX and CTC, capsule biosynthesis genes phyA and hyaE with AMX and CTC, phyB with CTC, and tight adherence gene tadB with AMX. Expression of dnaK, which is required for virulence in a number of bacterial pathogens [26], was decreased by CTC.
Fumarate/nitrate response regulator network in P. multocida in response to enrofloxacin visualized in Pathway Studio Figure 4 Fumarate/nitrate response regulator network in P. multocida in response to enrofloxacin visualized in Pathway Studio. P. multocida sub-MIC ENR response was marked by significant changes in fnr response regulated genes. Red nodes are genes with increased expression and green nodes are genes with decreased expression. The rest of the interacting nodes are shown in gray. These nodes had no significant changes in expression or were orthologs from additional gram-negative species in the molecular interaction database in Pathway Studio.

Conclusion
Our global transcriptome analysis of P. multocida response to three antibiotics with differing modes of action identified gene expression changes in a little over a quarter of the annotated open reading frames in the genome. The antibiotics had varying effects on various cellular and metabolic functions. Amoxicillin reduced the overall transcription rate as reflected by the majority of the identified significant changes in gene expression showing a downward trend in expression. An interesting aspect of the AMX response was the marked lack of expression of SOS response genes. Enrofloxacin, a DNA gyrase inhibitor, resulted in the significant overexpression of its target gene and induced a typical SOS response. Analysis of DE genes in the context of interacting protein networks facilitated the identification of coordinated regulation of gene expression across different COG functional categories. Sub-MIC of antibiotics influenced the expression of virulence factors in P. multocida. Our results form the framework for understanding the global effects of antibiotics on P. multocida, which will aid rational drug design for containing as well as treating infections caused by this pathogen across multiple species.

Methods
Bacterial culture and RNA isolation P. multocida Pm70 was cultivated in brain heart infusion broth (BHI) at 37°C with rotary aeration. Minimum inhibitory concentrations (MICs) of AMX, CTC and ENR for Pm70 are 0.5 g/ml, 4 g/ml, and 0.031 g/ml, respectively [10]. Growth kinetics of Pm70 in the presence of 1/ 4 MIC of the three antibiotics were previously described [10]. Stationary phase cultures of Pm70 were used to inoculate 50 ml BHI to an initial A 600 of 0.05; antibiotic treated cultures contained 1/4 MIC of AMX, CTC, or ENR, and control cultures were grown without antibiotics. For each treatment, cultures were grown in triplicate to midlog phase (A 600 of 0.8). RNAprotect (Qiagen) was added to culture samples, and bacteria were harvested by centrifugation (10,000 × g, 10 min, 4°C). Pellets were stored at -80°C.
Total RNA from three biological replicate cultures (for untreated control culture as well as antibiotic treated cultures) was isolated using Qiagen RNeasy kit using the manufacturer's protocol. The quality and concentration of RNA was determined by Agilent Bioanalyzer.  Raw intensity data (single channel) were log transformed, and normalization was done by quantile normalization [28] and robust multichip analysis (RMA) algorithm [29].

Microarrays and data analysis
The experimental design and all microarray data have been deposited in the NCBI Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/, accession number GSE12779). ANOVA followed by step-down Bonferroni correction for multiple testing was conducted to identify genes with significant changes in expression relative to untreated control (p  0.05) in any group. Significant changes between control vs specific antibiotic treatments were identified by performing Dunnett's test [30]. The differences in gene expression in response to antibiotics were calculated as the ratio of treatment vs control intensities.

Real time RT-PCR
The majority of significant changes in gene expression identified in this study were less than one-fold. Therefore, to validate changes in gene expression, we performed duplex real time quantitative RT-PCR (qPCR) of six genes: gyrb, impA, rpsg, hpkR, recQ, and murA using gapdh as the internal standard. The RNA template used for qPCR was the same RNA that was used for microarray hybridiza-tions. The probes and primers for duplex RT-PCR were designed using Beacon Designer software (Additional file 4, supplemental table 3). Reactions were performed using the Platinum ® SYBR ® Green One-Step qRT-PCR Kit (Invitrogen Corporation, Carlsbad, CA). Amplification and detection of specific products were done using the iCycler iQ Real-Time PCR Detection System (Bio-Rad Laboratories, Inc., Hercules, CA). Regression analysis of the C t values of the test RT-PCRs was used to generate standard curves. The mean C t value for the GAPDH mRNA-specific reactions was used to normalize the test values. Regression analysis of the gene expression trends determined by qPCR and microarray analysis was performed in Microsoft Excel.

Functional analysis of differentially expressed genes
Initial functional analysis of genes with significant changes in expression was done using clusters of orthologous groups (COG) classification [31]. For each COG category in the genome, the percent increase or decrease in expression of genes belonging to that category for each antibiotic treatment was calculated. The COG functional categories described for P. multocida genes [6] along with the description for proteins encoded by the genes at Uniprot [32] were used to identify overall themes in antibiotic response. Specific functions of P. multocida proteins were deduced from information available for E. coli orthologs at the EcoCyc database [33]. Network analysis of DE genes was done in Pathway Studio (Ariadne Genomics) as described earlier [11]. Briefly, we built interaction networks in Pathway Studio with proteins of interest including the upstream regulators and downstream targets. In the interaction networks, different colors were used for nodes to indicate significant increase (red), decrease (green), or no change (pink) in gene expression in response to sub-MIC of antibiotic. Entities in the interaction map that were not present in the Pm70 genome were shown in gray color.