Systems Biology Analysis of Temporal In vivo Brucella melitensis and Bovine Transcriptomes Predicts host:Pathogen Protein–Protein Interactions

To date, fewer than 200 gene-products have been identified as Brucella virulence factors, and most were characterized individually without considering how they are temporally and coordinately expressed or secreted during the infection process. Here, we describe and analyze the in vivo temporal transcriptional profile of Brucella melitensis during the initial 4 h interaction with cattle. Pathway analysis revealed an activation of the “Two component system” providing evidence that the in vivo Brucella sense and actively regulate their metabolism through the transition to an intracellular lifestyle. Contrarily, other Brucella pathways involved in virulence such as “ABC transporters” and “T4SS system” were repressed suggesting a silencing strategy to avoid stimulation of the host innate immune response very early in the infection process. Also, three flagellum-encoded loci (BMEII0150-0168, BMEII1080-1089, and BMEII1105-1114), the “flagellar assembly” pathway and the cell components “bacterial-type flagellum hook” and “bacterial-type flagellum” were repressed in the tissue-associated B. melitensis, while RopE1 sigma factor, a flagellar repressor, was activated throughout the experiment. These results support the idea that Brucella employ a stealthy strategy at the onset of the infection of susceptible hosts. Further, through systems-level in silico host:pathogen protein–protein interactions simulation and correlation of pathogen gene expression with the host gene perturbations, we identified unanticipated interactions such as VirB11::MAPK8IP1; BtaE::NFKBIA, and 22 kDa OMP precursor::BAD and MAP2K3. These findings are suggestive of new virulence factors and mechanisms responsible for Brucella evasion of the host's protective immune response and the capability to maintain a dormant state. The predicted protein–protein interactions and the points of disruption provide novel insights that will stimulate advanced hypothesis-driven approaches toward revealing a clearer understanding of new virulence factors and mechanisms influencing the pathogenesis of brucellosis.

To date, fewer than 200 gene-products have been identified as Brucella virulence factors, and most were characterized individually without considering how they are temporally and coordinately expressed or secreted during the infection process. Here, we describe and analyze the in vivo temporal transcriptional profile of Brucella melitensis during the initial 4 h interaction with cattle. Pathway analysis revealed an activation of the "Two component system" providing evidence that the in vivo Brucella sense and actively regulate their metabolism through the transition to an intracellular lifestyle. Contrarily, other Brucella pathways involved in virulence such as "ABC transporters" and "T4SS system" were repressed suggesting a silencing strategy to avoid stimulation of the host innate immune response very early in the infection process. Also, three flagellum-encoded loci (BMEII0150-0168, BMEII1080-1089, and BMEII1105-1114, the "flagellar assembly" pathway and the cell components "bacterial-type flagellum hook" and "bacterial-type flagellum" were repressed in the tissue-associated B. melitensis, while RopE1 sigma factor, a flagellar repressor, was activated throughout the experiment. These results support the idea that Brucella employ a stealthy strategy at the onset of the infection of susceptible hosts. Further, through systems-level in silico host:pathogen protein-protein interactions simulation and correlation of pathogen gene expression with the host gene perturbations, we identified unanticipated interactions such as VirB11::MAPK8IP1; BtaE::NFKBIA, and 22 kDa OMP precursor::BAD and MAP2K3. These findings are suggestive of new virulence factors and mechanisms responsible for Brucella evasion of the host's protective immune response and the capability to maintain a dormant state. The predicted protein-protein interactions and the points of disruption provide novel insights that will stimulate advanced hypothesis-driven approaches toward revealing a clearer understanding of new virulence factors and mechanisms influencing the pathogenesis of brucellosis.

INTRODUCTION
Brucella, an aerobic non-motile Gram-negative coccobacillus, is the etiological agent of brucellosis, a worldwide anthropozoonotic infectious disease that causes chronic infections with persistent or recurrent bacteremia in susceptible hosts, and mid-to late gestation abortion in pregnant animals. At present, there are 11 recognized species within the genus Brucella based on preferential host specificity (O'Callaghan and Whatmore, 2011;Whatmore et al., 2014). Goats and sheep are the preferred hosts for Brucella melitensis (Alton, 1990), although this pathogen also infects cattle depending on specific epidemiological conditions (Kalher, 2000;Banai, 2010;Alvarez et al., 2011;Liu et al., 2012;Wareth et al., 2014). B. melitensis is also the most pathogenic and the most frequent causative agent of brucellosis in humans (Traxler et al., 2013). Clinically, human brucellosis can be an incapacitating disease that results in intermittent fever, chills, sweats, weakness, myalgia, osteoarticular complications, endocarditis, depression, and anorexia with low mortality (Dean et al., 2012).
The predominant route for B. melitensis penetration after natural exposure is the alimentary tract (Adams, 2002). Usually B. melitensis enters through the oral mucosa and colonizes the lymph nodes that drain the facial area (Carpenter, 1924;von Bargen et al., 2015); although several studies have isolated Brucella from different sections of the alimentary tract and feces revealing the possibility that brucellae invade in multiple sites of the gastrointestinal tract (Carpenter, 1924;Davis et al., 1988). We previously found that under experimental conditions, B. melitensis was able to invade the bovine host through the domed epithelium of jejuno-ileal Peyer's patches followed by rapid systemic dissemination (Rossetti et al., 2013). The calf ligated jejuno-ileal loop model has been demonstrated to be a very useful model to study in vivo natural host:infectious agent molecular and morphological initial interactions (Santos et al., 2002;Khare et al., 2009Khare et al., , 2012Winter et al., 2010;Lawhon et al., 2011;Rossetti et al., 2013), a subject that has not been broadly studied in brucellosis.
Brucella lack several classical bacterial virulence factors such as exotoxins, cytolysins, a capsule, fimbriae, plasmids, resistant strains, lysogenic phages, antigenic variation, or endotoxic lipopolysaccharide among others (Moreno and Moriyón, 2002). Brucella has developed a stealthy strategy to avoid being recognized and successfully infect hosts. Succinctly stated, Brucella circumvents strong innate immune responses, obstructs the direct action of bactericidal substances, resists destruction by professional phagocytes and maintains the host cells alive to establish long lasting infections (Barquero-Calvo et al., 2007). Although, significant advances have been made lately (de Figueiredo et al., 2015), the molecular pathogenesis of brucellosis is still incompletely understood. To date, fewer than 200 gene products have been identified as Brucella virulence factors (He, 2012), and very few related to adhesion and invasion function have been characterized. For instance, the heat shock protein 60 (Hsp60) family proteins on the surface of Brucella abortus have been found to bind macrophage and intestinal M cells cellular prion protein (PrP c ) before internalization (Watarai et al., 2003;Nakato et al., 2012). The SP41 (UgpB) protein that has significant homology with the glycerol-3-phosphatebinding ABC transporter protein interacts with cellular sialic acid residues, facilitating efficient host invasion (Castaneda-Roldán et al., 2006). Other recently characterized proteins in the brucellae genome, such as, BmaC (Posadas et al., 2012), BtaE (Ruiz-Ranwez et al., 2013a, BtaF (Ruiz-Ranwez et al., 2013b) interact with different components of the extracellular matrix, while novel adhesion-encoding regions inv (Alva-Perez et al., 2014), or BigA  have been demonstrated to promote adhesion and invasion, but their target host molecules have not been identified yet.
A well-recognized Brucella virulence factor is the twocomponent response regulator, BvrR/BvrS, that modulates the host cell cytoskeleton upon Brucella invasion (Sola-Landa et al., 1998) and regulates the Brucella OMP expression (Guzmán-Verri et al., 2002). Dysfunction of BvrR/S response regulator system induces mutant strains with reduced invasiveness and failure to replicate and survive intracellularly. Brucella lipopolysaccharide is also a confirmed virulence factor , that prevents complement-mediated bacterial killing (Allen et al., 1998;Tumurkhuu et al., 2006), provides resistance against antimicrobial peptides such as defensins, lysozyme, and lactorerrin (Martínez de Tejada et al., 1995) and inhibits cell death (Pei and Ficht, 2004;Pei et al., 2006). Additionally, Brucella LPS masks recognition of the pathogen-associated molecular patterns (PAMPs) from immune-receptor recognition, and as a consequence, impedes, or attenuates proinflammatory responses and immune system activation (Forestier et al., 2000;Jiménez de Bagués et al., 2004). Simultaneously, the type four secretion system (T4SS) is also a key virulence factor for Brucella intracellular survival , persistent infection in mice and induction of the host immune response (Rolan and Tsolis, 2007;Roux et al., 2007). This virulence factor, encoded by a virB operon, is induced by early phagosome acidification after phagocytosis (Boschiroli et al., 2002;Celli et al., 2003) to translocate effector molecules directly into the host cell cytoplasm. Additional investigations have identified several of these effector proteins, although most of their functions remain undefined (de Jong et al., , 2013de Barsy et al., 2011;Marchesini et al., 2011Marchesini et al., , 2016Salcedo et al., 2013;Ke et al., 2015;Del Giudice et al., 2016). The secretion systems and secretomes of Brucella were recently computationally analyzed, resulting in the prediction of 29 host-pathogen specific interactions between cattle and B. abortus and 36 host-pathogen interactions between sheep and B. melitensis proteins (Sankarasubramanian et al., 2016). The two-component system RegB/A, including the aceA encoding isocitrate lysase, has been found to play a critical role in the persistence and in vivo pathogenicity of Brucella suis (Abdou et al., 2017).
An additional virulence mechanism used by Brucella to survive intracellularly is the periplasmic compound cyclic B-1,2 glucan, that interferes with cellular trafficking and maturation of the Brucella-containing vacuole by disrupting cholesterol-rich lipid rafts present on phagosomal membranes and preventing the phagosome-lysosome fusion (Arellano-Reynoso et al., 2005;Martirosyan et al., 2012). Other virulence elements reported to sustain a chronic infection include: phosphatidylcholine, a phospholipid compound abundant in eukaryotic cell membranes that facilitate Brucella avoidance of host recognition Conde-Alvarez et al., 2006); PrpA, a prolineracemase family compound that elicits B lymphocyte polyclonal activation (Spera et al., 2006); BtpA and BtpB, Brucella TIRcontaining effector proteins that suppress innate immunity and modulate host inflammatory responses during infection (Salcedo et al., 2008(Salcedo et al., , 2013; MucR, a transcriptional regulator involved in Brucella metabolism, cell wall/envelope biogenesis, replication, type IV secretion system, quorum sensing system, and stress tolerance (Dong et al., 2013) and a flagellar appendage, required for virulence in a mouse infection model (Fretin et al., 2005).
Most of these virulence factors were characterized individually without considering how they are temporally and coordinately expressed or secreted during the infection process. Recently, several experiments have been performed to more fully understand the sequential expression and coordinated regulation of the infection process (Kohler et al., 2002;Al Dahouk et al., 2008;Rambow-Larsen et al., 2008;Lamontagne et al., 2009;Rossetti et al., 2009;Wang et al., 2009;Viadas et al., 2010;Weeks et al., 2010;Hanna et al., 2013;Tian et al., 2013). These experiments using in vitro culture media, infected cell cultures, or infected mice were successful for generating initial hypotheses to enhance the understanding of the pathogenesis of brucellosis.
Here, we describe an integrative approach of experimentation and computation to analyze the in vivo temporal transcriptional profile of B. melitensis during the first 4 h of the interaction with a naturally susceptible host, using the established calf jejuno-ileal loop model of infection. We then performed a system-level analysis by applying both a traditional statistical differential analysis to determine significance of B. melitensis gene expression and a pathway and gene ontology (GO) analysis that employed a dynamic Bayesian network (DBN) technique (Khare et al., 2016) to identify perturbations trends over time. The fundamental concept of systems biology is to: (1) perturb a system-(time-course B. melitensis infected bovine Peyer's patch), (2) measure systems-wide responses-(B. melitensis and bovine transcriptomes), and (3) integrate measured responses into a model-(host:pathogen::Bovine:B. melitensis interactome model) to understand the observations and iteratively predict novel interactions and perturbations. The system-level analyses aided understanding of the strategies exploited by B. melitensis to invade, survive, and replicate intracellularly; and to identify perturbations of major genes modulating critical cellular pathways in the pathogenesis of brucellosis. Further, through systems-level in silico hostpathogen protein-protein interactions (PPIs) simulation (see File S1), we were able to make inferred predictions of interactions of close apposition with specific B. melitensis expressed genes/proteins to plausible host (bovine) pathway points of disruption or perturbations. The predicted PPIs and the points of disruption provide novel insights that will stimulate advanced iterative hypothesis-driven approaches toward revealing a clearer understanding of new virulence factors and mechanisms contributing to the evasion of the host's protective immune responses.

Infection Model
The in vivo infection model for Brucella was described previously (Rossetti et al., 2013). Briefly, five bovine jejuno-ileal segments from four calves were inoculated intraluminally with 3 ml of a suspension containing 1 × 10 9 CFU of B. melitensis 16 M/ml (total of 3 × 10 9 CFU) at late-log growth phase cultured in F12K medium (ATCC R ) supplemented with 10% heatinactivated fetal bovine serum (HI-FBS) (ATCC R ). One infected segment was removed at every time point (0.25, 0.5, 1, 2, 4 h post-inoculation from each of the four calves.), and six to ten 6 mm biopsy punches were collected from each segment.

RNA Isolation, Labeling, and Hybridization
RNA isolation, labeling, and hybridization procedures were performed as described in previous experiments (Rossetti et al., 2010(Rossetti et al., , 2011b. Total RNA from B. melitensis-infected bovine Peyer's patches was extracted according to the TRI-Reagent manufacturer's instructions. Tissue-associated B. melitensis total RNA was initially enriched (MICROBEnrich R , Ambion) and then amplified from 30 µg of total RNA from B. melitensisinfected bovine Peyer's patches (Rossetti et al., 2010). Briefly, the enriched RNA was precipitated in 100% ethanol at −20 • C, washed and re-suspended in 25 µl of DEPC-treated water (Ambion). Immediately, the total amount of RNA was linearly amplified in a 3 step-protocol. First, RNA was reverse transcribed to cDNA using B. melitensis genome direct primers (BmGDPs), T7 promoter-template switching primer (T7-TS) (Sigma Genosys, The Woodland, TX) and Moloney Murine Leukemia Virus Reverse Transcriptase (Clontech, Palo Alto, CA). In the next step, the second-strand cDNA was synthesized and purified (Qiagen, Valencia, CA), followed by concentration in a speed-vac with no heat. In the last step, the in vitro transcription, was performed using the double-stranded cDNA as the template and T7 polymerase (Ambion). Then, 10 µg of each experimental sample (n = 44, i.e., 4 were in vitrogrown cultures of B. melitensis at late-log phase of growth; 20 were enriched and amplified B. melitensis RNA from total RNA from infected bovine Peyer's patches; and an additional 20 were from total RNA of infected bovine Peyer's patches) were reverse transcribed overnight to amino-allyl cDNA using 1.5 ug of B. melitensis genomic directed primers (BmGDPs) (Rossetti et al., 2010), labeled with Cy3 (Amersham Pharmacia Biosciences, Piscataway, NJ), mixed with 0.5 µg of Cy5 labeled B. melitensis gDNA, and applied to a custom 3.2K B. melitensis oligoarray (Weeks et al., 2010). Since the enrichment procedure does not eliminate host RNA, total RNA from B. melitensisinfected bovine Peyer's patches were also reverse transcribed, labeled and hybridized on B. melitensis oligo microarray, due to a potential concern that eukaryotic RNA present in enriched and amplified samples could possibly overlap with sequences of the B. melitensis transcripts and cross hybridized with probes on B. melitensis oligo microarrays, resulting in falsely detected pathogen genes. The isolation and labeling of B. melitensis gDNA has been described in detail elsewhere

Data Acquisition, Normalization, and Microarray Data Analysis
Immediately after washing, the dried slides were scanned using a commercial laser scanner (GenePix 4100; Axon Instruments Inc., Foster City, CA). Scans were performed using the autoscan feature with the percentage of saturated pixels set at 0.03%. The genes represented on the arrays were adjusted for background and normalized to internal controls using image analysis software (GenePixPro 6.0; Axon Instruments Inc.). Genes with fluorescent signal values below background were disregarded in all analyses. Arrays were initially normalized against B. melitensis genomic DNA, and the resulting data were analyzed and modeled using an integrated platform termed the BioSignature Discovery System (BioSignatureDS R ) (Seralogix, LLC, Austin, TX; http://www. seralogix.com) explained in detail elsewhere (Lawhon et al., 2011;Rossetti et al., 2013). The tissue-associated B. melitensis gene expression at every time point (0.25-4 h) was compared to the gene expression of the inoculum (i.e., in vitro-grown cultures of B. melitensis at late-log phase of growth, cultured in F12K medium with 10% HI-FBS; n = 4). Significantly expressed genes were determined with the z-test (p < 0.025) (enhanced with Bayesian methods of variance estimation) after subtracting those genes also expressed at statistically significant levels when total RNA of B. melitensis-infected bovine Peyer's patches was compared to the gene expression of the inoculum. BiosignatureDS tools for statistical z-score gene thresholding, Brucella pathway and gene ontology (GO) perturbation scoring (scored using Bayesian Information Criterion and transformed to z-score), and mechanistic gene identification were used for the comprehensive analysis performed in this study. A specialized application was developed to implement algorithms that integrate multiple sources of prior biological knowledge (PBK) into the inference of host-pathogen protein-protein interaction (PPIs) prediction (see File S1 for complete details). Briefly, we adopted three algorithmic methods for the identification of candidate interaction points for use in network learning between the host and the pathogen from in vivo gene expression data. These algorithmic methods were: (1) a sequencesimilarity interaction transference procedure; (2) structural protein domain-based algorithm; and (3) a functional geneontology-based algorithm. Gene candidates for inclusion in our interaction prediction process were selected based on interpretation of pathway and GO analyses conducted by our Dynamic Bayesian Network methodology. The B. melitensis gene transcriptome was employed and ≈600 bovine host genes selected from 12 perturbed and immune response relevant pathways and 10 GO terms to form gene sets representing two "unconnected" system models for starting the "interactome model" network learning process. Those algorithms yielded 348, 68, and 295 potential host-pathogen PPIs, respectively, that comprised the set of potential interactions at the interface of the pathogen and host systems. These potential interactions were then included into the Bayesian host-pathogen network structure learning algorithm. The method employs model structures to initialize learning with biologically relevant structures and utilize actual time-course co-expressed gene and other "omic" data from pathogen and host to search for a set of structures in which the data best fit. Microarray data and metadata are deposited in the Gene Expression Omnibus at the National Center for Biotechnology Information (http://www.ncbi.nlm.nih.gov/geo/) Accession #GSE89053. For microarray results validation, six randomly selected genes with consistently differential expressed from 15 min to 4 h postinfection (p.i.) by microarray results, were analyzed at every time point by quantitative RT-PCR (qRT-PCR) following the protocol described elsewhere (Rossetti et al., 2011b). Primers (Sigma Genosys) of tested genes were designed by Primer Express Software v2.0 (Applied Biosystems) ( Table S1). For each gene tested, the individual calculated threshold cycles (Ct) were averaged among each condition and normalized to the Ct of the 16S rRNA gene from the same cDNA samples before calculating the fold change using the Ct method (Livak and Schmittgen, 2001). 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 quantitation. Statistical significance was determined by Student's t-test and expression differences considered significant when P < 0.05. As gene expression by microarray and qRT-PCR were based on zscore and fold-change, respectively, array data were considered valid if the fold change of each gene tested by qRT-PCR was expressed in the same direction as determined by microarray analysis.

B. melitensis Transcriptome Is Perturbed at the Onset of the Infection Process
We previously reported that the number of B. melitensis 16M organisms after intraluminal inoculation increased from 15 min to 4 h p.i. (Rossetti et al., 2013) in this model. To study pathogen alterations in gene expression, pathway, and GO perturbations during the initial infection process, Brucella RNAs extracted from infected bovine Peyer's patches at different times p.i. were hybridized on B. melitensis microarrays and analyzed. As expected, the traditional z-score analysis (|2.24|, 97.5% confidence) identified a total of 2,356 different B. melitensis genes (1,221 up-regulated vs. 1,135 down-regulated) differentially expressed (DE) at least once over the 4 h time course p.i., compared to the in vitro grown control ( Table S2). As opposed to the 1 h-peak of the host gene expression after infection (Rossetti et al., 2013), the total number of perturbed Brucella genes is rather constant in the first 4 h p.i. (15 min: 1,899 genes, 30 min: 1,937 genes, 1 h: 1,968 genes, 2 h: 1,909 genes and 4 h: 1,912 genes; Figure 1). The combined analysis of these results, i.e., the host and pathogen transcriptional profiles during bovine Peyer patch infection, clearly indicate that both host and Brucella gene expression responses are markedly perturbed at a very early time post-interaction. This is in concordance with other results (He et al., 2006;Rossetti et al., 2011bRossetti et al., , 2012, which corroborates an initial transcriptionally perturbed period followed by a more quiescent one. A group of 1,740 genes (55% of B. melitensis genome) was markedly perturbed in the same direction in at least 4 of 5 time points (Table S3). These genes were considered as the core set of genes associated with the adaptive changes of B. melitensis during the early in vivo bovine Peyer's patch infection, and therefore important in understanding key events in the early modulation of host response. From this set of 1,740 Differentially Expressed (DE) genes, 925 (53%) were activated and 815 (47%) were repressed compared with the in vitro grown culture. Interestingly, genes from the core set located on chromosome I were mainly activated (over 1,174 DE genes: 752 were up-and 422 were down-regulated), while genes located on chromosome II were mainly repressed in higher numbers (566 total DE genes: 173 were up-and 393 were down-regulated). Chromosome I encodes the majority of the core metabolic machinery for transcription, translation and protein synthesis, and Chromosome II is overrepresented in genes involved in pathways for utilization of specific substrates (membrane FIGURE 1 | Graphic representation of B. melitensis genes differentially expressed (DE) throughout the experiment. Blue bars represent genes activated; light red bars represent genes repressed. For differential analysis, the in vivo infected loop gene expression is compared to the in vitro log growth phase inoculum as the control. transport, central intermediary and energy metabolism, and regulation; Paulsen et al., 2002). Altogether, these results suggest that Brucella may restrain metabolic functions while inducing transcriptomic modifications to adapt from an extracellular to an intracellular lifestyle. These results are largely in concordance with previous publications (Lamontagne et al., 2009;Rossetti et al., 2011a,b) even though our study analyzes the complexity of in vivo invasion process in comparison with pathogen gene expression in other in vitro (i.e., one cell line) models of invasion.
Microarray gene expression data were validated by qRT-PCR. Six randomly selected Brucella genes, determined to be significantly affected throughout the first 4 h p.i., were chosen for verification at every time point (i.e., 30 data points). As shown by the representative examples in Figure S1, gene expression changes were consistent between microarray and qRT-PCR for genes with increased expression or genes with decreased expression relative to the control.

Major Brucella Virulence Factors Are Down Regulated at the Onset of the Infection
Within 15-30 min of in vivo exposure, B. melitensis adhere and immediately penetrate through the intestinal mucosa and Peyer's patch and rapidly disseminate through systemic circulation (Rossetti et al., 2013). Early stage host gene expression of Syndecan 2, Integrin alpha L and Integrin beta 2 genes coincide with initial Brucella adhesion which is coupled with simultaneous repression of two intestinal barrier-related pathways (Tight Junction and Trefoil Factors Initiated Mucosal Healing), subverting mucosal epithelial barrier function and facilitating Brucella transepithelial migration (Rossetti et al., 2013). To elucidate Brucella virulence mechanisms responsible for this host molecular response, we expanded our analysis on pathogen pathways (Table 1) and GO alterations (Tables S4-S8). Note that pathway molecular interactions and annotations are based on those provided by the Kyoto Encyclopedia of Genes and Genomes (KEGG; Kanehisa et al., 2017). Pathways and Gene Ontology groups are comprised of gene sets which may be either activated or repressed in some combination over time. The Bayesian scoring method computes the log-likelihood of the in vivo expressed data and measures its goodness-of-fit to a model trained with control data (the in vitro inoculum expression data). In this manner, it is possible to determine if a pathway or GO group is activated or repressed. In our computational system biology approach, if the sum of the individual gene scores within a pathway/GO group is positive then the pathway/GO score is considered to be activated. Otherwise, if the sum is negative, the pathway/GO score is considered repressed and assigned a negative score value. Table 1 shows the results of the pathway analysis scoring listed by pathway categories and sorted by activated or repressed state on the 15 min p.i. column. Specific gene expression scores within these pathways are provided in Table S7. Early in the infection process, the pathway category "Environmental Information Processing" has several important pathways involved in B. melitensis pathogenicity which are   * Indicates that current KEGG Database only includes the 'map' reference pathway. ** Indicates the "bme" pathway was combined with another pathway in the current KEGG database. The pathways are organized by category and then sorted by activated or repressed states on the T15 minute p.i. column. Pathway scoring employed the in vivo gene expression for the experimental treatment condition and the in vitro gene expression from the log growth phase of the inoculum as the control condition. The "bme" pathway molecular interactions were originally based on the 2009 version of Kyoto Encyclopedia of Genes and Genomes (KEGG) database and their KEGG pathway name designators (Kanehisa et al., 2017). Some pathway names indicated with a * or ** required updating to be consistent with the current online KEGG database.
repressed at 15 min p.i. that include "Type IV secretion system, " "Type III secretion system, " Two-component system, " and ABC transporters. It is interesting to note that the Twocomponent regulatory systems (TCRSs) and Type III secretion system pathways reverse to an activated state at 30 min. p.i. In the cellular processes category, the "Flagellar assembly" and "Bacterial chemotaxis" pathways are also repressed across all time point p.i. The TCRSs are signal transduction mechanisms that allow microorganisms to sense and respond to changes in environmental conditions. Bioinformatic analysis of Brucella genomes has identified 15 predicted bona fide TCRS pairs (Lavin et al., 2010). Several of these systems have been characterized in Brucella species, such as BvrSR (Sola-Landa et al., 1998), FeuQP (Dorrell et al., 1998), NtrBC (Dorrell et al., 1999), NtrXY (Foulongne et al., 2000;Carrica et al., 2012), PrlSR (Mirabella et al., 2012), the flagellar master regulator FtcR (Leonard et al., 2007), the blue-light-activated LOV HKs (Swartz et al., 2007), RegA (Abdou et al., 2017), and CenR (Zhang et al., 2009). Other TCRSs have been identified by transpositional mutagenesis during global screening for virulence factors (Lestrate et al., 2000;Wu et al., 2006) but remain uncharacterized. Our pathway analysis revealed that the TCRSs were initially repressed at 15 min. p.i. and were then activated for the remaining time points (Table 1) providing evidence that in vivo Brucella sense and actively regulate their metabolism through the transition to intracellular lifestyle. Changes in expression between 15 min to 30 min p.i. by the TCRSs genes dctM, glnG, glnL, phoB, phoQ, citE, and divJ resulted in the TCRSs pathway transitioning from a repressed state to an activated state with the exception of aceA which was activated early and then repressed. These genes went from a strongly repressed state to an insignificant expressed state.
Contrarily, other Brucella pathways involved in virulence such as "ABC transporters" and "T4SS system" were continuously repressed suggesting a silencing strategy to avoid stimulation of the host's innate immune response very early in the infection process ( Table 1). The highly repressed genes associated with the ABC transporters repression included BMEII0196, BMDII0861, PBMII0120, BMEI1138, proW, ybbP, pstB, potB, afuB, rbsC and several others. The T4SS system repression was induced by the repressed genes virB4, virB5, virB6, and virB9 (Table S2). In vitro studies have demonstrated that T4SS is not required for cellular invasion, and its expression begins 15 min after phagocytosis and maximizes at 5 h p.i. (Sieira et al., 2004). It has been shown to be indispensable for intracellular survival of Brucella Sieira et al., 2000;Delrue et al., 2005). Under our in vivo experimental conditions, expression of genes from the virB operon was repressed as was confirmed by qRT-PCR [ Figure S1: BMEII0033 (virB9)]. In addition, the transcriptional regulator vjbR (BMEII1116) that positively regulates the expression of B. melitensis virB operon (Delrue et al., 2005) was not differentially expressed in our microarray results (Table S2). These data show that the T4SS was repressed during the first 4 h p.i. of in vivo infection. Collectively, our results, in addition to those reported earlier (Roux et al., 2007;den Hartigh et al., 2008) that failed to detect differences in the number of B. abortus and B. melitensis WT and virB mutant recovered from mice spleens in the first 3 days p.i., suggest that the virB operon may not play a major role in the initial in vivo Brucella pathogenesis. There are likely in vivo environmental signals that modulate the expression of the virB operon differently than reported in in vitro systems of infection. Identification of the host molecule targets of the T4SS will help characterize its expression based on the host cellular response discussed further in the next section "Systems biology in vivo interactome modeling results."

Systems Biology In vivo Interactome Modeling Results
The simultaneous collection of host and pathogen gene expression data of the bovine host ileal loop infected with B. melitensis WT (Rossetti et al., 2013), provided us with a unique opportunity to examine temporal host pathway perturbations concurrent with those of the pathogen. A computational approach based on DBN machine learning was employed to infer protein-protein interactions (PPIs) and to create a novel in silico host-pathogen interactome model (File S1). To identify plausible PPIs, a specialized application was developed to implement algorithms that integrate multiple sources of PBK such as from KEGG, BIOCARTA, NCBI, PIBASE, and Brucella proteomic analyses (Delvecchio et al., 2002;Wagner et al., 2002;Connolly et al., 2006;Mol et al., 2016), into the inference of host-pathogen protein interactions. Such interactions aid in the identification of mechanisms of host invasion and evasion through manipulation of the host's immune response system. Our application employed Bayesian networks (BNs) (Friedman et al., 2000;Hartemink et al., 2001) that were expanded by others to include PBK (Imoto et al., 2004;Werhli and Husmeier, 2007). We employed methods similar to Werhli for learning PPIs from expression data and PBK (Werhli and Husmeier, 2007). We adopted three algorithmic methods for the identification of candidate interaction points for use in network learning between the host and pathogen from in vivo gene expression data. Through either: (1) protein binding domain, (2) sequence similarity, or (3) Gene Ontology-based functional algorithms of pathogen gene expression with the host gene perturbations, we identified potential B. melitensis interactions with host pathways.
The PPI analysis resulted in identifying the virB gene that encodes the T4SS proteins to have plausible interaction with a number of genes in the host's immune response pathways ( Table 2 and Table S8). For example, the significantly perturbed virB11 gene has a high protein domain binding prediction with a negative correlation with the host gene/protein PIK3R2 expression. PIK3R2 is a key regulatory protein in several key pathways including: mTOR signaling, T-cell and B-cell receptor, Integrin-mediated cell adhesion, regulation of actin cytoskeleton, apoptosis, and Toll-like receptor. The host gene PIK3R2 remained activated for all time points p.i. Interestingly, in our host pathway analysis of "Regulation of Actin Cytoskeleton, " the genes ABI2, PPN1, and ARPC5, downstream from PIK3R2, were repressed suggesting a mechanism of host pathway disruption or highjacking. The VirB11 also had a strong protein domain binding prediction with negative correlation with the host genes MAPK8IP1/2 of the MAPK signaling pathway. The MAPK8IP1/2 host genes were strongly repressed 15 min p.i. and insignificantly expressed thereafter. Down-stream of MAPK8IP1/2, the transcription factor JUND and nuclear factor of activated T-cells, cytoplasmic 3, NFATC2 are both strongly repressed.
Brucella flagellum is a virulence factor transiently expressed during vegetative growth and required for persistent infection, but not for internalization in vivo (Fretin et al., 2005). In agreement, our results during the first 4 h of infection showed a repression of the three flagellum-encoded loci (BMEII0150-0168, BMEII1080-1089, and BMEII1105-1114 (Table S3), with corresponding repression of the "flagellar assembly" pathway ( Table 1), and the cell components "bacterial-type flagellum hook" and "bacterial-type flagellum" (Table S5) in tissueassociated B. melitensis compared with in vitro-grown cultures. In addition, RopE1 sigma factor (BMEI0371), a flagellar repressor (Ferooz et al., 2011), was activated throughout the experiment (Table S3). We further examined potential interaction of the flagella-associated genes with the host. Three highly correlated PPIs were identified which included the flagella genes BMEI0324, BMEII1085 (flgA), and BMEII1113 (fliG-2). The ORF BMEI0324 had strong binding sequence similarity and positive correlation to the host expressed JUN (jun oncogene) which is part of the highly perturbed host pathways: Toll-like Receptor, ErbB Signaling, BRC Signaling, B-cell Signaling, T-cell Signaling, Epithelial Cell Signaling, WNT Signaling, and MAPK Signaling. The flagellum gene flgA also had strong binding sequence similarity and positive gene expression correlation to host CASP2 gene while having a reversed (negative) correlation with the host activated CASP3 gene. Interestingly, the lowly expressed CASP2 is only associated with the highly perturbed MAPK signaling pathway, while CASP3 has several pathway associations that include: Apoptosis, Epithelial Signaling, Natural Killer Cell, and MAPK Signaling. The third flagellum gene fliG-2 had strong binding sequence similarity to the host MAP4K1 gene. The highly activated MAP4K1 is associated with only the host MAPK Signaling pathway. Such interactions may be novel virulence candidates that facilitate circumventing the host immune response. An example of interaction is the pathogen gene flgA with the host gene/protein Casp4/6/7/9 which are MAPK pathway genes. Accordingly, down-stream of Casp genes are the repressed RAC1/2/3 genes of the Rho family of GTPases. RAC1 has been implicated in various downstream cellular functions, including, but not limited to, cellular plasticity, migration and invasion, cellular adhesions, cell proliferation, and apoptosis.
Host cells identify specific pathogen-associated molecular pattern (PAMP) motifs present in the bacteria by patternrecognition receptors (PRRs), such as Toll-like Receptors (TLRs). These receptors are key to establishing an important network between the innate and adaptive immune systems. TLR5 is the cellular receptor for extracellular flagellin, a major structural protein of Gram-negative flagella. Binding of flagellin to the extracellular domain of TLR5 rapidly induces a signal cascade that culminates in the production of proinflammatory mediators such as cytokines, chemokines, and costimulatory molecules (Honko and Mizel, 2005). Therefore, the absence of flagellum apparatus during extracellular life while inside the host suggests the Brucella strategy is to avoid triggering a host immune response and an initiation of a Brucella persistence mechanism . However, our previous analysis showed that TLR5 pathway is activated in B. melitensis-infected bovine Peyer's patches during the first hour p.i. (Rossetti et al., 2013), which may have been associated with remnants of flagella in the in vitro-growth culture media intraluminally inoculated. More detailed analysis of this pathway showed that down-stream of TLR5 there were several strongly repressed genes including PIK3C2B, PIK3R4, STAT1, AKT3, RAC3, IL6, and TICAM3. This may suggest that the pathogen is manipulating important signaling processes by some other mechanism. PPI analysis indicated that virB genes have predicted interactions with STAT1, PIK3C2B, and IL6 which may also be circumventing the TLR5 response to flagellin stimulation and preventing the host from mounting an effective immune response. The complexity of a complete system-level host and B. melitensis interaction model (G interactome ) is illustrated in Figure 2A. This Bayesian network model is comprised of approximately 528 nodes (genes) and 987 arcs that connect gene nodes (relationships). Of the 987 arcs, 101 arcs were learned for host-pathogen points of interaction which are highlighted by the orange colored arcs. The number of host genes were limited to a selected set of perturbed host pathways which included MAPK signaling, ErbB signaling, mTOR signaling, WNT signaling, VEGF signaling, Toll-like Receptor signaling, GnRH signaling, Tight junction, Phosphatidylinositol signaling, Notch signaling, Natural killer cell mediated cytotoxicity, and Apoptosis. The intent for using only perturbed pathways was to look for plausible points of pathogen interactions which could influence the hosts immune response. Although this model is visually complex, the model allows for the computational extraction of potentially important mechanisms of interaction. Of the 101 arcs, the prioritization of these interactions can be analyzed according to the most likely to least likely in the following order: "protein domain, " "sequence similarity, " "GO Functionality, " and "Correlated Data". The creation of the interactome model employed only PPI relationships based on protein domain or sequence similarity. For example, Figure 2B illustrates the simplification for the interaction between the pathogen's Type IV secretion system and a known cell surface protein, BtaE (Ruiz-Ranwez et al., 2013b), with the host's Toll-like receptor pathway. This demonstrates how predictive information can be employed to interpret host-pathogen responses. The thicker orange arcs are the connections from the pathogen to the host. From this type of analysis, it is possible to understand the state of host's gene expression down-stream from the potential points of interaction/disruption in any of the pathways showing potential manipulation by the pathogen. Table 2 lists a selected subset of predicted interactions representing the pathogen-host pairs based solely on "protein domain" and "sequence similarity" prediction. A complete listing of predicted PPIs based on protein domain or sequence similarity is provided in Table S8. Note that our computational approach did predict interactions of the B. suis gene BtaE with several host proteins listed in Table 2. BtaE belongs to the type II (trimeric) autotransporter family and is an orthologue of B. melitensis BMEI1873. BtaE has been shown to have an active role in host cell adhesion and binding with components of the extracellular matrix such as fibronectin, collagen, and vitronectin (Ruiz-Ranwez et al., 2013b). The interactome model and its predicted PPI list is the analysis output to be employed for further in vitro validation and model refinements. The resulting B. melitensis PPI gene set may represent important new virulence factors with the potential to disrupt or hijack the host immune response. Table 2 also lists the perturbed host pathways in which the host gene PPI is associated, intentionally unfiltered conceptually for subcellular locations so that all PPI are presented. Little is known about the complete secretome of B. melitensis during in vivo host invasion and proliferation although the secretion systems and secretomes of Brucella were recently computationally analyzed which predicted 29 host-pathogen specific interactions between cattle and B. abortus and 36 host-pathogen interactions between sheep and B. melitensis proteins (Sankarasubramanian et al., 2016). The PPI computational approach employed evidence of host pathway perturbation and gene expression disruption as possible indicators of pathogen interaction. If there was plausible potential for a PPI based on binding domain or sequence similarity to known protein interactions between the pathogen and host protein, then it was included in the PPI in list ( Table 2 and Table S8). The list of PPIs can be prioritized based on the normalized correlation weights. The larger the normalized weight indicates stronger likelihood of a relationship between the pathogen and host genes. Note that the normalized correlation weight is an output of structure learning and is employed in the acceptance or rejection of an arc (edge) in the final Bayesian network. Arc weight is dependent on the number of incoming arcs to a node and other factors and should not be confused as a true correlation measurement between two gene expression values.
It is thought that the unfiltered PPI predictions could, in future experiments, employ such information as normalized correlation weights and cellular localization to help prioritize the selection of which PPI to be experimentally validated. Further, it is proposed that the evidence driven computational approach (in vivo host-pathogen responses and machine learning) for predicting bacterial and host cell protein interactions will narrow the focus on likely PPI candidates and will greatly enhance our capacity to design hypothesis-driven experimental approaches to discover which Brucella proteins directly participate in host interactions.
In conclusion, the in silico interactome modeling offers informative insights leading toward new hypotheses regarding host-pathogen mechanisms of invasion and evasion. This modeling infers that B. melitensis has multiple points of host interaction that occur at the early stage post infection. A number of important innate immune response pathways appear to be potential targets of disruption by invading B. melitensis, such as Regulation of Actin Cytoskeleton, mTOR Signaling, MAPK, and Toll-like Receptor Signaling appear to be likely targets of pathogen manipulation that warrant further exploratory research and verification of PPIs. As we have reported and discussed here, identifying interactive host:pathogen PPIs is often the initial step to establish functional significance according to the principle of "guilty by association" (Schauer and Stingl, 2009) that may drive future research to a higher level of understanding of the molecular pathogenesis of brucellosis, thereby facilitating the design of novel immunotherapeutic drugs and vaccines.

ETHICS STATEMENT
This study was carried out in accordance with the recommendations of the Texas A&M University Institutional Animal Care and Research Advisory Committee. The protocol (AUP#2003-178) was approved by the Texas A&M University Institutional Animal Care and Research Advisory Committee.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fmicb. 2017.01275/full#supplementary-material Figure S1 | Validation of Brucella melitensis microarray results by quantitative real time PCR. Six randomly selected B. melitensis ORFs that were consistently perturbed in microarray results in the first 4 h p.i. as compared to the inoculum as validated by quantitative RT-PCR. Fold-change was normalized to the expression of B. melitensis 16s rRNA and calculated using the ∆∆C t method. All tested genes at every time point had expression altered in the same direction as microarray. Open bars represent fold-change by microarray analysis and black bars represent fold-change by qRT-PCR.
Table S1 | Primers for Real Time PCR analysis of genes in B. melitensis samples. Table S2 | Bayesian z-score of B. melitensis genes in tissue-associated B. melitensis from 15 min to 4 h post-infection. Genes with z-score >|2.24| were considered differentially expressed. Positive numbers in the body of the table indicate activated genes and negative (−) numbers indicate repressed genes. The tissue-associated B. melitensis gene expression at every time point was compared to the gene expression of the inoculum (i.e., in vitro-grown cultures of B. melitensis at late-log phase of growth). Measured time points were 15 (T15), 30 (T30), 60 (T60), 120 (T120), and 240 (T240) min p.i. Table S3 | Core set of B. melitensis genes differentially expressed (>|2.24|) in at least four of five time points in the first 4 h post-infection of bovine Peyer's patch. Genes with z-score >|2.24| were considered differentially expressed. Positive numbers in the body of the table indicate activated genes, negative (−) numbers indicate repressed genes. The tissue-associated B. melitensis gene expression at every time point was compared to the gene expression of the inoculum (i.e. in vitro-grown cultures of B. melitensis at late-log phase of growth). Measured time points were 15 (T15), 30 (T30), 60 (T60), 120 (T120), and 240 (T240) min p.i. Table S4 | Significantly perturbed Biological Processes (BP) (Bayesian z-score >|2.24|) of tissue-associated B. melitensis during the first 4 h post-bovine Peyer's patch infection. Measured time points were 15 (T15), 30 (T30), 60 (T60), 120 (T120), and 240 (T240) min p.i. The Bayesian z-scores represent the degree of perturbation of the group of BP genes in tissue-associated B. melitensis versus the control. Positive z-scores represent activation of the BP (BP score is dominated by more activated genes), while the negative (−) z-score represents BP repression (BP score is dominated repressed genes). Table S5 | Significantly perturbed Cellular Components (CC) (Bayesian z-score >|2.24|) of tissue-associated B. melitensis during the first 4 h post-bovine Peyer's patch infection. Measured time points were 15 (T15), 30 (T30), 60 (T60), 120 (T120), and 240 (T240) minp.i. The Bayesian z-scores represent the degree of perturbation of the group of CC genes in tissue-associated B. melitensis versus the control. Positive z-scores represent activation of the CC (CC score is dominated by activated regulated genes), while the negative (−) z-score represents CC repression (CC score is dominated by repressed genes). Table S6 | Significantly perturbed Molecular Functions (MF) (Bayesian z-score >|2.24|) of tissue-associated B. melitensis during the first 4 h post-bovine Peyer's patch infection. Measured time points were 15 (T15), 30 (T30), 60 (T60), 120 (T120), and 240 (T240) min p.i. The Bayesian z-scores represent the degree of perturbation of the group of MF genes in tissue-associated B. melitensis versus the control. Positive z-scores represent activation of the MF (MF score is dominated by more activated genes), while the negative (−) -score represents MF repression (MF score is dominated by repressed genes). Table S7 | Dynamic Bayesian pathway analysis scores along with the associated individual genes and their Bayesian Scores by time point post-inoculation of tissue-associated B. melitensis during the first 4 h post-bovine Peyer's patch infection. Measured time points were 15 (T15), 30 (T30), 60 (T60), 120 (T120), and 240 (T240) min p.i. Positive z-scores represent activation of genes while negative (−) -scores represent repressed genes). Table S8 | Comprehensive list of predicted host:pathogen tissue-associated B. melitensis:bovine protein-protein interactions (PPI). This list includes only those PPIs which were learned by either known protein domain binding or sequence similarity to known binding domains employing our Bayesian methods for PPI prediction.