The alteration of N6-methyladenosine (m6A) modification at the transcriptome-wide level in response of heat stress in bovine mammary epithelial cells

Background Heat stress has a substantial negative economic impact on the dairy industry. N6-methyladenosine (m6A) is the most common internal RNA modification in eukaryotes and plays a key role in regulating heat stress response in animals. In dairy cows, however, this modification remains largely unexplored. Therefore, we examined the effects of heat stress on the m6A modification and gene expression in bovine mammary epithelial cells to elucidate the mechanism of heat stress response. In this study, Mammary alveolar cells-large T antigen (MAC-T) cells were incubated at 37 °C (non-heat stress group, NH) and 40 °C (heat stress group, H) for 2 hours, respectively. HSP70, HSF1, BAX and CASP3 were up regulated in H group compared with those in the NH group. Results Methylated RNA immunoprecipitation sequencing (MeRIP-seq) and RNA sequencing (RNA-seq) were conducted to identify m6A peaks and to produce gene expression data of MAC-T cells in the H and NH groups. In total, we identified 17,927 m6A peaks within 9355 genes in the H group, and 18,974 peaks within 9660 genes in the NH groups using MeRIP-seq. Compared with the NH group, 3005 significantly differentially enriched m6A peaks were identified, among which 1131 were up-regulated and 1874 were down-regulated. In addition, 1502 significantly differentially expressed genes were identified using RNA-seq, among which 796 were up-regulated and 706 were down-regulated in the H group compared to the NH group. Furthermore, 199 differentially expressed and synchronously differentially methylated genes were identified by conjoint analysis of the MeRIP-seq and RNA-seq data, which were subsequently divided into four groups: 47 hyper-up, 53 hyper-down, 59 hypo-up and 40 hypo-down genes. In addition, GO enrichment and KEGG analyses were used to analyzed the potential functions of the genes in each section. Conclusion The comparisons of m6A modification patterns and conjoint analyses of m6A modification and gene expression profiles suggest that m6A modification plays a critical role in the heat stress response by regulating gene expression. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-022-09067-6.

Heat stress is a major limiting factor in dairy production because of its adverse effects on productivity, health, reproduction and overall welfare of dairy animals [22,23]. The economic loss due to heat stress is estimated at over US $1.2 billion per year due to the reduction in total milk yield and the reduced quality of milk [22,24]. The long-term increase in global temperatures aggravates the negative impacts of heat stress. And genetic antagonism between heat tolerance and milk production is an additional complicating factor [25]. Dairy cattle have higher metabolic turnover than beef cattle, and the continuous selection for higher milk production results in decreased heat tolerance [26]. Thus, there is growing interest in understanding the molecular mechanisms underlying heat stress response, and better balance of thermotolerance and high milk production traits. Milk is produced in the mammary glands, thus the effect of heat stress on gene expression is particularly important with regard to these glands. Recent large-scale transcriptome sequencing has shown that heat stress significantly changes the expression profiles of coding and noncoding RNAs in mammary glands, including mRNAs, microRNAs, long noncoding RNAs and circular RNAs [27][28][29][30][31]. The most extensively studied heat stress related genes are heat shock protein (HSP) genes, including HSP110 (HSPH), HSP90 (HSPC), HSP70 (HSPA), HSP40 (DNAJ), and the small HSP (HSPB) family [32]. Epigenetic changes, including DNA methylation and histone acetylation, play important roles in regulating heat stress response [33,34]. Recently, it has been confirmed that changes in m6A modifications are closely related to heat stress response as well. In sheep, the general characteristics of m6A modification and methylation profiles altered during heat stress [35]. Further, m6A modification has been reported to be involved in the response of turbot to heat stress [36]. However, the respective mechanisms have not been investigated in dairy animals. In addition, heat stress markedly increase methylation of the 5′ untranslated region (5'UTR) of HSP70 mRNA, which facilitates initiation of cap-independent translation of HSP70 [18]. However, the role of m6A RNA methylation in the regulation of heat stress in dairy cattle remains unclear.
This study was conducted to investigate the effects of heat stress on m6A modification and gene expression in mammary epithelial cells to elucidate the mechanism of the heat stress response. Mammary alveolar cells-large T antigen (MAC-T) cell are an established mammary epithelial cell line that retains the phenotypic characteristics of bovine mammary epithelial cells. We incubated MAC-T cells at 37 °C (non-heat stress group, NH) or 40 °C (heat stress group, H) for 2 hours. Subsequently, we investigated the whole transcriptome m6A profile and gene expression profile in the H and NH groups using methylated RNA immunoprecipitation sequencing (MeRIP-seq) and RNA sequencing (RNA-seq), respectively. Different m6A modification and gene expression patterns were systematically analyzed to improve our understanding of the role of m6A modification in regulating the heat stress response in dairy cattle.

Gene expression of mammary epithelial cells after heat stress
MAC-T cells were cultured at 37 °C and were then assigned to two treatments, i.e., the NH treatment with, continued culturing at cultured at 37 °C, and the H treatment, with culturing at 40 °C for 2 hours. Molecular indicators of heat stress were detected by quantitative real-time PCR (qRT-PCR). The relative expression of HSP70, HSF1, BAX and CASP3 was up-regulated 21.4, 6.1, 6.5, and 4.8 fold, respectively, in H group, compared to that of the NH group (Fig. 1).

General features of m6A methylation in heat stressed MAC-T cells
To investigate the role of m6A modification due to heat stress, MeRIP-seq was performed. More than 20 M reads were generated per sample (Table 1). Over 96% of the clean reads were mapped to the reference genome, with over 94% of the clean reads being unique ( Table 2). The quality control and mapping statistics of MeRIP-seq data from all samples are displayed in Tables 1 and 2, which   indicated that all samples were of high quality that was sufficient for subsequent analyses.
To investigate the distribution of m6A modifications in transcripts, we examined the metagene profiles of m6A peaks in the H and NH groups. As shown in Fig. 2a, m6A peaks were markedly enriched near the stop codon. Two further distinct peaks were observed in the start codon and 3′ UTR. The m6A density in the H group was lower than that in the NH group. To systematically assess the enrichment, the transcript was divided into five regions: transcription start site (TSS), 5′ UTR, coding sequence (CDS), 3′ UTR and the stop codon (Stop). Each m6A peak was assigned to one of the five segments, and the enrichment proportion was calculated. In both groups, more than 40% of m6A peaks were in 3′ UTR segment, more than 30% were in the CDS region, and less than 2% were at the TSS ( Fig. 2b and c). The proportion of m6A peaks enriched in the TSS, 5′ UTR, CDS, and Stop segments was higher in the H group than in the NH group. However, the trend in the 3′ UTR was the opposite. The m6A consensus motif was analyzed using MEME software [37]. The top five m6A consensus motifs were in the H and NH groups are shown in Fig. 2d and e. The consensus motif "RRACH" was among the top three in the H group and among the top five in the NH group. In addition, we characterized the pattern of m6A modification of mRNA. And a unique m6A peak occurred in more than 40% of the genes. Approximately 90% of genes contained one to three m6A peaks. Only approximately 10% of the genes contained more than three peaks per gene (Fig. 2f ).

Differential m6A modification between H and NH groups
In total, we identified 17,927 m6A peaks in the H group, which were associated with 9356 genes, and 18,974 peaks in the NH group, which were associated with 9660 genes (Supplementary Tables 1 and 2). Overall, 13,776 m6A peaks overlapped between the two groups, and 3273 and 4214 peaks were unique to the H and NH group, respectively (Fig. 3a). The top 20 significantly differentially expressed m6A peaks in the H and NH groups are shown in Table 3. Moreover, 8590 m6A-containing genes overlapped between the two groups, and 766 and 1070 unique genes were unique to the in H and NH groups (Fig. 3b). As shown in Fig. 3c, we analyzed all m6A peaks and 3005 were identified as significantly differentially enriched in the H group, compared with the NH group, of which 1131 peaks were upregulated and 1874 were down regulated (P < 0.05, fold change > 1.5). Overall, the results of MeRIP-seq data analyses showed an altered landscape of m6A methylation under heat stressed condition.
To further assess the role of m6A modification in MAC-T cells in response to heat stress, we analyzed the biological functions of differential m6A peak-related genes. All differential m6A peak-related genes were subjected to Gene Ontology (GO, http:// www. geneo ntolo gy. org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG, http:// www. kegg. jp/) pathway analyses. The top 30 significantly enriched GO items, ranked by P-value, are shown in Fig. 3d. In terms of biological processes, m6A-containing genes were mainly involved in metabolic processes, such as organic substance metabolic process, primary metabolic process, nitrogen compound metabolic process and cellular metabolic process. Cellular components were mainly related to intracellular part. In terms of molecular function, they were significantly correlated with heterocyclic compound binding, organic cyclic compound binding, nucleic acid binding, ion binding and zinc ion binding. The results of the KEGG analysis showed that the m6A-containing genes were significantly associated with ubiquitin mediated proteolysis, and the mammalian target of rapamycin (mTOR) signaling pathway and Hippo signaling pathway (Fig. 3e). Moreover, most genes were associated with the metabolic pathways contained.

Conjoint analysis of MeRIP-seq and RNA-seq data
To explore the potential role of m6A modification in gene expression, we examined the transcriptome profiles of H and NH groups using RNA-seq. More genes were expressed in the H group (Fig. 4a); Compared with the NH group, 796 genes were up-regulated and 706 were down-regulated in the H group (Fig. 4b). The top 20 differentially expressed genes are listed in Table 4. Subsequently, the RNA-seq and MeRIP-seq data were combined to elucidate the relationship between m6A modification and gene expression. 199 differentially expressed genes and synchronously differentially methylated genes were observed, which are referred to here as diff-diff genes All diff-diff genes were divided into four groups, including 47 hypermethylated and up-regulated genes (hyper-up), 53 hypermethylated and down-regulated genes (hyper-down), 59 hypomethylated and upregulated genes (hypo-up) and 40 hypomethylated and down-regulated genes (hypo-down) (Fig. 4c). All diff-diff genes were subjected to GO and KEGG analyses. The top 30 GO terms are shown in Fig. 4d. Diff-diff genes were mainly enriched with regard to biosynthetic and metabolic processes, such as organic substance biosynthetic process, cellular macromolecule biosynthetic process and cellular nitrogen compound metabolic process. The results of the KEGG analysis showed that diff-diff genes were related to steroid biosynthesis, fructose and mannose metabolism, and the HIF-1 signaling pathway (Fig. 4e).  showing the overlap of m6A containing genes between H and NH groups. c Volcano plots of the significantly different m6A peaks. d GO annotation of differential m6A peak-related genes. e KEGG enrichment analysis of differential m6A peaks related genes m6A modification of HSP genes HSP40, HSP70, HSP90, and HSP110 were significantly up-regulated ( Table 5). The most notable m6A modification peak was observed in the vicinity of the start codon in HSP genes (Fig. 5). However, only DNAJA1 was synchronously differentially methylated between NH and H conditions (Table 5).

Discussion
Improving our understanding of the molecular mechanisms underlying heat stress response is critical for addressing the problem of heat stress in the dairy industry. m6A modification and its related enzymes have emerged as regulators of the cellular stress response [38]. Extensive research on m6A modification and cellular stress has been conducted on humans and model organisms [18,21,39,40]. In Drosophila, m6A modification reduces the brain's biological response to heat stress in a METTL3-dependent manner [41]. Stress granules can assemble under stress conditions to regulate mRNA translation and degradation. m6A genes are enriched in stress granules, and YTHDF proteins play essential roles in their formation in human cells [40]. However, the role of m6A modification in the heat stress response of dairy cows has not yet been elucidated. To our knowledge, this is the first transcriptome-wide study to report an m6A map of bovine mammary epithelial cells under heat stress.
We produced an m6A modification map of MAC-T cells using MeRIP-seq. According to these data, there were 1.9 peaks per gene in the H group and 2.0 peaks per gene in the NH group. The m6A modification sites were more abundant than those in sheep. As sheep genes were previously reported to contain on average 0.32 and 0.42 m6A modification sites under heat stress and control conditions, respectively [35]. Compared with the NH group, we identified 3005 significantly differentially enriched m6A peaks that were related to 2620 genes in the H group. Of these 2620 genes, 73 were zinc finger protein genes, suggesting that m6A modification may be involved in the regulation of transcription by RNA polymerase II. We found no significant difference between the H and NH groups in terms of the distribution of m6A peaks (Fig. 2a), whereas previous study found that heat shock led to a marked enrichment of m6A in the 5'UTR [21]. This discrepancy may be attributed to the use of different cell lines. In addition, we found more expressed genes but fewer m6A peaks in the H group, implying that m6A methylation may be a critical impetus for gene expression (Figs. 3 and 4a).
HSP genes are induced by cellular stress to protect cells against apoptosis and provide thermotolerance [42][43][44]. HSPA1A and HSPA6 are the most heat-shockinducible HSP members [45]. Consistent with previous studies, the top three up-regulated genes in the H group were HSPA6, HSPA1A and HSPH1. The top  [46]. In a previous study on mouse embryonic fibroblast cells, the HSP gene HSPA1A showed elevated m6A modifications in the 5'UTR in response to heat stress [18]. In the present study, members of HSP40, HSP70, HSP90 and HSP110 were found to be m6A modified according to MeRIP-seq data. We identified significantly different m6A peaks in HSP genes including HSPD1, HSPA2, and HSPA12A; however, these genes were not significantly differentially expressed. Meanwhile, we identified 796 up-regulated genes by RNA-seq analysis, including HSP genes such as HSPA6, HSPA1A, and HSPH1 (Table 5). HSP genes, especially those of the HSP70 and HSP40 families, were significantly up-regulated at the mRNA level, which is in accordance with previous results at the protein level [47]. Intriguingly, however, no significantly different m6A peak was found in these HSP genes except for DNAJA1. DNAJA1 was hypo-up gene in our conjoint analysis of MeRIP-seq and RNA-seq data. DNAJA1, HSPA1A, DNAJB1, and HSPH1 are genes encoding components of the mammalian protein disaggregation/refolding machine, which dissociates and refolds aggregated proteins [48]. All these genes were significantly up-regulated under heat shock condition. This indicates that cells enhance their capacity to combat misfolded proteins caused by high temperatures by increasing protein disaggregation/refolding.m6A modification plays a key role in regulating mRNA splicing, stability, localization and translation efficiency [18][19][20][21].
In line with previous findings, we found that the differential m6A peaks were related to biological progress items such as RNA metabolic process in GO analysis and RNA transport pathway in KEGG analysis. Mammary gland growth, development and function are critical events that affect the lactation and reproductive performance of dairy cows. m6A genes are enriched for key signaling pathways such as Metabolic pathway, Hippo, Wnt, transforming growth factor β (TGFβ) and mTOR signaling pathways, and our KEGG analysis results were consistent with the results of previous studies [41]. Metabolic pathway was the most enriched pathway with the regard to m6A containing genes. Metabolic activity, especially carbohydrate and lipid metabolism, is decreased in mammary gland [49]. Our results suggest that m6A modification is involved in metabolic regulation under heat stress condition. The Wnt signaling pathway is associated with multiple stages of mammary   [50,51]. This implies that heat stress may affect mammary gland development through m6A regulated genes. The mTOR signaling pathway is relevant for lactation, and it can be activated by amino acids to promote milk synthesis in mammary epithelial cells [52]. A previous study found that heat shock can alter the expression levels of genes associated with the mTOR signaling pathway, thus the amino acid metabolism and concentration in bovine mammary epithelial cells were proposed to regulate the mTOR signaling pathway to adapt to heat stress [53]. However, the specific underlying mechanisms remain to be elucidated. In the current study, the key molecules of the mTOR signaling pathway, such as EIF4EBP1, EIF4E and TSC1, were enriched according to the KEGG analysis of the differentially m6A peak-related genes. EIF4EBP1 was hypodown gene in the conjoint analysis of MeRIP-seq and RNA-seq data. This suggests that mammary epithelial cells may affect the mTOR signaling pathway through m6A modification-mediated gene regulation.
Overall, we presented transcriptome-wide profiling of the m6A modification under heat stress and normal conditions and postulated a molecular mechanism that is shown in Fig. 6. HSPs, including HSP40, HSP70, HSP90, and HSP110, were up-regulated immediately after heat shock. Higher HSP expression has been associated with increased mechanistic target of rapamycin complex 1 (mTOCR1) activity [54]. mTORC1 enhances the translation of WTAP and thus activates the m6A writer complex [55]. Thereafter, the global m6A modification is increased. This can partially explain why our RNA-seq data showed that the mRNA expression levels of key m6A regulators did not change under heat stress. The m6A modification of genes associated with TGFβ, Wnt, and Hippo signaling pathways was influenced, thereby modifying cell signaling. Finally, cell proliferation was inhibited and mammary gland development was blocked during heat stress. Meanwhile, the expression of milk protein-encoding genes was reduced. These results reveal the potential regulatory mechanism of the m6A modification on milk production and milk quality under heat stress.

Conclusion
In summary, this study provides a resource for defining the role of m6A modification in the heat stress response of dairy cattle and reveals the potential genetic regulatory mechanism of m6A modification.

Cell culture
MAC-T cells, which is purchased from the ATCC Corporation (Beijing, China), were cultured in Dulbecco's modified Eagle' s medium (DMEM)/nutrient mixture F12 (Ham) medium (Gibco, Grand Island, NY, USA) supplemented with 10% fetal bovine serum (Gibco), 200 U/ mL penicillin and streptomycin. Cells were cultured in T-75 cell flasks in an incubator at 37 °C and under 5% CO 2 . After 72 hours, cells were assigned to two groups, with four replicates, each. The NH group cells were continuously cultured at 37 °C and the H group at 40 °C for 2 hours.

RNA extraction
Total RNA was isolated from MAC-T cells of the H and NH groups using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), in accordance with the manufacturer's instructions. RNA integrity and concentration were determined using an Agilent 2100 bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and a SimpliNano spectrophotometer (GE Healthcare, Freiburg, Germany), respectively.

qRT-PCR
RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, Waltham, MA, USA) was used to reverse-transcribe the RNA samples. qRT-PCR was performed using GoTaq qPCR Master Mix (Promega, Madison, WI, USA) according to the manufacturer's instructions. Thermocycling included 94 °C for 5 min and 35 cycles of 95 °C for 30 s, 58 °C for 15 s, and 72 °C for 30 s, and was carried out on a Bio-Rad IQ5 system (Bio-Rad, Hercules, CA, USA). β-actin was used as an internal control, and the relative gene expression level was calculated using the 2 -ΔΔCT method of. Primers are listed in Table 6.

MeRIP-seq and data analysis
Fragmented mRNA was incubated with anti-m6A polyclonal antibody (Synaptic Systems, Coventry, UK) in IP buffer for 2 h at 4 °C. Then, immunoprecipitated mRNA or Input were used for library construction using NEB-Next Ultra II Directional RNA Library Prep Kit for Illumina (New England Biolabs, Ipswitch, MA, USA). After library quality control, sequencing was performed on an Illumina NovaSeq platform (Illumina, San Diego, CA, USA) according to standard protocols. Raw data in fastq format were processed using fastp (version 0.19.11) [56]. High-quality clean data were obtained after Q30 quality control and removal of reads containing adapter. Then, clean reads were aligned to the cattle reference genome using BWA software (v0.7.12) [57]. The R package exomePeak (v2.16.0) was used for m6A peak identification of the IP and input samples [58]. Differential peak calling was performed using R package exomePeak (v2. 16.0) at P < 0.05 and fold change≥1.5. The m6A-enriched motifs in each group were identified by HOMER (v4.9.1). Enrichment analyses were carried out using GO and KEGG for the peak related genes [59]. GO enrichment analysis was carried out by the GOseq R package, in which gene length bias was corrected. GO terms with corrected P value< 0.05 were considered significantly enriched. KOBAS software (v 3.0) was used to test the statistical enrichment of genes in KEGG pathways.

mRNA-seq
mRNA was purified from the total RNA using magnetic beads with an attached poly-T oligo. cDNA libraries were prepared using the Ultra RNA Library Prep Kit for Illumina (New England Biolabs), according to the manufacturer's instructions. After cluster generation, the cDNA libraries were sequenced on an Illumina NovaSeq platform (Illumina), and 150 bp paired-end reads were generated. Raw data in fastq format were first processed using in-house Perl scripts, and clean data were obtained after removing reads containing adapter, poly-N, and low-quality reads. The read counts for each sample were adjusted and differential expression analysis of the H and NH groups was performed using edgeR package (v3.22.5). The P values were adjusted using the Benjamini and Hochberg method. Corrected P-value < 0.05 and fold change ≥1.5 was considered to indicate significantly differential expression.