Breed and adaptive response modulate bovine peripheral blood cells’ transcriptome
Journal of Animal Science and Biotechnology volume 8, Article number: 11 (2017)
Adaptive response includes a variety of physiological modifications to face changes in external or internal conditions and adapt to a new situation. The acute phase proteins (APPs) are reactants synthesized against environmental stimuli like stress, infection, inflammation.
To delineate the differences in molecular constituents of adaptive response to the environment we performed the whole-blood transcriptome analysis in Italian Holstein (IH) and Italian Simmental (IS) breeds. For this, 663 IH and IS cows from six commercial farms were clustered according to the blood level of APPs. Ten extreme individuals (five APP+ and APP- variants) from each farm were selected for the RNA-seq using the Illumina sequencing technology. Differentially expressed (DE) genes were analyzed using dynamic impact approach (DIA) and DAVID annotation clustering. Milk production data were statistically elaborated to assess the association of APP+ and APP- gene expression patterns with variations in milk parameters.
The overall de novo assembly of cDNA sequence data generated 13,665 genes expressed in bovine blood cells. Comparative genomic analysis revealed 1,152 DE genes in the comparison of all APP+ vs. all APP- variants; 531 and 217 DE genes specific for IH and IS comparison respectively. In all comparisons overexpressed genes were more represented than underexpressed ones. DAVID analysis revealed 369 DE genes across breeds, 173 and 73 DE genes in IH and IS comparison respectively. Among the most impacted pathways for both breeds were vitamin B6 metabolism, folate biosynthesis, nitrogen metabolism and linoleic acid metabolism.
Both DIA and DAVID approaches produced a high number of significantly impacted genes and pathways with a narrow connection to adaptive response in cows with high level of blood APPs. A similar variation in gene expression and impacted pathways between APP+ and APP- variants was found between two studied breeds. Such similarity was also confirmed by annotation clustering of the DE genes. However, IH breed showed higher and more differentiated impacts compared to IS breed and such particular features in the IH adaptive response could be explained by its higher metabolic activity. Variations of milk production data were significantly associated with APP+ and APP- gene expression patterns.
In the context of adaptation, stress response is an important neurobehavioral and physiological reaction and it is essential for the survival of living organisms. In response to a stressor, the body orchestrates changes in brain activity followed by the secretion of “stress mediators”, including cytokines, metabolic hormones and corticosteroids .
The body’s response during the first stage of stress is known as fight-or-flight response. It includes the activation of sympathetic nervous system and the stimulation of the production of adrenaline and noradrenaline by adrenal glands. These molecules increase the heart rate and the glycemia and modify blood distribution to supply greater levels of glucose to organs where they are needed, like brain and skeletal muscles. Shortly after, the hypothalamic-pituitary-adrenal (HPA) axis is activated and releases corticosteroids (in particular adrenal glucocorticoids). In turn, these produce a negative feedback onto immune cells and suppress further synthesis and release of cytokines, thereby protecting the host from the detrimental consequences of an overactive immune response (e.g., tissue damage, autoimmunity, septic shock) .
The long-term activation of the stress-response mechanisms may also cause irreversible damages, like cardiovascular diseases, immunosuppression, dysfunction of digestive and reproductive systems, type-II diabetes mellitus, impairment of thyroid function, weakening and loss of body lean mass [3, 4]. Such pre-pathological or pathological consequences seriously affect not only the efficiency of animal production and the quality of the product, but undoubtedly reduce animal welfare.
During the acute phase reaction (APR), the body mounts a multifactorial response trying to remove or replace damaged tissues and one of the mechanisms involved is the secretion of the so-called acute phase proteins (APPs). The concentration of some APPs increases several fold during the APR, while others, including albumin, decreases as the liver switches the production of proteins towards the synthesis of proteins required to deal with the damage [5, 6].
In ruminants, APPs are very sensitive factors that allow the early and precise detection of inflammation . The most frequently investigated proteins in cattle are: haptoglobin (Hp), serum amyloid A (SAA), fibrinogen (Fb), ceruloplasmin, α 1-antitrypsin and α 1-acid glycoprotein (α1-AGP) [5, 8–10]. It is possible that the synthesis of APPs in cattle is influenced by cortisol [11, 12], which is the key effector molecule of the HPA axis and is recognized as the physiological response to stress [13–16].
Stress response mechanisms in cattle are still not well understood and the research is complicated by individual differences in stress response . Today, the investigation of how dairy cattle adapt to intensive production is particularly important, since the animal welfare is a growing public concern and stressed animals are less efficient, producing less than predicted by their genetic potential mostly due to a higher environmental impact.
Next-generation high-throughput RNA sequencing technology (RNA-seq) is a recently-developed method for discovering, profiling, and quantifying RNA transcripts. Such approach is used to analyze the continually changing cellular transcriptome and might help identifying gene patterns involved in adaptive response. Applicability of RNA-seq for transcriptome analysis of whole blood samples was already confirmed by many research groups [18–20]. Among the most distinct advantages of RNA-seq over prior methods for mapping and quantifying the transcriptome are unbiased whole-transcriptome profiling, higher sensitivity and accurate estimation of lowly expressed transcripts in peripheral whole blood with or without globin depletion .
Up to date RNA-seq technique was highly applied for the assessment of changes in blood transcript abundance in response to stress events, pathogenic processes, and specific physiological and metabolic statuses in dairy cattle [19–22]. However, no study has comprehensively evaluated the adaptive response on molecular changes in dairy cattle whole blood cell transcriptome as an indicator of immune activity without the visible environmental perturbations.
In this context, we used a whole-transcriptome analysis to understand if and how differential gene expression contributes to such a complex phenomenon as adaptive response. Therefore, in the present research, the transcriptome of blood cells was analyzed in selected bovines belonging to Italian Holstein (IH) and Italian Simmental (IS) breeds from six commercial farms in Friuli-Venezia Giulia region, Italy. Cows were clustered for blood APPs, plasma Zn, milk cortisol and somatic cell count (SCC) in milk. The analysis included RNA isolation from blood , sequencing by RNA-seq with Illumina pipeline [24–27] and the use of the normalized data for the identification of genes expression of which is significantly associated to the adaptive response to the undefined stress conditions. Genes and metabolic pathways were further analyzed using the dynamic impact approach (DIA) and DAVID online software tool .
Animals and management
A total of 663 IH and IS cows from six commercial farms in Friuli-Venezia Giulia region of Italy were included in the experiment. All animals were kept under the same feeding and management conditions and were in the stage of lactation. Farm veterinary practitioner confirmed that all sampled animals passed the preliminary veterinary checkup, were clinically healthy and were not under any treatments for at least 1 month before the collection day. Composition of herds and characteristics of animals included in the analysis are reported in the Tables 1 and 2.
Farmers and farm veterinary practitioners gave an oral informed consent to the study and had a copy of all the data obtained from the laboratory analyses. All farms involved in the present study adhere to a high standard of veterinary care based on best practice manual under the supervision of the official veterinary service. Sample collection was approved by the Bioethics Committee of the University of Udine.
Cows were housed in a free stall barn with cubicle design and automated milking parlour. They were milked twice a day, at approximately 12 h interval. Cows had free access to water and were fed ad libitum twice a day a total mixed ration (TMR) based on corn silage and formulated to cover nutrient requirements . TMR was administered after each milking. To ensure that no dietary variations occurred during the time window of the study, the ration formulation and the offered amount were recorded using registrations of the TMR mixed feeder. In Table 3 are summarized composition, chemical properties and nutritional values of the diet in the six commercial farms.
Milk and blood sampling and assays
Milk was sampled on the day of the official record. Coccygeal vein blood samples were collected just before the morning milking and prior the feeding process. The same collection protocol was used across all farms. Blood was collected in PAXgene Blood RNA Tubes (PreAnalytiX GmbH, Switzerland), frozen 4 h after the collection and stored at -80 °C until the RNA isolation.
Prior to RNA isolation, blood samples were thawed at +4 °C for at least 12 h. RNA was isolated according to PAXgene Blood RNA Kit (PreAnalytiX GmbH, Switzerland) protocol.
Blood biochemical parameters, i.e., total protein, albumin, urea, glucose, creatinine, total bilirubin, cholesterol, AST/GOT (aspartate transaminase/glutamic oxaloacetic transaminase), gGT (gamma-glutamyl transpeptidase), zinc, ceruloplasmin, haptoglobin and paraoxonase were assayed as described in Sgorlon et al. . Milk composition data, i.e., fat and total protein percentage, casein, urea, SCC and milk cortisol, were obtained from milk samples collected the same day of blood sampling and are described in Sgorlon et al. .
Clustering of animals
To identify animals differing in their adaptive response to the environment, cows were clustered according to the level of acute phase proteins and molecules (total protein, albumin, zinc, ceruloplasmin, haptoglobin, paraoxonase, milk cortisol and SCC). To control the differences in adaptive response between breeds and to correct it for the potential effect of the environment the clustering was performed separately for each farm. For this the principal component analysis (PCA) using the correlation matrix in SPSS package was applied. According to the first two principal components the ten extreme individuals (five “plus” [APP+] and five “minus” [APP-] variants) from each farm were selected for gene expression analysis.
RNA quality control and sequencing
First RNA was quantified and quality controlled by NanoDrop ND-1000 Spectrophotometer analysis (Thermo Fisher Scientific Inc., United States). Further RNAs with the highest quality was assigned the RNA integrity number (RIN) score by the Agilent 2100 Bioanalyzer (Agilent Technologies, United States) . Finally samples with RIN ≥ 7 were selected for sequencing .
Previously was reported that high-throughput sequencing by RNA-Seq is highly reproducible within a large dynamic range of detection and provides an accurate estimation of RNA concentration in peripheral whole blood . Thus, the experimental globin depletion from RNA samples was avoided as it could significantly reduce the amount and quality of isolated RNA and biological samples in our occasion were not possible to replenish.
The 60RNA samples (30 APP+ and 30 APP- variants) were sequenced by RNA-seq technology with the Illumina pipeline [24–27]. Reads obtained from the sequencing were aligned against Bos taurus UMD 3.1 reference genome assembly .
Raw counts produced by RNA-seq were normalized with the DeSeq2 software [27, 35]. To identify differentially expressed genes APP+ and APP- cows were compared either ignoring or considering their breed of origin: i) all APP+ vs. all APP-; ii) IH APP+ vs. IH APP-; iii) IS APP+ vs. IS APP-.
For each comparison, normalized RNA-seq data were analyzed with DeSeq2 software to calculate differential expression values (as log2 of the fold change) and raw P-values. To identify the significant genes raw P-values were corrected with the false discovery rate (FDR) method , using the cutoff of 0.05.
Dynamic impact approach analysis
Gene expression data were also analyzed by the “dynamic impact approach” (DIA) developed by Bionaz and colleagues  for the transcriptome analysis.
DIA produces the list of the most impacted pathways integrating information coming from the dataset of the whole list of genes (regardless of their significance), differential expression values, FDR correction factor and raw p-value calculated by the DeSeq2 software. Graphically, the output is well demonstrated through two types of bars: the Impact bar indicating entity of the impact (colored in blue), and the Flux bar showing direction of the impact (red color represents the overexpression of the pathway, green color represents the under-expression of the pathway and yellow color indicates the absence in expression differences).
Significant genes were submitted to the Database for Annotation, Visualization and Integrated Discovery (DAVID) to perform a serial annotation clustering . This pipeline allowed us to form series of clusters with genes grouped according to their biological function. P-values automatically associated to each cluster were corrected by the Benjamini-Hochberg method. Clusters were considered significant if corrected P-values were lower than 0.05. Thereafter, significant genes from different clusters were grouped in a single list and further checked across comparisons to find out genes in common. The clustering procedure was applied separately for each comparison.
Statistical analysis of milk production data
To investigate differences in milk composition across farms and APP groups, a mixed model analysis of variance (ANOVA) with nested design and Fisher’s least significant difference (LSD) test was applied for each breed separately. Data were analyzed with the SPSS package using the following statistical model:
Y ijk : dependent variable;
m: general mean;
Farm i : fixed effect for the farm, with i ranging from 1 to 3;
APP(Farm) ji : nested effect for the APP group of animals, with j ranging from 1 to 2 within the i th farm;
e ijk : residual error.
Principal component analysis
The result of PCA is plotted in Fig. 1. The total variability explained by the first two components for each PCA separately was on average 50% (with a range from 43 to 54%) for the 6 commercial farms. The first component explained on average 31% of variability (range of 27-36%) and the second component explained on average 18% of variability (range of 15-25%). More accurate information about % of variability explained by PCA is reported in Table 4. Variables related to ceruloplasmin, haptoglobin and SCC were among the most important as they were highly correlated with the 1st PC within each farm (PCA). Their loadings correlation values ranged from 0.4-0.8 for each farm. Variable related to the total proteins was also very important as it was highly correlated (around 0.7) to the PC1 in the 4 out of 6 tested farms. Characteristics of groups of animals chosen for the final transcriptome analysis in the selected commercial farms are reported in Additional file 1.
Statistical analysis of milk production data
Significant differences between plus and minus APP groups were observed for milk urea in IS (P ≤ 0.001). Other parameters did not show significant differences (Table 5).
Between farms, IS cows showed significant differences in milk protein percentage (P ≤ 0.001), milk yield (P ≤ 0.01) and percentage of caseins (P ≤ 0.01). For IH animals the statistically significant differences between farms were observed only for milk fat percentage (P < 0.05).
Unlike IS animals, APP+ IH animals demonstrated a marked, even if not significant decrease in milk yield. Milk urea in APP+ animals showed a marked decrease in absolute values in both breeds; however in IS breed the decrease reached the significant level (P < 0.001).
Alignment of RNA-seq data to the UMD 3.1 bovine reference genome identified 13,665 genes expressed in bovine blood cells. A total of 1,152 significant differentially expressed genes (P < 0.05) were identified in the comparison of all APP+ vs. all APP- animals across breeds; 531 in comparison of IH APP+ vs. APP- and 217 in comparison of IS APP+ vs. APP-. The number of shared and unique transcripts within each comparison is indicated on the Venn diagram (Fig. 2).
The higher number of significant genes obtained in the global comparison among APP+ and APP- cows may be explained by the procession of data from all sampled animals, hence each gene had expression data from both IF and IS cows. This fact increased the level of significance of number of genes in the comparison of all APP+ vs. all APP- variants and showed a less robust significance within intra-breed comparisons.
This analysis also allowed to evaluate the number of over- and underexpressed genes in the list of significant DE genes. Since the expression rate was indicated as the log2 of the fold change, we assumed that genes with an expression rate greater than 1 were overexpressed and those with the expression rate lower than 1 were underexpressed. The three comparisons, despite the great diversity in the number of significant genes, showed similar ratios between over- and underexpressed genes: in each case, the number of overexpressed genes was far greater than the number of underexpressed ones (Fig. 3). In the global comparison between APP+ and APP- cows the upregulated genes were about 2-fold higher than the downregulated: 763 overexpressed and 389 underexpressed genes. In IF comparison overexpressed genes were about 3-fold higher than underexpressed ones: 396 overexpressed and 135 underexpressed genes. In IS comparison overexpressed genes were about 4-fold more numerous than underexpressed ones: 174 overexpressed and 43 underexpressed genes.
These data highlights that the adaptive response affects blood transcriptome principally by increasing the expression of a high number of genes, while the downregulation is a mechanism with much lower extent.
Dynamic impact approach
The ten most impacted KEGG pathways were identified by DIA for the comparisons with or without breed consideration (Fig. 4).
The three most impacted pathways in the comparison APP+ vs. APP- across breeds and within IH (Figs. 4a and b) are vitamin B6 metabolism, folate biosynthesis and nitrogen metabolism. In IS (Fig. 4c) these pathways are within the first five, in particular nitrogen metabolism is the second, folate biosynthesis the third and vitamin B6 metabolism the fifth most impacted pathway. Other pathways found in all three comparisons are linoleic acid metabolism, drug metabolism-other enzymes and African trypanosomiasis. Some other pathways are present in one comparison or in two out of three comparisons. Amoebiasis and dorso-ventral axis formation are present in the comparison across breeds and in IH; drug metabolism-cytochrome P450 is present in the comparison across breeds and in IS; NOD-like receptor signaling pathway is present only in the comparison across breeds; taurine and hypotaurine metabolism and caffeine metabolism pathways are present only in the IH comparison; retinol metabolism, steroid hormone biosynthesis and metabolism of xenobiotics by cytochrome P450 are present only in the IS comparison.
The complete list of all impacted pathways in each comparison is presented in Additional file 2.
Gene annotation clusters
Function terms of each cluster identified by DAVID within each comparison were predicted by Gene ontology (GO) (Table 6). After elimination of repeated gene terms, 369 significant differentially expressed genes across breeds, 173 in IH and 73 in IS remained included in significant clusters.
A total of 24 genes were differentially expressed in all 3 comparisons. Gene names and differential expression values for each of the 24 genes within each group are listed in Table 7.
The aim of the present study was to investigate the impact of stress response on gene expression patterns in peripheral blood cells of lactating cows. The analysis of transcriptome variation was performed after the peak of lactation as the transition period is the most challenging in dairy cows and can interfere with the metabolic imbalance of animals . Considering that animals on each commercial farm were kept under the same environmental conditions and that the farm factor was considered in the statistical model, the influence of management on animal adaptive response should have been minimized. Hence, the different levels of plasma APPs are likely to result from individual animal response to subclinical inflammatory/infective events or other stresses, since cows did not show visible clinical signs or symptoms of the presence of functional disorders.
Stress response is a very complex phenomenon as it can affect overall physiology through different mechanisms, like activation of sympathetic nervous system with the release of catecholamines, activation of HPA axis and non-circadian production of glucocorticoids [1, 4, 39], onset of an acute phase response . Activation of these mechanisms may cause harmful and sometimes irreversible effects on many body systems. Stress may affect circulating glucocorticoids [4, 40, 41] with consequences on female reproductive system , immune system, osteoblastogenesis and bone metabolism [43–45], muscle production , metabolism of nutrients [47–49], functioning of the thyroid gland  and growth hormone axis . Stress research is therefore complicated by these complex and diverse mechanisms and by individual and interspecies differences in stress response . Understanding the biological basis of stress response in livestock is important for improving animal welfare in intensive production systems. In addition of being a growing public concern, animal welfare is important for production efficiency and influence both farm economy and environmental footprint. Whole-transcriptome analysis is crucial to understand how stress influences gene expression to elicit the complex phenomenon of adaptive response. Here we investigated differential expression in blood samples obtained from cows with high and low levels of positive APPs as proxy of stress status and identified stress response-related genes pathways in white blood cells.
Impacted pathways by DIA
Once fed a list of differentially expressed genes, DIA exploits an online sheet of the Kyoto Encyclopedia of Genes and Genomes (KEGG) database  to detect significantly impacted pathways. It calculates the entity and direction of the impact and whether the pathway is entirely overexpressed, underexpressed or if the expression is not altered at all.
The most impacted pathways in APP+ cows in across and within breed analyses are presented in the Fig. 4. Significant genes in significant KEGG pathways have been analyzed in detail to understand gene and pathway function, since the names of KEGG pathways reported in the DIA output files are sometimes misleading.
The pathway of vitamin B6 metabolism (KEGG bta00750), was among the three most significant ones (rank 1 across breeds and in IH and rank 5 in IS). In this pathway we found two significant genes, PDXK (pyridoxal kinase) and AOX1 (aldehyde oxidase 1). The latter was not significant in IS comparison, but it shows xa similar expression pattern. Pyridoxal kinase is involved in the ATP-dependent phosphorylation of pyridoxal, pyridoxamine and pyridoxine to pyridoxal-5-phosphate (PLP), pyridoxamine-5-phosphate and pyridoxine-5-phosphate, respectively. There is a requirement for ubiquitous expression of piridoxal kinase in mammalian tissues as PLP, before entering a cell, must be dephosphorylated and after diffusing through cell membranes it is converted back to the active cofactor by cytosolic pyridoxal kinase . PLP is a very important enzymatic cofactor, as it participates to all transamination reactions and, in some cases, to decarboxylation, deamination and racemization of amino acids , catalyzes the rate-limiting step in glycogenolysis . The impact of PLP on amino acid metabolism has direct consequences also on protein synthesis, whose physiological level is altered during an adaptive response. This evidence suggests that overexpression of PDXK gene is fundamental for metabolism regulation during adaptive phenomenon, as the cofactor affects protein and energy metabolism, that are likely increased during adaptation, and this is true not only for peripheral white blood cells, but for the vast majority of tissues in an organism. Aldehyde oxidase 1 is an enzyme mainly found in liver, shows broad substrate specificity, including pyridoxal, and catalyzes the oxidation of several endogenous and exogenous aldehydes, with the production of hydrogen peroxide and superoxide ion. Aldehyde oxidase isolated from polymorphonuclear leukocytes showed a more narrow substrate specificity: for example, the enzyme found in leukocytes is inactive on xanthine . This enzyme has a role in oxidative stress and regulation of reactive oxygen species (ROS) homeostasis . The catalysis requires the presence of flavin adenine dinucleotide (FAD) and a molybdopterin cofactor (MoCo) [57–59]. The enzyme also has a role in nitric oxide (NO) biosynthesis . NO has various functions in the organism; white blood cells, mainly macrophages, secrete it as a chemical defense against bacteria and to induce vasodilation . The involvement of aldehyde oxidase 1 in several oxidative metabolisms, in oxidative stress and in NO signaling may explain the overexpression of AOX1 gene in APP+ cows.
Interesting significant genes emerged from the analysis of folate biosynthesis pathway (KEGG bta00790; rank 2 across breeds and in IH, rank 3 in IS). These are ALPL (alkaline phosphatase liver/bone/kidney) and MOCS1 (molybdenum cofactor synthesis 1). The ALPL gene is mainly expressed in neutrophils and monocytes . The role of alkaline phosphatase is fundamental as a high number of signal transduction cascades are involved in adaptive processes. Phosphorylation and dephosphorylation of signal proteins is the key determinant in the phenomenon that regulates the transduction and the amplification of a stimulus from the cell membrane receptor to the nucleus, where the modulation of gene expression occurs. Elevations in plasma alkaline phosphatase, whose sources include neutrophils and monocytes, can be also related to pathological conditions . MOCS1 encodes for a protein involved in the biological activation of molybdenum and it is highly expressed by peripheral white blood cells. Participating in the production of MoCo, it is indirectly involved in the catalytic activity of several enzymes, including aldehyde oxidase and xanthine oxidase . MOCS1 was significant only in APP+ vs. APP- comparison across breeds and the involvement of the gene in a number of metabolic oxidative pathways is likely the reason for its significant overexpression in the APP+ bovines.
In nitrogen metabolism pathway (KEGG bta00910; rank 3 across breed and in IH, rank 2 in IS), GLUL (glutamate-ammonia ligase) and CA4 (carbonic anhydrase IV) genes were significantly differentially expressed. Glutamate-ammonia ligase is a PLP-dependent enzyme that produces glutamine from glutamate and free NH3. Glutamine is a common metabolite in many amino acid, purine and pyrimidine biosynthetic pathways, so this enzyme has a major role in protein and nucleic acid metabolism. It is also involved in acid-base homeostasis, cell signaling, cell proliferation and biosynthesis of γ-aminobutyric acid (GABA) . Carbonic anhydrase IV is an important Zn-dependent enzyme present in several tissues and, among leukocytes, it is expressed principally by eosinophils and neutrophils. It has a main role in the control of acid-base balance in blood and other tissues . Particularly, this isoform exists in the form of a glycophosphatidylinositol (GPI)-anchored protein and plays an important role in maintaining an appropriate cellular environment for the reactions that occur during adaptive responses. Both these genes are overexpressed in APP+ cows.
The pathway for linoleic acid metabolism (KEGG bta00591; rank 4 across breeds, rank 6 in IH, rank 1 in IS) includes a number of significant genes involved in inflammatory response and metabolism of drugs and xenobiotics. The phospholipase A2 genes, PLA2G4F and PLA2G4A, selectively hydrolyze membrane phospholipids. The first one has high selectivity for phosphatidylethanolamine, hydrolyzing the ester bond in sn-2 position, and has a role in mitogen-associated protein kinase (MAPK) and Ras signaling pathways . The latter pathway leads to the production of free arachidonic acid, which is further converted in eicosanoids, involved in inflammatory response, and lysophospholipids, that are precursors of platelet-activating factor (PAF). Hence, this enzyme has a role in inflammatory response and hemodynamics, and is also involved in MAPK and G protein-coupled receptor (GPCR) signaling pathways . Some significant genes in this pathway belong to cytochrome P450 superfamily (in detail: CYP2E1, CYP3A4, CYP3A5). These genes are also involved in metabolism of steroid hormones, drugs and carcinogens, playing a role in steroid-mediated physiological responses, activation and metabolic drug clearance and carcinogenesis [68–70].
In the pathway of drug-metabolizing enzymes (KEGG bta00983; rank 5 across breeds, rank 4 in IH and IS) there are several significant genes with a role in the metabolism of nucleotides, suggesting their role in an adaptive response. Xanthine dehydrogenase (XDH) is a paralog of AOX1, which can operate either as a dehydrogenase or as an oxidase. Xanthine dehydrogenase is involved in metabolism of hypoxanthine and xanthine and in the generation of ROS . Recently, its role in recruiting macrophages through inflammasome activation has been investigated . Cytidine deaminase (CDA) preserves pyrimidine pool by irreversibly deaminating cytidine and deoxycytidine to uridine and deoxyuridine, respectively. It is also involved in antibody diversification . Uridine phosphorylase (UPP1) reversibly cleaves ribose-1-phosphate and deoxyribose-1-phosphate from uridine and deoxyuridine, releasing free uracil . Another significant enzyme involved in pyrimidine metabolism is dihydropyrimidine dehydrogenase (DPYD) .
In African trypanosomiasis (KEGG bta05143; rank 6 across breeds, rank 5 in IH, rank 9 in IS) and amoebiasis (KEGG bta05146; rank 7 across breeds, rank 9 in IH) pathways we found a large number of significant genes directly involved in the onset of a stress response or inflammation. Among these, we found several proinflammatory cytokine genes, as IL12B, IL18, IL1B and TNFα. These cytokines are produced by different types of immune cells involved in growth, differentiation, chemotaxis and proliferation of white cells during an inflammatory event [76–80]. NFkB1 and RELA, the nuclear transcription factor genes forming the same protein complex, were also among affected ones. NFkB1 is activated by cytokines, free radicals, UV ray and pathogens’ products and it is involved in regulation of inflammation-mediated pathways  and in regulation of the expression of number of genes involved in cell adhesion and migration across vascular endothelium, like vascular cell adhesion molecule (VCAM1), laminin genes (LAMA4, LAMC1), fibronectin (FN1) and integrin genes (ITGB2, ITGAM). Some of the latter genes are strongly induced by cytokine signaling, so they have a primary importance in leukocyte adhesion and in cell signal transduction [82–86]. Significant genes in the African trypanosomiasis and amoebiasis pathways include also a number of cell surface receptors involved in inflammation signaling cascades, regulation of cell physiology during these events and in functionality of activated leukocytes. Among these genes there are Fas cell surface death receptor (FAS), toll-like receptors (TLR2, TLR4), CD14 molecule and complement protein genes (C8A) [87–90]. These receptors transducer signals through different molecules, including myeloid differentiation factor (MYD88) . Moreover, some significant genes in these pathways include enzymes that regulate the production of second messengers, like phosphatidylinositide-3-kinases (PIK3R2, PIK3CD)  and phospholipases C (PLCB3, PLCB4) , and proteases that regulate cell protein homeostasis, limiting tissue damage produced by overexpressed proteolytic enzymes, like serpin peptidases (SERPINB1, SERPINB3) . In IS breed amoebiasis pathway was found to be not significantly impacted between APP+ and APP- animals.
NOD-like receptor signaling pathway (KEGG bta04621) was found to be significant only in the APP+ vs. APP- across breed comparison (rank 8). Significant genes belonging to the pathway are mostly involved in regulation of cell cycle and pro-apoptotic signaling, interacting with the nuclear factor NFkB1 to activate it. Pro-apoptotic genes include the nucleotide-binding oligomerization domain containing receptors (NOD1, NOD2)  and genes of the caspase recruitment domain family (CARD6, PYCARD, NLRP3), which also participate to the formation of inflammasomes [96–98]. In this pathway there are also two significantly overexpressed molecular chaperones, HSP90AA1 and HSP90B1.
Dorso-ventral axis formation (KEGG bta04320) pathway was significantly impacted in APP+ vs. APP- across breeds (rank 9) and in IH (rank 8). However, some genes belonging to this pathway were significant also in IS breed. In this pathway there are two significant transcription factors, ETV6 and ETS2, which are oncogenes with a role in hematopoiesis and apoptosis [99, 100].
The cytochrome P450-mediated metabolism of drugs (KEGG bta00982) is among the significantly impacted pathways across breeds (rank 10) and in IS (rank 7). Important genes involved in oxidative stress and detoxification of oxidation by-products are included in this pathway, such as membrane-bound microsomal glutathione S-transferase (MGST1), with a role in the development of inflammation and in cellular defense , aldehyde dehydrogenase (ALDH3B1)  and monoamine oxidase A (MAOA) . Monoamine oxidase A has a role in the metabolism of serotonin: this molecule has been shown to be synthesized, released and degraded also by T lymphocytes ). Among the significant genes, we found also some terms significantly overexpressed in other pathways, as AOX1 and CYP2E1.
Metabolism of taurine and hypotaurine (KEGG bta00430) was strongly impacted only in IH (rank 7) and the only significant gene detected was a member of γ-glutamyltransferase family (GGT5). This gene is involved in metabolism of glutathione and leukotrienes and plays a role in oxidative stress and inflammatory response .
The pathway of caffeine metabolism (KEGG bta00232) showed a high impact only in IH (rank 10). It includes the gene XDH, the function of which has been previously discussed.
Retinol metabolism (KEGG bta00830) was highly impacted only in IS (rank 6). Important differentially expressed genes in this pathway are AOX1, CYP3A4 and CYP3A5 and their role in adaptive response was described previously.
Another pathway impacted only in IS included steroid hormone biosynthesis pathway (KEGG bta00140; rank 8). It included CYP3A4, CYP3A5 and CYP2E1 genes, involved in the biosynthesis of different types of steroids, including corticosteroids, which have a relevant role in adaptive response.
The same scenario occurred also in the pathway relative to metabolism of xenobiotics mediated by cytochrome P450 (KEGG bta00980; rank 10 in IS). It included genes MGST1, CYP2E1 and ALDH3B1, which are significant also in other pathways.
The same significant gene plays different roles in different pathways, for example in the metabolism of different compounds or in the signaling of a number of signal transduction cascades. Thus, to obtain a comprehensive analysis and to confirm the significance of the relevant genes, the most impacted pathways were explored by DAVID.
DAVID annotation clusters
Annotation clusters produced by DAVID online tool confirmed the significance of a large number of genes involved in adaptive response. Genes were clustered according to their common structural characteristics or molecular function. Table 6 summarizes the significant clusters, the number of total genes per cluster and the number of genes in common with at least another cluster. In Table 7 are listed common genes between all three comparisons. A comprehensive list of genes included in each cluster is available in Additional file 3.
The highest number of common genes among clusters was included in C5, CH4 and CS8, which are clusters of glycoproteins, based on the GO terms. Twenty-three out of 24 shared genes are included in these clusters. The only one excluded is the GNAL gene. We can assume that genes present in these clusters are grouped only according to structural features, as they have various molecular functions, but are all classified as glycoproteins. GNAL is included in the cluster that reports purine nucleotide-interacting proteins.
In these clusters we can find some genes that were significant also in DIA output, i.e., CA4, ALPL and IL10. Other shared genes were not included in the most impacted pathways, but they have an important role in the regulation of adaptive response. For example, the receptor of insulin-like growing factor II (IGF2R) is overexpressed in all three comparisons, possibly because of its role in intracellular trafficking of lysosomal enzymes and degradation of IGF-II .
Angiotensin I-converting enzyme 2 (ACE2) has a direct effect on cardiac and renal functions . Its overexpression may be another consequence of hyperactivation of HPA axis, as glucocorticoids have a direct effect on blood pressure and cardiovascular system, and of an increased uremia, as this condition may trigger the recruitment of pro-atherogenic, ACE2-expressing monocytes .
As most of cytokines, present in all three comparisons, interleukin 10 (IL10) is significantly overexpressed. Such overexpression can be explained by its important role in lymphocytes differentiation and proliferation and in the production of antibodies by B cells. In the same list of overexpressed genes are present receptors for cytokines, like interleukin 2 receptor alpha (IL2RA), involved in the intracellular signaling pathways.
The only one underexpressed cytokine revealed among three comparisons is the pro inflammatory cytokine IL34.
Defensins are small antimicrobial, cytotoxic peptides produced by neutrophils . β-defensins (DEFB7 and DEFB10 in Table 7), belonging to one of the three existing groups of defensins, were overexpressed in all comparisons. Their overexpression could be associated to the presence of a bacterial infection in mammary gland  and this result is in agreement with the highly impacted pathway relative to S. aureus infection (see Additional file 2).
Another interesting overexpressed gene is the gene encoding matrix metalloproteinase 9 (MMP9). Its overexpression can be associated with an augmented production of hematopoietic stem cells  that give rise to all red and white blood cells and platelets. MMP9 has recently been investigated for its role in promoting the secretion of pro-inflammatory cytokines and the migration of T cells towards inflammation sites. Moreover, the protein increases the permeability of brain-blood barrier to cytokines, showing an important involvement in neuroinflammation .
Summarizing the results, a lot of highly impacted pathways were found in common between IH and IS breeds, indicating a similar variation in the gene expression under environmental adaptation. This observation was also confirmed by annotation clustering of the significantly expressed genes: in fact, some of the significantly clustered genes are present also in the most impacted pathways, with similar levels of expression.
Differences between breeds were observed only on the level of individual genes. Therefore, considering the overall variation in stress-inducing factors and stress-related gene expression, the general patterns can be considered very similar in the two breeds investigated.
The obtained data can be considered reliable even if the further validation analysis of the differential gene expression in the population could improve the relevance of the conclusions. This point can be developed in a future research.
Analysis of milk production data
The differences in milk production data between APP- and APP+ cows are reported in Table 5. Significant association with high APP level was observed only for milk urea in IS, which showed a marked decrease in APP+ animals. A marked decrease, close to being statistically significant, was observed also in IH milk yield and milk urea.
The decrease in milk yield can be associated to the onset of a stress response as energy resources are driven towards other organs and other more important physiological processes to guarantee the animal’s survival. The decrease in milk urea can be linked to the overall increased protein synthesis. In dairy cows milk urea reflects the catabolism of protein by the ruminant tissues and within the rumen by bacteria. The decrease of rumen ammonia may indicate that animals increase the protein synthesis during an adaptive response . Indeed, high levels of APPs and at the same time high expression of gene products of several pathways were observed in APP+ animals.
This transcriptomic study in lactating dairy cows from IH and IS breeds allowed to assess that the onset of an adaptive response to the environment involves a large number of pathways that are regulated in a similar way in the two breeds, with marginal differences in significantly over- or underexpressed single genes. The altered expression of these genes is statistically associated with the variations of only certain milk production parameters between groups of cows with activated and not activated stress response mechanisms.
Considering all the results we can conclude that the two studied breeds have similar patterns but vary in the degree of activation of metabolic and physiological mechanisms of adaptation to the environment. IH showed a higher rate of significant genes and impacted pathways, demonstrating a higher metabolic activity, respect to the IS breed.
Acute phase reaction
Acute phase protein
Dynamic impact approach
Principal component analysis
RNA integrity number
Somatic cell count
Kinlein SA, Wilson CD, Karatsoreos IN. Dysregulated hypothalamic-pituitary-adrenal axis function contributes to altered endocrine and neurobehavioral responses to acute stress. Front Psychiatry. 2015;6:31.
Silverman MN, Pearce BD, Biron CA, Miller AH. Immune modulation of the hypothalamic-pituitary-adrenal (HPA) axis during viral infection. Viral Immunol. 2005;18(1):41–78.
Selye HA. Syndrome produced by diverse nocuous agents. Nature. 1936;138:132.
Elenkov IJ, Webster EL, Torpy DJ, Chrousos GP. Stress, corticotropin-releasing hormone, glucocorticoids, and the immune/inflammatory response: acute and chronic effects. Ann NY Acad Sci. 1999;876:1–11.
Eckersall PD, Conner JG. Bovine and canine acute phase proteins. Vet Res Commun. 1988;12(2-3):169–78.
Haeryfar SM, Berczi I. The thymus and the acute phase response. Cell Mol Biol (Noisy-le-grand). 2001;47(1):145–56.
Kent J. Acute phase proteins: their use in veterinary diagnosis. Br Vet J. 1992;148(4):279–82.
Conner JG, Eckersall PD, Doherty M, Douglas TA. Acute phase response and mastitis in the cow. Res Vet Sci. 1986;41(1):126–8.
Horadagoda A, Eckersall PD, Alsemgeest SP, Gibbs HA. Purification and quantitative measurement of bovine serum amyloid-A. Res Vet Sci. 1993;55(3):317–25.
Dowling A, Hodgson JC, Schock A, Donachie W, Eckersall PD, McKendrick IJ. Experimental induction of pneumonic pasteurellosis in calves by intratracheal infection with Pasteurella multocida biotype A:3. Res Vet Sci. 2002;73(1):37–44.
Alsemgeest SP, Taverne MA, Boosman R, van der Weyden BC, Gruys E. Peripartum acute-phase protein serum amyloid-A concentration in plasma of cows and fetuses. Am J Vet Res. 1993;54(1):164–7.
Jawor P, Stefaniak T. Acute phase proteins in cattle. In: Veas F, editor. Acute phase proteins as early non-specific biomarkers of human and veterinary diseases. Rijeka: InTech d.o.o; 2011. p. 381–408.
Amadori M, Stefanon B, Sgorlon S, Farinacci M. Immune system response to stress factors. Ital J Anim Sci. 2009;8(1):287–99.
Cray C, Zaias J, Altman NH. Acute phase response in animals: a review. Comp Med. 2009;59(6):517–26.
Ranabir S, Reetu K. Stress and hormones. Indian J Endocrinol Metab. 2011;15(1):18–22.
DairyCare COST Action FA1308. 2014. http://www.dairycareaction.org. Accessed 21 June 2016
Hing S, Narayan E, Thompson RCA, Godfrey S. A review of factors influencing the stress response in Australian marsupials. Conserv Physiol. 2014;2(1):cou027.
Solano-Aguilar G, Molokin A, Botelho C, Fiorino A-M, Vinyard B, Li R, et al. Transcriptomic Profile of Whole Blood Cells from Elderly Subjects Fed Probiotic Bacteria Lactobacillus rhamnosus GG ATCC 53103 (LGG) in a Phase I Open Label Study. PLoS One. 2016;11:2.
McLoughlin KE, Nalpas NC, Rue-Albrecht K, Browne JA, Magee DA, Killick KE, et al. RNA-seq transcriptional profiling of peripheral blood leukocytes from cattle infected with mycobacterium bovis. Front Immunol. 2014;5:396.
Shin H, Shannon CP, Fishbane N, Ruan J, Zhou M, Balshaw R, et al. Variation in RNA-Seq transcriptome profiles of peripheral whole blood from healthy individuals with and without globin depletion. PLoS One. 2014;9:3.
Weikard R, Demasius W, Hadlich F, Kühn C. Different blood cell-derived transcriptome signatures in cows exposed to vaccination pre- or postpartum. PLoS One. 2015;10:8.
O’Loughlin A, Lynn DJ, McGee M, Doyle S, McCabe M, Earley B. Transcriptomic analysis of the stress response to weaning at housing in bovine leukocytes using RNA-seq technology. BMC Genomics. 2012;13:250.
Chai V, Vassilakos A, Lee Y, Wright JA, Young AH. Optimization of the PAXgene blood RNA extraction system for gene expression analysis of clinical samples. J Clin Lab Anal. 2005;19(5):182–8.
Mardis ER. Next-generation DNA, sequencing methods. Annu Rev Genomics Hum Genet. 2008;9:387–402.
Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10(1):57–63.
Wickramasinghe S, Cánovas A, Rincón G, Medrano JF. RNA-sequencing: a tool to explore new frontiers in animal genetics. Livest Sci. 2014;166:206–16.
Han Y, Gao S, Muegge K, Zhang W, Zhou B. Advanced applications of RNA sequencing and challenges. Bioinform Biol Insights. 2015;9 Suppl 1:29–46.
Bionaz M, Periasamy K, Rodriguez-Zas SL, Hurley WL, Loor JJ. A novel dynamic impact approach (DIA) for functional analysis of time-course omics studies: validation using the bovine mammary transcriptome. PloS One. 2012;7(3):e32455.
Institut National de la Recherche Agronomique. Introduction. Feeding standards for ruminants. In: Jarrige R, editor. Ruminant Nutrition. Recommended Allowances and Feed Tables. London: Eurotext; 1989. p. 15–22.
Sgorlon S, Fanzago M, Sandri M, Gaspardo B, Stefanon B. Association of index of welfare and metabolism with the genetic merit of Holstein and Simmental cows after the peak of lactation. Ital J Anim Sci. 2015;14(3):368–73.
Sgorlon S, Fanzago M, Guiatti D, Gabai G, Stradaioli G, Stefanon B. Factors affecting milk cortisol in mid lactating dairy cows. BMC Vet Res. 2015;11:259.
Schroeder A, Mueller O, Stocker S, Salowski R, Leiber M, Gassmann M, et al. The RIN: an RNA integrity number for assigning integrity values to RNA measurements. BMC Mol Biol. 2006;7:3.
Sgorlon S, Colitti M, Asquini E, Ferrarini A, Pallavicini A, Stefanon B. Administration of botanicals with the diet regulates gene expression in peripheral blood cells of Sarda sheep during ACTH challenge. Domest Anim Endocrinol. 2012;43(3):213–26.
Elsik GC, Unni DR, Diesh CM, Tayal A, Emery ML, Nguyen HN, et al. Bovine Genome Database: new tools for gleaning function from the Bos taurus genome. Nucleic Acids Res. 2015;44(D1):D834–9.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B Stat Methodol. 1995;57(1):289–300.
Huang DW, Sherman BT, Tan Q, Collins JR, Alvord WG, Roayaei J, et al. The DAVID gene functional classification tool: a novel biological module-centric algorithm to functionally analyze large gene lists. Genome Biol. 2007;8(9):R183.
Sundrum A. Metabolic disorders in the transition period indicate that the dairy cows’ ability to adapt is overstressed. Animals (Basel). 2015;5(4):978–1020.
Papadimitriou A, Priftis KN. Regulation of the hypothalamic-pituitary-adrenal axis. Neuroimmunomodulation. 2009;16(5):265–71.
Buckingham JC, Loxley HD, Christian HC, Philip JG. Activation of HPA axis by immune insults: roles and interactions of cytokines, eicosanoids, glucocorticoids. Pharmacol Biochem Behav. 1996;54(1):285–98.
Kitajima T, Ariizumi K, Bergstresser PR, Takashima A. A novel mechanism of glucocorticoid-induced immune suppression: the inhibition of T cell-mediated terminal maturation of a murine dendritic cell line. J Clin Invest. 1996;98(1):142–7.
Chrousos GP, Torpy DJ, Gold PW. Interactions between the hypothalamic-pituitary-adrenal axis and the female reproductive system: clinical implications. Ann Intern Med. 1998;129(3):229–40.
Weinstein RS, Jilka RL, Parfitt AM, Manolagas SC. Inhibition of osteoblastogenesis and promotion of apoptosis of osteoblasts and osteocytes by glucocorticoids. Potential mechanism of their deleterious effects on bone. J Clin Invest. 1998;102(2):274–82.
Doga M, Bonadonna S, Giustina A. Glucocorticoids and bone: cellular, metabolic and endocrine effects. Hormones (Athens). 2004;3(3):184–90.
Tamura Y, Okinaga H, Takami H. Glucocorticoid-induced osteoporosis. Biomed Pharmacother. 2004;58(9):500–4.
Schakman O, Gilson H, Kalista S, Thissen JP. Mechanism of muscle atrophy induced by glucocorticoids. Horm Res. 2009;72(1):36–41.
Tsigos C, Young RJ, White A. Diabetic neuropathy is associated with increased activity of the hypothalamic-pituitary-adrenal axis. J Clin Endocrinol Metab. 1993;76(3):554–8.
Anagnostis P, Athyros VG, Tziomalos K, Karagiannis A, Mikhailidis DP. Clinical review: The pathogenetic role of cortisol in the metabolic syndrome: a hypothesis. J Clin Endocrinol Metab. 2009;94(8):2692–701.
Rose AJ, Vegiopoulos A, Herzig S. Role of glucocorticoids and the glucocorticoid receptor in metabolism: insights from genetic manipulations. J Steroid Biochem Mol Biol. 2010;122(1-3):10–20.
Joseph-Bravo P, Jaimes-Hoy L, Charli JL. Regulation of TRH neurons and energy homeostasis-related signals under stress. J Endocrinol. 2015;224(3):R139–59.
Nicolaides NC, Kyratzi E, Lamprokostopoulou A, Chrousos GP, Charmandari E. Stress, the stress system and the role of glucocorticoids. Neuroimmunomodulation. 2015;22(1-2):6–19.
Kaneisha M, Goto S. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 2000;28(1):27–30.
Hanna MC, Turner AJ, Kirkness EF. Human pyridoxal kinase. cDNA cloning, expression, and modulation by ligands of the benzodiazepine receptor. J Biol Chem. 1997;272(16):10756–60.
Toney MD. Reaction specificity in pyridoxal phosphate enzymes. Arch Biochem Biophys. 2005;433(1):279–87.
Livanova NB, Chebotareva NA, Eronina TB, Kurganov BI. Pyridoxal 5’-phosphate as a catalytic and conformational cofactor of muscle glycogen phosphorylase b. Biochemistry Mosc. 2002;67(10):1089–98.
Beedham C. Similarity in the aldehyde oxidases from guinea-pig liver and polymorphonuclear leucocytes. J Pharm Pharmacol. 1986;38:57–8.
Kundu TK, Hille R, Velayutham M, Zweier JL. Characterization of superoxide production from aldehyde oxidase: an important source of oxidants in biological tissues. Arch Biochem Biophys. 2007;460(1):113–21.
Sahi J, Khan KK, Black CB. Aldehyde oxidase activity and inhibition in hepatocytes and cytosolic fractions from mouse, rat, monkey and human. Drug Metab Lett. 2008;2(3):176–83.
Garattini E, Fratelli M, Terao M. The mammalian aldehyde oxidase gene family. Hum Genomics. 2009;4(2):119–30.
Li H, Kundu TK, Zweier JL. Characterization of the magnitude and mechanism of aldehyde oxidase-mediated nitric oxide production from nitrite. J Biol Chem. 2009;284(49):33850–8.
Lundberg JO, Weitzberg E, Gladwin MT. The nitrate-nitrite-nitric oxide pathway in physiology and therapeutics. Nat Rev Drug Discov. 2008;7:156–67.
Izumi M, Ishikawa J, Takeshita A, Maekawa M. Increased serum alkaline phosphatase activity originating from neutrophilic leukocytes. Clin Chem. 2005;51(9):1751–2.
Hänzelmann P, Hernández HL, Menzel C, García-Serres R, Huynh BH, Johnson MK, et al. Characterization of MOCS1A, an oxygen-sensitive iron-sulfur protein involved in human molybdenum cofactor biosynthesis. J Biol Chem. 2004;279(33):34721–32.
Labow BI, Souba WW, Abcouwer SF. Mechanisms governing the expression of the enzymes of glutamine metabolism – glutaminase and glutamine synthetase. J Nutr. 2001;131Suppl 9:2467S–74S. discussion 2486S-7S.
Wen T, Mingler MK, Wahl B, Khorki ME, Pabst O, Zimmermann N, et al. Carbonic anhydrase IV is expressed on IL-5-activated murine eosinophils. J Immunol. 2014;192(12):5481–9.
Ohto T, Uozumi N, Hirabashi T, Shimizu T. Identification of novel cytosolic phospholipase A2s, murine cPLA2δ, ε, and ζ, which form a gene cluster with cPLA2β. J Biol Chem. 2005;280(26):247576–83.
Burke JE, Dennis EA. Phospholipase A2 structure/function, mechanism, and signaling. J Lipid Res. 2009;50(Suppl):S237–42.
Schulz-Utermoehl T, Mountfield RJ, Bywater RP, Madsen K, Jørgensen PN, Hansen KT. Structure-function analysis of human CYP3A4 using a specific proinhibitory antipeptide antibody. Drug Metab Dispos. 2000;28(7):718–25.
Patki KC, Von Moltke LL, Greenblatt DJ. In vitro metabolism of midazolam, triazolam, nifedipine, and testosterone by human liver microsomes and recombinant cytochromes p450: role of cyp3a4 and cyp3a5. Drug Metab Dispos. 2003;31(7):938–44.
Hanioka N, Yamamoto M, Iwabu H, Jinno H, Tanaka-Kagawa T, Naito S, et al. Functional characterization of human and cynomolgus monkey cytochrome P450 2E1 enzymes. Life Sci. 2007;81(19-20):1436–45.
Harrison R. Structure and function of xanthine oxidoreductase: where are we now? Free Radic Biol Med. 2002;33(6):774–97.
Ives A, Nomura J, Martinon F, Roger T, LeRoy D, Miner JN, et al. Xanthine oxidoreductase regulates macrophage IL1β secretion upon NLRP3 inflammasome activation. Nat Commun. 2015;6:6555.
Chung SJ, Fromme JC, Verdine GL. Structure of human cytidine deaminase bound to a potent inhibitor. J Med Chem. 2005;48(3):658–60.
Temmink OH, de Bruin M, Laan AC, Turksma AW, Cricca S, Masterson AJ, et al. The role of thymidine phosphorylase and uridine phosphorylase in (fluoro) pyrimidine metabolism in peripheral blood mononuclear cells. Int J Biochem Cell Biol. 2006;38(10):1759–65.
Van Kuilenburg AB, Meinsma R, Beke E, Bobba B, Boffi P, Enns GM, et al. Identification of three novel mutations in the dihydropyrimidine dehydrogenase gene associated with altered pre-mRNA splicing or protein function. Biol Chem. 2005;386(4):319–24.
Boraschi D, Tagliabue A. Human interleukin 1: structure-function relationship. Ann Ist Super Sanita. 1990;26(3-4):273–82.
Oshimi K, Hoshino S. Function and molecular structure of IL-12. Nihon Rinsho. 1992;50(8):1840–4.
Burdin N, Rousset F, Banchereau J. B-cell-derived IL-10: production and function. Methods. 1997;11(1):98–111.
Biet F, Locht C, Kremer L. Immunoregulatory functions of interleukin 18 and its role in defense against bacterial pathogens. J Mol Med. 2002;80(3):147–62.
Mukai Y, Shibata H, Nakamura T, Yoshioka Y, Abe Y, Nomura T, et al. Structure-function relationship of tumor necrosis factor (TNF) and its receptor interaction based on 3D structural analysis of a fully active TNFR1-selective TNF mutant. J Mol Biol. 2009;385(4):1221–9.
Pereira SG, Oakley F. Nuclear factor-kappaB1: regulation and function. Int J Biochem Cell Biol. 2008;40(8):1425–30.
Ponce ML, Nomizu M, Delgado MC, Kuratomi Y, Hoffmann MP, Powell S, et al. Identification of endothelial cell binding sites on the laminin γ1 chain. Circ Res. 1999;84(6):688–94.
Dib K, Andersson T. β2 integrin signaling in leukocytes. Front Biosci. 2000;5:D438–51.
Cook-Mills JM. VCAM-1 signals during lymphocyte migration: role of reactive oxygen species. Mol Immunol. 2002;39(9):499–508.
Sottile J, Hocking DC. Fibronectin polymerization regulates the composition and stability of extracellular matrix fibrils and cell-matrix adhesions. Mol Biol Cell. 2002;13(10):3546–59.
DeHahn KC, Gonzales M, Gonzalez AM, Hopkinson SB, Chandel NS, Brunelle JK, et al. The α4 laminin subunit regulates endothelial cell survival. Exp Cell Res. 2004;294(1):281–9.
Kawakami A, Eguchi K, Matsuoka N, Tsuboi M, Koji T, Urayama S, et al. Expression and function of Fas and Fas ligand on peripheral blood lymphocytes in normal subjects. J Lab Clin Med. 1998;132(5):404–13.
Viriyakosol S, Mathison JC, Tobias PS, Kirkland TN. Structure-function analysis of CD14 as a soluble receptor for lipopolysaccharide. J Biol Chem. 2000;275(5):3144–9.
Sabroe I, Dower SK, Whyte MK. The role of Toll-like receptors in the regulation of neutrophil migration, activation, and apoptosis. Clin Infect Dis. 2005;41 Suppl 7:S421–6.
Hadders MA, Beringer DX, Gros P. Structure of C8α-MACPFF reveals mechanism of membrane attack in complement immune defense. Science. 2007;317(5844):1552–4.
Qiu Y, Shen Y, Li X, Ding C, Ma Z. Molecular cloning and functional characterization of a novel isoform of chicken myeloid differentiation factor 88 (MyD88). Dev Comp Immunol. 2008;32(12):1522–30.
Koyasu S. The role of PI3K in immune cells. Nat Immunol. 2003;4(4):313–9.
Ting AE, Pagano RE. Detection of a phosphatidylinositol-specific phospholipase C at the surface of Swiss 3 T3 cells and its potential role in the regulation of cell growth. J Biol Chem. 1990;265(10):5337–40.
Christensen S, Sottrup-Jensen L. Characterization of two serpins from bovine plasma and milk. Biochem J. 1994;303 Pt.2:383–90.
Ekman AK, Cardell LO. The expression and function of Nod-like receptors in neutrophils. Immunology. 2010;130(1):55–63.
Dufner A, Pownall S, Mak TW. Caspase recruitment domain protein 6 is a microtubule-interacting protein that positively modulates NF-κB activation. Proc Natl Acad Sci USA. 2006;103(4):988–93.
Zhou R, Yazdi AS, Menu P, Tschopp J. A role for mitochondria in NLRP3inflammasome activation. Nature. 2011;469(7329):221–5.
Proell M, Gerlic M, Mace PD, Reed JC, Riedi SJ. The CARD plays a critical role in ASC foci formation and inflammasome signaling. Biochem J. 2013;449(3):613–21.
Wang LC, Swat W, Fujiwara Y, Davidson L, Visvader J, Kuo F, et al. The TEL/ETV6 gene is required specifically for hematopoiesis in the bone marrow. Genes Dev. 1998;12(15):2392–402.
Dwyer J, Li H, Xu D, Liu JP. Transcriptional regulation of telomerase activity: roles of the Ets transcription factor family. Ann NY Acad Sci. 2007;1114:36–47.
Siritantikorn A, Johansson K, Ahlen K, Rinaldi R, Suthiphongchai T, Wilairat P, et al. Protection of cells from oxidative stress by microsomal glutathione transferase 1. Biochem Biophys Res Commun. 2007;355(2):592–6.
Marchitti SA, Orlicky DJ, Vasiliou V. Expression and initial characterization of human ALDH3B1. Biochem Biophys Res Commun. 2007;356(3):792–8.
Medvedev AE, Gorkin VZ. The role of monoamine oxidase in the regulation of mitochondrial energy functions. Vopr Med Khim. 1991;37(5):2–6.
Chen Y, Leon-Ponte M, Pingle SC, O’Connell PJ, Ahern GP. T lymphocytes possess the machinery for 5-HT synthesis, storage, degradation and release. Acta Physiol (Oxf). 2015;213(4):860–7.
Heisterkamp N, Groffen J, Warburton D, Sneddon TP. The human gamma-glutamyltransferase gene family. Hum Genet. 2008;123(4):321–32.
Morgan DO, Edman JC, Standring DN, Fried VA, Smith MC, Roth RA, et al. Insulin-like growth factor II receptor as a multifunctional binding protein. Nature. 1987;329(6137):301–7.
Ingelfinger JR. Angiotensin-converting enzyme 2: implications for blood pressure and kidney disease. Curr Opin Nephrol Hypertens. 2009;18(1):79–84.
Trojanowicz B, Ulrich C, Kohler F, Bode V, Seibert E, Fiedler R, et al. Monocytic angiotensin-converting enzyme 2 relates to atherosclerosis in patients with chronic kidney disease. Nephrol Dial Transplant. 2016;0:1–12.
Selsted ME, Tang YQ, Morris WL, McGuire PA, Novotny MJ, Smith W, et al. Purification, primary structures, and antibacterial activities of β-defensins, a new family of antimicrobial peptides from bovine neutrophils. J Biol Chem. 1993;268(9):6641–8.
Kościuczuk EM, Lisowski P, Jarczak J, Krżyzewski J, Zwierzchowski L, Bagnicka E. Expression patterns of β-defensin and cathelicidin genes in parenchyma of bovine mammary gland infected with coagulase-positive or coagulase-negative Staphylococci. BMC Vet Res. 2014;10:246.
Janowska-Wieczorek A, Marquez LA, Nabholtz JM, Cabuhat ML, Montaño J, Chang H, et al. Growth factors and cytokines upregulate gelatinase expression in bone marrow CD34+ cells and their transmigration through reconstituted basement membrane. Blood. 1999;93(10):3379–90.
Vafadari B, Salamian A, Kaczmarek L. MMP-9 in translation: from molecule to brain physiology, pathology, and therapy. J Neurochem. 2016;139 Suppl 2:91–114.
Schrödl W, Büchler R, Wendler S, Reinhold P, Muckova P, Reindl J, et al. Acute phase proteins as promising biomarkers: perspectives and limitations for human and veterinary medicine. Proteomics Clin Appl. 2016;00:1–16.
This work was conceived within the EU COST Action FA1308 “DairyCare”. We are grateful to Prof. Bruno Stefanon for the helpful comments and discussion points leading to refinement of our concepts.
Project was approved by independent peer reviewers and on the basis of these evaluations the research was funded by the Italian Ministry of Education, University and Research (PRIN GEN2PHEN). The grant covered all the expensiveness of the study.
Availability of data and materials
The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request. The DIA software is not publically available, but can be requested from the developer. DAVID and KEGG databases are free and available online.
NP and TM performed RNA extraction from blood, elaboration of sequencing data applying DIA and DAVID tools, statistical elaboration of results and drafted the manuscript. MD and FP participated in the design of the study, in animal sampling and helped to draft the manuscript. DL elaborated raw RNA-Seq data using DeSeq2. PAM contributed to the final revision of the manuscript. AM participated in sample collection and helped in manuscript revision. SS participated in the design of the study and its coordination and also in animal sampling and biochemical analysis of blood and milk samples. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Consent for publication
Characteristics of animals chosen for the transcriptome analysis in the selected commercial farms. F: farm, IH: Italian Holstein, IS: Italian Simmental, DIM: days in milk, BCS: body condition score, SCC: somatic cell count, values are expressed as mean ± SD. (XLSX 11 kb)
Complete set of 10% most impacted pathways produced by DIA analysis in each comparison. Three images are presented. The first image (named A) collects the 10% most impacted pathways for APP+ vs. APP- comparison. The second image (named B) lists the 10% most impacted pathways for IH APP+ vs. IH APP- comparison. The third image (named C) lists the 10% most impacted pathways for IS APP+ vs. IS APP- comparison. (DOCX 211 kb)
Complete list of genes included in DAVID Annotation Clusters. Each sheet of the.xlsx file includes one of the Annotation Clusters listed in Table 6 of the present article. In each Annotation Cluster the included genes are listed. The genes found to be common in each comparison are bolded. (XLSX 57 kb)
About this article
Cite this article
Pošćić, N., Montanari, T., D’Andrea, M. et al. Breed and adaptive response modulate bovine peripheral blood cells’ transcriptome. J Animal Sci Biotechnol 8, 11 (2017). https://doi.org/10.1186/s40104-017-0143-y
- Acute phase proteins
- Adaptive response
- Dynamic impact approach (DIA)
- Hypothalamic-pituitary-adrenal (HPA) axis
- Stress response