Molecular regulation of whole genome DNA methylation in heat stress response of dairy cows

Abstract

Background

Heat stress seriously affects the production and health of dairy cows and is a key factor limiting the sustainable development of the dairy industry. DNA methylation serves as an important epigenetic regulatory mechanism closely associated with an animal's response to heat stress. However, the specific molecular mechanism of DNA methylation in cows' heat stress response is not fully understood.

Results

In this study, whole genome bisulfite sequencing analysis of blood identified 49861 specific differentially methylated regions corresponding to 7613 differentially methylated genes between spring and summer dairy cows. Among them, 4069 the promoter region of differentially methylated genes were significantly enriched in key biological pathways such as substance transport, reactive oxygen species metabolism, signal transduction, and energy metabolism. By integrating the expression data of 4069 promoter differentially methylated genes, 157 genes were further screened, and their DNA methylation levels were negatively correlated with gene expression. The changes in DNLZ, GNAS, and SMAD5 genes were most significant, and network analysis showed that DNLZ gene has high connectivity in the protein–protein interaction network, indicating its potential key function in heat stress response. Experimental verification shows that under heat stress conditions, the methylation level of CpG islands in the promoter region of DNLZ gene significantly increases, and its methylation level is significantly negatively correlated with gene expression level. The Dual-luciferase reporter assays using constructs containing the DNLZ promoter reporter gene experiment further confirms that promoter methylation significantly inhibits DNLZ transcriptional activity, and the higher the degree of methylation, the stronger the inhibitory effect.

Conclusions

The research results provide new insights into the mechanism of heat stress-related DNA methylation in dairy cows, clarify the key roles of genes such as DNLZ, and provide potential target genes and epigenetic markers for the cultivation of heat-resistant dairy cows.

Background

Dairy cows are important economic animals that provide meat and dairy products for humans [1, 2]. Due to their large body size, thick fur, and underdeveloped sweat glands, dairy cows are cold-resistant but not heat-resistant, making them especially susceptible to the effects of high summer temperatures [3, 4]. To maintain a constant body temperature, dairy cows regulate their thermal balance through physiological mechanisms, such as respiratory rate, as well as behavioural adjustments, including resting patterns. These regulations help balance the heat generated by metabolic processes (including maintenance, exercise, growth, lactation, pregnancy, and feeding) with the heat dissipated into the environment through conduction, convection, and evaporation [5, 6]. Based on reports, the optimal living temperature range for dairy cows is 5–25 °C. When the ambient temperature exceeds the upper limit, dairy cows may experience heat stress, resulting in an increase in core body temperature and negative impacts on their health [7]. As a systemic response of the body to high temperature stimuli, heat stress can disrupt the physiological homeostasis of dairy cows through multi-level mechanisms, thereby impairing their production performance. Firstly, heat stress activates the hypothalamic pituitary adrenal axis, leading to an increase in levels of glucocorticoids (such as cortisol), which in turn triggers immune suppression and metabolic disorders [8]; At the same time, high temperature induces mitochondrial dysfunction and accumulation of reactive oxygen species (ROS), triggering oxidative stress responses, causing protein denaturation, lipid peroxidation, and DNA damage, affecting cellular homeostasis and function [9]. In addition, heat stress elicits a reprogramming of energy metabolism, encompassing attenuated gluconeogenesis, suppressed adipose tissue lipolysis, and compromised milk component synthesis [10]. These changes collectively lead to a decrease in milk production and degradation of dairy quality, significantly weakening their production performance and economic value [11]. Therefore, in the context of global climate change and the growth of intensive dairy farming, heat stress has become one of the most common stress factors in dairy production [12]. According to statistics, in the United States alone, lactating dairy cows suffer economic losses of up to $1.5 billion annually due to heat stress, accounting for approximately 60% of the total heat stress losses in the entire livestock industry [13]. In response to this challenge, researchers have increasingly focused on its molecular regulatory mechanisms in recent years. The role of epigenetic modifications in the regulation of heat stress-related gene expression has gradually emerged, providing new research directions for mechanism interpretation and the establishment of intervention strategies [14].

DNA methylation, as a key epigenetic modification, is characterized by its epigenetic transmission, distinct distribution patterns, spatiotemporal specificity, and reversibility. It plays a pivotal role in regulating gene expression and is involved in various biological processes throughout an organism's development and adaptation to environmental changes [15, 16]. By regulating gene silencing or expression, DNA methylation participates in multiple biological processes including stress response and environmental adaptation, regulation of tissue-specific gene expression, genomic imprinting, cell differentiation and development [17, 18]. During heat stress, DNA methylation is important in regulating gene expression, particularly in tissue-specific responses [15]. For example, in Apostichopus japonicus, increased CpG methylation in promoter regions was associated with transcriptional repression of key heat shock proteins, including hsp70, which plays a critical role in cellular protection against thermal stress [19]. High methylation of CpG islands in the promoter region typically inhibits the transcription of downstream heat shock protein (HSPs) genes, thereby reducing their expression, and potentially impairing the cell's ability to respond to heat stress [20, 21]. In addition, the distribution and variation of differentially methylated regions (DMRs) may lead to selective expression or silencing of specific genes under heat stress conditions [22, 23]. This dynamic regulation of methylation allows cells to quickly respond and adapt to elevated temperature, maintaining protein homeostasis and functional integrity by activating protective mechanisms such as HSPs. Meanwhile, allelic differentially methylated regions (aDMRs) or clustered S-ASMs (sequence-dependent allele-specific methylation) may cause preferential expression or silencing of a certain allele during the heat stress response, thereby influencing overall genomic efficiency in managing stress [24, 25]. As an important epigenetic mechanism regulating gene expression, DNA methylation is gradually becoming a key entry point for analyzing individual differences in heat stress adaptation. In the field of molecular breeding, epigenetic markers such as DNA methylation can not only serve as important biomarkers for identifying heat stress resistant phenotypes, helping to screen germplasm resources with excellent heat adaptability, but also achieve precise intervention in key gene expression by regulating their status, providing a new strategy and technical path for cultivating heat resistant and high-yield dairy cattle lines [26]. Although the important role of DNA methylation in gene expression regulation has been widely recognized, its specific function in the expression of heat stress-related genes in dairy cows still needs to be further explored. Therefore, further exploration of the mechanism of DNA methylation in heat stress in dairy cows is expected to provide new research directions and potential solutions for improving the health and production performance of dairy cows.

The aim of this study is to reveal the changes in DNA methylation in Chinese Holstein cows under heat stress conditions. By constructing a genome-wide methylation map of the blood of Chinese Holstein cows, the regulatory effect of heat stress on DNA methylation was systematically analyzed. We hypothesize that heat stress induces hypermethylation in the promoter region of DNLZ, resulting in transcriptional repression. DNLZ is a gene closely associated with mitochondrial function and oxidative stress regulation, and its downregulation may help reduce energy metabolism and the production of ROS. We propose that this epigenetic modification represents a potential adaptive mechanism in dairy cows to cope with heat stress. The results of this study can not only provide data support for the epigenetic modification mechanism of dairy cows under heat stress conditions, but also provide potential candidate genes and epigenetic markers for molecular breeding of heat-resistant traits.

Methods

Sample collection

In order to study the physiological stress response of Holstein cows in different seasons, 15 healthy Chinese Holstein cows with the same lactation period and parity were selected, and no genetic relationship, from the Beijing Sanyuan Qinglian Treasure Island Ranch as experimental subjects (Additional file 2: Table S1). The feeding conditions and nutritional levels of these dairy cows strictly comply with ranch management requirements, and no animals were euthanized in the study. The study collected 30 blood samples of approximately 16 ml each from the tail vein of these 15 dairy cows in April 2017 (spring, non-stress period, 51 > THI > 58) and July 2017 (summer, heat stress period, 86 > THI > 77) for DNA and RNA extraction and analysis. Blood samples were collected without anesthesia, in compliance with current animal welfare and research laws and regulations in China (approval number: SS-QX-2014–06). These 15 individuals were then randomly divided into three groups, each containing five DNA samples from five dairy cows. Mix the DNA of each individual in equal proportions to prepare a composite DNA sample. Rename these mixed samples, with spring samples labeled as Spr1, Spr2, and Spr3, and summer samples labeled as Sum1, Sum2, and Sum3, and then use these 6 mixed DNA samples for whole genome methylation sequencing (WGBS). Due to the relatively small sample size (n = 15), DNA samples from five dairy cows per group were pooled to ensure sufficient sequencing depth and reduce technical noise. This pooling strategy was chosen to detect consistent methylation changes induced by heat stress at the population level. However, the pooling strategy itself has certain limitations, as it may obscure biologically meaningful inter-individual differences in methylation and reduce the ability to detect extreme phenotypes or low-frequency methylation events. To address this limitation, key DMGs identified from WGBS were further validated at the individual level using bisulfite-PCR (BSP) and gene expression analysis (RT-qPCR), allowing for confirmation of heat stress-induced methylation changes.

Whole-genome Bisulfite library preparation and sequencing

Genomic DNA extraction was performed using the Blood DNA Midi Kit (Chengdu Forensic Biotechnology Co., Ltd., China). The samples were subjected to bisulfite conversion using the DNA Methylation-Gold™ Kit. The sequencing libraries were constructed using the Illumina TruSeq DNA Methylation Library Preparation Kit. 0.5% λDNA was added to the library as quality control. WGBS sequencing was performed using ultrasonic fragmentation technology, and the fragmented DNA samples were subsequently screened and purified using LC Beads. The samples were incubated in a PCR machine under the following conditions: 37 °C for 15 min, followed by denaturation at 95 °C for 2 min. The double-strand extension reaction was then carried out with a program set at 98 °C for 1 min, 62 °C for 2 min, and 65 °C for 5 min. Purification was performed using LC Beads 1.2 to complete the 3'adapter ligation reaction. The 5'adapter ligation reaction was then performed, with incubation at 25 °C for 15 min, followed by purification of the product using LC Beads 1.0. During the library amplification stage, amplification primers with index sequences were added, with the PCR conditions set to an initial denaturation at 98 °C for 30 s, followed by an appropriate number of amplification cycles. The final amplified product was purified using LC Beads 0.8. During the analysis, the bisulfite conversion rate was calculated using lambda DNA, and the library construction was completed by combining fragment selection and PCR amplification (insert fragment range: 200–500 bp). Finally, the HiSeq 4000 platform with a 150 bp paired end sequencing strategy (150PE) was used for sequencing, with a sequencing depth of 30 ×.

WGBS data quality control and alignment

The raw data from high-throughput sequencing is stored in FASTQ format as the result file after base calling of image files, which often contains low-quality data. To ensure the accuracy of subsequent bioinformatics analyses, quality control and evaluation of the raw data are required to obtain clean data. The sequencing quality, adaptor content, and base distribution of the sequencing data were evaluated using FastQC (v0.11.0) with default parameters. Based on the evaluation results, Fastp (v0.21.0) [27] was used to filter the raw data. The quality control principles were as follows: (1) filter reads containing adaptors; (2) filter reads with a length less than 50 bp; (3) filter reads with unknown base (N) numbers greater than 10; and (4) filter low-quality reads (where more than 20% of the total reads have a quality value of Q ≤ 10). The clean data obtained after quality control was re-evaluated using Fastp to determine the effectiveness of data quality control.

Sequence alignment between the clean data after quality control and the Bos taurus reference genome (ARS-UCD1.2) was performed using Bismark (v0.20.1) [28] in combination with Bowtie2. The Bos taurus genome sequence was downloaded, extracted, and indexed prior to alignment. After alignment, Bismark was used to perform deduplication, statistical analysis, and extraction of methylation sites. The generated SAM files were then converted into BAM format using Samtools for storage and downstream analysis. Finally, Circos (v0.63) was used to generate a genome-wide methylation distribution plot.

Differentially methylated regions analysis and annotation

The identification and annotation of DMRs were conducted using the R package MethylKit (v1.1.7) [29]. DMRs were identified using a sliding window approach with a 1000 bp window and a 500 bp overlap method. The methylation level of each region was calculated as the ratio of methylated cytosines to the total cytosines detected in that region. Statistical significance was assessed using a t-test, and p-values were adjusted for multiple testing using the Benjamini-Hochberg (BH) false discovery rate (FDR) correction, setting the FDR threshold at < 0.05. To ensure biological relevance, DMRs are annotated according to their position in the genome, and only DMRs located in functional regions such as promoters, gene bodies, and enhancers are included in subsequent analysis. The identified DMRs were annotated with genomic location using ChIPseeker (v1.20.0) [30], thereby enabling the identification of differentially methylated genes (DMGs). The methylation level was calculated as the ratio of the number of reads supporting methylation to the total number of reads (i.e., the number of reads supporting methylation plus the number of reads not supporting methylation).

DMGs were compared against functional databases, including Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG), to obtain detailed gene annotations and conduct functional enrichment analysis. GO and KEGG enrichment analyses were conducted using the DAVID (https://david.ncifcrf.gov/) and KOBAS (http://kobas.cbi.pku.edu.cn/) online databases, with the significance threshold set at p < 0.05. Additionally, Bioconductor packages in R, including ggplot2, ComplexHeatmap, and circlize, were used for data visualization, enabling statistical plotting, clustering heatmaps, and circular genome representation, respectively. For protein–protein interaction network analysis, the STRING database (https://www.string-db.org/) was utilized.

CpG Island and transcription factor prediction

The sequence of the DNLZ gene was obtained from the NCBI database (ID: 100,848,156). Transcription factor binding regions were predicted using the PROMO online software (http://alggen.lsi.upc.es/cgi-bin/promo_v3/promo/promoinit.cgi?dirDB=TF_8.3). The CpG islands in the promoter region of the DNLZ gene were predicted using the Methprimer online software (http://www.urogene.org/cgi-bin/methprimer/methprimer.cgi) [31]. The identification criteria for CpG islands included a length greater than 300 bp, GC content greater than 50.0%, and an observed/expected value ratio (O/E) greater than or equal to 0.6.

Verification of gene expression

All RT-qPCR primers were designed with NCBI Primer Blast Online Tool (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) and synthesized by Shenggong Bioengineering (Shanghai) Co., Ltd, the detailed information was shown in Table 1. Verify the expression levels of key genes using RT-qPCR and set negative and positive controls in real-time fluorescence quantitative PCR experiments to ensure the effectiveness of the experiment. The detection of internal reference genes and target genes is carried out using at least three technical replicates, and the reaction is performed in a 96 well plate. The reaction system consists of 2 × SuperReal Premix Plus 5.0 μL, 0.3 μL of 10 μ M forward and reverse primers, 1.0 μL of cDNA, and water added to 10 μL. The amplification conditions are: pre denaturation at 95℃ for 15 min, followed by denaturation at 95℃ for 10 s, annealing at 60℃ for 20 s, and extension at 72℃ for 30 s, for a total of 40 cycles. After completion, perform a final extension at 72℃ for 32 s. Addition operations for all samples are performed on ice to ensure the accuracy of the reaction. The relative expression levels of target genes were calculated using the 2^-ΔCt method. The mean ΔCt value of each group was used to determine relative expression, and intergroup comparisons were conducted using an independent two-tailed Student’s t-test. Statistical analyses were performed using GraphPad Prism software (version 9.5), with a P < 0.05 considered statistically significant.

Table 1 Primer information of RT-qPCR

Bisulfite conversion and nested PCR

Using nested PCR to amplify target fragments of DNA modified with bisulfite: Design BSP nested PCR primers (product length < 700 bp, primer 23–30 bp, original sequence containing ≥ 4 Cs, avoiding CpG sites as much as possible, no A at the 3'end, and no ≥ 3 consecutive G/Cs), use the modified DNA as a template for the first round of PCR amplification, and then dilute the product moderately to serve as the template for the second round of PCR. The PCR reaction is carried out according to the preset system and conditions, as detailed in Additional file 2: Table S2. The amplified products were examined by 1% agarose gel electrophoresis, and the qualified products were sent to Biotechnology (Shanghai) Engineering Co., Ltd. for product purification and sequencing. The sequencing results were compared with the transformed reference sequence, and the methylation level of C at the CPG site was calculated based on the peak heights of C and T = C peak height/(C peak height + T peak).

Construction of promoter luciferase reporter plasmid

The DNLZ gene sequence (ID: 100848156) was retrieved from the NCBI database. Based on the WGBS and BSP results, a 648 bp promoter region spanning from − 596 to + 51 bp (relative to the transcription start site, ending at the start codon) was synthesized, with unique MluI and HindIII restriction sites added at both ends. The fragment was cloned into the pGL3-basic vector to construct the pGL3-DNLZ luciferase reporter plasmid, and pRL-TK (Renilla luciferase) was used as the internal control plasmid. The synthesis of target fragments and construction of recombinant plasmids were completed by Beijing Haichuang Science and Technology Co., Ltd.

Construction of PGL3-DNLZ recombinant plasmid and validation of its methylation modification

Double digestion identification was performed using the unique MIUI and HindIII cleavage sites on the PGL3-DNLZ plasmid to verify the successful construction of the target recombinant plasmid. Large scale cultivation of Escherichia coli transformed into PGL3 basic, PGL3-DNLZ, and pRL TK using LB liquid medium supplemented with ampicillin. Extract plasmids from 100 mL of bacterial solution using the EndoFree Midi Plasmid Kit (DP108). Adopting M.Sss I methyltransferase converts all CpG background C in DNA sequence to methylated C. The purified plasmid is recovered using TIANquick Midi Purification Kit (DP204), and then subjected to single enzyme digestion using HpaII enzyme to verify successful methylation.

Cell transfection and dual luciferase activity detection

Instantaneous co-transfection using HEK-293 T cells and Lipofectamine 2000: Inoculate 1 × 10 ^ 5 cells/well in a 24 well plate, discard the original medium when cultured to 80% density, replace with 400 μ L Opti MEM I for 3 h, and then dilute 2 μ L Lipofectamine 2000 with 500 ng plasmid (experimental plasmid 450 ng + internal reference pRL TK 50 ng) in each well with 50 μ L Opti MEM I and incubate at room temperature for 5 min. Mix the two at room temperature for 20 min, slowly drop them into the well, and incubate at 37℃ and 5% CO2 for 4–6 h before replacing with DMEM containing 10% FBS completely. Culture medium, detect luciferase activity 24–48 h later. During detection, cells were lysed with 1 × PLB, centrifuged to obtain the supernatant, and LAR II was added to measure firefly luciferase. Stop&Glo was then added to measure sea kidney luciferase, and the relative luciferase activity was calculated based on the ratio of sea kidney luciferase to firefly luciferase.

Results

The methylation patterns across the whole genome in spring and summer dairy cows

Among the 15 dairy cows tested, 421213695 and 448911171 cytosine sites were detected in the spring and summer groups, respectively, with a bisulfite conversion rate of 99.78% (Additional file 2: Table S3). Among them, 17636035 and 18238922 loci were methylated, and the proportion of methylated cytosine (mC) was about 4% (Additional file 2: Table S4). The average methylation levels in the whole genome of dairy cows in spring group were 75.64% for mCpG, 0.30% for mCHG, and 0.28% for mCHH, while the individuals in summer group exhibited 75.68% for mCpG, 0.31% for mCHG, and 0.28% for mCHH (Fig. 1A). The methylation levels between the spring and summer groups were not significant, indicating that heat stress had a relatively small impact on the overall methylation levels of dairy cows. Further calculation of the distribution ratios of different methylation types revealed that CpG methylation had the highest proportion, exceeding 90% in both groups, while CHH and CHG methylation accounted for approximately 6% and 2%, respectively (Fig. 1B and C).

Fig. 1

The Genome-wide methylation patterns in dairy cows (A). The average ratio of DNA methylation types of dairy cows in spring (B) and summer (C). The whole genome methylation map of dairy cows in spring (D) and summer (E): the outermost circle of the whole genome methylation map is the chromosome length scale; the CpG (purple), CHG (blue), and CHH (green) methylation types of chromosome intervals are arranged from the outside to the inside (darker colors indicate higher methylation levels); the innermost circle represents the number of genes in the corresponding interval (darker colors indicate more genes)

Based on the whole genome methylation map and chromosome methylation distribution map of spring or summer dairy cows (Fig. 1D and E, Additional file 1: Figure S1), the methylation levels of each chromosome both in spring and summer dairy cows were similar. Among the CpG methylation types, the methylation level of chromosome 24 was the highest while the methylation level of chromosome 27 was the lowest. In both CHH and CHG types, the methylation level of chromosome 13 was the highest while the methylation level of chromosome X was the lowest.

Identification and annotation of differentially methylated regions in spring and summer dairy cows

The study identified and analyzed the genomic distribution characteristics of DMRs in dairy cows under heat stress. A total of 49,861 DMRs were identified in the summer heat stressed group compared to the spring non-stressed group. These DMRs are predominantly hypomethylated in terms of methylation levels, with 34,784 hypomethylated DMRs, which is 2.3 times more than the 15,077 hypermethylated DMRs (Fig. 2A, Additional file 2: Table S5). Further annotation and statistical analysis of the locations of DMRs in the genome revealed that the majority of DMRs are located in the intergenic region (30,961, 62.1%), followed by the exon region (8,260, 16.6%) and intron region (5,854, 11.7%), with the lowest number of DMRs involving the gene promoter region (4,789, TSS upstream 2 kb, 9.6%). Further, most DMRs (63.42%) in gene promoter regions showed an hypermethylated, which is opposite to the methylation patterns observed in other genomic regions. This suggests that heat stress may lead to an increase in methylation levels in the promoter regions of genes in dairy cows.

Fig. 2

DMRs and DMGs in the genomes of spring and summer dairy cows: A Distribution statistics of DMRs across different genomic regions; B Volcano plot of DMRs, with blue representing hypomethylated DMRs and red representing hypermethylated DMRs in the summer group relative to the spring group; C Cluster heatmap of DMR methylation for the top 100 promoter DMGs, Color changes represent different degrees of methylation; D PPI network analysis of DMGs in the heat shock protein (HSP) family. Colored nodes indicate core proteins and their direct interactors, white nodes represent indirectly linked proteins. Lines denote functional associations, including experimentally validated interactions and bioinformatics-based predictions

These identified DMRs between the spring and summer groups correspond to 7,613 genes, known as DMGs, of which, 4,069 DMGs involve differential methylation in the promoter region. DMGs include 214 non-coding RNAs, 36 coding rRNAs (34 coding 5S rRNA, 2 coding 5.8S rRNA), and 3090 protein coding genes. Figure 2C shows the heatmap of the top 100 DMRs in the promoter region. The top ten protein coding genes are RTL5, CDC42BPB, DNLZ, TNNC2, FZD5, ZNF613, GLB1L3, GPR27, KCNV1, and HSPA1L. Additionally, among these 4,069 promoter DMGs, 19 heat shock protein family genes [such as HSPA1 A, HSPA1L, DNAJC16, etc.; see Additional file 2: Table S6] showed significant changes in promoter methylation levels (FDR < 0.05). Figure 2D shows the protein–protein interaction (PPI) network analysis results of these 19 genes, with HSPA1 A, HSPA1L, and HSPA6 proteins being key nodes that interact with multiple other proteins.

GO enrichment and pathway analysis of promoter region DMGs in spring and summer dairy cows

Promoter methylation plays a crucial role in the regulation of gene expression, we focused on the impact of differential methylation in promoter regions in current study. Annotation and functional enrichment analysis of the 4069 DMGs in promoter regions revealed that these DMGs are involved in 3341 GO terms and 70 significantly enriched KEGG pathways. Figure 3A shows the top 50 significantly enriched GO terms (P < 0.05), and the predominantly enriched in the GO terms were"regulation of transport"and"positive regulation of transport"containing 69 and 39 genes, respectively. Also, terms such as"reactive oxygen species metabolic process","regulation of GTPase activity"and"regulation of ATPase activity"are closely related to heat stress. The significant enrichment of the term"regulation of DNA binding in transcription regulatory region"suggests that DNA methylation may respond to heat stress by influencing gene transcription level. Figure 3B shows the 50 significantly enriched KEGG pathways (P < 0.05), among which"oxidative phosphorylation","longevity regulation","cortisol synthesis and secretion"and"AMPK signaling pathway"are highly related to heat stress. Further, the heat shock protein family such as HSPA1 A, which were identified as the promoter DMGs in the summer heat stressed group compared to the spring non-stressed group (Fig. 2D), is directly enriched in the"thermogenesis"pathway (Fig. 3B), further supporting their association with heat stress.

Fig. 3

GO and KEGG pathway enrichment analysis of DMGs in promoter region: A GO term enrichment analysis of DMGs. The y-axis represents the pathway names, and the x-axis represents the Count. B KEGG pathway enrichment analysis of DMGs. The y-axis represents the pathway name, and the x-axis represents the gene number. The color of the dots corresponds to different ranges of P-values, darker color represent smaller P-values, and the size of the dots represents the number of genes

Identification of key DMRs and DMGs in promoter in dairy cows under heat stress

To further investigate the relationship between promoter methylation and gene expression under heat stress. This study combined the 4,069 promoter DMGs from WGBS with 1,997 differentially expressed genes (DEGs, FDR < 0.05 and |log2 FC|> 1) from the results of transcriptome sequencing in our team (Fig. 4A). Ultimately, 264 common genes between promoter DMGs and DEGs, including GNAS and DNLZ, were identified as differentially methylated and expressed genes (DMEGs) (Fig. 4B). Among these 264 DMEGs, 157 exhibited methylation changes that were inversely related to their gene expression levels [see Additional file 2: Table S7]. Specifically, the methylation levels of 97 genes, including GNAS and DNLZ, increased while their expression levels decreased. Conversely, the methylation levels of 60 genes, including SMAD5, decreased while their expression levels increased compared to the spring group. Notably, the DNLZ gene demonstrated the smallest differential FDR value (4.5 × 10⁻²⁰) among the157 DMEGs. The SMAD5 gene is involved in the most heat stress-related GO terms, such as metal ion binding. Transcriptome sequencing data showed that the GNAS gene exhibited the most significant downregulation in the summer group. Figure 4C, D illustrated the specific methylation and expression patterns of GNAS, DNLZ, and SMAD5.

Fig. 4

Identification of key DMRs and DMGs in promoter in dairy cows under heat stress (A) The volcano plot of DEGs; B Venn diagram of promoter DMGs and DEGs; C The methylation levels of DMRs in the promoters of GNAS, DNLZ, and SMAD5; D The expression levels of GNAS, DNLZ, and SMAD5 using TPM (number of transcripts per million reads). E PPI Network interaction of DMEGs with negative regulation of expression by promoter methylation. Note: In (C) and (D), arrows indicate the direction of change, upward arrows indicate an increase, and downward arrows indicate a decrease in summer expression; *** FDR < 0.001

The results of the PPI network interaction analysis involving 157 DMEGs are shown in Fig. 4E. RPS18, MRPL36, and ATP5IF1 are identified as key interaction nodes. DNLZ, along with genes such as PHPT1 and ATP8, participates in mitochondrial functions and is involved in GO terms like DNA binding transcription factor activity and protein stability (Fig. 3A). Notably, DNLZ exhibited significant connectivity in the PPI network, unlike the GNAS, DNLZ, and SMAD5 genes. Given the importance of methylation differences, inter gene network connectivity strength, and the downregulation of gene expression, DNLZ was selected as one of key DMRs and DMGs associated with heat stress for further studies to understand its role in this study.

The DNLZ expression was suppressed as its promoter methylation increased in dairy cows under heat stress

Methprimer software was used to predict the CpG islands in the DMR of DNLZ gene (−596 bp and + 404 bp). A 395 bp CpG island is shown between −436 bp and −41 bp in the promoter region of the DNLZ gene (Fig. 5A). Due to primer design limitations and the clarity of sequencing peaks, we only analysed the methylation levels in the region from −323 bp to −61 bp (designated as DNLZ1). PROMO prediction shows that the DNLZ1 region contains 39 CG sites and potential binding sites for important transcription factors such as HIF1, TFIID, and HES1, and the 17 and 21 st CG located in HIF1 transcription factor binding region are significantly methylated in summer dairy cows (P < 0.05, Fig. 5B, D). These results suggest that the expression of the DNLZ gene may be regulated by DNA methylation.

Fig. 5

The DNLZ expression was suppressed as its promoter methylation increased in dairy cows under heat stress. A Predicted CpG island in the DNLZ DMR sequences; B Nucleotide sequence of DNLZ1(part of DNLZ promoter), with 1 to 39 indicating the cytosine numbering of CG sites, and single underlines representing predicted transcription factor binding site; C Detection of mRNA expression levels of DNLZ, MeCP2, DNMT1, DNMT3 A, and DNMT3B. D The methylation level of each CG site in DNLZ1. E The overall methylation level of CpGs in DNLZ1 in dairy cows. Note: * P < 0.05, ** P < 0.01, *** P < 0.001

To verify the changes in expression level of DNLZ gene in heat stressed dairy cows, blood samples were collected from 24 dairy cows in spring and summer, and RNA was extracted from the white blood cell layer. RT-qPCR was used to detect the gene transcription levels. Also, the expression of methyltransferases DNMT1, DNMT3 A, DNMT3B, and the methyl-binding protein MeCP2 were detected in 24 dairy cows in spring and summer (Fig. 5C). Compared with the spring group, the expression level of the DNLZ gene was significantly downregulated in summer group (P < 0.01), consistent with the trend observed in the transcriptome sequencing results (Fig. 4D), indicating that heat stress affects the transcription of the DNLZ gene. Furthermore, the expression levels of MeCP2, DNMT3 A, and DNMT1 were significantly increased after heat stress (P < 0.05), while the change in DNMT3B was not significant but showed an upward trend. This suggests that heat stress may increase the DNA methylation levels of genes by upregulate the expression of methyltransferases and methyl-binding proteins. At the same time, we obtained the cytosine methylation status of the 38 CpG sites in DNLZ1 with BSP method (Fig. 5D, E). The results showed that heat stress significantly increased the methylation levels of the 38 CG sites in the DNLZ1 region, consistent with the trend of WGBS detection results (Fig. 4D).

The DNLZ expression was suppressed as its promoter methylation increased in Mac-T cells under heat stress

To verify the methylation changes of DNLZ1 region at the cellular level under heat stress, this study used Mac-T cells as a model and subjected them to 39℃ for different times to detect the methylation levels at each time point. The results showed that after heat treatment at 39℃ for 48 and 72 h, the overall methylation level of the CG locus of the DNLZ1 region significantly increased (P < 0.05) (Fig. 6A), while the mRNA expression level significantly decreased (Fig. 6B). In the same time, four methylation related enzymes were upregulated to varying degrees after heat stress treatment of 48 and 72 h (Fig. 6C, D). Overall, heat stress significantly increased the promoter methylation level of DNLZ gene and inhibited its expression.

Fig. 6

The DNLZ expression was suppressed with increased methylation of CpG island in its promoter in Mac-T cells under heat stress. A The overall methylation level of CpGs in DNLZ1 in Mac-T cells. B The expression level of DNLZ gene in Mac-T cells after 24 and 48 h of heat stress treatment. C and D The expression levels of four methylases in Mac-T cells after 48 and 72 h of heat stress treatment. * Indicates P < 0.05, ** indicates P < 0.01

The methylated promoter of DNLZ gene suppressed its expression in vitro

To further clarify the regulatory effect of methylation in DNLZ gene promoter on gene expression, this study constructed a PGL3-DNLZ promoter (ID: 100848156; −596 ~ + 51 bp) luciferase reporter plasmid and conducted in vitro methylation treatment and co transfection experiments (Fig. 7A). The constructed plasmid was identified by MluI and HindIII double digestion, and the digestion product was analysed by agarose gel electrophoresis. The band was clear, and the fragment size was in the range of 600–1000 bp, which was consistent with the theoretical expectation [see Additional file 1: Figure S2]. The PGL3-DNLZ luciferase reporter plasmid synthesized in vitro was treated with M.Sss I (Spiroplasma CpG-Specific DNA Methyltransferase Spiroplasma I) to induce CpG site methylation, and the fully methylated plasmid was named Met-PGL3-DNLZ. To verify the degree of methylation, different enzyme concentrations and reaction times were grouped, including: 0U + 0 h, 2U + 0.5 h, 4U + 0.5 h, 2U + 1 h, 4U + 1 h, 2U + 2 h and 4U + 2 h. The results showed that groups 4U + 0.5 h, 2U + 1 h, and 2U + 2 h exhibited partial methylation, groups 4U + 1 h and 4U + 2 h achieved complete methylation, while group 2U + 0.5 h showed minimal methylation (Fig. 7B) [see Additional file 1: Figure S3].

Fig. 7

The methylated promoter of DNLZ gene suppressed its expression in vitro. A Schematic diagram of PGL3-DNLZ promoter luciferase reporter plasmid structure. B Partial methylation identification, where Lane 1 ~ Lane 7 are single enzyme cleavage products of DNLZ plasmids and HpaII under different conditions: Lane 1 is 0U + 0 h, Lane 2 is 2U + 0.5 h, Lane 3 is 4U + 0.5 h, Lane 4 is 2U + 1 h, Lane 5 is 4U + 1 h, Lane 6 is 2U + 2 h, Lane 7 is 4U + 2 h, Lane 8 is plasmid, and M is DNA marker. C The relative activity detection results of dual luciferase after complete methylation of PGL3-DNLZ plasmid. D Detection results of dual luciferase relative activity of PGL3-DNLZ plasmid at different M.Sss I enzyme concentrations and treatment times. Note: The complete gel image can be found in Additional file 1: Figure S3. This image represents the full experimental process, with no parts cropped before the electrophoresis experiment

To further analyze the effect of different methylation levels on the transcriptional activity of gene promoter, this study co transfected HEK-293 T cells with methylated and unmethylated PGL3-DNLZ plasmids and pRL-TK reference plasmid and detected the relative transcriptional activity through dual luciferase reporter assay. The results showed that the relative luciferase activity of unmethylated PGL3-DNLZ was significantly higher than that of empty vector PGL3 basic (P < 0.001) (Fig. 7C), indicating the successful construction of the DNLZ promoter luciferase reporter plasmid. The relative luciferase activity of fully methylated Met-PGL3-DNLZ was significantly downregulated by 96.91% (P < 0.001) compared to the unmethylated group, indicating that promoter methylation significantly inhibits gene transcription expression. Further analysis of the effect of different degrees of methylation on gene transcriptional activity (Fig. 7D) showed that the average activity of PGL3-DNLZ luciferase decreased by 40.13% in partially methylated groups (such as 2U + 0.5 h group). However, as the enzyme dosage increased and the reaction time prolonged, the degree of methylation deepened, and gene transcriptional activity gradually decreased. Under the 4U treatment condition, transcriptional activity basically disappeared. In summary, these results indicates that the methylation level of the DNLZ gene promoter has a significant regulatory effect on its transcriptional expression, and the degree of methylation is positively correlated with the intensity of transcriptional inhibition. In a highly methylated state, gene transcriptional expression is significantly reduced or even completely inhibited.

Discussion

With the development of global warming and intensive farming, heat stress has become one of the major challenges facing the dairy industry [4, 32]. The ability of dairy cows to cope with heat stress determines the stability of their milk production. Therefore, maintaining the temperature regulation and metabolic balance of dairy cows is undoubtedly crucial [3]. Although DNA methylation has been extensively studied in mammals, the systematic landscape of DNA methylation during heat stress in dairy cows is still largely unknown [33]. This study generated high-quality WGBS data from the blood of fifteen Holstein cows, producing a total of approximately 3.52 billion reads. Based on the WGBS sequencing data, about 4% of the C sites in the entire dairy cow genome were methylated. The average methylation levels in CpG sites exceeded 75% (75.64% in the spring group and 75.68% in the summer group), while the methylation levels in CHG and CHH sites were lower (0.30–0.31% for CHG and 0.28% for CHH). Our results are consistent with the DNA methylation pattern commonly observed in mammals, where CpG sites are predominantly methylated, while non-CpG methylation remains at low levels [34,35,36]. Similar patterns were observed in pigs, with CpG methylation levels ranging from 67.62% (normal group) to 67.54% (heat stress group), and CHG/CHH methylation levels below 1% [37]. The intergenic region had the highest methylation levels while the promoter region had the lowest ones, which indicates that heat stress impact methylation levels differently across various regions of the genome [38]. Given the complexity and size of the mammalian genome, various functional regions of genes may respond individually or synergistically to changes in the organism's phenotype [39, 40]. This process involves complex mechanisms that require further research and exploration.

Previous studies have shown that heat stress significantly affects the metabolic pathways of poultry, pigs, and ruminants, leading to reorganization of energy metabolism, such as the regulation of carbon metabolism and glycolysis pathways [41,42,43]. During this process, heat stress can trigger oxidative stress, resulting in increased levels of ROS in the body, which can damage cellular function [44, 45]. In this study, we used blood samples from dairy cows in the summer and spring groups to identify 1,997 DEGs and 7,613 DMGs. Given that promoters are key genomic elements regulating gene expression, their methylation levels are often associated with transcriptional repression [46, 47]. After annotating and conducting functional enrichment analysis on 4,069 DMGs in the promoter region, we found that these DMGs were significantly enriched in the GO terms"transport regulation"and"forward transport regulation”. This finding is consistent with research on Apostichopus japonicus, which demonstrated that heat stress induces significant changes in DNA methylation, particularly in genes associated with transport processes [48]. This phenomenon occurs because organisms rapidly adapt to heat stress by regulating substance transport genes through DNA methylation, thereby maintaining cellular homeostasis, optimizing energy metabolism, and protecting cells from high-temperature damage [26]. Meanwhile, DNA methylation, as a relatively stable epigenetic marker, partially environmentally induced methylation patterns can be preserved and transmitted to offspring during meiosis, providing organisms with a rapid adaptation strategy that does not require changes to the genome sequence [49]. For example, when parents experience heat stress, their offspring may exhibit enhanced physiological adaptability, indicating that DNA methylation plays an important role in environmental adaptation and evolutionary processes [49]. This regulatory mechanism optimizes resource allocation and improves the survival ability of organisms in harsh environments, which is one of the important adaptive strategies formed during the evolutionary process [50]. In addition, GO analysis showed that terms such as"reactive oxygen species metabolism process","GTPase activity regulation"and"ATPase activity regulation"were closely related to heat stress responses of dairy cows. The process of reactive oxygen species metabolism is vital as it helps cells manage oxidative stress induced by heat. This is achieved by regulating intracellular ROS levels, controlling both the generation and clearance of ROS, and maintaining redox balance. These mechanisms are essential for preventing oxidative damage and ensuring cellular integrity during stress conditions [51, 52]. While the regulation of GTPase activity is crucial to signal transduction and cell cycle control. Proper GTPase function ensures that cells can effectively respond to and navigate stressful environments, enabling appropriate pathways to activate protective mechanisms [53]. ATPase activity regulation may optimize energy metabolism, participate in processes such as maintaining cell membrane potential, facilitating ion transport, and aiding in protein folding. These functions are essential for ensuring that cells have adequate energy reserves to sustain their physiological functions and promote survival under stress [54, 55]. Additionally, the significant enrichment of terms such as"DNA binding regulation in transcriptional regulatory regions"suggests that cells may modify gene expression through the activation of specific transcription factors, such as heat shock factor (HSF) [56]. This factor can induce the expression of genes belonging to the HSP family. The HSPs play a critical role in refolding or repairing damaged proteins, thereby mitigating protein denaturation and aggregation that can occur due to heat stress [57]. In our analysis, significant changes in DNA methylation levels were observed in the promoter regions of 19 HSP family genes (FDR < 0.05), further indicating the important role of DNA methylation in heat stress response. DNA methylation affects gene expression levels by regulating the activity of gene promoter regions. Under heat stress conditions, this epigenetic mechanism may help cells adapt to environmental changes and enhance their survival ability by regulating the expression of heat shock proteins [58]. Overall, our results indicate that, compared to the spring dairy cows, the summer dairy cows experienced substantial alterations in reactive oxygen species metabolism, signaling pathways, energy metabolism, and the overall heat stress response. These findings are consistent with previous research, which has documented similar changes in response to heat stress in various organisms [59, 60].

By integrating DNA methylation and transcriptome data, we identified a total of 264 overlapping DMGs and DEGs. Among these, 157 DMEGs exhibited changes in methylation that were opposite to the direction of gene expression, whereas 107 showed changes in the same direction. This observation indicates that the relationship between promoter DNA methylation and gene expression regulation under heat stress in dairy cows is complex and cannot be simply categorized as a positive or negative regulatory interaction. This complexity may stem from the multi-level effects of DNA methylation on gene expression: on the one hand, methylation changes may regulate mRNA transcription levels by affecting the accessibility of promoter regions; On the other hand, methylation may indirectly affect protein levels, such as regulating miRNA expression, thereby altering protein degradation or stability after translation [61]. Factors such as chromatin accessibility and other epigenetic modifications, including histone modifications, may also influence this relationship [62]. Despite this complexity, DNA methylation remains a predominant mechanism of negative gene expression regulation and continues to be a significant focus of genomic regulation studies [37]. This negative correlation was further validated across multiple vertebrate species, including humans, frogs, and pufferfish, utilizing WGBS and RNA sequencing technologies [63]. Given this context, our study emphasizes the 157 DMEGs that exhibit differential promoter methylation with opposite directions of methylation and expression changes. Among these genes, we observed a significant increase in the methylation level of the DNLZ gene promoter following heat stress. This suggests that DNLZ plays a pivotal regulatory role in the response to heat stress in dairy cows and may serve as a potential molecular marker. The DNLZ gene located on chromosome 11 is associated with mitochondrial function [64] and crucial for the assembly and functional maintenance of mitochondrial complex I [65, 66]. Research has shown that heat stress often correlates with increased oxidative stress, particularly elevated ROS levels [67]. The excessive accumulation of ROS can trigger lipid peroxidation, protein oxidative modification, and DNA damage in cell membranes, thereby impairing mitochondrial function, inducing cell apoptosis or necrosis, and exacerbating inflammatory reactions [68, 69]. The electron leakage of the mitochondrial electron transfer chain (ETC) is considered the main source of heat stress-induced ROS overproduction [70]. In poultry models, Mujahid et al. found that acute heat stress significantly increased the level of superoxide (O ₂⁻) in chicken skeletal muscle mitochondria, indicating that heat stress directly affects mitochondrial redox homeostasis [71]. In addition, Hall et al. used electron paramagnetic resonance spectroscopy (EPR) to detect a significant increase in semiquinone free radical content in the portal vein blood of heat stressed rats, which further confirms the abnormal accumulation of ROS in high-temperature environments [72]. In our study, we found a significant increase in the overall methylation level of CpG island within the promoter of the DNLZ gene, both at the individual and cellular levels. After methylation treatment of the PGL3-DNLZ promoter luciferase reporter plasmid, we found that the methylation level of the DNLZ gene significantly affected its expression level. We hypothesize that this downregulation mechanism may represent an adaptive response of cells to oxidative stress induced by heat stress. Specifically, the downregulation of DNLZ may impair mitochondrial functionality by reducing the activity and solubility of the mitochondrial chaperone HSPA9, thereby disrupting proper protein folding and mitochondrial homeostasis [73]. As HSPA9 is essential for the import and maturation of mitochondrial proteins involved in oxidative phosphorylation, DNLZ deficiency may indirectly affect the efficiency of the electron transport chain (ETC), leading to decreased ATP production and impaired energy metabolism [74]. Moreover, reduced ETC activity may limit electron leakage and subsequently attenuate excessive ROS generation [75]. Under heat stress conditions, abnormal ROS accumulation can induce lipid peroxidation, protein oxidation, and DNA damage, ultimately triggering apoptosis or necrosis [76]. Therefore, by downregulating DNLZ, cells can reduce energy metabolism rate, decrease ATP consumption, and effectively inhibit ROS generation, thereby reducing oxidative stress levels and enhancing tolerance to the adverse effects of heat stress. From a biological perspective, downregulation of DNLZ expression may be a regulatory mechanism by which cells actively reduce energy metabolism and ROS load under heat stress conditions, helping to alleviate oxidative damage and maintain cell survival. Further analysis shows that the key transcription factor HIF1 can bind to the CG site in the DNLZ promoter region, and the methylation level in this region significantly increases after heat stress, especially at sites 17 and 21, showing more significant changes. Due to the abundance of CpG islands in HIF1 binding sites, methylation may directly hinder HIF1's recognition of HREs by altering DNA spatial conformation, or competitively inhibit HIF1 binding by recruiting methylation binding proteins (such as MeCP2), thereby weakening its transcriptional activation of DNLZ. As a key hypoxia response factor, HIF1 can promote angiogenesis and improve oxygen transport efficiency, enhancing the survival ability of cells in low oxygen environments [77]. Based on this, we speculate that the increase in methylation levels at CG sites may inhibit the binding of HIF1 to the DNLZ promoter region, leading to downregulation of DNLZ gene expression. This epigenetic regulatory mechanism may represent a protective adaptive strategy for cells to cope with excessive accumulation of ROS induced by heat stress, by inhibiting DNLZ expression to reduce the activity of mitochondrial complex I, thereby reducing electron leakage and ROS production, ultimately enhancing the survival ability of cells under heat stress conditions.

This study systematically revealed the important role of DNA methylation in the heat stress response of dairy cows, but there are still certain limitations. Firstly, limited by experimental conditions, the sample size is small, and the WGBS stage adopts a mixed sample strategy, which helps to improve sequencing depth and reduce technical errors, but may mask the true epigenetic variations between individuals and reduce the detection accuracy of differentially methylated regions. Secondly, this study only used blood as the sampling tissue. Although blood is representative as a peripheral tissue, its methylation characteristics may not fully reflect the tissue-specific regulatory mechanisms in key target tissues such as the liver, breast, or hypothalamus. In addition, although the luciferase reporter assay and BSP sequencing have verified the regulatory effect of DNA methylation on gene expression, the relevant conclusions are mainly based on in vitro models and it is still difficult to fully simulate the complex physiological state in dairy cows. Due to limitations in sample acquisition, ethical approval, and experimental techniques, in vivo functional validation remains a bottleneck in current research. In the future, animal models and primary cells can be combined to further explore the epigenetic regulatory mechanisms of key genes in heat stress, providing theoretical basis and technical support for molecular breeding of heat-resistant high-quality dairy cows.

Conclusions

This study reveals the regulatory mechanism of whole genome DNA methylation in the heat stress response of dairy cows. Through WGBS analysis, 49861 DMRs and 7613 DMGs were identified, significantly enriched in key pathways such as substance transport, oxidative stress metabolism, signal transduction, and energy homeostasis. Among them, the promoter region of the DNLZ gene showed a significant increase in methylation levels under heat stress, leading to suppressed gene expression. In vitro experiments have validated the inhibitory effect of DNLZ promoter methylation on transcriptional activity, indicating that its methylation affects gene expression and thus participates in the response of dairy cows to heat stress. This study provides new molecular targets for cultivating heat-resistant dairy cows and potential epigenetic markers for precision breeding and management strategies, providing new basis for the development of molecular markers for heat-resistant traits. These achievements not only provide theoretical support for analyzing the heat adaptation mechanism of dairy cows, but also provide potential epigenetic markers and molecular intervention targets for precision breeding and management strategies, which are expected to accelerate the cultivation of high-yield and heat-resistant high-quality dairy cow lines.