Bacterial and viral pathogen-associated molecular patterns induce divergent early transcriptomic landscapes in a bovine macrophage cell line

Background Pathogens stimulate immune functions of macrophages. Macrophages are a key sentinel cell regulating the response to pathogenic ligands and orchestrating the direction of the immune response. Our study aimed at investigating the early transcriptomic changes of bovine macrophages (Bomacs) in response to stimulation with CpG DNA or polyI:C, representing bacterial and viral ligands respectively, and performed transcriptomics by RNA sequencing (RNASeq). KEGG, GO and IPA analytical tools were used to reconstruct pathways, networks and to map out molecular and cellular functions of differentially expressed genes (DE) in stimulated cells. Results A one-way ANOVA analysis of RNASeq data revealed significant differences between the CpG DNA and polyI:C-stimulated Bomac. Of the 13,740 genes mapped to the bovine genome, 2245 had p-value ≤0.05, deemed as DE. At 6 h post stimulation of Bomac, poly(I:C) induced a very different transcriptomic profile from that induced by CpG DNA. Whereas, 347 genes were upregulated and 210 downregulated in response to CpG DNA, poly(I:C) upregulated 761 genes and downregulated 414 genes. The topmost DE genes in poly(I:C)-stimulated cells had thousand-fold changes with highly significant p-values, whereas in CpG DNA stimulated cells had 2–5-fold changes with less stringent p-values. The highest DE genes in both stimulations belonged to the TNF superfamily, TNFSF18 (CpG) and TNFSF10 (poly(I:C)) and in both cases the lowest downregulated gene was CYP1A1. CpG DNA highly induced canonical pathways that are unrelated to immune response in Bomac. CpG DNA influenced expression of genes involved in molecular and cellular functions in free radical scavenging. By contrast, poly(I:C) highly induced exclusively canonical pathways directly related to antiviral immune functions mediated by interferon signalling genes. The transcriptomic profile after poly(I:C)-stimulation was consistent with induction of TLR3 signalling. Conclusion CpG DNA and poly(I:C) induce different early transcriptional landscapes in Bomac, but each is suited to a specific function of macrophages during interaction with pathogens. Poly(I:C) influenced antiviral response genes, whereas CpG DNA influenced genes important for phagocytic processes. Poly(I:C) was more potent in setting the inflammatory landscape desirable for an efficient immune response against virus infection. Electronic supplementary material The online version of this article (10.1186/s12864-018-5411-5) contains supplementary material, which is available to authorized users.


Background
Gene expression profiles generated from a wide variety of cell types treated with different pathogens or substances that mimic pathogens can yield valuable insights into the mechanisms of host-pathogen interaction. Moreover, gene expression profiles permit identification of specific responses to microbial stimuli that can lead to rationale design of therapeutic approaches or next-generation vaccines. With the completion of the bovine genome [1], it is possible to analyze transcriptional profiles or interactomes to understand the host-pathogen relationships that lead to cellular responses. Therefore, highly sophisticated techniques such as RNASeq have become useful to study whole transcriptomic landscapes of cells during infection or treatment with therapeutic agents. There is now abundant evidence that the innate immune response immensely influences the development of an infection into disease. Among other cells at the forefront of innate immunity, macrophages (MΦ) play a major role in containing primary infections. For instance, infection of mice with herpes simplex 1 virus followed by induction of M2 macrophages through CSF-1 DNA was associated with reduced virus replication in the eye, reduced latency and reduced levels of CD4, CD8, IFN-γ and PD-1 transcripts in the trigeminal ganglion [2]. Some reports show that MΦ alter their gene expression profile upon infection, e.g., Mycobacterium bovis infection is associated with the repression of host gene expression in MΦ [3,4]. Therefore, reaction of MΦ to various pathogens is variable and is not yet completely understood. Lewandowska-Sabat et al. [5] have reported the early phase transcriptional program of bovine monocyte-derived MΦ infected with Staphylococcus aureus and show that S. aureus induces both, alternative and classical MΦ activation pathways. They concluded that activation of MΦ through the alternative pathway possibly contributes to intracellular persistence of S. aureus during mastitis in dairy cattle. Infection of an epithelial cell-MΦ co-culture with Mycobacterium avium subspecies paratuberculosis (MAP) revealed a number of metabolic, DNA repair and virulence genes that are worthy to investigate for new drug targets [6]. In particular, this study revealed a novel iron assimilation system for carboxymycobactin. Another RNASeq study of MAP infection [7] of monocyte-derived MΦ showed expression of genes that account for protective host immunity and those that might support MAP survival and proliferation in MΦ.
Antigen presenting cells (APCs), such as MΦ, express pattern recognition receptors (PRRs) such as Toll-like receptors (TLRs), which are used for detecting pathogen-associated molecular patterns (PAMPs). PRR signal through intermediate molecular adaptors to activate transcription factors that drive gene transcription and expression of pro-inflammatory cytokines responsible for antimicrobial activity. CpG DNA and synthetic dsRNA (poly(I:C)) are classically used as model ligands to represent specific pathogens such as bacteria and dsRNA viruses, and are potent adjuvants of immune response in many animal species. Therefore, it is imperative that the influence of adjuvants on the responding cells be thoroughly understood. During stimulation with dsRNA, the adaptor TRIF (TIR-domain-containing adapter-inducing interferon-β) is recruited to the Toll/ Il-1 Receptor (TIR) domain of activated, dimeric TLR3, which leads to stimulation of the transcription factors IRF3, NF-κB and AP-1 [8,9]. By contrast, the chaperon protein UNC93B mediates transport of TLR9 to the endosome where the binding with DNA molecules induces the formation of the signalling complex that includes MyD88, TRAF6 and IRAK4. The resulting response activates the NF-κB pathway leading to the production of pro-inflammatory cytokines. In addition, activation of IRF7 upon TLR9 activation rather leads to the production of type-I and -III interferon (IFN) which depend on numerous co-activators such as IRAK1, TRAF3, OPN, IKKα and PI3Kδ [10].
Bomac cells are often used but the responsiveness to CpG DNA and poly(I:C) has not been previously assessed using high throughput genomics techniques. To shed more light on the initial transcriptional response of an innate immune cell upon encounter with different PAMPs, we studied the influence of CpG DNA and poly(I:C) on a bovine MΦ cell line. Thus, the purpose of this study was to examine differential expression of genes in bovine macrophages under the influence TLR3 and TLR9 agonists in order to better understand the early phase of transcriptional changes induced by CpG DNA or poly(I:C) in the Bomac macrophage cell line. Our results show that poly(I:C) induces a transcriptional profile rather resembling an antiviral immune response, whereas CpG DNA induces a varied transcriptional profile with a highly marked function of free radical scavenging.

Cells
The Bomac cell line is a transformed bovine peritoneal macrophage cell line described by Stabel and Stabel [11]. The original Bomac cell line present in the laboratory of one of the authors [12] was found to be contaminated with bovine viral diarrhea virus (BVDV), probably present already in the original cell type. The cell line was subsequently cured of BVDV (Marti, S. and Schweizer, M., manuscript in preparation; [13], and maintained in IMDM supplemented with Glutamax, Hepes, 10% fetal calf serum, non-essential amino acids, MEM vitamins, penicillin and streptomycin, and 2-mercaptoethanol. The cured cell line used in these studies was obtained from Dr. Matthias Schweizer. Cells were cultured at 37°C in 5% CO 2 in the atmosphere and subcultured every two days. Additionally, the cells were tested for Mycoplasma contamination using the PlasmoTest™ Mycoplasma Detection Kit from Invivogen, and found to be negative (Additional file 1). Further, a conventional RT-PCR was performed to detect BVDV in the cured cells. We did not detect and BDVD sequences (Additional file 1).

Flow cytometry
Because the Bomac cells have been reported not to be a homogeneous population, and have lost certain functions, we sought to determine the phenotype of the cells that we were working with, by flow cytometry. Bomac cells were stained with anti-CD14 (cat. nr MCA2678GA), CD16 (cat. nr MCA5665GA), CD11b (cat. nr MCA1425F), CD172a (cat. nr MCA6079), CD44 (cat. nr MCA2433F), MHC II (cat. nr MCA2445PE) from BioRad, and CD40 (cat. nr ABIN2480301), CD68 (cat. nr ABIN2472322), CD71 (cat. nr ABIN2560526), CCR2 (ABIN2787667) from Antibodies Online. Another set of cells were stimulated with p(I:C) as described in the Cell stimulation section, and then similarly stained for flow cytometry. At least 100,000 events were acquired on a BD FACS Calibur flow cytometer, and later analyzed in FlowJo. A total of six determinations were performed and means compared using the descriptive statistical analysis tool in Microsoft Excel.

Cell stimulation
Cells were first cultured in 6 well plates at a density of 10 4 cells per well for 24 h at the time they reached confluence. The maintenance medium was then removed and replaced with medium containing CpG DNA or poly(I:C) at 10 μg/ml or 50 μg/ml, respectively. A separate culture plate was set with unstimulated cells to serve as control. We used the class B CpG oligonucleotide (ODN 2007) reactive with bovine TLR9, whereas the poly(I:C) was a low molecular weight synthetic analog of dsRNA specific for TLR3. Both reagents were obtained from InvivoGen (San Diego, USA). In order to analyze the early phase of transcriptomic changes, cells were incubated in the presence of PAMPs for only 6 h prior to RNA isolation. This time point was chosen based on ELISA assays measuring the earliest appearance of IL-1β, IL-6 and IFNβ protein in the supernatants following stimulation with CpG DNA or poly(I:C). For that purpose, Bomac cells, plated in 6-well plates, were stimulated with the individual PAMPs and incubated for 2, 4, 6, 8 and 12 h prior to collecting the cell culture supernatants and quantification of the cytokines by ELISA using reagents purchased from BioSource according to the manufacturers' instructions. The first cytokines to be detected were IL-6 and IFNβ at 8 and 12 h post  stimulation, whereas IL-1β protein was not detectable  at time points earlier than 18 h (data not shown). Therefore, we selected 6 h as the earliest time point to track early transcriptomic changes in Bomac cells stimulated with PAMPs.

RNA isolation
After the 6 h incubation period, the maintenance medium was removed from the wells containing stimulated cells followed by two washes with cold phosphate-buffered saline. Cells were lyzed with Qiazol directly in the wells and lysates were transferred to DNase/RNase free microfuge tubes for further processing. RNA isolation was carried out as described [14] and according to the manufacturers' instructions. The Qiagen miRNeasy Mini kit was used for RNA isolation. RNA concentration was measured using the Implen NanoPhotometer (Implen GmbH, Munich, Germany). RNA integrity was monitored by agarose RNA electrophoresis (Additional file 2). RNA samples were dried in GenTegra-RNA tubes (GenTegra, USA) using a VentriVap DNA concentrator (Labconco, USA) for 55 min. at 25-27°C and sent for sequencing. Each stimulation was repeated in three separate experiments and, thus, 9 RNA samples (control, CpG DNA-stimulated, poly(-I:C)-stimulated cells) were submitted for RNASeq.

mRNA -sequencing and analysis
The samples were sequenced at the Roy J. Carver Biotechnology Center at the University of Illinois, USA. Ribosomal RNA was removed with the Ribozero Human/mouse/rat kit (Illumina, USA). Strand-specific single-end libraries were prepared using the TruSeq Stranded mRNA-Seq Sample Prep kit (Illumina, USA). The libraries were quantitated by qPCR, pooled in equimolar amounts and sequenced on one lane of a HiSeq 4000 (Illumina; sequencing kit v 1), generating over 373 million 100 bp single-end reads. Fastq files were generated and demultiplexed per sample with the bcl2fastq v2.17.1.14 Conversion Software (Illumina), which also trims Illumina adapters from the reads and removes any resulting sequences shorter than 35 bp. All bases across the reads showed quality scores greater than Q30 (FASTQC v 0.11.2) and, thus, quality trimming was not performed. The Bos taurus UMD 3.1.1 reference genome was downloaded from NCBI along with Annotation release 105 gene models containing 35,315 Entrez Gene IDs. The gene models were converted from gff3 to gtf format using the rtracklayer package [15] (v 1.36.3) in R [16] (v 3.4.0) while extracting the Entrez Gene IDs to an attribute named "gene_id". Reads were aligned to the genome using STAR (version 2.5.2b) [17] with parameters --sjdbGTFfeatureExon exon, −-sjdbGTFtagExonPar-entGene gene_id and --sjdbOverhang 99. Read counts per gene were generated using featureCounts [18] from the subread package [19] (v 1.5.0). The datasets generated for this study can be found in the NCBI, GEO accession number GSE106843.
The read counts were put into R [16] (v 3.4.0) for data pre-processing and statistical analysis using packages from Bioconductor [20] as indicated below. The samples had 23-36 million reads assigned to genes (Table 1), thus any gene without 0.5 Count Per Million (CPM) reads in at least 3 samples were filtered out. 13,740 genes passed this filter and were analyzed using edgeR [21] (v 3.18.1) in the quasi-likelihood pipeline [22,23] that also accounted for the total library size for each sample and an extra TMM normalization factor [24] for any biases due to changes in total RNA composition of the samples. A one-way ANOVA test was employed from the model, along with pairwise comparisons between the two stimulation groups. Adjustment for multiple testing was done using the False Discovery Rate (FDR) method [25].
Gene Ontology terms for the B. taurus genes were analyzed from Bioconductor's org.Bt.eg.db (v 3.4.1) database while the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway mappings were performed directly from KEGG (http://www.genome.jp/kegg/) on July 12, 2017 using Bioconductor's KEGGREST package (v 1.16.0). Overrepresentation testing for GO and KEGG terms was done using the goana() and kegga() functions from edgeR [21], which correct for any trends due to expression level [26]. Separate tests were done for the upand downregulated genes for each of the three pairwise comparisons between stimulation groups. Comparison of GO terms' raw p-values between all gene sets from all pairwise comparisons was done using heat maps of -log10 (p-values) made separately for biological processes (BP), molecular functions (MF) and cellular components (CC) categories using terms with raw p-values <1e-5 in any comparison for BP and < 1e-4 for MF and CC, due to their having fewer terms. A similar heat map was made for KEGG terms with raw p-values <1e-5 in any comparison, and also pathway maps showing both upand downregulated genes were made for each pairwise comparison for any pathway with either up or downregulation with a p-value <1e-3 using the pathview package [27] (v 1.16.1), (data not shown). Figure 1 shows the mean fluorescence results from a single colour staining of Bomac cells with antibodies against several surface markers of bovine macrophages. These results largely confirm what others [12] have already shown, with the exception of CD16. We observed expression of CD16, CD172a, CD40 and CD44. Only CD16 and CD40 (p ≤ 0.05) were upregulated following stimulation of the cells with p(I:C). Additionally, we scrutinized the phagocytic function of Bomac. As shown in Additional file 1, the Bomac can phogocytose Staphylococcus aureus labelled with FITC, albeit at very low level compared to freshly isolated bovine monocytes. However, the capability to phagocytose is not homogeneous, as majority of the cells did not engulf bacteria.

Categorization and abundance of cDNA
The statistics of categorization and abundance of sequence reads generated from 9 cDNA libraries used for differential gene expression analysis are shown in Table 1. A one-way ANOVA analysis revealed differences between the 2 test groups. 2245 genes out of 13,740 expressed in the data set had FDR p-values ≤0.05 (Fig. 2). The hierarchical clustering dendrogram of Euclidian distances between the genes (rows) was cut at height = 4, resulting in ten clusters. In comparison to control experiments, the number of genes with pairwise FDR p ≤ 0.05 for CpG DNA-stimulated cells was 210 downregulated and 347 upregulated. In poly(I:C)-stimulated cells 414 Table 1 Statistics of categorization and abundance of sequence reads generated from 9 cDNA libraries used for differential gene expression analysis genes were downregulated while 761 were upregulated. When the 2 treatment groups were compared, 911 genes were found to be downregulated whereas 1240 genes were upregulated (Fig. 3). Tables 2 and 3 show the top 10 upregulated and downregulated genes at 6 h post stimulation in both CpG DNA vs control and poly(I:C) vs control ranked by fold-change.

Functional characterization of differentially expressed genes KEGG pathways
There are 312 different KEGG pathways, of which 82 were highly overrepresented (raw p-value <1e-6) in at least 1 of the 3 pairwise gene sets (6 in total). Figure 4 shows a heat map comparing the 6 gene sets across these 82 pathways. Both poly(I:C) vs Con and poly(I:C) vs CpG DNA show similar up and downregulated pathways, indicating that poly(I:C) treatment induces a very different transcriptional profile in Bomac cells at 6 h post-stimulation than CpG DNA. In Table 4 (see also Additional file 3) an expansion is made on these data to show the top 5 KEGG pathways in each data set and the top 5 immune related pathways. In CpG DNA-stimulated cells, the pathways included axon guidance, focal adhesion, and protein digestion and absorption. Among immune related pathways, PI3K-Akt signalling and Cytokine-cytokine receptor signalling were overrepresented (Table 4Ab). Poly(I:C) stimulation induced upregulation of a number of genes in many different pathways. Upregulated pathways included various disease and immune-related pathways while various signalling pathways were downregulated. Among the highly represented ones were the NOD-like receptor signaling pathway (39 genes), RIG-I-like receptor signaling pathway (22 genes), and all pathways indicating activation of mechanisms associated with response to viral pathogens. An example of this can be observed in the three viral pathways, herpes simplex, influenza, and measles virus infection (Table 4Ba). Comparison of CpG DNA and poly(I:C) stimulations showed marked domination of pathways related to response to infection, induced by poly(I:C) stimulation (data not shown).

Gene ontology
Gene ontology (GO) analysis was applied to these data to gain insight into the represented gene categories among the differentially expressed genes. Overall, 9490 terms were identified in the BP, 2629 terms in MF and 1273 terms in CC. A stringent raw p < 1e-5 was applied across the 6 gene sets, which reduced the number of terms to 76 BP terms, 13 MF terms and 11 CC terms. Due to a high number of GO terms in BP satisfying the condition raw p < 1e-5, the Panther functional classification tool was used to show only the main categories of GO terms. Figure 5 shows a pie chart representation of the Panther GO-Slim Bioprocess. In the bioprocesses, mainly cellular and metabolic processes were highly represented in all three pairwise comparisons. Immune system process upregulated genes represented 8.

Ingenuity pathway analysis (IPA)
Qiagen's Ingenuity Pathway Analysis was employed to further scrutinize the DE genes https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis [28]. The genes were mapped to human/mouse/rat genome for  To show relevant relationships between molecules in the 2 datasets, core analysis in IPA was performed to search for networks. In total, 11 networks were identified in the CpG DNA dataset and 81 in the poly(I:C) dataset. The top 5 in each dataset are shown in Table 6 (Fig. 7b). Further, the molecular and cellular function were closely examined in the datasets. Table 7 shows the top 5 significant functional networks in the two datasets. The topmost function in CpG DNA-stimulated cells was the free radical scavenging (p-value 1.45 × 10 − 8 ) containing 21 genes (e.g., NR4A1, + 3.2-fold change FDR p-value 6 × 10 − 3 , DDIT4 + 1.67-fold change FDR p-value, 5.9 × 10 − 4 , PLAUR + 1.76-fold change, FDR p-value 2.1 × 10 − 5 ). Figure 8a shows a network of molecules responsible  for formation of reactive oxygen species, generation of reactive oxygen species, synthesis of reactive oxygen species and production of reactive oxygen species in macrophages. The cellular function and maintenance was the most significantly represented function group in poly(-I:C)-stimulated cells (  (Fig. 8b). Upstream regulator analysis was performed to identify upstream regulators that might have influenced the gene expression changes observed in the two datasets. 1334 and 1721 molecules in CpG DNA and poly(I:C)-stimulated cells, respectively, were identified of which the top 5 regulators for each data set are listed in Table 8. The top most regulator in the CpG DNA dataset was the tumor necrosis factor (TNF) (upregulated, FDR p-value 2.3 × 10 − 10 ) and it affected 40 target molecules, while the poly(I:C) dataset revealed IRF7 (upregulated, FDR p-value 1.38 × 10 − 57 ) as the topmost regulator with 57 target molecules (Additional file 3). Mechanistic networks involving TNF and IRF7 to show a plausible set of connected upstream regulators that have contributed to the gene expression changes observed in the CpG DNA and poly(I:C) datasets are shown in Fig. 9.
Finally, the regulator effects in both datasets were addressed. The most significant regulator effects predicted for CpG DNA and poly(I:C)-stimulated cells are shown in Table 9. Three regulators Cg, HIF1A, miR-3202 (miR-NAs w/seed GGAAGGG) in CpG DNA-stimulated cells had the highest consistency score and regulated 5 molecules (FOXF1, MIST1R, NR4A1, PLAU, PLUAR) responsible for cell movement of embryonic cells (Fig. 10a). All these genes were significantly upregulated in the dataset. One of the regulators was a miRNA that apparently regulates FOXF1, a transcription regulator. Poly(I:C)-stimulated cells revealed IFN-α as the main regulator with the highest consistency score, regulating at least 13 targets (ADAR, APOBEC3B, CCL5, EIF2AK2, IFIT1, IFITM1, ISG15, ISG20, MX1, OAS1, RNASEL, RSAD2, ZC3HAV1) (Fig. 10b), all of them participating in inhibition of replication of viral replicons. Similar IPA analysis was done on the dataset comparing poly(I:C) to CpG

Discussion
RNASeq and systems biology tools were applied to unravel the transcriptomic landscape in vitro MΦ stimulated with CpG DNA or poly(I:C) as a model for bacterial and viral pathogen-associated molecular patterns, respectively. Using the data derived from gene expression analysis, IPA, KEGG and GO were used to reconstruct cardinal networks, pathways, and cellular function groups of genes differentially expressed in the Bomac cell line stimulated with CpG DNA or poly(I:C).
It should be noted that in IPA analyses, the genes were    The latter part of this observation was consistent with observations by others using monocyte-derived MΦ infected with MAP [4,7]. Further, unlike CpG DNA-stimulated cells, poly(I:C) stimulation of Bomac cells produced an expression profile that may be expected following infection with RNA viruses such as influenza or measles virus. Responses to these two viruses involve activation of transcription factor genes such as STAT1, STAT2 and NFκB1, which were upregulated in the Bomac cells stimulated with poly(I:C). Most canonical pathways represented within the DE genes in CpG DNA-stimulated cells were not related to immune functions, whereas the majority of canonical pathways in poly(I:C)-stimulated cells were related to immune functions. Stimulation of Bomac cells with CpG DNA highly affected the genes that regulate the phagocytic function of MΦ in contrast to stimulation with poly(I:C). This also confirms the importance to work with bovine cells in the absence of contaminations with ruminant pestiviruses, as BVDV was reported to inhibit the phagocytic activity of MΦs and neutrophils [12,29]. However, immunological studies with these cells should be cautious, as these cells appear to poorly represent a lineage-specific phenotype [12]. Common transcriptomic profiles between CpG DNAstimulated and poly(I:C)-stimulated cells were: (i) upstream and regulator effects, both stimulations induced  . TNFSF18 a ligand for TNFRSF18 expressed on the cell surface of regulatory T cells [30]. TNFSF18 is also released upon activation of DCs with CpG DNA [31]. In that report, plasmacytoid in this dataset. It is expressed at significant levels in most tissues including MΦs [32], and activates cysteine-type endopeptidase activity involved in apoptosis of transformed or tumor cells but does not kill normal cells [33]. However, our result is in contrast to that obtained by Casey et al. [7] who found that TNFSF18 was downregulated at 6 h post infection of monocyte-derived bovine MΦs with Mycobacterium avium subspecies paratuberculosis. Although the analytical technique used was the same as ours, in their case, cells were derived from live animals compared to our bovine MΦ cell line. The most downregulated gene in both treatments was CYP1A1 (cytochrome P450, subfamily I (aromatic compound-inducible), polypeptide 1, transcript variant X2). We could not find a specific role described for CYP1A1 among the immune functions of MΦs, but assume that the gene product of CYP1A1 could be involved in other cellular processes. According to a review by Androutsopoulos et al., [34] CYP1A1 plays an important role in the detoxication of environmental carcinogens, as well as in the metabolic activation of dietary compounds with cancer preventative activity. The KEGG pathways overrepresentation results were very similar to those generated in the IPA pathways analyses, both reporting a marked polarization of the response of MΦ to poly(I:C) towards antiviral activity.
In contrast to stimulation of Bomac with CpG DNA, cells stimulated with poly(I:C) had the highest number of DE genes suited to antiviral response of MΦs. IFIT2 (> 2300-fold change, FDR p-value 1.5 × 10 − 13 ) induced by interferon and regulated by IRF1 (+ 4.8-fold change, FDR p-value 5.13 × 10 − 16 ) gives rise to gene products that have antiviral functions [35]. IRF1 further regulates RSAD2 (> 2258.2-fold change, FDR p-value 5.2 × 10 − 25 ), which has been reported to be highly expressed during viral infections [36,37]. IFI44L and IFI44 were upregulated (fold change + 2010.7 and + 1254.3, FDR p value 4.8 × 10 − 13 and 2.6 × 10 − 11 , respectively) and play a role in immune responses against viruses, although the precise mechanism has not been described. The remaining highly expressed genes OAS1, ISG15, MX1, MX2 and OAS2 are well characterized antiviral response genes in many animal species [35,[38][39][40]. A study in mouse central nervous system tissue by Pomeranz et al. [41] showed a very similar profile to ours of highly expressed antiviral genes after infection with an RNA virus. In regard to highly expressed genes in the CpG DNA-stimulated cells, genes such as SLC26A10, IGFLR1, GSDMB, LBH, MYOZ2 and RRAD have no clearly described role in immune responses of MΦs. GSDMB appears to be associated with asthma and certain types of autoimmune diseases such as, type 1 diabetes, inflammatory bowel disease, and rheumatoid arthritis [42,43]. Ekwall et al. identified LBH as a candidate gene in rheumatoid arthritis and it was expressed in the synovial lining layer in such patients [44]. NR4A1 is a nuclear transcription factor that belongs to the nur77 signalling in T lymphocytes and calcium-induced T lymphocyte apoptosis canonical pathways. These pathways were both represented as top canonical pathways in the CpG DNA dataset. NR4A1 is expressed in MΦs in response to pro-inflammatory stimuli [45]. PLAU participates in the coagulation system canonical pathway represented as the 5th ranking canonical pathway in the CpG DNA data set. Expression of PLAU appears to be regulated by TGFB1 [46]. Such a pattern of gene expression provides evidence for the advantage of whole genome analysis because it allows the detection of other genes highly expressed that might not directly relate to the immune function of MΦ but rather account for other functions of cells during an immune response. This provides an opportunity to further study functions of such genes in order to fully understand regulation of immune processes in which bovine MΦs may be engaged. An interesting finding was the occurrence of PRL (encoding prolactin) as a regulator in poly(I:C)-stimulated cells. Although prolactin is widely known for its essential role in lactation, it appears to have important roles in immune processes, indicating that MΦ may have a prolactin regulated function during reproduction [47,48]. Close examination of the canonical pathways based on IPA analysis of poly(I:C) dataset shows that majority of these pathways are directly involved in immune response of MΦs. In the interferon pathway, at least five genes were interferon-inducible (IFIT3, IFI35, IFIT1, IFI6, IFITM1), and two genes were interferon regulatory genes (IRF1, IRF9) and the remaining 8 genes (SOCS1, OAS1, MX1, PSMB8, TAP1, ISG15, STAT2, STAT1) have activity highly regulated by interferon. Signaling through TLR3 upon stimulation with synthetic dsRNA (poly(I:C)) leads to production of IFNβ [49]. It is likely that IFNβ influenced the expression of the genes mentioned above, particularly that supernatants from Bomac cells cultured in the presence of poly(I:C) contained a detectable amount of IFNβ measured at 8 h of stimulation (data not shown) compared to cells stimulated with CpG DNA. Sivori et al. [50] studied CpG DNA and poly(I:C) in regard to their ability to stimulate NK cells through TLRs and found that both could stimulate NK cells to Fig. 9 Mechanistic networks generated from differentially expressed genes in datasets derived from Bomac cells stimulated with CpG DNA or poly(I:C) for 6 h. a CpG DNA, TNF upstream regulator, b poly(I:C), IRF7 upstream regulator. The mechanistic networks were created with the use of IPA (QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis   acquire cytotoxic activity against tumor cells. However, poly(I:C) induced more robust responses than CpG DNA, because poly(I:C)-stimulated NK cells produced more IFNγ and TNFα and were highly cytotoxic compared to CpG DNA-stimulated NK cells. The fact that pro-inflammatory cytokines are released was shown in the study by Casey et al. [7] involving bovine MΦs infected with MAP where IL-1 was produced as early as 10 min under their experimental conditions. The canonical pathways and networks, molecular and cellular functions as identified by IPA analysis of DE genes differed greatly between CpG DNA and poly(I:C) stimulated cells. In particular, CpG DNA promoted expression of genes involved in free radical scavenging, a group of molecules responsible for primary function of MΦs. Among the 21 genes in this group of molecules, 13 (ABCA1, AGER, APLN, BNIP3, DDIT3, DDIT4, HK2, HSPB6, NR4A1, PLAU, PLAUR, S100A6, SREBF1) were significantly upregulated with p-values ranging from 5.8 × 10 − 4 to 1.45 × 10 − 8 . The downregulated genes in this function group could possibly be regarded as control mechanisms in the processes involving reactive oxygen species. In human MΦs AGER encodes a nucleic acid receptor that promotes inflammatory responses to DNA [51]. BNIP3 and DDIT4 are directly involved in generation of reactive oxygen species [52,53] and formation of reactive oxygen species [54]. The remaining molecular and cellular functional groups in both stimulation groups mainly concerned cell development, morphology, survival, signalling or death. Further, on the molecular functions, careful examination of GO terms revealed that many top functions included chemokine activity. At least seven CXCL chemokines were found in both CpG DNA and poly(I:C) datasets. However, upregulation of CXCL5, CXCL8, CXCL16 was observed only in poly(I:C)-stimulated Bomac cells and had an unaltered expression in CpG DNA-stimulated Bomac cells with the exception of CXCL5, which was also upregulated in CpG DNA-stimulated cells. CXCL5 has been reported to activate MΦ and also increases the expression of ABCA1 (mRNA upregulated in this dataset) [55].
Although we show marked differences in the transcriptomic landscape of MΦs treated with a viral or a bacterial PAMP, we cannot generalize the differences or compare them to in vivo conditions because, first, the CpG DNA and poly(I:C) used are highly purified reagents much different from live pathogens; second, we have used a transformed cell line that may not accurately represent the transcriptomic changes that occur in vivo (compare [12]). However, the results obtained in this study are largely in line with what is known about MΦ stimulation through TLR and further provide an insight into the early phase of transcriptomic alterations in MΦs during initial interaction with agonists through their pattern recognition receptors.

Conclusion
Very early following stimulation of Bomac, poly(I:C), more than CpG DNA, triggered signals that exerted a transcriptional profile suited to an antiviral response, whereas CpG DNA influenced genes important for the phagocytic processes. Besides influencing the genes related to immune function, CpG DNA highly affected many genes in metabolic pathways and other non-immune biological processes. Although both, poly(I:C) and CpG DNA, are proposed vaccine adjuvants, poly(I:C) appears to be more potent in setting up the inflammatory landscape required to induce an efficient immune response.