Methylome Patterns of Cattle Adaptation to Heat Stress.

Heat stress has a detrimental impact on cattle health, welfare and productivity by affecting gene expression, metabolism and immune response, but little is known on the epigenetic mechanisms mediating the effect of temperature at the cellular and organism level. In this study, we investigated genome-wide DNA methylation in blood samples collected from 5 bulls of the heat stress resilient Nellore breed and 5 bulls of the Angus that are more heat stress susceptible, exposed to the sun and high temperature-high humidity during the summer season of the Brazilian South-East region. The methylomes were analyzed during and after the exposure by Reduced Representation Bisulfite Sequencing, which provided genome-wide single-base resolution methylation profiles. Significant methylation changes between stressful and recovery periods were observed in 819 genes. Among these, 351 were only seen in Angus, 366 were specific to Nellore, and 102 showed significant changes in methylation patterns in both breeds. KEGG and Gene Ontology (GO) enrichment analyses showed that responses were breed-specific. Interestingly, in Nellore significant genes and pathways were mainly involved in stress responses and cellular defense and were under methylated during heat stress, whereas in Angus the response was less focused. These preliminary results suggest that heat challenge induces changes in methylation patterns in specific loci, which should be further scrutinized to assess their role in heat tolerance.

Introduction

Rapid and unpredictable climate change and extreme climatic events (floods, drought, extreme temperatures) are increasing in frequency, which have an impact on agricultural productivity and affect food security. These changes are not only limited to tropical and subtropical regions but also affect more temperate regions (Rajaud and Noblet-Ducoudré, 2017). In the most vulnerable areas, management and breeding of livestock should address heat stress and decreasing water availability, while maintaining fitness and productivity in order to remain viable (Rojas-Downing et al., 2017). To satisfy the demand for animal products, which is increasing at a faster rate than human population growth, it is necessary for the livestock sector to respond to rapid changes in the climate.

Heat stressed animals are not able to adequately dissipate the excess of endogenous and exogenous heat to maintain the body thermal balance (Bernabucci et al., 2014). Heat stress is a primary stressor for high-producing dairy cows which have a high metabolic rate with associated endothermy, which results in heat-induced depression of milk production and growth (Cowley et al., 2015). Therefore, heat stress results in economic losses due to reduced production and reproduction performance (Aguilar et al., 2010; Nardone et al., 2010; Mehla et al., 2014; Biffani et al., 2016; Macciotta et al., 2017). Physiological changes are observed affecting the ability of the animal to cope with heat stress. For example, rectal temperature and respiration rate increase (Dikmen et al., 2015; Perano et al., 2015). Some of these reponses have a genetic component, moderate heritability and have been associated with genomic variants, e.g., several candidate genes have been identified for rectal temperature (Dikmen et al., 2013, 2015). Traditional phenotype-based selection for adaptation and fitness, while maintaining productivity, is difficult and slow as these traits are difficult to measure. For this reason, genomic information could accelerate the genetic progress for adaptation traits. SNP markers associated with heat tolerance have been detected in dairy cattle (Hayes et al., 2009; Dikmen et al., 2012, 2013), although few causative genes have been identified (reviewed in Collier et al., 2008). One of these, the prolactin receptor gene (PRLR), carries mutations having a major effect on heat tolerance (Littlejohn et al., 2014; Porto-Neto et al., 2018). This variant causes the SLICK phenotype which is easily detected as it confers animals a short and sleek coat (Olson et al., 2003; Dikmen et al., 2008). Holstein cattle into which a SLICK haplotype was introduced had superior thermoregulatory ability compared with wild-type Holsteins (Dikmen et al., 2014).

Little is known about the biology of adaptation, and even less about mechanisms that mediate changes in gene expression and metabolism in animals subjected to environmental stress. Technologies developed to sequence and characterize the human genome have paved the way for the routine sequencing of the genomes of agriculturally important species. In addition to the genome sequence, epigenetic features control gene expression and have attracted much attention in the last few years. At least four molecular systems are involved in the control of gene expression, including DNA methylation (Bernstein et al., 2007), non-coding RNAs (Mitchell et al., 2009), histone post-translational modifications (Strahl and Allis, 2000), and chromatin remodeling (Cosma et al., 2010). DNA methylation was the first identified and is currently the better-studied epigenetic regulatory mechanism of gene expression. In mammalian cells, DNA methylation mostly occurs at the cytosine of a CpG dinucleotide, and CpG islands, which are enriched for CpG dinucleotides, are frequently located in promoter regions at the 5′ of coding sequences.

There is still limited knowledge on the relationship between epigenetic patterns and phenotypic variation in livestock. Nevertheless, epigenetic effects on some trait have been described, e.g., growth is affected by epigenetic imprinting of IGF2 in pig (Van Laere et al., 2003), and epigenetic regulation of callipyge affects development in sheep (Georges et al., 2003; Vuocolo et al., 2007). Recent studies have characterized pig and chicken methylomes (Li et al., 2012; Nätt et al., 2012). An atlas of the porcine methylome identified that differentially methylated regions contain ∼80% of the known or candidate human obesity-related genes, 72% of which mapping in QTL regions that affect fatness and pork quality (Li et al., 2012). Nätt et al. (2012) suggested that the heritability of methylation in the chicken genome might play a role in adaptation and selection since they observed an over-representation of differentially methylated genes in selective sweep regions associated with domestication. This could mean that novel methylation patterns have been acquired by domestic animals during their selection history. In dairy cows, recent work has shed light on the role of epigenetics in the regulation of milk synthesis and mammary development (Singh et al., 2010). In the lactating cow, remethylation of the typically hypo-methylated casein gene promoter drastically reduces casein mRNA expression and milk protein synthesis during acute udder infection.

The association between DNA methylation patterns and heat stress has been suggested. For example exposure of male guinea pigs to chronic heat stress has been shown to alter DNA methylation patterns in the liver of both the F0 and F1 generations (Weyrich et al., 2016).

In the current paper, epigenomic responses to heat stress were investigated in indicine (Bos indicus) and taurine (Bos taurus) cattle. We hypothesize that, in addition to morphological traits determined by genetics (coat color, body size and conformation, large ears, loose skin, sweat gland size), the adaptation to heat and the better ability to cope with deleterious effects of a high temperature of the indicine Nellore compared to the taurine Angus breed is also associated with changes in DNA methylation patterns of specific regulatory genes and metabolic pathways. In the study reported here, cattle of these two breeds were exposed to the sun during the high temperature and high humidity season of the Brazilian South-East region and the physiological and epigenetic responses assessed.

Materials and Methods

Animals Treatment and Experimental Design

A total of 25 Nellore (indicine heat-tolerant breed) and 25 Angus (taurine heat-susceptible breed) bullocks were included in the investigation. Animals were all healthy young males of about 15 months of age at the time of the investigation. We selected half-sibs animals within each breed to minimize genetic variation. Angus bullocks were purchased at 7 months of age in Uruguaiana (Rio Grande do Sul state, Brazil) and arrived at UNESP Aracatuba (Sao Paulo state, Brazil) experimental station (located at −21.186244 latitude and −50.439053 longitude) on June 13, 2015. Seven-month old Nellore bullocks were purchased in Dourados (Mato Grosso do Sul state, Brazil), kept at Agua Branca farm in Birigui (Sao Paulo, Brazil) until 10 months of age and moved to UNESP on September 2, 2015. Animals were kept at the experimental station of the Veterinary Faculty, Universidade Estadual Paulista – UNESP- Araçatuba Campus. A full description of the diet fed during the study is given in Supplementary File 1, while a schematic representation of experimental design is given in the Figure 1.

FIGURE 1

Schematic representation of the experimental design.

. On October 3, 2015, the steers were randomly assigned to 2 groups (1 Nellore and 1 Angus groups) of 12 or 13 animals and kept in two 200 square meter paddocks of which 100 square meters was covered by a shadowing net (80% sunblock) with regular access to pasture until December 3 (60 days). This period will be further referred to as pre-challenge adaptation period (PRE).

The experiment started on December 4, 2015, by removing the shadowing net from paddocks. Animal groups were therefore kept without shade until February 3, 2016 (56 days). This period is subsequently referred to as the experimental challenge period (EXP).

On February 4th shadowing nets were replaced such that all animals were kept with shade available and were allowed access to pasture until slaughter at the end of June 2016. This period will be further referred to as post challenge recovery period (POST).

Nellore adaptation period terminated in November 2015, when the season was already hot and humid and animals were under moderate stres. Sampling was therefore planned at the beginning of February 2016, at peak heat stress period, and in the middle of June 2016, during the cool season, after full recovery from heat stress.

Environmental Data Collection

Environmental data were collected hourly from June 1, 2015 to October 10, 2016 from CETESB meteorological station located in the UNESP campus, nearby the experimental site. Experimental variables collected were Temperature (T, in centigrades) and Humidity (H, in percentage). Missing data were imputed using the median among −10/ + 10 day data at the same time.

DNA Isolation, Library Preparation and Sequencing

Blood (10 ml) was collected from the same 5 randomly selected Nellore (heat stress resilient) and 5 randomly selected Angus (heat stress susceptible) belonging to sun-exposed group, in two periods: in February 2016 (stressful/challenge period) and June 2016 (recovery period). We used QIAamp DNA Blood Midi Kit (Qiagen) procedures to extract DNA from whole blood. One μg input was used in the MSP1 digest by overnight incubation at 37°C, following manufacturer instruction. Libraries were prepared using the TruSeq® DNA PCR-Free Library Preparation Kit (Illumina) including a step of bisulfite treatment. Subsequently, ligated products corresponding to DNA fragments 150–400-bp long were converted with EpiTectBisulfite Kits (Qiagen) and finally, PCR amplified with KAPA HiFi Uracil + (KapaBiosystems) to obtain RRBS libraries. DNA with a known methylation level was used as a spike control, and all conversion rates were >99%, ranging from 99.11 to 99.53% (Supplementary Table 1). Twenty Reduced Representation Bisulfite Sequencing (RRBS) libraries were used for cluster generation and subsequent sequencing on Illumina HiSeq 2500 PE 2 × 50 bp (NXT-Dx Ghent, Belgium).[1]

Environmental Indexes

From environmental data we obtained THI (Temperature Humidity Index), which considers temperature (T) and humidity (H) in a single index: (1.8 ∗ T + 32) − [(0.55 − 0.0055 ∗ H) ∗ (1.8 ∗ T - 26)] (Thom, 1959). Values above 75 were considered as moderate stress, above 85 as severe stress and over 90 as extreme stress (McDowell et al., 1976).

Data Analysis

Preliminary quality control of raw reads was carried out with FastQC v0.11.9.[2] The FASTQ sequence reads were generated using the Illumina Casava pipeline 1.8.2. The quality and adapter trimming of Illumina raw sequenze was performed with Trim Galore v0.6.1[3] using a two-step approach, which allowed us to remove two additional bases containing a cytosine, which were artificially introduced in the end-repair step during the library preparation. Bismark software (version 0.22.3) (Krueger and Andrews, 2011) was used to align each bisulfite-treated read to the bovine reference genome (ARS-UCD1.2) with option-N 1 (maximum number of mismatches allowed). The reference genome was first transformed into a bisulfite-converted version (C-to-T and G-to-A conversions) and then indexed using bowtie2 software (Langmead and Salzberg, 2012). Sequence reads were also transformed into fully bisulfite-converted versions (C-to-T and G-to-A conversions) before they were aligned to similarly converted versions of the genome in a directional manner. Sequence reads that produced the best unique alignment from the two alignment processes (original top and bottom strand) were then compared to the normal genomic sequence, and the methylation state of all cytosine positions in the reads was inferred using the Bismark_methylation_extractor function. Read duplicates were marked and removed using Picard Tools v2.23.[4]

Identification of Differentially Methylated Regions (DMRs) and Genes (DMGs)

DMRs were identified within the statistical environment R using the library methylKit (Akalin et al., 2012), which applies a sliding-window approach with a window of 1,000 bp and a step size of 500 bp. We first filtered out bases with less than 10 reads or more than the 99.9th percentile of coverage distribution. This step was done to ensure a good quality of the data and great confidence in the methylation percentage. Coverage values were then normalized by default and bases were merged in order to retain the ones that were covered in all samples. A logistic regression was then implemented to obtain p values. Specifically, in any given region methylKit models the methylation proportion P, for sample i = 1,…,n (where n is the number of animals) by the following logistic regression model:

where T denotes the treatment indicator for sample i, T = 1 if sample i is in the treatment group (stress period) and T = 0 if sample i is in control group (recovery period). The parameter β0 denotes the log odds of the control group and β1 the log oddsratio between the treatment and control group. Therefore, independent tests for all the regions of interest are run against the null hypothesis H0: β1 = 0. The null hypothesis is rejected when the logodds (and hence the methylation proportions) are different between the treatment and the control group and the region therefore classified as a DMR. P values were then adjusted to q values using the Sliding Linear Model (SLIM) method (Wang et al., 2011). We excluded covariates and over dispersion correction from the model since they showed not to provide significant changes to the results. DMRs were considered statistically significant if the q value was lower than 0.05, read coverage greater than 10 in all samples and methylation difference between groups at least 15%. When a DMR and a specific gene region overlapped, the corresponding gene was selected as the DMR-related gene, namely a differentially methylated gene (DMG). Gene features within DMGs were annotated by matching information available in the General Feature Format (GFF) file downloaded from the Ensembl database, release 101.[5] The promoter region was defined as an area 2-kb upstream of the transcription start sites.

Gene Ontology and KEGG Enrichment Analysis of DMR-Related Genes

Gene Ontology (GO) and KEGG pathway enrichment analysis of DMGs were performed by clusterProfiler, an ontology-based R package able to automate the process of biological term classification and enrichment analysis of gene clusters and to provide a visualization module for displaying analysis results (Yu et al., 2012). GO terms and KEGG with p values of less than 0.01 and q values of less than 0.1 were considered significantly enriched by DMR-related genes.

Results

Heat Stress Based on Environmental Parameters

The environmental conditions were analyzed and Temperature Humidity indices (THI), were calculated to evaluate the severity of heat challenge before during and after the experimental period (Figure 2).

FIGURE 2

Ridgeline plot showing distribution of daily average THI recorded at CETESB meteorological station located in the UNESP campus from October 2015 to September 2016.

Global Mapping of DNA Methylation

A mean of 60 million raw reads were generated for Nellore samples and 67 million for Angus samples during challenge period and 89 million and 80 million, respectively during the recovery period. After removing low-quality reads, 57 and 64 million of clean reads for the experimental period and 85 and 76 million for the recovery period for Nellore and Angus respectively remained. In the two periods the average value of uniquely mapped reads covering the cattle genome corresponds to 43.18 and 43.82% in Nellore and 44.5 and 46.1% in Angus (Table 1 and Supplementary Table 1).

TABLE 1

Data generated by Reduced Representation Bisulfite Sequencing (RRBS).

	Nellore		Angus
	Challenge	Recovery	Challenge	Recovery
Raw reads	60,091,512	89,179,323	67,901,537	80,179,472
Clean reads	57,699,800	85,587,160	64,763,022	76,551,213
Total mapped reads	25,980,704	38,923,431	30,123,254	36,899,326
Average mapping rate (%)	43.18%	43.82%	44.5%	46.1%

DNA Methylation Patterns

We observed an average genome-wide level of 66.26% CG, 0.44% CHG, and 0.32% CHH methylation (H = A, T or C) in Nellore during challenge and 66.8% CG, 0.5% CHG, and 0.34% CHH methylation during recovery. Values in Angus were 67.7% CG, 0.3% CHG, and 0.2% CHH methylation in the experimental period and 67.14% CG, 0.3% CHG, and 0.22% CHH methylation in the post challenge period (Table 2).

TABLE 2

Genome-wide percentage of methylated CG, CHG, and CHH in Nellore and Angus breed during challenge and recovery period.

	Nellore		Angus
	Challenge	Recovery	Challenge	Recovery
mCG percent (%)	66.26%	66.8%	67.7%	67.14%
mCHG percent (%)	0.44%	0.5%	0.3%	0.3%
mCHH percent (%)	0.32%	0.34%	0.2%	0.22%

Over 66% of the CpG sites were methylated in both breeds during challenge and recovery periods. Heat stress effects did not significantly change DNA methylation of either the CpG and non-CpG sites in both breeds at a genome-wide level (Table 2).

DNA Methylation Level in Different Genome Features

We investigated the DNA methylation level in five different genome features including promoters, 5′UTRs, exons, introns and 3′UTRs, using a sliding window approach. In both breeds methylation level of CG dinucleotides showed a similar pattern during challenge and recovery periods. Promoter regions exhibited a higher variability compared with other features, with the maximum distance between the first and the third quartile. Exon, intron and 3′UTR regions showed the highest methylation level, while the 5′UTR had the lowest average level of methylation (Supplementary Figure 1).

Differentially Methylated Regions (DMRs) and Related Genes (DMGs)

Principal component analysis (PCA) based on genome-wide DNA methylation clearly distinguished Angus from Nellore but not the two different periods (stress and recovery) (Supplementary Figure 2). Differential methylation analysis identified a total of 4662 DMRs. There were 2695 DMRs between Angus and Nellore that did not change between stress and recovery periods and therefore are related to breed differences rather than a response to heat. It is interesting to note that 62% of these were hyper-methylated and 38% hypo-methylated in Angus compared with Nellore, suggesting an overall higher level of metylation in Angus vs Nellore. A total of 1967 windows were differentially methylated in the stress period compared with the recovery period (Supplementary Table 2), 857 in Angus, 896 in Nellore of these only 214 were seen in both breeds, indicating that epigenetic signatures related to heat stress response and recovery are mostly breed specific (Figure 3).

FIGURE 3

Venn diagram of differentially methylated windows in Angus and Nellore (C, Challenge, R, Recovery).

Heat map analysis indicates that in Angus, a large proportion of methylated regions were maintained in both periods possibly indicating difficulty in recovering after a severe stress. Conversely, Nellore seems to respond to stress with a demethylation of several key genomic regions during challenge (Figure 4).

FIGURE 4

Heatmap of methylation level (y axis) in Angus and Nellore populations (x axis) within 1967 DMRs. The blue and azure bars correspond to methylation level of Angus samples during challenge and recovery season and the orange and yellow bars correspond to methylation level of Nellore samples during challenge and recovery season respectively.

The blue and azure bars correspond to methylation level of Angus samples during challenge and recovery season and the orange and yellow bars correspond to methylation level of Nellore samples during challenge and recovery season respectively.

We annotated all 1967 DMRs using the genomic location and the matching annotation information on the genome structure of the cattle genome. The 857 DMRs containing windows identified in Angus overlapped with 351 annotated genes (192 hyper-methylated and 159 hypo-methylated during stress), while in Nellore the 896 windows overlapped with 366 genes (85 hyper-methylated and 281 hypo-methylated during stress). Among these genes 39 were common between Angus and Nellore even though they mapped to different genomic regions. Furthermore the 214 DMRs that were shared between breeds overlapped with 102 genes (among these, 61% were hyper or hypo-methylated in both breeds, while 39% show an opposite trend). In both breeds more than 50% of the differentially methylated regions were located in introns, 31% in the exons, withless than 5% in each of 3′,5′ UTR and promoter regions.

GO and KEGG Analysis of DMGs

Pathway enriched analysis was carried out to identify the metabolic and the immune system pathways affected and identify those that are relevant to a heat stress response. For each breed all DMGs detected were mapped to terms within the GO and the KEGG databases. Both breeds had a similar proportion of DMGs within GO pathways enriched: cellular process (Angus 184; 61.5%; Nellore 188; 62.7%); single-organism process (Angus 73; 24.4%; Nellore 68; 22.7%) and metabolic process (Angus 42; 14%; Nellore 44; 14.7%). KEGG analysis identified 77 pathways in Angus and 78 in Nellore. Among theme, 49 overlap between the two breeds, but a comparison between the top 30 pathways (Figures 5, 6) showed only few cellular functions (a single one in GO enrichment analysis, phosphoric diester hydrolase activity) in common between breeds. It is also interesting to note that shared pathways are enriched only minimally by the same DMGs (32% in Angus and 30% in Nellore), and more than half of these presents opposite behavior, having higher methylation in the stress period and lower in the recovery period in one breed and the opposite in the other. This reflects again the specific response that we previously observed from the Venn diagram analysis.

FIGURE 5

KEGG analysis of DMGs in Angus (left side) and Nellore (right side) in challenge and recovery periods. Peak hight is proportial to the number of differentially methylated genes, color to p value. Negative values indicate hypo-methylation, positive values hyper-methylation during heat stress.

FIGURE 6

Gene Ontology analysis of DMGs in Angus (left side) and Nellore (right side) in challenge and recovery periods. Peak hight is proportial to the number of differentially methylated genes, color to p value. Negative values indicate hypo-methylation, positive values hyper-methylation during heat stress.

Discussion

The tropical environment of the Brazilian Southeast region is characterized by high temperature and humidity and episodes of feed and water scarcity. Nellore cattle have greater resistance to such conditions, surviving, breeding and being productive in the tropics (Valente et al., 2015). In comparison Angus have difficulties adapting to extreme climatic conditions, failing to thrive with tropical heat and humidity, and are susceptible to pathogen challenges. Animal homeostasis in different environments is maintained by modulating metabolic processes to counteract the deleterious effects of stressful conditions (Chovatiya and Medzhitov, 2014). Epigenetic signatures, such as DNA methylation, act as a mediator between changing environmental conditions and animal metabolic response to maintain homeostasis. The present work investigated: (i) changes in DNA methylation that occur in cattle under heat stress, compared with non-stressful conditions and; (ii) differences in methylation in response to heat stress observed in tolerant indicine Nellore versus non-tolerant taurine Angus cattle.

We analyzed genome-wide DNA methylation using the RRBS method. We extracted DNA from blood, as an easily accessible and highly informative indicator of animal response to environmental challenges. The main source of DNA was therefore obtained from white cells, even if we can’t exclude a possible co-extraction of some circulating free DNA, that under stressful conditions may increase its usually low concentration in plasma (Cree et al., 2017; Henriksen et al., 2020). It has been shown that methylation levels at CpG sites in leukocytes are predictors of variability of CpG methylation or may at least partially be considered as proxies of epigenetic processes in other tissues (Ma et al., 2014; Braun et al., 2019; Del Corvo et al., 2020).

Key Pathways Related to the Immune System and Metabolic Processes in Angus and Nellore

In the present study several KEGG gene-sets with differential methyation were identified both in Nellore and Angus breeds in response to heat stress. In Nellore we detected metabolic pathways, circadian entrainment, cGMP-PKG, cAMP and PI3K-Akt signaling pathways that were enriched in differentially methylated genes. Several of these pathways are related to the immune system and inflammatory response. Circadian entrainment pathway has emerged as an internal regulator with the potential to impact disease (Scheiermann et al., 2013). Furthermore, a previous study demonstrated that signaling pathways (cGMP-PKG, cAMP and PI3K-Akt) might form a network with genes related to immuno regulatory mechanisms (Zaccolo and Movsesian, 2007; Wehbi and Taskén, 2016).

In the current study we observed gene enrichment in KEGG metabolic pathways in blood that are known to be involved in heat stress response in other tissues (Rhoads et al., 2010). In particular, pathways for hormone synthesis and/or secretion (cortisol, aldosterone, insulin, renin, growth hormone, GnRH; Figure 6) were all enriched in hypo-methylated genes. Blood cells have been previosly shown to be indicative of methylation status of other tissues in human (Ma et al., 2014). This characteristic is to be further confirmed, however, it would be remarkable, as blood is an easily accessible tissue.

In Angus patterns of response were detected from pathway analysis, that would not have been detectable in a single gene analysis. Several KEGG pathways were enriched in DMGs during heat exposure, including signaling networks, such as MAPK, mTOR, ErbB, PI3k/AKT and pathways related to hepatitis C and other viral infections that can be stimulated in response to heat stress (Danaher et al., 1999). A study of rat skeletal muscle demonstrated that genes in mTOR and PI3k/AKT pathways are activated by a temperature increase (Yoshihara et al., 2013). The ErbB signaling pathway plays a role upstream of MAPK signaling, which affects lipid metabolism (Zhang et al., 2017). Heat stress affects the regulation of lipolysis and the rate-limiting enzymes of lipogenesis (Faylon et al., 2015). In addition, the MAPK pathway is a downstream target of ErbB receptors, which are mediators of cell proliferation, differentiation, apoptosis, and cell motility (Holbro and Hynes, 2004).

Furthermore phosphoric diester hydrolase activity, the only GO pathway shared by both breeds as shown by enrichment analysis, seems to be activated by differentially methylated genes in liver and mammary gland of bull calves in another study as a consequence of maternal conditions of heat stress (Skibiel et al., 2018).

Key Differentially Methylated Genes Related to the Immune System and Inflammatory Response in Angus and Nellore

We identified several genes involved in the biological processes significant for the immune system and inflammatory response. Among them fatty acid elongase 5 (ELOVL5), F2R like thrombin or trypsin receptor 3 (F2RL3), fatty acid desaturase 1 (FADS1), mitogen-activated protein kinase kinasekinase 1 (MAP3K1), phosphodiesterase 5A (PDE5A), netrin 1 (NTN1), RAS p21 protein activator 3 (RASA3) and peroxiredoxin 1 (PRDX) were identified in Nellore. Genes linked to immune defense were in DMRs between stress and recovery periods also in Angus. These included autophagy related 16 like 2 (ATG16L2), calcium voltage-gated channel subunit alpha1 C (CACNA1C)and growth arrest and DNA damage inducible alpha (GADD45A). These genes are expressed in blood cells and have been shown to be differentially methylated also in the brain and are linked to nervous system functions.

In the present study, DMGs identified in Nellore were mainly hypo-methylated in the heat stress period, including in prometer regions, suggesting a pro-active response of this breed to the challenge, to maintain homeostasis.

For example, F2R like thrombin or trypsin receptor 3 (F2RL3) encodes the thrombin protease-activated receptor-4 (PAR-4), which is expressed on the surface of various body tissues, including circulating leukocytes (Xu et al., 1998; Vergnolle et al., 2002). The activation of PAR-4 is involved in leukocyte recruitment, as well as in the regulation of vascular endothelial cell activity (Vergnolle et al., 2002; Kataoka et al., 2003; Leger et al., 2006; Gomides et al., 2012). These pathophysiological events are considered to be the early steps of the vascular inflammatory reaction (Vergnolle et al., 2002; Steinhoff et al., 2005; Leger et al., 2006). Sibbons et al. (2018) showed that in peripheral blood mononuclear cells (PBMCs), the up-regulation of enzymes coded by the fatty acid elongase 5 (ELOVL5) and fatty acid desaturase 1 (FADS1) genes increases the synthesis of n-3 PUFA. Polyunsaturated fatty acids (PUFAs) play key roles in the immune response by acting as substrates for the synthesis of lipid second messengers, including eicosanoids, implicated in cell activation, and for cell membrane biosynthesis. In the present study, these two genes were hypo-methylated during the challenge period, possibly indicating an increase in gene expression. Interestingly, supplementation of PUFAs to the diet can protect animals from the deleterious effects of heat stress on bull sperms by improving quality and motility parameters of fresh semen (Gholami et al., 2011).

A mitogen-activated protein kinase (MAP3K-1) which is member of family of serine/threonine kinases (Pham et al., 2013) was in a region that was hypo-methylated in Nellore during the heat stress period. MAPK signaling is important for T lymphocyte development, homeostasis, and effector responses (von Boehmer, 1990; Kronenberg and Gapin, 2002). Among these, natural killer T cells (NKT) constitute a unique T cell population of the immune system (Kronenberg and Gapin, 2007). In mice, it was demonstrated that a deletion of MAP3K-1 leads to an increased NKT cell infiltration into the liver and a higher degree of liver damage (Suddason et al., 2016). It is therefore not surprising that in Nellore this gene is involved in heat response, since it has been widely documented in cattle that acute thermal stress decreases the number of helper T cells and increases the number of NK cells (Nagai and Iriki, 2001).

Phosphodiesterase 5A (PDE5A) encodes a cGMP-specific phosphodiesterase involved in cell signaling, protein binding, phosphorylation, cell activation and cGMP binding was found to be hypo-methylated in Nellore during heat stress. A hypo-methylated window located within an intronic region of this gene has been described in the above mentioned work of Skibiel et al. (2018), where DNA methylation of bull calves and heifers has been measured after gestation under maternal conditions of heat stress. In the present study another intronic region of the same gene was also hypo-methylated.

Both RAS p21 protein activator 3 (RASA3) and netrin 1 (NTN1) were hypo-methylated in Nellore during heat stress. RASA3 encodes for an enzyme member of the GAP1 family of GTPase-activating proteins. This gene is highly expressed in peripheral blood mononuclear cells during activation of the immune system cells and the inflammatory reaction (Wu et al., 2018). NTN1 is a neuroimmune guidance protein expressed in vascular endothelial cells that regulates inflammatory cell recruitment (Ly et al., 2005; Boneschansker et al., 2016). In human, a significant increase in core body temperature during heat stress is associated with a worsening of cognitive functions of people coping with mental illness (Hämäläinen et al., 2012). Differential methylation of both NTN1 and RASA3 has been reported in blood mononuclear cells of patients affected by the neurodegenerative disease, multiple sclerosis (MS), a disease known to be influenced by high body temperature (Flensner et al., 2011). This suggests that epigenetic regulation of NTN1 and RASA3 may be involved in homoeostatic changes of specific nervous system functions due to heat stress. Multiple sclerosis is associated with systemic sclerosis (SSc) (Pelidou et al., 2007). A recent methylome and transcriptome study carried out in blood cells of patients affected by systemic sclerosis detected five DMGs that were hypo-methylated in heat challenged Nellore in the present study: inositol polyphosphate-5-phosphatase A (INPP5A), Spi-1 proto-oncogene (SPI1), transcription factor 3 (TCF3), in addition to F2RL3 and MAP3K1 (Zhu et al., 2018). This suggest a link between the increase of body temperature and changes in the nervous system, through regulation of inflammatory processes.

Heat stress induced oxidative stress (Altan et al., 2003; Mujahid et al., 2005; Lin et al., 2006) can damage cellular macromolecules and interfere with cell signaling pathways. Peroxiredoxin1 (PRDX1) encodes a member of the peroxiredoxin family, which acts as a peroxidase as well as a chaperone to protect proteins from oxidative damage (Nassour et al., 2016) and plays key roles in innate immunity and inflammation (Knoops et al., 2016). The effect of heat stress on in vitro bovine granulosa cells resulted in incresed expression of PRDX1 in response to oxidative stress and to re-establish cellular homeostasis (Alemu et al., 2018). PRDX1 was hypo-methylation within intronic regions during the challenge period in Nellore.

The level of methylation of Angus samples was higher than those of Nellore, more than half of the significant DMGs were hyper-methylated in the promoter regions during heat stress and therefore gene expression of these genes was likely to be reduced in response to the challenge. This may be associated with the greater vulnerability of this breed to heat stress compared to the Nellore (Bagath et al., 2019). Some of the DMGs were associated with the immune system and were hyper-methylated during heat stress. Autophagy related 16 like 2 (ATG16L2) showed opposite behavior in the two breeds, having higher methylation in the stress period and lower in the recovery period in Angus and the opposite in Nellore. ATG16L2 codes for a protein which plays a role in autophagy (Ishibashi et al., 2011). It was observed that heat stress in pig populations leads to a suppression of autophagy and autophagosomal degradation, which results in the persistence of damaged mitochondria in cells that in turn leads to a dysfunctional intracellular environment (Brownstein et al., 2017). These authors reported a decrease in the autophagy related gene transcripts, which is consistent with our results that show a hyper-methylation in the gene body of ATG16L2. This gene has also been shown to be differentially methylated in blood mononuclear cells of patients affected by multiple sclerosis compared with controls (Kulakova et al., 2016).

Calcium voltage-gated channel subunit alpha1 C (CACNA1C), hyper-methylated in Angus during the heat stress period, belongs to a gene family coding for calcium channels. These channels transport positively charged calcium ions that play a key role in the cell’s ability to generate and transmit electrical signals (Liao et al., 2005). Calcium signaling is essential for T cell activation, tolerance of self-antigens, and homeostasis (Lewis, 2001). In lymphocytes, it has been shown that a signaling cascade leads to an increase of intracellular free calcium that in turn activate immune receptors (Fracchia et al., 2013). White blood cells of mice exposed to acute stress conditions, such as sleep deprivation, progressively loose intracellular Ca2 + and deplete endoplasmic reticulum (ER) reserves, resulting in an impaired immune response (Lungato et al., 2012). Genetic variations in CACNA1C are associated with multiple forms of neuropsychiatric diseases linked to stress and anxiety (Lee et al., 2016).

Growth arrest and DNA damage inducible alpha (GADD45A) was hyper-methylated in Angus during the challenge. Under various stress types, GADD45A maintains genomic integrity in many cell types by surveillance of a DNA-damage (Zhan, 2005) and promoting cell death, cell cycle arrest, and DNA repair (Moskalev et al., 2012). Increased mRNA levels of GADD45A have been reported in hematopoietic stem cells (HSCs), that are required for the continuous regeneration of the blood and the immune system and regulate stress responses (Orkin and Zon, 2008).

Changes in MAPK signaling, a pathway significantly enriched in DMGs in Angus during heat stress, are correlated with defects in innate immune responses of neutrophils and macrophages lacking GADD45 family members (Salerno et al., 2010).

Conclusion

In summary, the present study describes DNA methylation profiles in blood samples of Nellore and Angus steers exposed to the sun during the high temperature season of the Brazilian summer. We detected regions that were differentially methylated that have been previously been associated with activation of immune responses and counteracting the effect of heat stress. Methylation analysis identified specific breed patterns with different methylation responses to stress and recovery. The Nellore appears to actively respond to high temperature by reducing methylation of several key genes and pathways associated with heat response. Angus, however, increases the methylation levels in genes in some of the few pathways that are common between breeds. Further research is also needed to explore the function of demethylation of non-CpG dinucleotides, to improve our knowledge about the biological significance of changes seen as a response to stress.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE161113, GSE161113.

Ethics Statement

The animal study was reviewed and approved by CEUA (Ethics Committee on the Use of Animals) of the Universidade Estadual Paulista (FOA n° 2014-01445). Written informed consent was obtained from the owners for the participation of their animals in this study.

Author Contributions

GN, JG, and PA-M conceived and designed the study. LZ, MM, AU, and GN performed the experiments. EC performed the DNA extraction, libraries preparation, and sequencing. MD and BL analyzed the data. MD, BL, YU, AS, GN, and PA-M prepared the manuscript. All authors reviewed and approved the final manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.