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⁹ CFU of B. melitensis 16 M/ml (total of 3 × 10⁹ CFU) at late-log growth phase cultured in F12K medium (ATCC®) supplemented with 10% heat-inactivated fetal bovine serum (HI-FBS) (ATCC®). 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. The mucosal layer of Peyer's patch was immediately dissected, macerated and homogenized in TRI-Reagent® (Ambion, Austin, TX). Subsequently, samples were appropriately contained and transported to an inspected and approved BSL-3 laboratory for immediate RNA extraction. Calves were euthanized with an intravenous bolus of sodium pentobarbital at the completion of the procedures. All animal experiments were approved by the Texas A&M University Institutional Animal Care and Research Advisory Committee (AUP#2003-178). Surgeries were performed under biosecurity laboratory III (BSL-3) conditions in CDC-approved isolation buildings at the Texas A&M University experimental farm (College Station, TX).

RNA isolation, labeling, and hybridization procedures were performed as described in previous experiments (Rossetti et al., 2010, 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®, Ambion) and then amplified from 30 μg of total RNA from B. melitensis-infected 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 vitro-grown 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. melitensis-infected 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 (Rossetti et al., 2009). Slides were hybridized at 45°C for approximately 20 h in a dark humid chamber (Corning, Corning, NY). Then, washed for 10 min at hybridization temperature with low stringency buffer [1X SSC, 0.2% SDS] followed by two 5-min washes with a higher stringency buffer [0.1X SSC, 0.2% SDS and 0.1X SSC] at room temperature with mild agitation, dried by centrifugation and immediately scanned.

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®) (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 sequence-similarity interaction transference procedure; (2) structural protein domain-based algorithm; and (3) a functional gene-ontology-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 post-infection (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 z-score 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.

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., 2011b, 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 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.

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 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 Two-component 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 (O'Callaghan et al., 1999; 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.”

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, down-stream 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 tissue-associated 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 pattern-recognition 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 (Terwagne et al., 2013). 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.

KD is a Chief Technology Officer in Seralogix, LLC. Seralogix is a bioinformatics research and services company commercializing computational systems biology software tools that are being sponsored by the National Institute of Allergy and Infectious Diseases and the National Human Genome Research Institute. KD participated in conducting certain genomic data processing involving pathway analyses and modeling that helped to provide a more system-level perspective of the host-pathogen interaction to the Texas A&M University researchers. Data were processed by KD utilizing Seralogix's proprietary computational pipeline for biological systems analysis. The relation between Seralogix and Texas A&M University, College of Veterinary Medicine and Biomedical Science is strictly on a collaborative (mutually beneficial) research basis with no financial arrangements, commitments or interests. KD's motivation is to see their computational tools produce results that contribute to the improved understanding of host response to pathogen invasions (an objective of his National Health Institute research grants). KD contributed to the interpretation of the analysis results provided to the Texas A&M University researchers. Seralogix has no ownership of the data, nor results produced by their tools. The other authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.