Whole-Blood Gene-Expression Profiles of Cows Infected with Mycobacterium avium subsp. paratuberculosis Reveal Changes in Immune Response and Lipid Metabolism

Abstract

Mycobacterium avium subsp. paratuberculosis (MAP) is the causative agent of Johne's disease, a chronic debilitating disease affecting ruminants worldwide. In the present study, we aimed to determine the major gene networks and pathways underlying the immune response to MAP infection using whole-blood cells, as well as provide the potential transcriptional markers for identifying the status of MAP infection. We analyzed the transcriptional profiles of whole-blood cells of cattle identified and grouped according to the presence of MAP-specific antibodies and the MAP shed by them. The grouping was based on the results obtained by ELISA and PCR analyses as follows: i) Test1 group: MAP-negative results obtained by ELISA and positive results obtained by PCR; ii) Test2 group: MAP-positive results obtained by ELISA and negative results obtained by PCR; iii) Test3 group: MAP-positive results obtained by ELISA and positive results obtained by PCR; iv) uninfected control: MAP-negative results obtained both by ELISA and PCR analysis. The results showed down-regulated production and metabolism of reactive oxygen species in the Test1 group, activation of pathways related to the host-defense response against MAP (LXR/RXR activation and complement system) in the Test2 and Test3 groups, and anti-inflammatory response (activation of IL-10 signaling pathway) only in the Test3 group. Our data indicate a balanced response that serves the immune-limiting mechanism while the host-defense responses are progressing.

Introduction

Paratuberculosis (Johne’s disease) is a chronic and debilitating disease of ruminants, caused by infection with Mycobacterium avium subsp. paratuberculosis (MAP) [22,34]. The disease is characterized by persistent diarrhea, progressive wasting, and eventual death in ruminants, for which there is no available treatment [7,9,21]. The prevalence of MAP was 3-15% in herds in most regions, but can exceed 50% in dairy herds in some areas in major dairy-producing countries such as the United States and Europe [6,12,26,34]. There is a particular interest in controlling MAP infection owing to the huge economic losses in the dairy industry caused by this disease and its possible relationship with Crohn’s disease in humans [26]. The net economic losses due to Johne’s disease in the US dairy industry have been estimated to be 200-1500 million dollars per year, due to decreased milk production, reproductive dysfunctions, poor feed conversion, shortened production age, and increased susceptibility to other diseases [9,11].

According to the symptomatic assessment and quantification of the shed bacteria, Johne’s disease can be divided into four stages of progression: silent, subclinical, clinical, and advanced cellular infection. During the silent stage, the cattle do not shed any detectable amounts of bacteria and show no symptoms [33]. In the subclinical stage, the animals shed a small amount of MAP in their feces and milk, thereby contaminating the surrounding habitat and spreading MAP throughout the herd [33]. Bacteria are shed at high levels in the clinical stage, and then in the advanced cellular infection stage, typical disease symptoms begin to appear [33]. The disease progresses from the silent to the subclinical stage without any observable symptoms, which leads to the spread of the bacteria to the entire herd. Hence, the accurate detection and culling of infected animals at the early stage of MAP infection are important for control of Johne’s disease.

Currently, functional genomic technologies are being used to investigate the cellular pathways and molecular mechanisms implicated in the host immune response to mycobacterial infection in order to understand the disease pathogenesis [15,16,18,21,35]. Based on genomic profiling data, some reports have identified potential molecules that are critical for host-pathogen interactions during the infection [3,30,35]. In addition, many groups have performed transcriptional profiling of macrophages or peripheral blood mononuclear cells from Mycobacterium-infected animals; thus, transcriptomic research for the hostimmune response against mycobacterial infection has already been recognized in immune cells. In the present study, we characterized the transcriptional profiles of whole-blood cells in cattle, which were identified and grouped according to the presence of MAP-specific antibodies and shed MAP bacteria. Whole-blood cells contain dynamic and interactive information in the body, such as changes in association with a disease process; therefore, transcriptional profiling from whole blood has the potential to identify genetic biomarkers for diseases [13,24,31].

In the current study, we performed the transcriptional profiling of peripheral blood from eight MAP-infected cows using microarray analysis. The resulting transcriptional profiles were characterized by functional, network, and canonical pathways using Ingenuity Systems Pathway Analysis (IPA) Knowledge Base. These results may provide novel insights into understanding the host-immune response during the different stages of progression of MAP infection. Therefore, this research can potentially discover molecules that may serve as reliable transcriptional markers of MAP infection. This will in turn help to improve traditional diagnostic methods.

 

Materials and Methods

Animals and Sample Source

All animal procedures were performed according to the guidelines of the Institutional Animal Care and Use Committee at Seoul National University. The cattle used in this study were Holstein cows, which were evaluated for disease status by a commercial ELISA (IDEXX Laboratories, Inc., Westbrook, ME, USA) for the detection of MAP-specific antibodies and for IS900 and ISMap02 by fecal PCR. The fecal PCR was performed as previously describ ed [4,5]. Eight cows were divided into four groups: control (n = 2) with MAP-negative results ob tained b y ELISA and PCR, Test1 (n = 2) with MAP-negative results obtained by ELISA and positive results obtained by PCR, Test2 (n = 2) with MAP-positive results obtained by ELISA and negative results obtained by PCR, and Test3 (n = 2) with MAP-positive results obtained by both ELISA and PCR. Blood samples from the cattle were immediately collected into PAXGene Blood RNA tubes (PreAnalytiX/Qiagen, Hilden, Germany).

Microarray Hybridization

The total RNA was isolated from blood of the cattle according to the manufacturer's instructions using a PAXGene blood RNA kit (PreAnalytiX/Qiagen). RNA quantity and quality were assessed by a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA). All samples fulfilled the criteria of a 260/280 ratio greater than 1.7, a 28S/18S rRNA ratio greater than 1.7, and RNA integrity numbers greater than 8.

Labeled cRNA was prepared from 5 µg of total RNA according to the manufacturer’s instructions using Agilent’s Quick Amp Labeling Kit (Agilent Technologies, Inc.). Following fragmentation, 1.65 µg of labeled cRNA was hybridized to the Agilent expression microarray according to the protocols provided by the manufacturer. The arrays were scanned using the Agilent Technologies G2600D SG12494263 (Agilent Technologies, Inc.). Array data export processing and analyses were performed using Agilent Feature Extraction v11.0.1.1 (Agilent Technologies, Inc.). RNA amplification, cDNA labeling, array hybridization, and scanning for the microarray were carried out by Macrogen Inc. (Seoul, Republic of Korea). Microarray experiments were repeated three times for each sample.

Analysis of Microarray Data

The Bovine Oligo Microarray Chip (Bovine 4X44K G2519F) from Agilent was used in this study. The array contains 43,803 bovine probes that were developed by clustering more than 450,000 mRNA and EST sequences of the bovine genome (btau 4.0 assembly). Raw data were extracted using the software provided by the manufacturer (Agilent Feature Extraction v11.0.1.1). The selected gene signal value was transformed by a logarithm and normalized by the quantile method. The statistical values for the expression data were determined using fold change and the local pooled-error test. Correction for the false discovery rate (FDR) was performed by the Benjamini and Hochberg multiple testing method with a FDR adjusted p-value cut-off at ≤0.05 and a foldchange cut-off of ≥1.5. All data analyses and visualization of differentially expressed genes were conducted using R 2.15.1 (http://www.r-project.org).

Functional Analysis of Gene Expression Data

Data were analyzed using Qiagen’s Ingenuity Pathway Analysis (IPA; Ingenuity Systems Inc., Redwood City, CA, USA) for biological processes, canonical pathways, and networks analysis. The DEGs with an adjusted p ≤ 0.05 and fold changes ≥1.5 or ≤-1.5 were filtered and uploaded into the IPA program. Each gene was mapped to its corresponding gene object in Ingenuity’s Knowledge Base. Go-ontology analysis was performed using IPA to determine the biological functions of the differentially expressed genes in the MAP-infected cattle. Each GO category was ranked based on the number of DEGs that fit in each functional group. Right-tailed Fisher’s exact test was adopted to measure the p-value for each of the biological functions. Canonical pathways, from the IPA library of canonical pathways, were investigated to identify major biological pathways associated with MAP infection in cattle. The significance of the association between the data set and the canonical pathway was determined based on two parameters: (i) a ratio of the number of genes from the data set that map to the pathway divided by the total number of genes that map to the canonical pathway and (ii) a p-value calculated using the Fisher's exact test determining the probability of the association between the genes in the data set and the canonical pathway being due to chance alone. In the case of upstream regulators, the predicted activation state and activation z-score were based on the foldchange values for genes in the input data set for which an experimentally observed causal relationship has been established. The ten upstream regulators were investigated based on an activation z-score, which, when below (inhibited) or above (activated) 2, was considered as significant. Additionally, canonical pathways and networks were presented as upstream and downstream effects of activation or inhibition of molecules predicted by the Molecule Activity Predictor tool from the IPA software.

Reverse Transcriptase and Quantitative Real-Time PCR Validation

The total RNA was reverse transcribed with random primers by using the QuantiTect Reverse Transcription Kit (Qiagen Inc., Valencia, CA, USA) following the manufacturer’s instructions. Eight genes selected with differential expression were analyzed by real-time quantitative RT-PCR in order to validate the microarray results (Table 3). Real-time quantitative RT-PCRs were performed with 2 µl of cDNA using the Rotor-Gene SYBR Green PCR kit (Qiagen Inc.) and Rotor-Gene Q real-time PCR cycler (Qiagen Inc.). Amplification was performed for 40 cycles at 95℃ for 20 sec, followed by 45 sec at 60℃, with fluorescence detected during the extension phase. The expression level was determined by the 2-∆∆Ct method using a housekeeping gene, beta-actin, as a reference. The expression level was compared with the control group to determine the expression-fold change of each gene. A Pearson correlation coefficient was calculated for logarithmically transformed data obtained with quantitative RT-PCR and microarray analysis, using the IBM Statistical Package for Social Sciences software (SPSS, ver. 21; IBM SPSS Inc., Chicago, IL, USA). Differences were considered significant if a value of p < 0.05 was obtained.

 

Results and Discussion

Gene Expression Profiling of Blood from Cattle Naturally Infected with MAP

Whole-blood samples from 275 cows were examined by serum enzyme-linked immunosorbent assays (ELISAs) for MAP-specific antibodies and by fecal PCR testing for IS900 and ISMap02. Among them, eight cows were selected and divided into four groups depending on the results of the ELISA and fecal PCR analysis: Test1 with ELISA-negative and PCR-positive results; Test2 with ELISA-positive and PCR-negative results; Test3 with ELISA-positive and PCR-positive results; and control with ELISA-negative and PCRnegative results. In particular, Test2 and Test3 animals had moderate (S/P ratio ≥55 and <100) and high levels (S/P ratio ≥200) of MAP-specific antibodies, respectively. Transcriptional profiles of Test1, Test2, and Test3 group animals were compared with the control, and the gene expression was considered significant when p was ≤0.05 and the fold-change of expression was ≥1.5 for both upand down-regulation. Scatter plots of normalized signal values (Fig. 1A) show that the number of differentially expressed genes (DEGs) increased from Test1 to Test2 to Test3 animals. Among the 12,029 genes analyzed, 127 (1.06%), 443 (3.68%), and 655 (5.45%) genes were differentially expressed in Test1, Test2, and Test3 animals, respectively. Among the total genes, 61, 354, and 510 genes were upregulated and 66, 89, and 145 genes were down-regulated in Test1, Test2, and Test3, respectively (Fig. 1B). In total, 32 and 8 genes were commonly up- and down-regulated among the groups, respectively. In particular, Test2 and Test3 animals had many genes in common, including 280 up-regulated and 57 down-regulated genes (Fig. 1C). A list of the DEGs and the relative fold changes in Test1, Test2, and Test3 animals is provided in Tables S1, S2, and S3, respectively. The top 10 up- and down-regulated genes were enriched with functional annotation in Tables 1 and 2, respectively. The most up-regulated genes play a role in immune response (DEFB10, KLK12, DEFB1, DAB2, LAP, LOC616942, BOLA, JSP.1, and DEFB5), metabolic processes (ANG2, PRSS2, KLK12, LAP, DEFB10, and DEFB5), localization (HBE4, HBM, and TTYH2), and cellular processes (IGFBP1 and CDH10). In the most up-regulated genes, beta-defensin genes that efficiently kill bacteria and can participate in the early elimination of bacilli were detected in Test1 (DEFB10, DEFB1, and LAP) and Test3 (DEFB5) animals [29]. Although DEFB4A (2.73-fold up-regulation; 4.12-fold up-regulation), DEFB5 (2.67-fold up-regulation; 4.37-fold up-regulation), and TNF-alpha (2.02-fold up-regulation; 2.33-fold up-regulation) genes were not included in the top gene expression list, they were increased in both Test2 and Test3, and are related to the induction of beta-defensin. These data suggest that the antimicrobial peptides that increased in the animals that shed MAP bacteria (Test1 and Test3 animals) and the defensins that increased in the animals with MAP-specific antibodies (Test2 and Test3 animals) were associated with the maintenance of latent infection, constituting a bridge between innate and adaptive immunity [29]. The most down-regulated genes were involved in immune system processes (BOLA-DQA5, KLRC1, LOC100140174, P TX3, ADRA1B, LOC100126815, GBP6, ULBP 3, ULBP27, and CD163L1) and metabolic processes (RGS20, ADRA1B, and MYO5B). The top up- and down-regulated genes overlapped between the Test2 and Test3 groups. In particular, these genes, such as BOLA, JSP.1, ULBP 3, ULBP 27, BOLA-DQA5, LOC100126815, KLRC1, LOC100140174, and LOC616942, were involved in antigen processing and presentation of antigen, indicating a downregulation of MHC class I and II genes in the Test1 group, up-regulation of MHC class I genes and down-regulation of MHC class II genes in Test2 animals, and up-regulation of MHC class I and II genes in the Test3 group. KLRC1 and LOC100140174, as killer cell lectin-like receptor subfamily C, are genes related to modulating antiviral, bacterial, protozoal, and cancer-associated responses in gamma delta T lymphocytes [28]. They were strongly down-regulated in the Test1 (4.8-fold down-regulation) and Test3 (1.78-fold down-regulation) groups, demonstrating down-regulation of natural killer (NK) cell activity and enhanced survival of the shed MAP bacteria [8]. Additionally, interferoninducible guanylate-binding protein (Gbp) 6, one of the Gbps (Gbp1, Gbp6, Gbp7, and Gbp10) that may play a role in the innate immune response by regulating autophagy formation against intracellular pathogens, was downregulated in Test2 and Test3 [17,23].

Fig. 1.Gene expression levels of cattle infected with Mycobacterium avium subsp. paratuberculosis. The animals were grouped on the basis of their Mycobacterium avium subspecies paratuberculosis (MAP)-infection status determined by the results of ELISA and PCR analyses. The groups were defined as follows: ELISA-negative and PCR-positive results (Test1), ELISA-positive and PCR-negative results (Test2), ELISA-positive and PCR-positive results (Test3), and ELISA- and PCR-negative results (controls). (A) Scatter plots comparing gene expression levels between MAP-infected and control groups. Red dots indicate an expression level change of |fold change| ≥ 2. Expression level was calculated using a base-2 logarithm of the normalized hybridization signals from each sample. (B) Numbers of genes with altered expression levels during this experiment. (C) Venn diagram showing the overlapping genes that were significantly up- or down-regulated in MAP-infected cattle.

Table 1.aThe overlapped genes among the experimental groups are indicated and presented as (fold-changes).

Table 2.aThe overlapped genes among the experimental groups are indicated and presented as (fold-changes).

Table 3.aOnly pathways with the number of associated genes ≥5 are shown. bGenes/Total = number of differentially expressed genes (1.5-fold change; p < 0.05) out of total of genes associated with the canonical pathway according to IPA.

Functional Analysis of DEG in MAP-Infected Cattle

In total, 113, 328, and 528 up- and down-regulated transcripts in Test1, Test2, and Test3, respectively, were mapped to molecules in the Ingenuity Knowledge Base. Like the previous results, such as those shown by the Venn diagram and the top 10 regulated genes, the DEGs in Test2 and Test3 animals were involved in the same categories of “molecular and cellular functions,” “physiological system development and functions,” and “diseases and disorders” (Fig. 2). Although the DEGs in the Test2 and Test3 groups were found to be involved in the same categories of biological functions, the number of associated genes in Test3 animals was greater than in Test2 animals. In the Test1 group, the related biological functions included a small number of genes and different categories, such as lipid metabolism and free radical scavenging, compared with the Test2 and Test3 groups.

Fig. 2.Enrichment of the top five biofunctions in cattle infected with Mycobacterium avium subspecies paratuberculosis. (A) Physiological system development and function, (B) molecular and cellular functions, and (C) diseases and disorders attributed to the differentially expressed genes in cattle infected with Mycobacterium avium subspecies paratuberculosis (MAP). The top five biofunctions were derived from Ingenuity Systems Pathway Analysis of differentially expressed genes in the MAP-infected cattle.

As shown in Fig. 3, the representative network for the identified genes in the Test1 group showed enrichment of factors associated with cellular movement, immune cell trafficking, cellular function, and maintenance; they included haptoglobin (HP, 2.7-fold down-regulation), myeloperoxidase (MPO, 1.8-fold down-regulation), C-C motif ligand 5 (CCL5, 2.1-fold down-regulation), pentraxin 3 (PTX3, 3.5-fold down-regulation), oxidized low density lipoprotein (lectin-like) receptor 1 (OLR1, 2.5-fold upregulation), interleukin 12B (IL12B, 2.3-fold up-regulation), and serpin peptidase inhibitor, clade F, member 1 (SERPINF1, 1.9-fold up-regulation). This network showed the relationships in the DEGs as well as the predicted upstream and downstream effects of activation or inhibition of molecules. In particular, 23 of the DEGs in Test1, including the genes in this network, were selected and the activation state of diseases or function annotation was predicted by IPA. The inhibition of functions in the “free radical scavenging” category, including metabolism of reactive oxygen species (activation z-score = -2.309, p-value = 7.31E-05), synthesis of reactive oxygen species (z-score = -2.384, p-value = 1.37E-05), and production of reactive oxygen species (z-score = -2.758, p-value = 1.23E-02), were significantly predicted and contained 16 genes: HP, MPO, CYGB, GZMH, GZMA, GZMB, CCL4, CCL5, SERPINF1, FBP1, RETN, FOS, CXCL12, CEBPE, ALOX5, and OLR1. The inhibition of functions in the “lipid metabolism” and “small molecule biochemistry” categories, including synthesis of lipids (z-score = -2 .392, p-value = 1.03E-02 ), synthesis of fatty acids (z-score = -2.165, p-value = 1.34E-02), and synthesis of eicosanoids (z-score = -2.162, p-value = 1.23E-02), was also predicted and associated with 13 genes: ME1, ANG, PTX3, PDK4, FOS, RETN, IGFBP7, ALOX15, CCL5, IL1R2, CXCL12, ALOX5, and OLR1. Myeloperoxidase (MPO), which plays an antimycobacterial role as one of the most abundant granule proteins in neutrophils, was down-regulated in both Test1 and Test3 animals, suggesting that the shed MAP bacteria interfered with the microbicidal mechanism [2].

Fig. 3.Representative network of the alerted genes in the Test1 group. Individual nodes represent proteins with relationships represented by edges. Nodes are colored on the basis of gene expression, where red indicates up-regulation, green indicates down-regulation, and white indicates that the gene/factor is not differentially expressed, but has a defined relationship with other genes in the network. Arrows indicate directional relationships.

Canonical Pathways Associated with Differentially Expressed Genes in MAP-Infected Cattle

The top canonical pathways for the DEGs in Test1, Test2, and Test3 animals are reported in Table 3. Since the canonical pathways in Test1 animals included a very low number of DEGs, we could not determine their relationship to MAP infection (data not shown). In addition, all four canonical pathways in Test2 animals were also seen in Test3 animals, but the DEGs associated with the pathways in Test3 animals were highly enriched compared with Test2 animals. The genes whose products play a role in the five canonical pathways were up-regulated in Test2 and Test3 animals. These genes were involved in multiple pathways related to the immune response (Liver X receptor/ Retinoid X receptor (LXR/RXR) activation, complement system, interleukin-10 signaling, dendritic cell maturation, and recognition of bacteria and viruses by patternrecognition receptors), and lipid metabolism (LXR/RXR activation) and were strongly activated in both Test2 and Test3. The LXR/RXR activation pathway was primarily identified in both Test2 and Test3 animals as ELISA-positive in the naturally infected cattle. As shown Fig. 4, LXR activation in macrophages showed predicted activation and expression of several genes involved with inflammatory mediators, including genes encoding inducible NO synthase (iNOS), cyclooxygenase 2 (COX 2), IL-1, IL-6, and matrix metallopeptidase 9 (MMP9), as well as those involved in cholesterol trafficking and efflux, including genes encoding ATP-binding cassette, sub-family A, member 1 (ABCA1), ATP-binding cassette, sub-family G, member 1 (ABCG1), ATP-binding cassette, sub-family G, member 4 (ABCG4), and apolipoprotein E (ApoE). Several reports have demonstrated that Mycobacterium in the nutrient-deficient macrophage phagosome uses lipids as its main carbon source; therefore, lipid metabolism is involved in virulence and is important for long-term survival of the bacilli in macrophages [10,25,32]. In addition, cholesterol has been reported as an essential factor for the phagocytosis of the bacterium by macrophages and the inhibition of phagosome maturation [14,25]. LXR/RXR activation induced the protective immune response against M. tuberculosis [25], and the pathway was also identified in MAP-infected bovine macrophages [20]. The naturally infected cattle that showed MAP-specific ELISA-positive results also showed induction of LXR/RXR activation in the present study, and it is thought to be due to host responses modulating both innate and acquired immune responses as well as lipid metabolism.

Fig. 4.LXR/RXR activation pathways in the Test3 group. Nodes are colored on the basis of gene expression and predicted activity of upstream and downstream molecules, where red indicates up-regulation, green indicates down-regulation, orange indicates predicted activation, and blue indicates predicted inhibition. The intensity of the node color indicates the degree of regulation. Genes shown as uncolored nodes were not differentially expressed in our experiment and were integrated into the computationally generated networks based on the evidence stored in the IPA Knowledge Base.

Test2 and Test3 animals also showed induction of the complement system pathway with up-regulation of the complement component 1, q (C1q), C2, C3, and C3-receptor genes (C3aR) (Fig. 5). According to several reports, the complement receptor is used by pathogenic mycobacteria to gain entry into macrophages; [30] and M. tuberculosis-specific antibodies increase complement activation through both classical and alternative pathways, thus enhancing phagocytosis of the antibody-opsonized bacteria by macrophages [19]. The complement system pathway was more strongly induced in Test3 animals than in Test2 animals, which is thought to be due to a higher MAP-antibody titer in the cattle of Test3 compared with those of Test2. Among the three complement pathways, activation of the alternative pathway was predominantly observed, as shown in Fig. 5.

Fig. 5.Complement system pathway in the Test3 group. Nodes are colored on the basis of gene expression and predicted activity of upstream and downstream molecules, where red indicates up-regulation, green indicates down-regulation, orange indicates predicted activation, and blue indicates predicted inhibition. The intensity of the node color indicates the degree of regulation. Genes shown as uncolored nodes were not differentially expressed in our experiment and were integrated into the computationally generated networks based on the evidence stored in the IPA knowledge base.

Among the signaling pathways, the IL-10 signaling pathway was seen to be induced in Test3 animals (Fig. 6). In order to adopt various immune-evasion strategies and sustain against long-term infections of M. tuberculosis, the macrophages actively suppress the host-protective immune responses by secreting high levels of immunosuppressive cytokine IL-10 [1]. Therefore, these data suggest that the activation of the IL-10 signaling pathway in Test3 animals was a balanced response that serves the immune-limiting mechanism while the host-defense responses are progressing.

Fig. 6.IL-10 signaling in the Test3 group. Nodes are colored on the basis of gene expression and predicted activity of upstream and downstream molecules, where red indicates up-regulation, green indicates down-regulation, orange indicates predicted activation, and blue indicates predicted inhibition. The intensity of the node color indicates the degree of regulation. Genes shown as uncolored nodes were not differentially expressed in our experiment and were integrated into the computationally generated networks based on the evidence stored in the IPA knowledge base.

Predicted Upstream Regulators of DEG in MAP-Infected Cattle

The inhibition or activation of upstream regulators was predicted based on the absolute activation z-scores obtained by IPA. The 10 upstream regulators were considered when comparing the Test1, Test2, and Test3 groups (Fig. 7A). The regulators showing a different activation state in each group were selected for further investigation as biomarker candidates for the diagnosis of MAP infection. STAT3, IL-27, ERK, and CD28, transcription factors showing a significant absolute activation z-score in Test1, and STAT3 showing significant activation in Test3, were visualized in their predicted networks with their target genes (Fig. 7). STAT3 is essential for anti-inflammatory responses, and the pathway mediated by STAT3 is known to be associated with blocking phagosome maturation in Mycobacterium-infected macrophage and developing an immunosuppressive niche in granulomas [20,27]. Among the target genes of STAT3, we found some genes showing different gene expression patterns in each group, suggesting that they could be utilized as candidate diagnostic biomarkers for MAP infection: CCL4, 2.27-fold down-regulated in Test1; CCL5, 2.26-fold down-regulated in Test1; HP, 2.79-fold down-regulated in Test1; SERPPINE1, 2.12-fold up-regulated in Test3; IL-10, 2.69-fold up-regulated in Test3; MMP-9, 2.03-fold up-regulated in Test3; HGF, 2.03-fold up-regulated in Test3. All genes are located in the extracellular space. However, further studies checking the gene expression in a larger number of cattle is needed to determine potential diagnostic biomarkers.

Fig. 7.Differentially expressed target genes of predicted upstream regulators in cattle infected with Mycobacterium avium subspecies paratuberculosis. IPA-derived upstream regulators were selected by absolute activation z-score (A). STAT3 is significantly involved in the control of target genes (p ≤ 0.05 for Test1 (B) and Test3 (C)). Genes in red are up-regulated and genes in green are down-regulated.

Validation of Microarray Data

Quantitative RT-PCR (qRT-PCR) validation of the microarray results was performed for seven putatively differentially expressed genes. The PGA5, JSP, DEFB4, BOLA-DQA5, and BOLA genes from Test1 animals, TRD@, JSP, BOLA-DQA5, and BOLA genes from Test2 animals, and EXOC3L4, JSP, BOLA-DQA5, and BOLA genes from Test3 animals were selected for qRT-PCR (Table 4). The same RNA samples used for the microarray were used for qRT-PCR. Log2 foldchange data of qRT-PCR and microarray analysis were analyzed for validation. The Pearson correlation coefficient between the two analyses was 0.769 (p < 0.01). The differential expression of all the selected genes was considered to be validated since similar values were obtained for the qRT-PCR and microarray data analysis (Fig. 8).

Table 4.Primers used for quantitative RT-PCR.

Fig. 8.Validation of microarray data via quantitative RT-PCR. The relative expression level was normalized by the 2-∆∆CT method in terms of the beta-actin expression level relative to the control group.