Tissue Tropism in Host Transcriptional Response to Members of the Bovine Respiratory Disease Complex

Bovine respiratory disease (BRD) is the most common infectious disease of beef and dairy cattle and is characterized by a complex infectious etiology that includes a variety of viral and bacterial pathogens. We examined the global changes in mRNA abundance in healthy lung and lung lesions and in the lymphoid tissues bronchial lymph node, retropharyngeal lymph node, nasopharyngeal lymph node and pharyngeal tonsil collected at the peak of clinical disease from beef cattle experimentally challenged with either bovine respiratory syncytial virus, infectious bovine rhinotracheitis, bovine viral diarrhea virus, Mannheimia haemolytica or Mycoplasma bovis. We identified signatures of tissue-specific transcriptional responses indicative of tropism in the coordination of host’s immune tissue responses to infection by viral or bacterial infections. Furthermore, our study shows that this tissue tropism in host transcriptional response to BRD pathogens results in the activation of different networks of response genes. The differential crosstalk among genes expressed in lymphoid tissues was predicted to be orchestrated by specific immune genes that act as ‘key players’ within expression networks. The results of this study serve as a basis for the development of innovative therapeutic strategies and for the selection of cattle with enhanced resistance to BRD.


Results
Tissue gene expression profiles. The RNA-Seq experiments generated, on average, ~50 million reads per sample. The overall read alignment rate ranged from 84.8% to 94.1% and concordant pair mapping rate ranged from 74.2% to 93.1% (Supplementary Table 1). Transcript abundance differences between challenged and control animals were estimated for samples derived from lung lesion (LNGL) and four lymphoid tissues: bronchial lymph node (BLN), retropharyngeal lymph node (RLN), nasopharyngeal lymph node (NLN) and pharyngeal tonsil (PGT). For the lung, we also determined transcript abundance differences between LNGL and a nearby sample devoid of lesions collected from the lung (LNGH) of the same animal, with all tissue samples being collected at the peak of clinical signs. The gene expression changes in each tissue due to experimental challenge by individual pathogens (Bovine Respiratory Syncytial Virus, BRSV; Bovine Viral Diarrhea Virus, BVDV; Bovine herpesvirus 1, BoHV-1; Mannheimia haemolytica, MANNHE; or Mycoplasma bovis, MYCO) revealed variation in transcriptional profiles that were clustered into only 3 groups by both hierarchical clustering and principal component analysis (Fig. 1). The first two principal components captured 22.7% and 13.9% of the variance, respectively, in gene expression across all of the samples. The expression changes in BLN in response to challenge by the different BRDC pathogens clustered into a group that was distinct from other groups. Similarly, the gene expression changes identified in lung lesion samples relative to either unchallenged control animal lung samples or apparently healthy lung samples from the same animal (and therefore exposed to the same challenge pathogen) Differentially expressed genes in response to BRDC pathogens. The RNA-Seq data analysis revealed differential gene expression in different tissues between control and challenged animals, as well as between healthy lung and lung lesion samples from the same animals. The expression values along with significance (p-value and q-value) for the tests of differential expression are provided in the Supplementary Dataset. Different numbers of differentially expressed (DE) genes were identified for each of the tissue × pathogen combinations with fewer genes DE in the animals challenged with the bacterial pathogens than those challenged with viruses (p < 0.00001; Table 1). Across pathogens, the RLN had, on average, the fewest and BLN had the greatest number of DE genes. However, variation in the number of DE genes across pathogens was higher for RLN and NLN compared to the other tissues suggesting a greater specificity of pathogen response for these tissues. The number of DE genes differed across the tissue × pathogen combinations (p < 0.00001) with far more than the expected numbers of genes DE in BLN and LNGL in response to M. haemolytica, NLN and PGT in response to BVDV, and RLN in response to BRSV (Table 1). Table 2 shows the number of genes found to be DE in all 31 combinations of tissues ranging from a single tissue to all 5 analyzed tissues. Among the multi-tissue response genes, there were many more than the expected number of genes DE in both PGT and RLN in response to BRSV, in BLN, NLN and RLN in response to BoHV-1, in NLN and PGT in response to BVDV, in BLN and LNGL in response to M. haemolytica, and LNGL and NLN in response to M. bovis. The proportion of genes that were DE in two or more tissues was similar for all three viral pathogens (0.49-0.52) and both bacterial challenges (0.20-0.25). However, the proportion for bacterial challenges combined is less than half that of the combined viral challenges (2 × 2 contingency table; p < 0.00001), suggesting that gene expression in response to bacterial infections is largely tissue-specific.
The complexity of the relationships between DE genes among tissues is shown in Fig. 2B in a chord diagram. This map provides a graphical representation of the mutual information (a non-linear measure of dependency) of expression changes among genes across all tissue × pathogen combinations. Thus, Fig. 2B represents the extent of relationships between host transcriptional responses in different tissues that are due to the different pathogens. Among the DE genes identified in LNGL relative to LNGH, ~85% (2742 of 3211) were also found to be DE in the comparison of LNGL to the uninfected lungs of the control animals (Fig. 3A). In comparing lung lesion to healthy lung, we observed that the M. haemolytica challenge elicited a host transcriptome response that was the least related to the host transcriptional responses to the other pathogens (Fig. 3B). This may be due to the relatively larger changes in gene expression that occurred in response to M. haemolytica than for the other pathogens. Furthermore, the correlations between gene expression changes between pairs of tissues within each pathogen challenge group again revealed the tissue specificity of immune response to the bacterial challenges relative to the viral challenges (Fig. 4). In particular, the sign of the correlation between gene expression responses between PGT and NLN appears to discriminate between bacterial and viral infections. This is consistent with the fact that the immune response in these two tissues failed to distinctly cluster across the challenge pathogens (Fig. 1).
The Ingenuity Pathway Analysis software indicated that the majority of the DE genes were involved in pathways related to antimicrobial response, mostly innate, but also adaptive. In the challenged animals we generally found the up-regulation of pathways for acute phase signaling, complement system, regulation of cytokine production, interleukin and interferon signaling, granulocyte and agranulocyte adhesion and diapedesis, as well as the predominant down-regulation of a few lipid and cholesterol metabolism-related pathways such as peroxisome proliferator-activated receptor (PPAR) signaling, liver X receptor (LXR)/retinoid X receptor (RXR) activation and farnesoid X receptor (FXR)/RXR activation and antioxidant action of vitamin C. Likely downstream effects on cellular and organismal biology, such as the stimulation of lymphocyte activating factor IL1B, were also predicted from the expression data. We regularly predicted immunological, inflammatory and respiratory diseases as well as inflammatory responses. These findings suggest that we captured transcriptional variation within these tissues that were indicative of the organismal changes induced by infection. Genes responsible for nonspecific defense mechanisms against all respiratory disease pathogens 11 such as those encoding mucins, pattern recognition receptors (PPRs), host defense peptides (such as defensins, lactotransferrin and secretory leukoprotease inhibitor), and matrix metallopeptidase family members were consistently found to be DE across all challenge groups and tissues as were genes with reactive oxygen and wound healing (coagulation factors, THBD and VWF) functions. Among the pathways most enriched for DE genes across all pathogens, we consistently observed the activation of acute phase response signaling. While for most pathogens Oncostatin M (OSM) and IL1 appeared to regulate the expression of acute phase proteins, TNFα was induced by the M. haemolytica challenge. We also observed an enrichment of DE genes in pathways related to IL8 signaling, leukotriene, zymosterol, estrogen and cholesterol biosynthesis, eicosanoid signaling, hypercytokinemia and inhibition of matrix metalloproteases but the status of stimulation (up-or down-) could not always be predicted.
We hypothesized that in addition to a general immune response, pathogen-specific immune responses might be elucidated by sequencing the global RNA profiles of tissues collected at the peak of clinical disease 12 . We found many DE genes that were exclusive to each challenge group that were related to specific immune responses. While genes that were exclusively found to be DE between BRSV challenged and control animals, regardless of tissue, were involved in oxidative phosphorylation, mitochondrial dysfunction and hypoxia, the genes found to be exclusively DE in response to BoHV-1 infection were primarily involved in pathways involved in geranylgeranyl diphosphate biosynthesis via mevalonate, glutamate receptor signaling, mevalonate pathway, serotonin receptor and TGFβ signaling. The genes uniquely DE in response to BVDV challenge, were primarily related to the Ephrin receptors, which may function as entry receptors for BVDV. Pathways related to inhibition of viral replication such as Eukaryotic Initiation Factor 2 (EIF2) and Cell Cycle: G1/S Checkpoint Regulation, were also found to  Table 2. Numbers of differentially expressed genes in common between tissues or specific to a single tissue.    Signatures of tissue-specific gene expression. We next sought to understand how individual tissues respond to BRDC pathogens and whether they were associated with specific patterns of expression changes in response to challenge. To accomplish this, we first identified genes that were DE in a tissue-specific manner or that were ubiquitously DE in all five tissues (Table 2). Principal component analyses of the expression of genes DE in all tissues (marked as ' ALL') or specific tissues (marked using tissue abbreviations) are in Fig. 5A for each pathogen. The genes that were DE in all tissues tended to cluster differently for the different pathogens. We also conducted a principal component analysis of expression of DE genes between lung lesions and healthy lung tissue for each pathogen (Fig. 5B). These plots together show that gene expression patterns tend to cluster differently by host tissues and challenge pathogens. For example, infection by M. haemolytica induces more diverse transcriptional changes in lung lesion relative to healthy lung than any of the other pathogens (Fig. 6). In Fig. 7A, we show the relative proportions of genes that were DE between healthy lung and lung lesion for different combinations of BRDC pathogens. In particular, there were 1,256 genes DE in response to M. haemolytica but not in response to any of the other pathogens, 835 genes that were DE only in response to BVDV and 299 genes were DE in response only to BRSV. As a proportion of the total number of DE genes for each pathogen, the proportion for M. haemolytica (69.4%) was greater than that for BVDV (64.4%, p < 0.003) or for BRSV (39.7%, p < 0.00001). This suggests that the host immunological response in lung was predominantly associated with gene expression in response to the M. haemolytica challenge. Moreover, following M. haemolytica challenge, the expression levels for the DE and non-DE genes shown in the principal component analysis differed more for M. haemolytica than for all of the other BRDC pathogens (Fig. 7B-F), a result that was also supported by the clustering shown in Fig. 6. These data further support that the host transcriptional response in lung lesion differs substantially for M. haemolytica relative to the other BRDC pathogens.
Tissue tropism is associated with differential gene networking. While individual tissues appear to possess different roles in mounting a host response to infection, an interaction between tissue transcriptional responses may be necessary to mount an appropriate immune response to a specific pathogen. Accordingly, we hypothesized that the correlation between gene expression changes between tissues could be used to predict a tissue cooperative host response. As the number of DE genes varied considerably between the challenge pathogens (see Table 1), we randomly sampled 1,000 genes from the set of all DE genes across all tissues for each pathogen and generated mutual information matrices for each pathogen. From the mutual information analysis, we observed a strong deviation in mean mutual information for the viral challenges (BRSV = 0.2344, BVDV = 0.1238 Networks and gene ontology for response genes. We annotated the functions of the DE genes for which the host response was tissue-specific or tissue-agnostic. For each pathogen, these two gene sets were analyzed for the over-representation of gene ontology (GO) terms using DAVID software. The results indicated that genes that were DE in multiple tissues were more likely to play roles in immune response and host defense mechanisms than were the genes that were DE in a tissue-specific manner ( Table 3). Genes that were DE in all five tissues had functions in TLR, complement and coagulation cascades, endocytosis, chemokine and cytokine signaling, leukocyte transendothelial migration, cell adhesion and MAPK signaling. These events represent an orchestrated immune defense that occurs in all tissues to combat the infection. On the other hand, genes that were DE in a tissue-specific manner were more likely to be associated with GO terms such as extracellular region, contractile fiber part, myofibril, myosin filament and others. These genes included specific complement receptors regulating proinflammatory immune responses that are likely related to microenvironment immune responses and appear to define the nature of tissue tropism. For example, CD70, a member of Tumor Necrosis Factor ligand superfamily, was exclusively DE in the comparison of LNGL to control lung, as were genes encoding mitogen-activated protein kinases such as MAP2K6, MAP2K3 and MAPK8. Moreover, surfactant genes encoding proteins secreted by type II alveolar macrophages that induce immune responses 13 were also exclusively differentially regulated between LNGL and control lung.
The comparison of GO terms for the tissue-specific or tissue-agnostic DE genes in response to bacterial or viral challenges is shown in Supplementary Figure 1. This figure shows that tissue-specific or tissue-agnostic host transcriptional response GO terms tend not to overlap between the bacterial and viral challenges and when they were in common, tended to be for the bacterial pathogens.

Identification of key immune genes.
We identified immune function-related genes by querying the DE genes to the innate immunity database InateDB (http://www.innatedb.com) for Bos taurus (Supplementary Table 2). The immune function related genes that were DE between lung lesion and apparently healthy lesion-free lung tissues are listed in Supplementary Table 3. The Maximum Relevance Minimum Redundancy network analysis of expression changes for these immune function-related genes revealed an extensive interaction between tissues in response to the challenges (e.g., for BoHV-1 shown in Supplementary Figure 2). Using degree centrality statistics to predict the key players within the gene expression networks, we predicted the top three key players among the immune function gene networks in response to each of the challenge pathogens ( Table 4). The predicted key players have major roles in the defense against bacterial and viral infections. For example, BPIFA1 is expressed in the upper airways and nasopharyngeal regions in human and encodes an antimicrobial protein with antibacterial activity 14 . Several of the predicted key players are members of the C1Q 'complement' system gene family and are involved in host-pathogen interactions including respiratory tract inflammations 15 . APCS encodes amyloid P component, serum that is associated with laryngeal amyloidosis in humans 16 , and AKIRIN2 has been described as a 'novel player' in the transcriptional control of innate immunity 17 . Pathway prediction of differentially expressed immunity genes. Finally, we sought to establish which pathways might be involved in host transcriptional responses to BRDC pathogens. The genes that were found to be responsive across tissues were tested to determine if any Bos taurus KEGG 18 pathway was significantly over-represented by DE genes. We found that several disease-associated pathways, including those involved in host response to viral infections, were significantly associated with the genes that were consistently DE across      (Table 5). We also observed the RIG-I-like receptor, NF-kappa B and NOD-like receptor signaling pathways to be significantly enriched for DE genes suggesting that these pathways play roles in a tissue cooperative host response to infection. Phagosome and complement and coagulation cascades were also significantly enriched for DE genes further suggesting their roles in response to early infection events that may be indicative of infection spreading from one tissue to another. In lung tissues, we analyzed lesion-associated DE genes for KEGG pathway enrichment, and found that pathways such as the ECM-receptor interaction, focal adhesion, PI3K-Akt signaling pathway, along with other specific metabolic pathways were significantly enriched ( Table 6). The greatest number of DE genes was associated with the PI3K-Akt pathway, and this pathway is upstream of important molecular cascades known to be associated with cell cycle, programmed cell death and p53 signaling (Fig. 9). In Fig. 9, we observe specific DE genes (shown in red text) such as PEPCK, CCND1, FASLG and MYB that are directly upstream of these events (cell cycle, apoptosis and p53 signaling). This suggests that these DE genes within the PI3K-Akt pathway may play important roles in eliciting downstream cascades in lung lesions.

Discussion
We performed a detailed analysis of tissue transcriptional responses to understand the involvement of lymphoid and lung tissues in the normal immune response to infection by BRDC pathogens at the peak of clinical signs of disease. Relative to spontaneous infections, the intranasal inoculations for viruses and intra-tracheal inoculations for bacteria may have affected pathogen spreading and the localization and extent of the lesions as well as potentially the extent of lymphoid tissue involvement. While we cannot exclude this possibility, the respiratory viruses are known to be infective by the inhalation route and we considered the tracheal administration of the bacterial pathogens to best represent the normal mode of lung infection and using these routes of inoculation we successfully induced clinical signs and lung pathology in the beef steers following each of the single pathogen experimental challenges 9 . Consequently, we also expect that tissues analyzed produced appropriate genetic responses to each of the pathogens.  Table 4. Predicted top 3 key players in the immune function gene networks in response to challenge by BRDC pathogens. Prediction was based on degree centrality estimation within each of the mutual information networks of differentially expressed immune function genes. The degree centrality score is a value that represents how well the model predicts the three genes to occupy the central nodes in the network.  Respiratory epithelial cells can sense viral pathogens through diverse molecules such as TLR3, TLR7, TLR9, RIGI, MDA5 and NLRP3, leading to the activation of type I interferon and the induction of a range of antiviral mechanisms 11,19,20 . While we detected the up-regulation of IFNγ, we detected no type I interferon response in the animals challenged with BVDV. A recent study has shown that BVDV infection induces IFNγ secretion during acute phase signaling and that lymphoid tissues serve both as possible sources of IFNγ and as target organs for its effects 21 . Production of nitric oxide and reactive oxygen species were consistently predicted to be activated in the tissues of the BRSV challenged animals. The up-regulation of nitritic oxide and reactive oxygen species either reflects the induction of these species by the virus to facilitate replication 22 , or are hallmarks of oxidative stress and the induction of apoptosis 23 . We also detected the differential regulation of cytokine production by IL17A and IL17F which down-regulate IL10, IL6 and CCL3 in LNGL. Furthermore, we detected the suppression of leukocyte extravasation signaling in PGT and the up-regulation of thrombin signaling in all of the tissues. Non-cytopathic BVDV strains such as strain 890 used for inoculation in this study are known to suppress proinflammatory cytokines and co-stimulatory molecules 24 .  Table 6. Pathways associated with genes differentially expressed between lesion and healthy lung tissue. The bovine tracheal antimicrobial peptide was found to be DE in all tissues except NLN. While this gene was predominantly up-regulated, it was found to be down-regulated in NLN for animals challenged with BoHV-1 and M. haemolytica. However, pathogen produced toxins are known to have the ability to repress host defensin transcription as a mechanism to combat the host immune response 25 . Consistent with the role of reactive oxygen species in the first line of defense against pathogens, genes with functions involving reactive oxygen species were also found to be DE. Reactive oxygen species produced by phagocytes may help in the removal of pathogens that are resistant to antimicrobial peptides; however, potentially at the expense of tissue damage. We appear to have captured most of the events that caused lung tissue damage. The leukotriene biosynthesis, γ-glutamyl cycle and eicosanoid signaling pathways were inferred as being the most strongly induced pathways in the lung lesions of the BVDV challenged animals. Genes encoding γ-glutamyl transpeptidases, which are present in these pathways, are known to cleave leukotriene 26 and examination of the genes induced in this pathway revealed that they were involved in converting leukotriene-D4 to leukotriene-E4. Leukotrienes C4, D4 and E4 are released from lung tissue exposed to allergens and are involved in immediate hypersensitivity reactions 27 . The leukotriene biosynthesis pathway was also induced by M. haemolytica infection, however, ALOX5 was identified as the up-regulated gene responsible for the conversion of leukotriene A4 into other leukotrienes in the tissues of animals from this challenge group.

KEGG
We previously analyzed pathogen-specific host responses in the bronchial lymph nodes of these steers 9 . Consistent with our current findings, the common pathways underlying response to several of the challenge organisms included pathways with major roles in innate immunity that are not tissue or pathogen-specific. However, we also identified several pathways that were enriched for DE genes in all analyzed tissues including BLN, such as Toll-like receptor signaling, chemokine signaling and granulocyte adhesion and diapedesis, as well as conserved predicted upstream regulators including IL1B thus supporting their roles in the immune response to the pathogens of the BRDC. Despite the finding of common pathways enriched for DE genes across challenge pathogens, different genes appear to be responsible for pathway activation. Whether this result reflects tissue-specific activation of pathways or a lack of power of single tissue analyses with four biological replicates is not clear. The majority of the DE genes identified in the current study were associated with functions in complement and coagulation cascades, endocytosis, chemokine and cytokine signaling, leukocyte transendothelial migration, cell adhesion and MAPK signaling. These events represent an orchestrated immune defense that occurs in all tissues to combat the infection. We observed numerous genes exclusively DE in the PGT that operate within pathways related to the ribosome and spliceosome as well as the post-transcriptional modification of RNA. In particular, we detected genes related to EIF2 signaling and the regulation of the eIF4 and p70S6K signaling pathways. Similarly, the genes that were exclusively DE in the RLN were related to transcriptional regulation pathways. While these responses to the viral infections were also observed in other tissues, an increase in the number of DE genes related to signaling functions was observed in both the PGT and RLN.
Besides identifying genes with expression changes in response to different BRDC pathogens, our study shows that the host response to respiratory diseases is modulated by the nature of cross-talk between response genes in different lymphoid tissues. The network of response genes is modulated by key immune function related genes that we predicted using 'key player' analysis. While key player analysis has previously been applied primarily to the study of social networks, our implementation in the analysis of gene expression networks provided a useful approach for characterizing the nature of immune function genes in BRD. Our results further show that tissue gene expression profiles in response to infection by BRDC pathogens can generally be clustered based on tissues rather than pathogens (Fig. 1) suggesting that different lymphoid tissues play signature roles in mounting the host's response to respiratory infections. This is necessary to trigger appropriate host immunity to defend against pathogen infections where clinical manifestation varies depending on strain of virus, host, immunity and other factors 28,29 . We show in this work that tissue cooperative expression changes in genes differ significantly between viral and bacterial infections. This suggests that tissue tropism of the host's transcriptional response is dependent upon the invading pathogen and may be dynamic to counter the evolution of pathogen virulence. Host-pathogen interactions are dependent upon the binding receptors of appropriate cells in the target tissues and tropism in host tissue gene expression may be a key determinant of how well a pathogen can successfully spread its infection from one tissue to another.

Materials and Methods
Animal ethics statement. This study was conducted under animal use protocol #16424 approved by the Institutional Animal Care and Use Committee of the University of California at Davis. The protocol adheres to the Federation of Animal Science Societies Guide for the Care and Use of Agricultural Animals in Research and Teaching.
Animal sampling and challenge. Crossbred steers were generated at the University of California Davis Sierra Field Station located in Brown's Valley, California by mating Angus bulls to advanced generation Angus-Hereford crossbred dams. Blood from steers that had not been vaccinated against any pathogen from the BRDC was tested for antibodies against each bacterial and viral pathogen and those found to be seronegative, or to have the lowest titers against each of the challenge pathogens, were selected at 6-8 months of age for the challenge experiment that was conducted at the University of California Davis. The steers were provided water ad libitum while maintained in pens and a fed a 65% concentrate starter diet. Pathogen challenges were performed sequentially on pen-grouped animals under strict biosecurity protocols. A period between challenges with different pathogens was used to prevent cross-infection 12 .
A pilot project was conducted during the summer of 2011 to determine the optimum pathogen doses to administer to the challenged animals and also the timing of clinical signs of peak infection to determine the euthanasia end-point for each challenge group. In the optimized challenge experiment conducted in summer of 2012, groups of N = 4 animals were individually challenged with BRSV, BVDV or BoHV-1 via a nebulizer or with M. haemolytica or M. bovis administered by an intratracheal tube inserted through the ventral nasal meatus until the end was approximately in the mid-trachea. Optimized pathogen doses per animal were: BRSV (1.6 × 10 5 / ml × 8.5 ml), BVDV (2.0 × 10 8 /ml × 8.5 ml), BoHV-1 (1.0 × 10 7 /ml × 8.5 ml), M. haemolytica (4.8 × 10 11 CFU) and M. bovis (7.0 × 10 10 CFU). Using a nebulizer, two of the control animals were aerosol administered with 8.5 ml of tissue culture media while a second pair of controls were inoculated with 8.5 ml of phosphate buffered saline. The intratracheal inoculation of phosphate buffered saline was followed by sufficient air to remove all saline from the tube. The aerosol administration of tissue culture media was designed to mimic one of the mechanisms of viral exposure. Strains were previously reported 12 .
Clinical scores were recorded daily on each animal from 3 days prior to challenge until the day of euthanasia using an index that was a combination of observations. None of the animals required analgesic or anesthetic treatment except topical lidocaine, which was applied to the nasal passages prior to the insertion of the intratracheal tube that was used to administer each challenge pathogen. Clinical signs peaked on days 7, 15, 6, 5, and 15 post-challenge for BRSV, BVDV, BoHV-1, M. haemolytica, and M. bovis, respectively, at which point the steers from each group were euthanized. These were not considered to be humane end-points. A certified veterinary pathologist removed tissues for examination and the samples collected for transcriptome analysis were immediately frozen in liquid nitrogen. Tissues remained frozen at −80 °C until RNA was isolated for analysis. One animal that was challenged with M. bovis was euthanized early and tissues were not collected. This left 23 animals from which samples of each tissue were harvested from control animals (N = 4) and those animals artificially RNA isolation and sequencing. RNA was extracted from a total of 50-100 mg of frozen tissue using RNeasy Plus Universal Mini Kits (Qiagen, Hilden, Germany). Briefly, 900 μl of QIAzol (similar to acidic buffered phenol) was added to the frozen tissue sample, which was immediately homogenized and incubated at room temperature for 5 min for the dissociation of the nucleoprotein complex, and then 100 μl of genomic DNA eliminator solution and 180 μl of chloroform were added. Samples were centrifuged at 12000 g for 15 min at 4 °C and after transferring the aqueous layer to a fresh tube, 600 μl of 100% ethanol was added and mixed by pipette. The mixture was transferred to an RNeasy spin column (based on silica-membrane technology) and centrifuged at 8000 g for 15 sec at room temperature and was then washed with 700 μl of Qiagen RWT buffer and 500 μl of Qiagen RPE buffer. The RNA was eluted with 100 μl RNase-free water and was placed at 4 °C. RNA purity and concentration were evaluated using a NanoDrop 1000 v1.3.2 (Thermo Fisher Scientific, Waltham, MA) and Qubit 3.0 Fluorometer (Thermo Fisher Scientific). The extent of RNA degradation was initially assessed by the electrophoresis of 1 μg of RNA on a 1.0% agarose gel. Finally, RNA quality was assessed for each sample using a Fragment Analyzer TM (Advanced Analytical Technologies, Inc, Ames, IA) with a standard sensitivity RNA analysis kit and we accepted samples with an RNA quality number of at least 8.0 for library construction.
A total of 10 μg of RNA from each tissue was processed using the TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA) to prepare samples for sequencing. Oligo dT magnetic beads were used to purify polyadenylated RNA from the total RNA which was then fragmented with divalent cations under elevated temperature. First strand cDNA was synthesized using random hexamer primers and then the second strand was synthesized. The double-stranded cDNA was end-repaired and the 3′ ends adenylated. Finally, to produce the sequencing library, universal adapters were ligated to the cDNA fragments that were then amplified by solid phase polymerase chain reaction. Each library was evaluated using a Fragment Analyzer and equimolar amounts were used to create 7 library pools, which were each sequenced (2 × 50 bp) on a single lane of a HiSeq. 2000. We produced an average of 50,778,115 reads per sample.
Processing of sequence reads. Adapter sequences were trimmed from the sequence reads as previously described 30 . Read alignment. Computation for this work was performed on the high performance computing infrastructure provided by Research Computing Support Services at the University of Missouri, Columbia MO. TopHat v2.0.6 31 was used to map the trimmed sequence reads to the NCBI Bos taurus virtual transcriptome build and the UMD3.1 reference assembly. TopHat first used Bowtie to align the sequence reads to the virtual transcriptome build. Reads that failed to map to the virtual transcriptome build were next mapped to the UMD3.1 reference assembly. Reads that mapped to UMD3.1 were converted to genomic mappings, spliced if required, and then merged with the transcriptome mapped reads. No more than 2 mismatches were allowed in the alignment of each read pair and the remaining TopHat parameters were left at default values.
Testing for differential expression. Cuffdiff 32 was used to estimate transcript abundance levels as Fragments Per Kilobase of exon per Million fragments mapped for each gene in each challenge group and to test transcripts for differential expression against controls. The p-values for each performed test were transformed to q-values using the Benjamini-Hochberg correction 33 to correct for multiple testing.
Cluster analysis, network inference and key player analysis. Expression data for all genes were subjected to hierarchical clustering among samples to generate a dendrogram. The Euclidean distance measure was used to calculate distance and data fitting was performed using the method of Ward 34 . The model-based cluster analysis was performed using the R 'mclust' package using a finite normal mixture modeling method 34 . The Bayesian Information Criterion, implemented in mclust, was chosen for model selection to cluster the genes based on their expression data Eigenvalue Decomposition Discriminant Analysis (EDDA) method 35 . To infer the extent of interactions between genes within a cluster, network analysis was performed using 'minet' 36 based on mutual information values that were calculated using the Spearman correlation estimator. The square symmetric matrix containing all mutual information values for all genes was used to generate a weighted adjacency matrix by the 'maximum relevance minimum redundancy' method 37 . The adjacency matrix was then used to identify and plot the co-expression networks. Key genes in each network were predicted from degree centrality scores (that predict how central the genes are relative to the overall network structure) using a 'key player' analysis approach, a method used in analyzing social networks 38 . Statistical tests. Chi-Square contingency table tests were performed using counts of numbers of DE genes for single and multi-tissue expression to examine associations between host tissue and BRDC pathogen. All of the principal component analyses and data visualization were performed using R. For the randomization test of DE genes in the network analysis, we selected a fixed number of genes (n = 1000) to create the expression networks. We repeated the process of sampling 1000 genes and calculating pairwise mutual information scores and used the average mutual information scores across replicates to perform Chi-Square tests to examine associations between pathogen and tissue combinations.
Annotation of differentially expressed genes. The Ingenuity Pathway Analysis software (http://www. ingenuity.com) was used to understand the biological processes regulated by the DE genes. The gene ontology and KEGG pathway analyses of DE genes were performed using DAVID Bioinformatics Resources 6.8 (https:// david.ncifcrf.gov). Innate immune genes of cattle curated at InnateDB (http://www.innatedb.com) were used to identify immune function genes found to be DE in the RNA-Seq data. Pathway painting was performed using KEGG mapper (http://www.genome.jp/kegg/mapper.html).