Heat stress and immune response phenotype affect DNA methylation in blood mononuclear cells from Holstein dairy cows

Heat stress negatively affects health and production in cows. Examining the cellular response to heat stress could reveal underlying protective molecular mechanisms associated with superior resilience and ultimately enable selection for more resilient cattle. This type of investigation is increasingly important as future predictions for the patterns of heat waves point to increases in frequency, severity, and duration. Cows identified as high immune responders based on High Immune Response technology (HIR) have lower disease occurrence compared to their average and low immune responder herd-mates. In this study, our goal was to identify epigenetic differences between high and low immune responder cows in response to heat stress. We examined genome-wide DNA methylation of blood mononuclear cells (BMCs) isolated from high and low cows, before and after in vitro heat stress. We identified differential methylation of promoter regions associated with a variety of biological processes including immune function, stress response, apoptosis, and cell signalling. The specific differentially methylated promoter regions differed between samples from high and low cows, and results revealed pathways associated with cellular protection during heat stress.


Sample collection and BMC isolation. Sample collection was approved by Animal Care Services at the
University of Guelph (Animal Utilization Protocol #3555) and all experiments were conducted in accordance with their rules and regulations. Blood was collected from healthy multiparous Holstein cows that had previously been classified as high (n = 3) or low (n = 3) immune responding animals based on their EBVs using the HIR technology. All blood collections were done during the winter months (November-March) to avoid previous exposure to heat and to ensure animals were not heat stressed. Whole blood was centrifuged at 1200×g for 20 min with the centrifuge break turned off. The buffy coat was isolated and diluted in a 1:2 ratio with PBS. The buffy coat-PBS solution was then layered over histopaque (Sigma-Aldrich, Oakville, ON) in a 1:1 ratio and spun at 1200×g for 10 min with the break off. The buffy coat was isolated and the volume was made to 25 ml with PBS. The cell suspension was spun at 150×g for 10 min with the brake on. Supernatant was discarded and cell pellet was resuspended in 3 ml PBS. The total number of cells per ml was determined using an ORFLO cell counter (ORFLO Technologies, Ketchum, ID). Cell viability was checked using trypan blue solution (Sigma-Aldrich, Oakville, ON.) via hemocytometer method.
Heat treatment. Cells were diluted in RPMI media (Thermo-Fisher Scientific, Ottawa, ON) containing fetal bovine serum (Sigma Aldrich, Oakville, ON) and Penicillin-Streptomycin (Sigma Aldrich, Oakville, ON) to a concentration of 10 5 viable cells/ml. Cells were then plated in three 1 ml replicates on a 24-well flat bottom cell culture plate (Sigma Aldrich, Oakville, ON) destined for TN treatment and three more 1 ml replicates on a second 24-well flat bottom plate for HS treatment. All plates were placed at 37 °C with 5% CO 2 overnight. On the following day HS plates were moved to an incubator that was set at 42 °C with 5% CO 2 for 4 h. TN plates remained at 37 °C for this period. Samples were collected immediately after completion of the 4-h heat challenge. For each plate, cell supernatant was collected from each well and the 3 wells for each cow were pooled. Adherent cells were lifted using Trypsin (Sigma Aldrich, Oakville, ON) and added to the pool before spinning at 300×g. Supernatant was discarded and the cell pellet was homogenized in 1 ml of Trizol reagent following the manufacturer's protocol (https:// tools. therm ofish er. com/ conte nt/ sfs/ manua ls/ trizol_ reage nt. pdf.), flash frozen, and stored at − 80 °C.
The heat treatment method used in this study was done to assess molecular changes in response to acute heat stress. Previous studies have indicated changes in the concentration of heat shock proteins, including HSP70 concentrations, following a heat stress treatment of one-two hours [44][45][46][47] . Furthermore, as part of a connected study, the concentration of HSP70 was measured in BMC cultures (n = 45), that included the animals that were used in this study, before and after the four-hour heat treatment using a commercial ELISA kit (Abclonal, Woburn, MA) following the manufacturer's instructions (Supplementary Figure 1). DNA isolation. DNA was isolated from one control and one heat stressed (HS) sample from each of three high and three low cows, making a total of 12 samples (Table 1). The DNA was extracted following the manufacturer's instructions for DNA extraction from Trizol Reagent (https:// tools. therm ofish er. com/ conte nt/ sfs/ manua ls/ trizol_ reage nt. pdf).
Library construction for DNA methylation analysis. Library 48 . Libraries were prepared by DNA end polishing and adaptor ligation followed by library amplification using indexed primers and library purification. Purified library DNA was eluted in 12 μl of water (Epigentek, NY).
Enhanced MethylSeq sequence alignment and statistical analysis. Libraries were sequenced on an Illumina HiSeq 4000 with 50 bp single reads. After quality control of raw reads with FASTQC 49 version 0.11.8, low quality and adapter trimming was performed using Trim Galore 50 using the following parameters: minimum Sanger Phred score of 20, 3′ adapter trimming, and removal of reads shorter than 20 bp (RRBS-specific: additional 2 bps from any adapter-containing reads are removed. This option is used to avoid the possible filledin cytosine after bisulfite treatment and restriction enzyme digestion).
Cleaned reads were aligned to the Bos Taurus reference genome (BosTau9) using Bismark 48 , version 0.203.0. Bismark uses Bowtie 51 , version 2.2.5 with the parameters "-directional" for targeted bisulfite sequencing libraries and "-pbat" for post-bisulfite prepared RRBS libraries. Samtools 52 , version 0.1.9, was used to sort the SAM files produced by Bismark and remove the duplicate reads resulting from PCR amplification. Methylation information was extracted from the final Bismark mapping result at the base resolution using the Bismark methylation extractor (version 0.203.0). A minimal read coverage of 5 and minimal quality score of 20 at each base position were applied. MethylKit 53 was used to summarize methylation across promoter regions and conduct differential methylation analysis using default parameters. Promoters were defined as genomic regions from − 2000 bp to the transcription start site. Samples were filtered by coverage (minimum 5), normalized, merged (only promoter regions that were covered by a minimum of two samples in the comparison were included), and subjected to DMP identification with a cut off of 15% methylation difference and a q-value < 0.05. False discovery rate and multiple hypotheses testing were controlled using a Sliding Linear Model 54 .
The promoter region was selected for differential methylation analysis because DNA methylation in the promoter region is known to regulate gene activity: increased DNA methylation in the promoter region has been reported to suppress transcription 38,55 . Therefore, hypomethylation of a promoter region was assumed to suggest increased gene transcription whereas hypermethylation was assumed to indicate decreased transcription. The identified differentially methylated promoters (DMPs) were annotated against the RefSeq genes. Venn diagrams were generated with VENNY (version 2.1) 56 to visualize overlapping gene sets. qiagen. com/ ingen uity) was used to perform metabolic pathway and gene network analyses using the gene sets associated with DMPs. Ensembl gene IDs were converted into human gene symbols with Biomart (http:// www. bioma rt. org) before performing analyses in IPA 57,58 . Only the non-annotated genes, with a percentage of identity with the human homolog higher than 80%, were annotated by this approach 58,59 . Additionally, the Regulator Effects tool in the IPA package was used to identify potential regulators of the heat stress response in samples from high and low cows 58,59 . Table 1. Experimental design. This study contained a total of 12 samples from each of three cows immune response phenotyped as high (H) or low (L). A control (C) and heat stressed (HS) sample were generated for each cow, as indicated in the treatment column.

Results
Overview of genome-wide DNA methylation profiling. In order to better understand the changes to DNA methylation resulting from heat stress and immune phenotype, we carried out Enhanced MethylSeq on BMCs from three high and three low immune responder dairy cows. Illumina HiSeq 4000 sequencing generated between 39.7 and 58.1 million reads per sample. After quality filtering, 72% to 81% of the reads were successfully aligned to the bovine reference genome (bosTau9). Of these reads, between 40 and 66% mapped uniquely across libraries. In total, we identified 9 to 22 million CpG sites per sample, of which 6.7% on average were methylated in the CpG context. Raw sequencing data and mapping statistics are summarized in Supplementary Table S1. After correcting for multiple tests, differential DNA methylation analyses detected a total of 172 significant DMPs (q-value < 0.05) across the four comparisons described below. Principle component analysis separated the samples according to individuals indicating that individual differences were stronger than the effect of heat treatment on DNA methylation patterns (Fig. 1a). Some DMPs were shared across the four groups, but most were specific to the comparison (Fig. 1b).
Differentially methylated promoters between immune phenotypes. Heat stressed samples. The comparison between HS samples from high cows and low immune responder cows (samples 2, 4, and 6 versus 8, 10, and 12, see Table 1) detected 47 significant DMPs, of which 34% were hypermethylated in high samples and 66% were hypomethylated in high samples. Five lncRNAs were among the 47 genes associated with the DMPs (see Supplementary Table S2).
The top hypermethylated promoter regions in HS samples from high immune responder cows were associated with the genes OPRD1 (opioid receptor delta 1) and NUDT15 (Nudix Hydrolase 15). The top hypomethylated promoter regions in samples from high immune responder cows were associated with the genes BCL2L12 (BCL2 Like 12) and CPEB1 (Cytoplasmic Polyadenylation Element Binding Protein 1). Additional genes associated with DMPs in this subset that were hypomethylated in samples from high cows include HSPB9 (Heat Shock Protein Family B (Small) Member 9), IL15 (Interleukin 15), and NDRG1 (N-Myc Downstream Regulated 1) (see Supplementary Table S2).
Control samples. The comparison between control samples from high immune responder cows and low immune responder cows (samples 1, 3, and 5 versus 7, 9, and 11, see Table 1) revealed 44 DMPs, of which 43.5% were hypermethylated in high samples and 56.5% were hypomethylated in high samples. The top hypermethylated promoter regions in control samples from high immune responder cows were associated with the gene SLC37A4 (Solute Carrier Family 37 Member 4) and an uncharacterized lncRNA. The top hypomethylated regions in high samples were associated with the genes ALDOA (Aldolase, Fructose-Bisphosphate A) and microRNA bta-mir2302. The promoter region of IL15 was also hypomethylated in high samples in this subset (see Supplementary Table S3).  www.nature.com/scientificreports/ Differentially methylated promoters between HS and control samples. High samples. The comparison between control and HS samples from high responder cows (samples 1, 3, and 5 versus 2, 4, and 6, see Table 1) revealed 18 DMPs, of which just over half were hypermethylated in HS samples (55.5%). One lncRNA is included in these 18 genes (see Supplementary Table S4). The top hypermethylated promoter regions in HS samples were associated with the genes ARSG (Arylsulfatase G), APC2 (APC Regulator of WNT Signalling Pathway 2), and an uncharacterized lncRNA. The top hypomethylated promoter regions in HS samples were associated with the genes IGF2 (Insulin Like Growth Factor 2) and CPEB1 (Cytoplasmic Polyadenylation Element Binding Protein 1). JMJD8 (Jumonji domain-containing protein 8) was also hypomethylated in HS samples (see Supplementary Table S4).
Low samples. The comparison between control and HS samples from low immune responder cows (samples 7, 9, and 11 versus 8, 10, and 12, see Table 1) detected 64 DMPs, of which 53% were hypermethylated HS samples and 47% were hypomethylated in HS samples. Among these 64 DMPs were six lncRNAs (see Supplementary  Table S5).
In HS samples, the top hypermethylated promoter regions were associated with the genes RHOT2 (Ras Homolog Family Member T2) and NCKIPSD (NCK Interacting Protein With SH4 Domain), and top hypomethylated promoters were associated with the genes MMTAG2 (Multiple Myeloma Tumor-Associated Protein 2 Homolog), CCDC40 (Coiled-Coil Domain Containing 40), and NGF (Nerve Growth Factor) (see Supplementary Table S5).
Metabolic pathways and gene networks connecting DMPs. Genes associated with DMPs were mapped to the IPA database for metabolic pathway and gene network exploration. Significantly enriched metabolic pathways were identified in all four sets of DMPs. However, this result should be interpreted with caution because pathways were generally represented by a small number of genes 58 (Fig. 2a). In control samples, the top canonical pathways associated with DMPs in high versus low samples included IL15 Production (p = 1.33E−02) (Fig. 2b).
The DMPs were also grouped in gene regulatory networks with the IPA software. The top scoring regulatory network connecting DMPs identified in high versus low heat stressed samples was Cell-to-cell Signalling and Interaction, Amino Acid Metabolism, Cell Death and Survival (Fig. 3a). The top network connecting DMPs in high versus low control samples was Neurological Disease, Skeletal and Muscular Disorders, Immunological Disease (Fig. 3b) The Regulator Effects tool in the IPA package was used to identify potential regulators of the heat stress response in samples from high and low cows. This tool integrates upstream regulator results with downstream effects results to build causal hypotheses that help to interpret what may be occurring upstream to cause particular phenotypic or functional outcomes. Three regulators were identified [NFkB (Nuclear Factor-kB), IL15, and  (Fig. 4).

Discussion
Heat stress (HS) negatively affects production 60,61 and suppresses immune 62,63 and reproductive function in bovines 63 . There is evidence that HS increases mastitis incidence 64 and other health problems including ketosis, liver lipidosis, and even mortality 64,65 . Studies investigating the molecular changes in response to heat stress have shown epigenetic and transcriptional modulation of genomic regions associated with immune function, stress response, metabolism, and cell signalling 42,43,66 . To our knowledge, this is the first attempt made to study genome-wide DNA methylation changes to promoter regions resulting from heat stress in a population of BMCs isolated from dairy cows of distinct immune response phenotypes.
Analysis of differential methylation of promoter regions by principle component analysis revealed that the most variability in promoter regions was between individuals. However, 172 significant DMPs (q-value < 0.05) were identified across the four comparisons discussed below. Indeed, results from this study suggest that high cows may also have increased expression of interleukin 15 (IL15): within the HS and control groups, the IL15 promoter region was significantly hypomethylated in samples from high cows compared to low cows (see Supplementary Tables S2 and S3). Indeed, increased DNA methylation within gene promoter regions has been reported to suppress transcription 39 , suggesting that hypomethylation of gene promoter regions could be associated with up-regulated gene expression. IL15 is important for proliferation, survival, and differentiation of natural killer (NK) cells and T cells, and is implicated in the signalling pathway categories for apoptosis, cellular immune response, and cytokine signalling. The IL15 gene plays a key role in the immune response to viral and bacterial www.nature.com/scientificreports/ infections by regulating cells of both the innate and adaptive immune system 67 . Previous work has shown that IL15 and its receptor, IL15Ra, support NK cell homeostasis under resting conditions, and mediate NK cell survival and differentiation into functional NK cells capable of killing virally infected cells 68 . As such, up-regulation of IL15 expression in BMCs could contribute to the superior immune response phenotype of high immune responder cows by enhancing the function of their immune cells. Further expression studies to investigate this potential difference in IL15 between high and low cows will be informative. The analysis of individual DMPs across the bovine genome revealed genes involved in a range of functions. Within the HS groups, signalling pathways associated with DMPs between high and low samples included cytokine communication, cell death and survival, and cell-to-cell signalling and interaction. Many of the DMPs identified in this comparison were associated with genes that are important for protecting cells against stress. For example, the gene BCL2L12 is an anti-apoptotic member of the Bcl2 protein family 69 . Out of all of the DMPs in heat stressed high immune responder versus low immune responder samples, the promoter region of BCL2L12 showed the lowest methylation density in high samples (Table S2). Higher expression of BCL1L12 in high immune responder cows could be related to more anti-apoptotic activity and better cell survival. Furthermore, HSP70 has been show to interact with and protect BCL2L12 from degradation 70 , supporting the protective profile of DMPs in HS samples from high cows. Indeed, heat shock proteins have been extensively studied for their role in protecting cells from heat and other forms of stress (reviewed in 71 ). These molecular chaperones have cryoprotective and protein re-folding functions that preserve protein structure and transport 8 and therefore, elevated expression of HSPs is a desired response to an increase in body temperature [9][10][11] . The cytoprotective effect of hsp70 is partly due to its ability to impede apoptosis 72 . As a means of promoting cell survival, the heat stress response and HSPs are also known to play a role in inflammatory signalling by regulating the production of inflammatory cytokines 72 .
Additional genes with protective functions whose promoters were hypomethylated in BMCs from high cows include heat shock protein B9 (HSPB9) and N-myc downstream-regulated gene 1 (NDRG1) ( Table S2). Expression of HSPB9 in response to heat stress has been shown in brain, liver, and muscle tissue from chicken 73    www.nature.com/scientificreports/ NDRG1 gene is an important stress response protein that responds to a variety of cellular stressors and has a putative function in suppression of tumour metastasis 74,75 . The NDRG1 gene is also induced by hypoxia and iron depletion 76 . The top scoring regulatory network connecting genes with DMPs between high and low HS samples was Cellto-cell Signalling and Interaction, Amino Acid Metabolism, Cell Death and Survival (Fig. 3a). The tumour necrosis factor (TNF) molecule occupies a central position in the network shown in Fig. 3a. The TNF gene has a broad range of functions, including cell proliferation and differentiation, inflammation, and apoptosis 77 . TNF is also part of the IL-1 and IL-6 fever cascade that acts on the hypothalamus during infection with certain pathogens 78 . This network may represent the dual effect of HSPs on cell survival-on the one hand, as discussed here, HSPs play a critical role in cell survival in response to heat stress. However, it has also been shown that HSPs sensitize cells to certain apoptotic stimuli, such as TNF. Ran et al. 79 showed that HSP70 enhances TNF-mediated apoptosis by binding to IκB kinase γ and impairing NF-κB survival signalling. However, a study by Imao et al. 80 demonstrated that repeated heat stress was needed to initiate this apoptotic pathway.
The top network connecting DMPs in high versus low control samples (Fig. 3b) was Immunological Disease, Neurological Disease, and Skeletal and Muscular Disorders. The network shown in Fig. 3b shows groups and complexes important for immune function occupying central positions and having direct connections with genes identified with DMPs. For example, Jun N-terminal kinases (JNKs) and p38 mitogen-activated protein kinases (MAPKs) have important roles in the cellular response to many types of stress, as well as regulating the activity of inflammatory mediators 80 . Also in a central position was NF-κB, which is critical for regulating immune function 81 . The Regulatory Effects tool in IPA also identified NF-κB, along with IL15, and CD40, as regulators of four genes (IL15RA, TRAF3, CD44, and CSF3) with DMPs between high and low HS samples (Fig. 4). The NF-κB complex functions as a transcription factor to regulate a broad range of biological processes including immune function, inflammation, and stress responses 79 . These regulators may represent some of the differences in immune function between high and low cows and further study could reveal molecular differences between immune response phenotypes.
The network shown in Fig. 3b also shows the V-akt murine thymoma viral oncogene homologue molecule (AKT) in a central position, having connections with genes identified with DMPs (for example, IL15 and BNIP3). The AKT gene is a main regulator of glucose homeostasis 82 , suggesting that high and low cows could have differences in this process. Indeed, the promoter region of the solute carrier family 36 member 4 gene (SLC37A4), better known as the glucose-6-phosphate transporter (G6PT), had the highest density of methylation (hypermethylated) in high samples compared to lows. This sugar-phosphate exchanger maintains glucose homeostasis and neutrophil energy homeostasis, and deficiency is responsible for glycogen storage disease type Ib 83 . Hypermethylation of the G6PT promoter region in BMCs from high cows could point to decreased expression compared to low cows under resting conditions. Further expression profile studies should be conducted to investigate this potential difference in energy metabolism between high and low cows.
Significant DMPs were also identified in the comparison between control and HS treatment for both high and low samples. However, the specific DMPs that were identified in this comparison differed between high and low groups, suggesting that the cellular response to HS is different in cows of distinct immune response phenotypes. In samples from high cows, the top up-regulated metabolic pathway associated with the DMPs identified between control and HS samples was Ceramide Biosynthesis (Fig. 2c); ceramide plays a role in mediating apoptosis in response to cytokines and environmental stress. The promoters with the lowest methylation density in HS samples were insulin-like growth factor (IGF2) and cytoplasmic polyadenylation element-binding protein (CPEB1). CPEB1 is a key factor in controlling mRNA translation. Xiaoping et al. 84 demonstrated that hypomethylation of the CPEB1 promoter resulted in overexpression, suggesting that hypomethylation of the promoter region in our study could also result increased expression.
The jumonji C domain-containing protein JMJD8 gene promoter was hypomethylated in HS compared to control samples from high cows. Jumonji proteins have been shown to regulate cellular processes by hydroxylating or demethylating histone and non-histone targets 85 . An important function of genes in the jumonji family is modulation of gene expression via histone post-translational modifications 86 . The JMJD8 gene is involved in angiogenesis and cellular metabolism 87 , and is also a positive regulator or TNF-induced NF-κB signalling pathway 88 . Mass spectrometry analysis revealed two HSP proteins bound to JMJD8: HSPA5 and HSP90B1 89 , suggesting that JMJD8 could play a role in the stress response. A previous study demonstrated a role for another jumonji family protein (JMJD1A) in regulating the response to cold stress (4 °C for one week) in mice by upregulating genes associated with the thermogenic response in brown adipose tissue 90 .
Eight DMPs were hypermethylated in HS samples from high cows, including APC2 and BNIP3. APC2 is a regulator of the WNT signalling pathway and BNIP3 plays a critical role in inducting autophagy during heat stress 91 and is also proapoptotic 92 . The hypermethylation of the BNIP3 gene promoter in HS samples from high cows suggests decreased mRNA expression compared to control samples, potentially increasing cell survival during heat stress.
The comparison between control and HS samples from low cows revealed 64 DMPs associated with a range of metabolic and signalling pathways. The promoter region with the lowest methylation density (hypomethylated) in HS samples was the multiple myeloma tumour-associated protein 2 (MMTAG2) gene. Luo et al. 93 showed that MMTAG2 (also known as C1orf35) is involved in cell growth by promoting cell cycle progression from G1 to S. Overexpression of MMTAG2 results in up-regulation of c-MYC and subsequent accelerated cell proliferation 94 , perhaps revealing a cell survival mechanism activated in low cows. The promoter regions of two histone deacetylases (HDAC10, HDAC4) were hypermethylated in HS samples from low cows. In general, acetylation of histones is permissive to gene expression because it opens up chromatin so that it is accessible to transcription factors, whereas deacetylation represses gene expression. Previous 97 demonstrated that heat-shock factor 1 (a transcription factor that activates HSPs), is a master regulator of chromatin acetylation across the genome in response to heat stress. Hypomethylation of the promoter regions of both HDAC4 and HDAC10 in HS samples suggests that expression of these genes was decreased upon heat stress. Hence, a mechanism of cell protection used by BMCs from low cows could be activating expression of HSPs via decreasing expression of histone deacetylases.

Conclusion
We identified significant DMPs in BMCs isolated from immune phenotyped Holstein dairy cows, identified as high or low immune responders, under both control and heat stress conditions. These results suggest that DNA methylation of promoter regions could contribute to variation in immune phenotypes in dairy cows, and also variation in the response to heat stress. The DMPs differed between samples from high and low cows, but our results revealed pathways that could provide protection to the cells during heat stress in both immune phenotypes. In samples from high cows, heat stress resulted in differential methylation of gene promoters associated with stress response and apoptosis prevention, whereas in low cows, heat stress affected promoter methylation of genes associated with cell proliferation and histone deacetylases. Our results also revealed potential differences between high and low cows under control conditions: hypomethylation of the IL15 promoter region in samples from high cows suggests higher expression of this cytokine in high cows compared to their low herd mates. Additionally, hypermethylation of the G6PT promoter in samples from high cows could point to lower mRNA expression compared to low cows, perhaps revealing a difference in energy metabolism. Future studies evaluating the transcription profile of genes in response to heat stress will provide more insight into the functional relevance of our findings. Furthermore, in addition to the changes in DNA methylation of promoter regions, it is likely that heat stress alters gene expression through other epigenetic modifications as well, such as histone modifications and microRNAs. More information about changes in the epigenomic landscape in response to heat stress as well as corresponding transcriptional changes would improve our understanding of the molecular differences between immune phenotypes in bovine. This knowledge will allow us to better understand the relationship between immune response phenotype and response to heat stress. Considering the moderately high heritability of HIR (0.35) 31 , immuno-genetics could provide a means to genetically select for improved response to heat stress. Identifying cattle with both genetically enhanced immunity and resilience to heat stress would provide a new tool to increase livestock efficiency in the context of climate change.

Data availability
The datasets generated and analysed during this study are available in the GEO repository with the accession number GSE163222. www.nature.com/scientificreports/