Role of SPI-1 Secreted Effectors in Acute Bovine Response to Salmonella enterica Serovar Typhimurium: A Systems Biology Analysis Approach

Salmonella enterica Serovar Typhimurium (S. Typhimurium) causes enterocolitis with diarrhea and polymorphonuclear cell (PMN) influx into the intestinal mucosa in humans and calves. The Salmonella Type III Secretion System (T3SS) encoded at Pathogenicity Island I translocates Salmonella effector proteins SipA, SopA, SopB, SopD, and SopE2 into epithelial cells and is required for induction of diarrhea. These effector proteins act together to induce intestinal fluid secretion and transcription of C-X-C chemokines, recruiting PMNs to the infection site. While individual molecular interactions of the effectors with cultured host cells have been characterized, their combined role in intestinal fluid secretion and inflammation is less understood. We hypothesized that comparison of the bovine intestinal mucosal response to wild type Salmonella and a SipA, SopABDE2 effector mutant relative to uninfected bovine ileum would reveal heretofore unidentified diarrhea-associated host cellular pathways. To determine the coordinated effects of these virulence factors, a bovine ligated ileal loop model was used to measure responses to wild type S. Typhimurium (WT) and a ΔsipA, sopABDE2 mutant (MUT) across 12 hours of infection using a bovine microarray. Data were analyzed using standard microarray analysis and a dynamic Bayesian network modeling approach (DBN). Both analytical methods confirmed increased expression of immune response genes to Salmonella infection and novel gene expression. Gene expression changes mapped to 219 molecular interaction pathways and 1620 gene ontology groups. Bayesian network modeling identified effects of infection on several interrelated signaling pathways including MAPK, Phosphatidylinositol, mTOR, Calcium, Toll-like Receptor, CCR3, Wnt, TGF-β, and Regulation of Actin Cytoskeleton and Apoptosis that were used to model of host-pathogen interactions. Comparison of WT and MUT demonstrated significantly different patterns of host response at early time points of infection (15 minutes, 30 minutes and one hour) within phosphatidylinositol, CCR3, Wnt, and TGF-β signaling pathways and the regulation of actin cytoskeleton pathway.


Introduction
Despite extensive efforts to understand and prevent Salmonella infection, the number of people infected with Salmonellae each year has not changed appreciably in the United States over the last two decades [1]. Approximately 1.4 million people in the United States are infected each year with Salmonellae. This infection is particularly devastating to children and immunocompromised adults. In the United States, children under the age of 4 are disproportionately affected by salmonellosis (72 cases per 100,000 children) [2]. A recent study estimates that there are annually 93.8 million cases of Salmonella gastroenteritis worldwide, resulting in 155,000 deaths [3]. In parts of Africa where HIV infection is widespread, Non-typhoidal Salmonella (NTS) infection is recognized as a major health problem and is the leading cause of pediatric bacteremia [4,5,6]. In severe cases, particularly in immunocompromised individuals, Salmonellae penetrate the intestinal mucosa and enter the bloodstream leading to systemic infection. Understanding the acute phase of intestinal infection is important to designing intervention strategies for children and immunocompromised adults where CD4 + T cell help, which is required to prevent disseminated infection, may be diminished or absent.
Calves are ideal models for human Salmonella infection, because they develop diarrhea with clinical and pathologic features similar to those in people [7,8,9,10]. In both people and calves, salmonellosis is characterized by influx of neutrophils into the intestine and fluid secretion into the intestinal lumen, accompanied by severe gastroenteritis with abdominal cramping and diarrhea. During the first few hours of intestinal infection, polymorphonuclear cells (PMNs) leave the vasculature and migrate across the mucosal barrier, and ultimately arrive at the lumen of the intestine [9]. Increased vascular permeability is associated with PMN migration and movement of fluid from the circulation into the intestinal lumen characterizes the first 8 to 12 hours of infection. This movement of fluid may form the basis for the diarrhea exhibited during enteric salmonellosis [9]. Influx of PMNs and fluid secretion into the intestinal lumen are hallmarks of Salmonella-induced enteritis in both humans and cattle. These naturally occurring clinical and pathological similarities between human and bovine infection make the calf an ideal animal to model the clinical signs of salmonellae infection in humans [8,11,12].
Salmonella enterica Serovar Typhimurium (S. Typhimurium) is one of the top Salmonella serotypes that cause disease in humans and cattle. Invasion of the intestinal epithelium is critical for the establishment of infection and occurs within the first 15 minutes of bacterial contact [9,13,14]. S. Typhimurium invades the nonphagocytic enterocytes of the host intestinal epithelium through manipulation of the host actin cytoskeleton by Salmonella effector proteins secreted through the Type III secretion system (T3SS-1) encoded at Salmonella Pathogenicity Island I [13,15,16,17,18,19]. Salmonella effector proteins bind to the actin cytoskeleton and induce the formation of cell membrane ruffles that facilitate Salmonella internalization [20,21,22,23,24,25,26,27,28,29,30,31].
Beyond invasion, the T3SS-1 delivered effector proteins are required for PMN migration and fluid secretion in the intestine. Previous studies have shown that PMN migration and fluid secretion in bovine intestinal loop model are the result of the coordinated actions of five effector proteins secreted by T3SS-1: SipA, SopA, SopB, SopD, and SopE2 [12,32]. Recent studies demonstrate the individual functions of these effector proteins: SipA binds actin directly; SopA is a HECT E3 ubiquitin ligase; SopB is a multi-functional phosphatidyinositol phosphatase with roles in actin cytoskeleton signaling, apoptosis, and trafficking of the Salmonella containing vacuole (SCV); SopE2 activates actin cytoskeleton rearrangements through Rho family GTPases; and SopD participates in invasion and has a role in membrane fission and macropinosome formation [29,33,34,35,36,37,38,39,40]. Individual loss of the genes encoding these proteins only partially decreases fluid secretion in the bovine intestine relative to that induced by wild type S. Typhimurium. Loss of all five genes together (DsipA, sopABDE2) significantly diminishes fluid secretion, PMN influx, and lesions in the intestine [12]. Further, SipA, SopA, SopB, SopD, and SopE/E2 have combined roles in C-X-C chemokine expression [12,32]. Loss of all five genes decreases host expression of the genes encoding C-X-C chemokines, GROa, GROc, interleukin 8 (IL-8), and CXCL6 (GCP-2) [32]. The CXC chemokines are chemotactic for PMNs and induce their migration from the vasculature to sites of infection. Therefore, the increased expression of these genes may in part explain PMN migration to the intestine. It is well known that T3SS-1 secreted effectors SipA, SopA, SopB, SopD and SopE2 are required for Salmonella induced bovine and human diarrhea, however, the host cellular pathways and pathway mechanistic genes critical to fluid secretion, PMN transmigration, enteritis and diarrhea are less well defined. We hypothesized that comparison of the response of bovine intestinal mucosa to wild type Salmonella and a SipA, SopABDE2 effector mutant relative to uninfected bovine ileum would reveal heretofore unidentified diarrhea-associated host cellular pathways. In addition, we also hypothesized that this comparison would lead to the discovery of specific pathway mechanistic genes. To test our hypothesis, we challenged bovine ligated ileal loops containing Peyer's patch with three conditions; wild type Salmonella Typhimurium, an isogenic SipA, SopA, SopB, SopD and SopE2 deletion mutant, or media as a negative control. We measured the host mucosal gene expression profile sequentially for 12 hours post-inoculation. We then compared the data from the host response to the three conditions by conventional and dynamic Bayesian analytical methods, confirming known cellular pathways associated with salmonellosis and identifying additional cellular pathways. Known and newly discovered host cellular pathways were merged to establish functional models of the host:pathogen interface and to identify pathway mechanistic genes for future study. In this study, we report the comparison of host responses to wild type Salmonella as compared to the SipA, SopABDE2 mutant demonstrating significantly different patterns and perturbations of host cellular pathways at early time points of infection and disease expression relative to the inflammatory response and diarrhea.

Bacterial Strains and Culture
Derivatives of Salmonella enterica serotype Typhimurium strain IR715 and ZA21, were used in this study. IR715, a derivative of ATCC 14028, is a spontaneous nalidixic acid resistant derivative referred to hereafter as ''wild type'' (WT) and ZA21 is a previously described DsipA sopABDE2 mutant (MUT) [12]. Cultures were grown in Luria-Bertani (LB) broth supplemented with the following antibiotics as appropriate at the indicated concentrations: nalidixic acid (50 mg/l), kanamycin (100 mg/l), chloramphenicol (30 mg/l), and tetracycline (20 mg/l). To prepare the bacterial inoculum, IR715 and ZA21 were cultured in LB broth for 18 hours at 37uC with shaking at 150 rpm in a shaking incubator (Model C24, New Brunswick Scientific, Edison, NJ). The resultant cultures were diluted 1:100 in sterile LB broth and incubated as before for an additional 4 hours. Bacteria in the exponential phase of growth were harvested by centrifugation for 15 min at 1,5006g and re-suspended in fresh LB broth to obtain a final concentration of approximately 0.75610 9 colony-forming units/ml (CFU/ml). The bacterial concentration of the inoculum was confirmed by spreading serial 10-fold dilutions on agar plates containing the appropriate antibiotics.

Animals
Four male Holstein calves, 4-6 weeks of age and weighing 45-55 kg, were used in each experiment under an approved animal use protocol in accordance with animal use policy under the supervision of the Texas A & M University Institutional Animal Care and Use Committee (IACUC). All animal work was approved under Animal Use Protocols 2004-029 and 2004-037. Calves were obtained from a local, commercial source and were unrelated. The calves were fed antibiotic-free milk replacer twice daily and water ad libitum. To minimize the possibility of interference from other enteropathogens on the host gene expression profile, calves were tested twice for fecal excretion of Salmonella spp. and Eimeria spp. oocysts and only animals that tested negative for these bacteria and parasites were used. Fecal specimens and rectal swabs were collected from calves two weeks prior to and immediately before the experiment and cultured for the presence of Salmonella serotypes. Swabs were placed in tetrathionate broth (BBL, Sparks, MD) overnight and subsequently plated onto XLT-4 agar plates (BBL). Fecal flotation with a hypertonic saline solution was used for detection of Eimeria spp. oocysts. All four calves were clinically healthy before the experiment and were culture negative for fecal excretion of Salmonella. We elected to use four calves because previously published work in our lab indicated that four calves would be sufficient to detect significant differences between the treatments used in this study [32].

Bovine Ligated Ileal Loop Surgeries
The bovine ligated ileal loop surgical procedure has been described previously [41]. Briefly, the calves were fasted for 12 hours prior to the surgery. Anesthesia was induced with propofol (Abbot Laboratories, Chicago, IL) followed by placement of an endotracheal tube and maintenance with isoflurane (Isoflo, Abbot Laboratories) for the duration of the experiment. A laparotomy was performed, the distal jejunum and ileum containing Peyer's patches were exposed, and 21 loops approximately 6 cm in length were ligated with umbilical tape, leaving 1 cm loops between them. Overnight bacterial cultures were grown as described above. The ileal loops were inoculated by intralumenal injection of 3 ml of a suspension of either S. Typhimurium strain IR715 (WT) or DsipA sopABDE2 mutant (MUT) in LB broth containing approximately 0.75610 9 CFU/ ml. Loops injected with sterile LB broth (3 ml) served as uninfected, negative controls (UI). Loops were alternately inoculated with LB broth, WT, or MUT starting at the loop closest to the ileocecal junction and were subsequently excised in the same order as inoculated. Three loops were removed at each time point, one loop for each condition. Each loop was inoculated with only one bacterial strain or the negative control. The loops were replaced into the abdominal cavity until collection by excision. At each of seven time points (15 min, 30 min, 1, 2, 4, 8, and 12 hours), three loops (one LB inoculated control loop, one wild type inoculated loop, and one mutant inoculated loop) were excised. From each of these loops, samples were collected for histopathology, bacteriology, electron microscopy, frozen sections, and RNA extraction.

Morphologic Analysis
For morphologic analysis by light microscopy, full cross-sections of each loop (including the Peyer's patch) were fixed in formalin, processed according to the standard procedures for paraffin embedding, sectioned at 5 mm thickness, stained with hematoxylin and eosin, and examined by light microscopy. All tissue sections collected contained Peyer's patch, which is continuous for approximately 1.2 meter proximal to the ileocecal junction in the calf. Previous work has demonstrated that ligated ileal loops constructed in the manner described have virtually identical morphology and cell composition [9].
Quantification of Tissue-associated S. Typhimurium Two-6 mm biopsy punches (0.1 g) from the anti-mesenteric side of the intestinal mucosa to maximize Peyer's patch collection were excised from every loop, washed three times in PBS to remove extracellular bacteria, homogenized, and diluted in 1 ml of PBS. Similar procedures were followed for the control loops to identify the presence of viable Salmonella as a consequence of leakage of the loops, cross contamination during material processing, or Salmonella migration via lymphatic or blood vessels. To determine the number of viable CFU of S. Typhimurium in tissues, lysates were serially diluted and cultured on LB agar with nalidixic acid (50 mg/l). Bacterial counts represent tissue-associated bacteria and do not distinguish between intracellular and extracellular bacteria. For RNA isolation, only the intestinal mucosa and submucosa were collected, but for bacteriology, the tissues were full thickness intestinal tissue containing the mucosa, submucosa, and smooth muscle layers.

Statistical Analysis of Bacteriology
For analysis of invasion percentages, the data were transformed logarithmically. The geometric means were determined, and the statistical significance of differences was calculated using two-tailed Student's t test.

RNA Extraction
Total RNA was extracted after dissection of the ileal loops obtained at 15 min and 30 min and at 1, 2, 4, 8, and 12 hours post-infection. A 6 mm biopsy punch was used to collect 8 tissue pieces of Peyer's patch from the excised loop. The mucosa and submucosa of the samples were immediately dissected, the tissue pieces finely minced with a sterile scalpel blade, and two biopsy punches (0.1 mg of tissue) transferred to 1 ml of Trizol reagent (Molecular Research Center, Cincinnati, OH). Tissues were immediately further disrupted with an RNase, DNase free plastic disposable pestle. The RNA extraction was performed using the recommended protocol from the manufacturer. The resultant RNA pellet was re-suspended in DEPC-treated water (Ambion, Austin, TX). Genomic DNA was removed by RNase-free DNase I treatment (DNA-free, Ambion, Austin, TX) according to the manufacturer's instructions, and samples were stored at 280uC until used.

RNA Sample Analysis
The quality of the RNA was initially assessed by measurement of the l260/280 ratio and agarose gel electrophoresis. RNA concentration was quantified by measuring absorbance at l260 nm using a NanoDropH ND-1000 (NanoDrop, Wilmington, DE), and the RNA quality was determined using a Nano-ChipH on an Agilent 2100 Bioanalyzer (Agilent, Palo Alto, CA). Reference RNAs (rat liver total RNA and rat liver mRNA) were included on each of the four chips used for analysis. The RNA analyzed in all samples was of good to excellent quality (results not shown). The average RNA integrity number (RIN) was 7.42 and the median was 7.7 (Range 6.1-8.8 with the exception of four samples that had RIN values between 2.5 and 3.4). Samples were characterized by distinct 18S and 28S rRNA peaks and RNA size distribution. Ten micrograms of RNA was used as starting material for each array.

Preparation of Bovine Reference RNA
A reference RNA was used to normalize signal intensity across all arrays [42,43,44,45]. This method was used rather than cohybridizing treatment and control samples to the same array to allow comparison of bovine intestinal response not only to Salmonella but also to other pathogens, including Brucella melitensis, Mycobacterium avium subsp. paratuberculosis. Each of these pathogens cause significant disease in ruminants and potentially humans. Comparison of intestinal host response to these pathogens was a first step in building a model of intestinal host pathway interactions with the intent of building an in silico interactome model. These studies and description of the interactome model will be described elsewhere. A reference RNA was chosen over genomic DNA as reference RNA is subjected to the same procedures, such as reverse transcription, as sample RNA thereby providing control for RNA or dye degradation during processing. Total RNA was isolated from Madin-Darby bovine kidney (MDBK) and bovine B lymphocyte (BL-3) cell lines (American Type Culture Collection, Manassas, VA) and fresh bovine cerebral cortex and cerebellum collected from surgery calves following euthanasia at the conclusion of each experiment. RNA from this combination of tissues was experimentally shown to produce cDNA that would hybridize to the majority of the open reading frames (ORFs) represented on the microarray [46]. Cell lines were grown in 150 cm 2 cell culture flasks with Eagle's minimum essential medium (E-MEM) (ATCC) supplemented with 10% heatinactivated fetal bovine serum (HI-FBS). Tissues were homogenized in ice-cold Tri-reagentH (Ambion). RNA concentration and quality from each sample was determined before and after pooling the samples. Total RNA isolated from three samples was pooled together in equal amounts, aliquoted and stored at 280uC until needed.

Construction and Annotation of Bovine cDNA Microarrays
A custom, bovine cDNA array consisting of 13,257 unique 70mer oligonucleotides representing 12,220 cattle ORFs was designed from normalized and subtracted cattle cDNA libraries. The oligos were printed in 150 mM phosphate buffer at 20 mM concentration in duplicate on aminosilane-coated glass slides at the W. M. Keck Center (University of Illinois at Urbana-Champaign). The oligos were annotated based on the GenBank accession number, when available. Further description of the microarray and its development has been previously published [47,48]. The gene array annotation was based on bovine Unigene Build #95.

Microarray Sample Preparation and Slide Hybridization
Microarrays were used to examine the transcriptional profiles of bovine intestinal epithelia (control and wild type or mutant infected) across seven time points (15 min, 30 min, 1, 2, 4, 8, and 12 hours). Experiments were performed in quadruplicate (four animals), generating a total of 84 slides with the array printed in duplicate on each slide. The labeling and hybridization procedures have been described previously [47]. Briefly, cDNA from bovine experimental samples (i.e. from infected and control loops) and cDNA generated from the bovine reference RNA sample were cohybridized to the previously described custom 13K bovine 70-mer oligoarray. The cDNA was reverse-transcribed using Superscript III reverse transcriptase (Invitrogen, Carlsbad, CA) and labeled with amino-allyl-UTP (Ambion, Austin, TX). The cDNA generated from the bovine reference RNA was labeled with Cy3 while the cDNA generated from the RNA samples isolated from the intestinal loops was labeled with Cy5. Cy3 and Cy5 dye esters were covalently linked to the amino-allyl group by incubating the samples with the dye esters in 0.1 M sodium carbonate buffer. Prior to hybridization, the microarrays were denatured by exposing to steam from boiling water for three seconds, UV cross-linked and then immersed in pre-hybridization buffer [5X sodium chloride, sodium citrate buffer (SSC), 0.1% sodium dodecyl sulfate (SDS) (Ambion, Austin, TX), 1% bovine serum albumin (BSA)] at 42uC for a minimum of 45 min. The arrays were then washed four times in RNase-, DNase-free, distilled water, immersed in 100% isopropanol for 10 seconds, and dried by centrifugation. Slides were hybridized at 42uC for approximately 40 hours in a dark, humid chamber (Corning, Corning, NY), washed for 10 min at 42uC with low stringency buffer [1 X SSC, 0.2% SDS] and then washed twice for 5 min each wash in a higher stringency buffer [0.1 X SSC, 0.2% SDS and 0.1 X SSC] at room temperature with agitation.

Microarray Data Acquisition and Analysis
Immediately after washing, the slides were scanned using a commercial laser scanner (GenePix 4100; Axon Instruments Inc., Foster City, CA). The spots representing genes on the arrays were adjusted for background and normalized to internal controls using image analysis software (GenePix Pro 4.0; Axon Instruments Inc.). Spots with fluorescent signal values below background were disregarded in all analyses. Samples were normalized against the bovine reference RNA signals across slides and within each slide (across duplicate spots). Regression analysis was performed to verify reproducibility of replicate arrays, all of which had an R 2 value of 0.70 or greater. This R 2 value may have resulted from poor RNA quality in a small number of the arrays (predominantly arrays generated from the 12 hour time point). Pairwise comparisons and Student's t test with Benjamini-Hochberg correction were performed using GeneSifter (VizX Labs, Seattle, WA), and Spotfire DecisionSite 9.0 (Spotfire, Inc., Somerville, MA). Genes were considered to be significantly altered if the average fold-change was at least 1.5, the adjusted p value was less than 0.05, the alteration was reproducible across replicates, and the magnitude of difference between conditions was greater than 50% of the fold-change obtained for any two replicate controls. Replicate spots, representing a single gene were also expected to ''agree.'' Comparative pathogenicity analysis and modeling was performed using an integrated platform termed the BioSignature Discovery System (BioSignatureDS TM ) (Seralogix, LLC, Austin, TX). This approach for genomic data analysis and modeling at the system biology level offers an integrated view of biological mechanisms and networks of interactions. Specifically for the analysis reported herein, the tools were used to: 1) determine significant gene modulations via a statistical inference z-score method employing a Bayesian estimate of variance [49] thresholding technique and fold-change; 2) conduct biological system level analysis employing dynamic Bayesian network models for scoring and ranking of metabolic pathways, signaling pathways and gene ontology (GO) groups; 3) conduct dynamic Bayesian candidate mechanistic gene analysis to identify genes within the network models that are most responsible for causing pathway and GO group perturbations; and 4) create a genetic network system model derived from the candidate mechanistic genes and their genetic interactions. Description of the computational techniques employed by BioSignatureDS TM for the analyses and modeling of the data reported herein is described in File S1. Individual gene zscores and fold-change data are provided in Tables S1 and S2 respectively.

Microarray Data Deposited in the Gene Expression Omnibus
All microarray data is MIAME compliant and the original images and extracted data were deposited in the Gene Expression Omnibus at the National Center for Biotechnology Information (http://www.ncbi.nlm.nih.gov/geo/), Accession #GSE13847.

Microarray Results Validation
Ten selected bovine genes with differential expression by microarray results were analyzed by quantitative RT-PCR (qRT-PCR). Two micrograms of RNA samples were reversetranscribed using TaqManH Reverse Transcription Reagents (Applied Biosystems). For relative quantification of target cDNA, samples were either run in individual tubes in SmartCycler II (Cepheid, Sunnyvale, CA) using one SmartMix bead (Cepheid) for 2-25 ml PCR reaction along with 20 ng of cDNA, 0.2X SYBR Green I dye (Invitrogen) or run using 20 ng of cDNA in Sybr Green PCR Master Mix (Applied Biosystems) on an ABI 7500 Real-time PCR System according to the manufacturer's instructions. In either case, Primer Express Software v2.0 (Applied Biosystems) was used to design forward and reverse primers that were used at a concentration of 0.3 mM (Table 1). For each gene tested, the individual calculated threshold cycles (C T ) were averaged among each condition and normalized to the C T of the bovine GAPDH, from the same cDNA samples before calculating the fold-change using the DDC T method [50]. For each primer pair, a negative control (water) and an RNA sample without reverse transcriptase (to determine genomic DNA contamination) were included as controls during cDNA quantification.

Salmonella Invasion and Inflammation in the Intestinal Epithelium
Four male, unrelated, Salmonella culture negative calves underwent ligated ileal loop surgery. Loops in each calf were inoculated with wild type S. Typhimurium (WT), a DsipA, sopABDE2 mutant (MUT), or sterile culture medium as a negative control (UI). Individual loops from each treatment were collected at seven time points (0.25, 0.5, 1, 2, 4, 8, and 12 hours) post inoculation. Each intestinal loop contained Peyer's patch which in the calf extends approximately 1.2 meters cranial from the ileocecal junction. The morphologic description of the ligated ileal loops and Salmonella infection in the calf intestine has been previously described [9]. Loops were consistent in cellular composition and response to infection throughout the length of intestine used. Prior to analyzing gene expression, bacterial invasion and uptake by M cells and phagocytes were assessed by bacterial counts measuring tissue-associated bacteria and analyzing morphologic changes in the tissue (Figure 1 and data not shown). There was a small but significant difference in bacterial numbers between WT and MUT infected loops at 30 minutes and 1 hour post inoculation (Figure 1). This difference may have resulted from the invasion defect that has been previously described for this mutant. There were no differences at later time points which could result from either extracellular bacteria or uptake of the mutant by M cells or other cells. Consistent with previous reports, the MUT inoculated loops had significant reductions in PMN migration into the intestinal mucosa and submucosa, inflammation score, and fluid secretion relative to the WT inoculated loops (data not shown) [32].
Transcriptional Profile of Bovine Host Peyer's Patches Inoculated with S. Typhimurium RNA was isolated from the mucosa and submucosa of the Peyer's patch of each intestinal loop. There were a total of 84 RNA samples that were used to generate indirectly labeled cDNA. The cDNAs were co-hybridized on a custom bovine microarray with indirectly labeled cDNA generated from a bovine reference RNA. The same bovine reference RNA was used throughout the experiment and was used for normalization of the microarrays. Gene expression changes associated with infection by WT and the MUT were compared with gene expression in the uninfected loops (UI) and with each other. Two methods of microarray data analysis were employed to determine gene expression changes: Student's t test with Benjamini & Hochberg (B&H) multiple hypothesis testing and a statistical inference z-score method employing a Bayesian estimate of variance [49,51]. Multiple hypothesis testing provides an index for defining the expected proportion of false positives (false discovery rate, FDR) among all significance tests to identify a set of ''candidate positives,'' of which a high proportion are likely to be true. This t test employed a B&H-corrected p value of ,0.05 (95% confidence) combined with a |fold-change| of 1.5 as the cutoff values for significance testing. Likewise, the Bayesian statistical inference method employed a |z-score| .1.96 (95% confidence) combined with a |fold-change| .2 as the threshold values for inference testing. Whereas significance testing is a case of post-hoc analysis, Bayesian inference generates probabilistic statements about hypotheses based on the microarray data (a priori estimation) [49]. Significance testing and Bayesian inference identified respectively 732 and 1046 genes that exhibited significant changes in expression in WT loops relative to UI loops (Table 2 and 3 and  detailed gene lists in Tables S4 and S5). In MUT loops there were 203 and 637 significantly altered genes as measured by significance testing and Bayesian statistical inference, respectively ( Table 2 and 3 and Tables S4 and S5). According to the Bayesian statistical inference analysis, there were 983 genes significantly altered (uniquely different) between WT and MUT (Table 3 and  Table S4). A surprising number of transcripts represented genes of unknown function. Within the 1046 altered genes identified by Bayesian inference, there were 444 transcripts of unknown function in WT loops (Table S4). Within the 637 altered genes identified by Bayesian inference, there were 300 transcripts of unknown function in MUT loops (Table S4). Significance of gene expression changes was consistent (i.e. overlapped) between significance testing and Bayesian inference in 348 genes in WT loops and only 11 genes in MUT loops (Table 4 and Table S5). We propose that the data normalization processing steps and the Bayesian variance estimation technique are the source for not having a higher number of overlapping altered genes between the two analysis methods. The Bayesian inference method employed: 1) a less stringent spot quality filtering technique; 2) a more sophisticated universal reference normalization method in conjunction with Lowess correction; and 3) a Bayesian variance estimator that infers a better prediction of the standard deviation for genes which have a low number of biological replicates. Overall, the Bayesian statistical inference method identified a larger number of altered genes than the t test alone [49].

Temporal Differences in Gene Expression
Temporally, WT infection induced a biphasic host response with a period of intense gene expression activity at one hour followed by a second phase at later time points of progressively  increased activity over time. In MUT loops there was no biphasic response. The number of significantly affected genes increased progressively over time with a pattern of primarily increased gene expression at later time points (Tables 2A and B). While there was substantial overlap between gene expression patterns exhibited by the WT loops and MUT loops, there were significant differences between the WT and MUT at individual time points (Tables S4  and S5). As shown by the representative examples in Figure 2, gene expression changes were consistent between microarray and qRT-PCR for both up-regulated and down-regulated genes. Microarray gene expression data were validated by qRT-PCR. Six genes determined to be significantly affected in the host response to infection by both methods of analysis (significance testing and Bayesian methods) were chosen for verification ( Figure 3). In addition to these genes, two genes identified by Bayesian inference as significantly regulated, the genes encoding apolipoprotein C3 and interleukin 6, were also tested by qRT-PCR and found to be significantly regulated by that method as well ( Figure 3G and H). This was true for all of the genes tested by qRT-PCR. Differences in gene expression due to infection, either with WT or MUT, that were detected by the microarray were consistent with differences detected by qRT-PCR both for genes with increased expression and genes with decreased expression relative to the negative control ( Figure 3 and data not shown). Additionally, 602 genes were represented on the microarray with more than one probe. In general, multiple probes for the same  Genes with significant z-score change and significant p value and 1.5-fold-change in response to infection with wild type Salmonella Typhimurium or an isogenic DsipA, sopABDE2 mutant. WT loops = wild type inoculated loops; MUT loops = DsipA, sopABDE2 mutant inoculated loops; Detailed gene list is provided in Table S5. For genes to be on this list, they had to be found significant based on both statistical analysis methods (t test and Bayesian statistical inference). Genes with common regulation are defined by the intersection of the sets of genes found significant for WT and for MUT conditions, whereas the genes that differ are the set of significant genes in WT that are not found in MUT condition. doi:10.1371/journal.pone.0026869.t004 gene showed consistent expression between the probes. Probes with inconsistent regulation are noted where significant.

Regulation of Chemokine Expression
Chemokines are signaling molecules responsible for recruitment of leukocytes. Previous work demonstrated increased expression of C-X-C chemokines during Salmonella infection in the bovine intestine [9,32,52,53]. For additional confirmation of microarray results, gene expression patterns from previously published work were compared to the expression patterns in the data presented here and were found to be consistent. In particular, increased expression of C-C motif and C-X-C motif chemokines including CCL2, CCL8, CXCL6 (GCP2), and IL8 (CXCL8), and their receptors including CCR1 and CCR3 was seen in both WT and MUT inoculated loops at 4, 8, and 12 hours (Table 5). Inclusion in this table required a |z-score| $2.24 at any one time and reflects 97.5% confidence in the data (note that no fold-change cutoff was applied). The magnitude of expression of these genes was greater in WT than in MUT, but there was no significant difference between these conditions. Expression of IL8, IL8RB, CXCL6, and CCL8 differed between WT and MUT at 8 hours between these conditions (Table 5). Increased gene expression of three additional chemokines, CXCL9, CXCL10, and CXCL11 was also seen in both WT and MUT infection but did not differ between the two conditions (Table 5). These genes encode chemokines induced by interferon-c (IFN-c) and Type I interferons (Interferon-a and Interferon-b). These chemokines are chemotactic for T cells through the chemokine receptor CXCR3. IFN-c is produced in response to Salmonella infection and is encoded by IFNG which was strongly induced in both WT and MUT loops in this study consistent with previously published results. Additional genes that were compared and consistent with previously published data in calves included TNF, IL1B, IL1RN, and IL17.

Dynamic Bayesian Modeling Analysis of Pathway Activation
To better understand the complex molecular interactions between host and pathogen, gene expression data were mapped to 219 molecular interaction pathways using the Dynamic Bayesian Gene Group Activation (DBGGA) method described more fully in File S1. The DBGGA method is similar to Bayesian network modeling methods employed by other investigators [54,55,56,57,58,59,60]. DBGGA creates a Dynamic Bayesian network (DBN) model for each pathway based on KEGG [61,62] while a naïve Bayesian classifier [63] type network model is created for each GO category. Employing DBGGA, significantly activated pathways under both WT and MUT conditions in comparison to uninfected controls were identified and a differential comparison between WT and MUT was made. As an example, the WT S. Typhimurium modeling result for the Tolllike Receptor Pathway is shown in Figure 4 where the discovered mechanistic genes are indicated by concentric rings around the node. In this figure, the red rings indicate mechanistic genes that are unique to the WT condition while the blue rings are mechanistic genes in common between the WT and MUT conditions. The graphical model is used to view the state at any time point of all measured genes upstream and downstream of each mechanistic gene candidate. The products of mechanistic genes control key regulatory points in pathways through a variety of methods including altering gene transcription or translation or through post-translation processes such as protein phosphorylation. Alteration of mechanistic genes induces substantial effects on biological processes. See File S1 for the full mechanistic gene definition found significant from DBGGA Bayesian z-score method.
There were 193 pathways that were significantly altered between WT and MUT infected loops (WT vs MUT) at various time points throughout the experiment (Table 6). Pathways within the categories of cell mobility, cell communication, cell growth and death, infectious diseases, immune system, membrane transport, and signal transduction were compared. There were 49 pathways represented in these categories ( Figure 5). Note that WT versus MUT compares the WT expression values to the MUT expression values and the MUT is considered the control condition. Within the 49 perturbed pathways (either activated or repressed) were several that have been extensively studied in Salmonella infection including Regulation of Actin Cytoskeleton, Apoptosis [64,65], and the Toll-like receptor (TLR) signaling pathway [66,67]. A temporal pattern of activation or repression for these pathways was generated ( Figure 5). These 49 pathways are proposed to represent target pathways for intervention in Salmonella-induced diarrhea. Because a temporal pattern of increased gene activity at one hour was previously noted in WT loops, differences between WT and MUT loops at early time points (15 minutes, 30 minutes, 1 hour) and late time points (2 hours, 4 hours, and 8 hours) were explored in greater depth (Tables 2 and 3). The 12 hour time point was not evaluated because of reduced signal quality in the microarrays from the 12 hour time points. The pathways that differed in significant regulation between early and late time points at a minimum of at least one time point are summarized in Table 6. Pathways within the categories of cell mobility, cell communication, cell growth and death, infectious diseases, immune system, membrane transport, and signal transduction were compared for differences between WT and MUT in early and late time points (Figure 6). Within this pathway sub-group, Apoptotic Signaling in Response to DNA Damage, Hedgehog Signaling and Notch Signaling differentiated WT from MUT loops at early time points. Multiple pathways differentiated the two conditions at late time points including the Complement and Coagulation Cascades, Cytokine Receptor Signaling, and Epithelial Cell Signaling in Helicobacter pylori (Table 6; Figure 6). At eight hours, there were very few pathways that were significantly different between WT and MUT. To explore differences between WT and MUT further, a set of the significantly affected pathways within the subgroup was investigated to identify key regulatory genes that differed in expression during early infection (Figure 7). Pathways with primary functions in metabolism were not considered in this group. This subset included 1) Adherens Junction; 2) Apoptosis; 3) BRC signaling; 4) CCR3 Signaling in Eosinophils; 5) Calcium Signaling; 6) ErbB Signaling; 7) MAPK Signaling; 8) Phosphatidylinositol Signaling System; 9) Regulation of Actin Cytoskeleton; 10) TGF-b Signaling Pathway; 11) Tight Junction; 12) Toll-like Receptor Signaling; 13) Wnt Signaling and; 14) mTOR Signaling pathways. The changes in activation status between WT and UI, MUT and UI, and WT and MUT over time can be seen in Figure 7. The inter-relationship between genes in these pathways is summarized in the host disease model in Figure 8. The method for implementing the disease (system) model is described in File S1.
All 14 of the pathways exhibited activation in the MUT loops at 2 hours, 4 hours, and 12 hours. This was not the case in WT and Tight Junctions were all activated. These differences between WT and MUT are consistent with differential regulation of host response through the five deleted effector proteins. These seven pathways are connected through regulation of the actin cytoskeleton and phosphatidylinositol signaling and despite early downregulation, by one hour all but BRC Signaling, Wnt Signaling and Phosphatidylinositol Signaling System were activated. Within the 14 selected pathways, there were 103 genes that were common key mechanistic genes within the pathways. Of these genes, 11 were identified as key intersection points having mechanistic functions in more than five of the 14 pathways (Table S7). These mechanistic genes are AKT2, JUN, MAPK1, NFKB1, NFKB2, PIK3C2G, PIK3C3, PIK3R2, PRKCG, RAC1, and RHOA. Based on Bayesian z-score, only AKT2, RHOA, and NFKB1 at 1 hour were appreciably regulated (|Bayesian z-score| .2.24). None of these genes were significantly modulated at later time points, however small alterations in expression of these key genes can have significant effects on other genes as reflected in the dynamic Bayesian model.
Although it is not possible to separate the effects of SipA, SopA, SopB, SopD and SopE2 in the data presented here, the known direct and indirect effects of SopE2 and SopB provide crucial links between these 11 mechanistic genes and 14 pathways and the differences between the early phase of WT and MUT infection. The role of SopB during Salmonella infection is complex. SopB is a phosphoinositol phosphate phosphatase required for activation of Akt, early trafficking of the Salmonella-containing vacuoles (SCV), and activates an SH3-containing guanine nucleotide exchange factor (SGEF) which is an exchange factor for RhoG [33,35,68] SopE and SopE2 catalyze nucleotide exchange on Rho GTPases including Rac1, Cdc42, and RhoG and activate MAPK and NFkb signaling [34]. These roles of SopE2 and SopB during the early stage of Salmonella infection explain the presence and absence of a biphasic response in gene expression in WT and MUT infection respectively. Rac1 and MAPK1, link multiple pathways including the Adherens Junction, BRC signaling, CCR3 Signaling in Eosinophils, MAPK Signaling, Regulation of Actin Cytoskeleton, TGF-b Signaling, Toll-like Receptor Signaling, mTOR Signaling and Wnt Signaling pathways. With the addition of Akt2, Apoptosis, ErbB Signaling, and Tight Junction pathways can be added to the model. Calcium Signaling can be added through the effects of NFKB1 and NFKB2. The Phosphatidylinositol Signaling System cannot be directly linked through these genes but SopB can hydrolyze inositol phosphates and phosphoinositides in vitro and is required for formation of phosphatidylinositol 3-phosphate on SCVs [29,33,35,69]. Within the Phosphatidylinositol Signaling System, PIK3C2G, PIK3C3, and PRKCG may represent attractive targets for future studies.

Discussion
Using microarray analysis, we were able to generate a temporal gene expression profile for 12,257 genes concurrently.
Analysis of a relatively large number of genes reduces investigator bias and allows the generation of a gene expression model. Selection of the genes printed on the microarray was independent of our expectations about the effects of Salmonella infection or the role of SipA, SopA, SopB, SopD or SopE2. We expected that this method would yield a dataset that would include genes with a known role in Salmonella infection as well as identification of expression changes in genes whose role in host response was unexpected. We identified alterations in 38 genes whose functions are as yet unknown. Undoubtedly, the function of these transcripts will be better understood as our understanding of the bovine genome continues to improve and may enhance our understanding of host response to Salmonella infection. The microarray data presented here were validated by confirmation of a portion of the altered genes using real time PCR as well as consistency of our results with previously published data. The limitations of this approach are that only changes in expression of genes spotted on the array are measured and that the intestinal tissue used is comprised of multiple cell types, such that it is not possible to differentiate gene expression from one cell type, such as PMNs, from another type, such as epithelial cells; however, we have employed lasercapture microdissection in other studies to evaluate specific gene expression in specific intestinal cell types. Additionally, the data from this study are limited to gene transcription and do not describe alterations in protein expression or state of activation or phosphorylation that may affect activity in the pathways discussed. Studies to address these issues are currently underway in our labs. Ideally, unbiased evaluation of gene expression would include expression data for all genes but at the time that this study was performed, deep-sequencing techniques were cost-prohibitive. Upcoming studies are in progress to utilize these methods and to further develop our understanding of bovine host response to Salmonella.
Limitations of our study included the use of four animals instead of a larger number of animals. Previous studies in our laboratory indicated that four animals would be sufficient to detect significant differences between wild type Salmonella and the sipA, sopABDE2 mutant [32]. Additionally, we wished to limit the number of animals sacrificed for the study and there were financial limitations on the number of animals and number of microarrays that we could perform. Recognizing that low number of replicates is a common problem in microarray studies, techniques to compensate for low replicates were employed. Specifically, an intensity dependent smoothing variance function was used in the Bayesian inference significance tests [49]. This method borrows strength across the large number of genes measured across all time points. The genes were ranked by intensity levels and a new standard deviation computed by finding the mean of the variances of the next 50 higher genes and the next 50 lower genes surrounding the gene of interest. This produced a Bayesian estimator of the standard deviation that is found to be a reliable approach for statistical testing compensating for the low replicates.
The primary result of this work was the establishment of a temporal gene expression profile and development of an in silico Table 6. Pathways that differed in regulation between early and late time points in WT vs MUT based on Bayesian z-score described in Table S2.  model for host response to both wild type S. Typhimurium and a DsipA sopABDE2 mutant. SipA, SopA, SopB, SopD, SopE2 were previously demonstrated to have a combined role in the generation of fluid, induction of cytokines and C-X-C chemokines, and PMN recruitment in bovine ligated ileal loops [9,12,32,70]. The major difference between WT and MUT infections was the presence of a biphasic response with increased activity at one hour post infection in the intestinal loops infected with WT. Using dynamic Bayesian network modeling (DBN) in conjunction with data from KEGG, the Gene Ontology and GenBank, patterns of pathway activation were established for both WT and MUT infection. The differences in gene activity seen in the WT but not MUT were studied in a subset of pathways selected from categories of cell mobility, cell communication, cell growth and death, immune system, membrane transport, and signal transduc-tion. These pathways could be linked through the interactions of two Salmonella effectors, SopB and SopE2, and the products of three bovine genes, RAC1, MAPK1, and AKT2. Although this subset of interactions does not explain the entire cascade of interactions observed in the data presented herein, they do serve to validate the model and the DBN approach Our original hypothesis and intent was to generate a profile of host gene expression associated with inflammation and diarrhea using two Salmonella strains, wild type and a mutant that generates reduced fluid secretion and inflammatory responses during intestinal infection. Although we have generated gene expression profiles for these two strains, we have not correlated activation or repression of these pathways with the production of fluid in the bovine intestine as a measure of diarrhea. In vitro studies employing gene silencing are underway to determine the roles of the 49 Figure 6. Temporal differences between pathways scores for wild type S. enterica serovar Typhimurium and the DsipA, sopABDE2 mutant. Pathway activation in wild type and mutant inoculated loops was compared and evaluated for temporal changes. Pathways activated in early time points (15 minutes, 30 minutes, and 60 minutes) were summarized (A) and pathways activated in late time points (120, 240, and 480 minutes) were summarized (B). Pathways were limited to the categories of cell mobility, cell communication, cell growth and death, infectious diseases, immune system, membrane transport, and signal transduction. The score magnitudes are shown as a gradient color from light to bright red for induced and from light to bright green for suppressed pathway activity. pathways identified as differing between the WT and MUT conditions. In addition to the effects on pathway regulation and mechanistic gene identification, this study identified effects on individual genes in response to infection. We identified increased expression of CXCL9, CXCL10, and CXCL11. This is consistent with previous work in both rhesus macaques and mice which demonstrated that IL17, IL12, and CXCL10 are increased in response to wild type S. Typhimurium and that this response is ablated in animals with diminished numbers of CD4+ T cells in the intestinal mucosa [71]. Identification of IFNc as a mechanistic gene in both WT and MUT loops supports a role for T cells and natural killer cells (both of which produce this powerful cytokine) which is independent of SPI-1 secreted effectors. Increased expression of CXCL9, CXCL10, and CXCL11 was observed in WT and MUT loops ( Table 5). It has been suggested from work in murine cell lines that TLR4 signaling induces expression of Type I Interferons (interferon-a and interferon-b) [72,73]. Increased expression of CXCL9, CXCL10, and CXCL11 suggests either that IFN-c induces the expression of these chemokines or that Type I interferons are increased during the first 12 hours of infection in the bovine intestine. Support for a role of Type I interferons includes studies in mice, which demonstrate that IFNab are increased by Salmonella infection [73,74]. Additional indirect evidence includes the increased expression of genes such as MX2 that is induced during viral infections by IFN-ab [75]. IFN-ab is critical for the function of natural killer cell activity [76,77,78]. Further study is warranted to define the source of expression of these chemokines and their role in recruitment of T cells and natural killer cells to the intestine during Salmonella infection.
In conclusion, this report catalogued the temporal pattern of intestinal gene expression in response to wild type S. Typhimurium in a natural Salmonella host, the calf, in response to S. Typhimurium infection. We have identified a temporal gene expression pattern comparing wild type S. Typhimurium and a strain attenuated for production of diarrhea and inflammation, and we have catalogued the differences observed due to the combined effects of SipA and SopABDE2 on host gene expression. Pathway analysis and the construction of a system level genetic interaction model of the host response have improved our ability to identify important pathways, genetic mechanisms, and novel genes previously not associated with S. Typhimurium infection. Further, since the genetic model defined by a dynamic Bayesian network can be trained and used for modeling the effects of silencing selected genes and inferring the possible host response, it will be an invaluable tool for: 1) selecting potential intervening gene targets; 2) guiding future experiments for validating host response models; and 3) identifying dynamic gene expression patterns for early diagnostics.

Supporting Information
Table S1 z-score for All Genes. The z-scores for the individual comparison of gene expression between wild type inoculated loops and control loops, mutant inoculated loops and control loops, and between wild type inoculated and mutant    Table 3 Summary of significant genes determined by Bayesian inference zscore using cutoff |z-score| $1.96 and |fold-change| $1.5 in response to infection with wild type Salmonella Typhimurium or an isogenic DsipA, sopABDE2 mutant. UI = Uninfected control loops inoculated with Luria-Bertani broth. WT = Wild type inoculated loop versus uninfected control loops. MUT = DsipA, sopABDE2 mutant inoculated loops versus uninfected control loops. (XLSX)  Table 4 Summary of significantly altered genes that overlap between the two statistical methods (Table 2 and 3). Genes with significant zscore changes and significant p value and 1.5-fold-change in response to infection with wild type Salmonella Typhimurium or an isogenic DsipA, sopABDE2 mutant. UI = Uninfected control loops inoculated with Luria-Bertani broth. WT = Wild type inoculated loop versus uninfected control loops. For genes to be on this list, they had to be found significant based on both statistical analysis methods (t test and Bayesian statistical inference). (XLSX) Table S6 z-score of pathways that differed significantly between loops inoculated with wild type Salmonella Typhimurium or an isogenic DsipA, sopABDE2 mutant over time. Pathways are based on Kyoto Encyclopedia of Genes and Genomes descriptions and pathway numbers. Symbol_Count refers to the total number of gene symbols present in a pathway. Observed_Count refers to the number of genes that were observed to be significantly altered at a minimum of one time point. Probe_Count refers to the number of spots on the array that were significantly altered at a minimum of one time point (time in minutes). (XLSX)