Skip to main content
Download PDF

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

BMC Genomics volume 23, Article number: 829 (2022) Cite this article

Abstract

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.

Peer Review reports

Introduction

N6-methyladenosine (m6A) is the most frequent internal RNA modification in eukaryotes, which was first discovered in the 1970s [1]. To date, over 160 types of chemical modifications have been identified in RNA molecules [2, 3], with methylation accounting for more than 60% [4]. m6A has been identified in various taxa, including yeast [5], plants [6], insects, mammals [7, 8], and some viruses [9]. m6A modification is present transcriptome-wide in more than 25% of all RNAs and preferably occurs in highly conserved regions with the consensus motif “RRACH” (R = G or A; H = A, C, or U) [10]. The methylation of “A” is highly dynamic and reversible [11]. It is catalyzed by “writers” (m6A methylases) including methyltransferase-like (METTL)3, METTL14, and Wilms’ tumor 1 associated protein (WTAP) [12, 13]. This modification can be reversed by “erasers” (m6A demethylases), including fat mass and obesity associated protein (FTO) and AlkB homolog H5 (ALKBH5) [14]. The m6A modification site is recognized by m6A “readers” (m6A binding proteins) such as YTH m6A-binding protein (YTHDF)1, YTHDF2, and YTHDF3; the heterogeneous nuclear ribonucleoprotein (HNRNP) A2/B1 and HNRNPC; and the YTH domain containing (YTHDC)1, YTHDC2, and IGF2BPs [15,16,17]. RNA modification through m6A may affect mRNA translation efficiency, stability, splicing, and nuclear export [4, 18,19,20,21]. Accumulating evidence indicates that m6A modification plays a vital role in cellular pathways and processes, such as cell differentiation, development, and metabolism.

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.

Results

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).

Fig. 1
figure 1

Effect of heat stress on HSP70, HSF1, BAX, and CASP3 expression in MAC-T cells. Expression of HSP70 (a), HSF1 (b), BAX (c), and CASP3 (d) was detected through qRT-PCR

Full size image

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.

Table 1 Statistics and quality control of data generated by MeRIP sequencing
Full size table
Table 2 Statistics of mapping data
Full size table

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).

Fig. 2
figure 2

m6A topological patterns in MAC-T cells. a Metagene plots showing the m6A peak density distribution across the transcripts in the NH and H groups. Gene segments distribution of m6A peaks in the H (b) and NH (c) groups. The top sequence motifs enriched within m6A peaks in the H (d) and NH (e) groups. (f) Proportion of genes containing different numbers of m6A peaks in the NH and H groups

Full size image

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.

Fig. 3
figure 3

Differential m6A peaks between H and NH groups. a Venn diagram showing the overlap of m6A peaks between H and NH groups. b Venn diagram 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

Full size image
Table 3 Top 20 significantly differentially expressed m6A peaks in H and NH groups
Full size table

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.geneontology.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 up-regulated 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).

Fig. 4
figure 4

Conjoint analysis of MeRIP-seq and RNA-seq data. a Venn diagram showing the numbers of genes expressed in the H and NH groups. b Volcano plots showing differentially expressed genes. c Four-quadrant diagram showing the distribution of genes with significant changes in both m6A modification and mRNA expression. d GO annotation of genes with significant changes in m6A modification and mRNA expression. e KEGG enrichment analysis of genes with significant changes in m6A modification and mRNA expression

Full size image
Table 4 Top 20 significantly differentially expressed genes in H and NH groups
Full size table

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).

Table 5 Expression and m6A modification of HSP genes
Full size table
Fig. 5
figure 5

The distribution of m6A peaks in HSP mRNA. Visualization of HSPA1A, HSPA5, HSPA8, and HSPH1 using Integrative Genomics Viewer (IGV) software. The blue box indicates an exon and the blue line indicates an intron

Full size image

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-shock-inducible HSP members [45]. Consistent with previous studies, the top three up-regulated genes in the H group were HSPA6, HSPA1A and HSPH1. The top five included DNAJB9, and the top eight was DNAJA4, which belong to HSP40 (DNAJ) family. These results indicate that the heat shock treatment of MAC-T cells was effective. Moreover, HSP genes, such as HSPA1A, HSPA1B, HSPB1, HSPA9, HSP90AA1, and HSPD1, were m6A methylated [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 development. In mammary gland, there are 10 Wnt proteins including Wnt 2,4, 5a, 5b, 6, 7b, 9a, 9b, 10b, and 16. Genes including Wnt10b, Wnt5a, LRP5, and LRP6 showed different m6A peaks under H and NH conditions and were enriched in the Wnt signaling pathway, according to KEGG analysis. LRP5 and LRP6 are co-receptors for Wnt proteins. Both of which are required to respond to Wnt10b, and Wnt10b overexpression leads to excessive branching and precocious alveolar development [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 hypo-down 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.

Fig. 6
figure 6

Effect of heat stress on cells through the m6A modification. The speculated molecular mechanism by which heat stress affects m6A is shown. The cell signaling pathways enriched in the KEGG analysis are indicated

Full size image

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.

Materials and methods

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% CO2. 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.

Table 6 Primer sequences of mRNA in qRT-PCR
Full size table

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 NEBNext 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.

Availability of data and materials

The datasets generated and analyzed during the current study are available in the NCBI Sequence Read Archive under the accession number PRJNA863056, http://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/bioproject/863056.

Abbreviations

m6A:

N6-methyladenosine

NH:

Non-heat stress group

H:

Heat stress group

TSS:

Transcription start site

5′ UTR:

5′ untranslated regions

CDS:

Coding sequence

3′ UTR:

3′ untranslated regions

Stop:

Stop codon

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

mTOR:

Mammalian target of rapamycin

mTORC1:

Mammalian target of rapamycin complex 1

TGFβ:

Transforming growth factor β

References

  1. Wei CM, Gershowitz A, Moss B. Methylated nucleotides block 5′ terminus of HeLa cell messenger RNA. Cell. 1975;4(4):379–86.

    Article  CAS  Google Scholar 

  2. Boccaletto P, Machnicka MA, Purta E, Piatkowski P, Baginski B, Wirecki TK, et al. MODOMICS: a database of RNA modification pathways. 2017 update. Nucleic Acids Res. 2018;46(D1):D303–7.

    Article  CAS  Google Scholar 

  3. Frye M, Harada BT, Behm M, He C. RNA modifications modulate gene expression during development. Science. 2018;361(6409):1346–9.

    Article  CAS  Google Scholar 

  4. Cantara WA, Crain PF, Rozenski J, McCloskey JA, Harris KA, Zhang X, et al. The RNA modification database, RNAMDB: 2011 update. Nucleic Acids Res. 2011;39(Database issue):D195–201.

    Article  CAS  Google Scholar 

  5. Clancy MJ, Shambaugh ME, Timpte CS, Bokar JA. Induction of sporulation in Saccharomyces cerevisiae leads to the formation of N6-methyladenosine in mRNA: a potential mechanism for the activity of the IME4 gene. Nucleic Acids Res. 2002;30(20):4509–18.

    Article  CAS  Google Scholar 

  6. Wang Z, Tang K, Zhang D, Wan Y, Wen Y, Lu Q, et al. High-throughput m6A-seq reveals RNA m6A methylation patterns in the chloroplast and mitochondria transcriptomes of Arabidopsis thaliana. PLoS One. 2017;12(11):e0185612.

    Article  Google Scholar 

  7. Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR. Comprehensive analysis of mRNA methylation reveals enrichment in 3′ UTRs and near stop codons. Cell. 2012;149(7):1635–46.

    Article  CAS  Google Scholar 

  8. Wei CM, Moss B. Nucleotide sequences at the N6-methyladenosine sites of HeLa cell messenger ribonucleic acid. Biochemistry. 1977;16(8):1672–6.

    Article  CAS  Google Scholar 

  9. Manners O, Baquero-Perez B, Whitehouse A. m(6)A: Widespread regulatory control in virus replication. Biochim Biophys Acta Gene Regul Mech. 2019;1862(3):370–81.

    Article  CAS  Google Scholar 

  10. Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar L, Osenberg S, et al. Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature. 2012;485(7397):201–6.

    Article  CAS  Google Scholar 

  11. Jia G, Fu Y, Zhao X, Dai Q, Zheng G, Yang Y, et al. N6-methyladenosine in nuclear RNA is a major substrate of the obesity-associated FTO. Nat Chem Biol. 2011;7(12):885–7.

    Article  CAS  Google Scholar 

  12. Wang Y, Li Y, Toth JI, Petroski MD, Zhang Z, Zhao JC. N6-methyladenosine modification destabilizes developmental regulators in embryonic stem cells. Nat Cell Biol. 2014;16(2):191–8.

    Article  CAS  Google Scholar 

  13. Ping XL, Sun BF, Wang L, Xiao W, Yang X, Wang WJ, et al. Mammalian WTAP is a regulatory subunit of the RNA N6-methyladenosine methyltransferase. Cell Res. 2014;24(2):177–89.

    Article  CAS  Google Scholar 

  14. Mauer J, Luo X, Blanjoie A, Jiao X, Grozhik AV, Patil DP, et al. Reversible methylation of m(6)Am in the 5′ cap controls mRNA stability. Nature. 2017;541(7637):371–5.

    Article  CAS  Google Scholar 

  15. Xu C, Wang X, Liu K, Roundtree IA, Tempel W, Li Y, et al. Structural basis for selective binding of m6A RNA by the YTHDC1 YTH domain. Nat Chem Biol. 2014;10(11):927–9.

    Article  CAS  Google Scholar 

  16. Wang X, Zhao BS, Roundtree IA, Lu Z, Han D, Ma H, et al. N(6)-methyladenosine modulates messenger RNA translation efficiency. Cell. 2015;161(6):1388–99.

    Article  CAS  Google Scholar 

  17. Zhao BS, Roundtree IA, He C. Publisher correction: post-transcriptional gene regulation by mRNA modifications. Nat Rev Mol Cell Biol. 2018;19(12):808.

    Article  CAS  Google Scholar 

  18. Zhou J, Wan J, Gao X, Zhang X, Jaffrey SR, Qian SB. Dynamic m(6)A mRNA methylation directs translational control of heat shock response. Nature. 2015;526(7574):591–4.

    Article  CAS  Google Scholar 

  19. Edupuganti RR, Geiger S, Lindeboom RGH, Shi H, Hsu PJ, Lu Z, et al. N(6)-methyladenosine (m(6)A) recruits and repels proteins to regulate mRNA homeostasis. Nat Struct Mol Biol. 2017;24(10):870–8.

    Article  CAS  Google Scholar 

  20. Wang X, Lu Z, Gomez A, Hon GC, Yue Y, Han D, et al. N6-methyladenosine-dependent regulation of messenger RNA stability. Nature. 2014;505(7481):117–20.

    Article  Google Scholar 

  21. Meyer KD, Patil DP, Zhou J, Zinoviev A, Skabkin MA, Elemento O, et al. 5′ UTR m(6)A promotes cap-independent translation. Cell. 2015;163(4):999–1010.

    Article  CAS  Google Scholar 

  22. Tao S, Orellana Rivas RM, Marins TN, Chen YC, Gao J, Bernard JK. Impact of heat stress on lactational performance of dairy cows. Theriogenology. 2020;150:437–44.

    Article  Google Scholar 

  23. Yan G, Liu K, Hao Z, Shi Z, Li H. The effects of cow-related factors on rectal temperature, respiration rate, and temperature-humidity index thresholds for lactating cows exposed to heat stress. J Therm Biol. 2021;100:103041.

    Article  Google Scholar 

  24. Wankar AK, Rindhe SN, Doijad NS. Heat stress in dairy animals and current milk production trends, economics, and future perspectives: the global scenario. Trop Anim Health Prod. 2021;53(1):70.

    Article  Google Scholar 

  25. Carabano MJ, Bachagha K, Ramon M, Diaz C. Modeling heat stress effect on Holstein cows under hot and dry conditions: selection tools. J Dairy Sci. 2014;97(12):7889–904.

    Article  CAS  Google Scholar 

  26. Ravagnolo O, Misztal I. Genetic component of heat stress in dairy cattle, parameter estimation. J Dairy Sci. 2000;83(9):2126–30.

    Article  CAS  Google Scholar 

  27. Li L, Sun Y, Wu J, Li X, Luo M, Wang G. The global effect of heat on gene expression in cultured bovine mammary epithelial cells. Cell Stress Chaperones. 2015;20(2):381–9.

    Article  CAS  Google Scholar 

  28. Qi Y, Zhang L, Guo Y, Wang J, Chu M, Zhang Y, Guo J, Li Q. Genome-wide identification and functional prediction of circular RNAs in response to heat stress in Chinese Holstein cows. Anim Biotechnol. 2022;33(6):1170-1180.

  29. Li Q, Yang C, Du J, Zhang B, He Y, Hu Q, et al. Characterization of miRNA profiles in the mammary tissue of dairy cattle in response to heat stress. BMC Genomics. 2018;19(1):975.

    Article  CAS  Google Scholar 

  30. Li Q, Qiao J, Zhang Z, Shang X, Chu Z, Fu Y, et al. Identification and analysis of differentially expressed long non-coding RNAs of Chinese Holstein cattle responses to heat stress. Anim Biotechnol. 2020;31(1):9–16.

    Article  Google Scholar 

  31. Wetzel-Gastal D, Feitor F, van Harten S, Sebastiana M, Sousa LMR, Cardoso LA. A genomic study on mammary gland acclimatization to tropical environment in the Holstein cattle. Trop Anim Health Prod. 2018;50(1):187–95.

    Article  CAS  Google Scholar 

  32. Kampinga HH, Hageman J, Vos MJ, Kubota H, Tanguay RM, Bruford EA, et al. Guidelines for the nomenclature of the human heat shock proteins. Cell Stress Chaperones. 2009;14(1):105–11.

    Article  CAS  Google Scholar 

  33. Marinova Z, Leng Y, Leeds P, Chuang DM. Histone deacetylase inhibition alters histone methylation associated with heat shock protein 70 promoter modifications in astrocytes and neurons. Neuropharmacology. 2011;60(7–8):1109–15.

    Article  CAS  Google Scholar 

  34. Miozzo F, Saberan-Djoneidi D, Mezger V. HSFs, stress sensors and sculptors of transcription compartments and epigenetic landscapes. J Mol Biol. 2015;427(24):3793–816.

    Article  CAS  Google Scholar 

  35. Lu Z, Liu J, Yuan C, Jin M, Quan K, Chu M, et al. m(6)A mRNA methylation analysis provides novel insights into heat stress responses in the liver tissue of sheep. Genomics. 2021;113(1 Pt 2):484–92.

    Article  CAS  Google Scholar 

  36. Huang Z, Guo X, Wang Q, Ma A, Zhao T, Qiao X, et al. Digital RNA-seq analysis of the cardiac transcriptome response to thermal stress in turbot Scophthalmus maximus. J Therm Biol. 2022;104:103141.

    Article  CAS  Google Scholar 

  37. Bailey TL, Boden M, Buske FA, Frith M, Grant CE, Clementi L, et al. MEME SUITE: tools for motif discovery and searching. Nucleic Acids Res. 2009;37(Web Server issue):W202–8.

    Article  CAS  Google Scholar 

  38. Wilkinson E, Cui YH, He YY: Context-dependent roles of RNA modifications in stress responses and diseases. Int J Mol Sci. 2021;22(4):1949. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms22041949.

  39. Yu F, Wei J, Cui X, Yu C, Ni W, Bungert J, et al. Post-translational modification of RNA m6A demethylase ALKBH5 regulates ROS-induced DNA damage response. Nucleic Acids Res. 2021;49(10):5779–97.

    Article  CAS  Google Scholar 

  40. Fu Y, Zhuang X. m(6)A-binding YTHDF proteins promote stress granule formation. Nat Chem Biol. 2020;16(9):955–63.

    Article  CAS  Google Scholar 

  41. Perlegos AE, Shields EJ, Shen H, Liu KF, Bonini NM. Mettl3-dependent m(6)A modification attenuates the brain stress response in drosophila. Nat Commun. 2022;13(1):5387.

    Article  CAS  Google Scholar 

  42. Beere HM, Green DR. Stress management - heat shock protein-70 and the regulation of apoptosis. Trends Cell Biol. 2001;11(1):6–10.

    Article  CAS  Google Scholar 

  43. Zhang J, Miano FN, Jiang T, Peng Y, Zhang W, Xiao H. Characterization of three heat shock protein genes in Pieris melete and their expression patterns in response to temperature stress and Pupal diapause. Insects. 2022;13(5): 430. https://0-doi-org.brum.beds.ac.uk/10.3390/insects13050430.

  44. Yadav P, Yadav B, Swain DK, Anand M, Yadav S, Madan AK. Differential expression of miRNAs and related mRNAs during heat stress in buffalo heifers. J Therm Biol. 2021;97:102904.

    Article  CAS  Google Scholar 

  45. Brocchieri L, Conway de Macario E, Macario AJ. hsp70 genes in the human genome: conservation and differentiation patterns predict a wide array of overlapping and specialized functions. BMC Evol Biol. 2008;8:19.

    Article  Google Scholar 

  46. Yu J, Li Y, Wang T, Zhong X. Modification of N6-methyladenosine RNA methylation on heat shock protein expression. PLoS One. 2018;13(6):e0198604.

    Article  Google Scholar 

  47. Miao W, Yang YY, Wang Y. Quantitative proteomic analysis revealed broad roles of N(6)-Methyladenosine in heat shock response. J Proteome Res. 2021;20(7):3611–20.

    Article  CAS  Google Scholar 

  48. Deane CA, Brown IR. Components of a mammalian protein disaggregation/refolding machine are targeted to nuclear speckles following thermal stress in differentiated human neuronal cells. Cell Stress Chaperones. 2017;22(2):191–200.

    Article  CAS  Google Scholar 

  49. Gao ST, Ma L, Zhou Z, Zhou ZK, Baumgard LH, Jiang D, et al. Heat stress negatively affects the transcriptome related to overall metabolism and milk protein synthesis in mammary tissue of midlactating dairy cows. Physiol Genomics. 2019;51(8):400–9.

    Article  Google Scholar 

  50. Goel S, Chin EN, Fakhraldeen SA, Berry SM, Beebe DJ, Alexander CM. Both LRP5 and LRP6 receptors are required to respond to physiological Wnt ligands in mammary epithelial cells and fibroblasts. J Biol Chem. 2012;287(20):16454–66.

    Article  CAS  Google Scholar 

  51. Chu EY, Hens J, Andl T, Kairo A, Yamaguchi TP, Brisken C, et al. Canonical WNT signaling promotes mammary placode development and is essential for initiation of mammary gland morphogenesis. Development. 2004;131(19):4819–29.

    Article  CAS  Google Scholar 

  52. Qi H, Wang L, Zhang M, Wang Z, Gao X, Li M. Methionine and leucine induce ARID1A degradation to promote mTOR expression and milk synthesis in mammary epithelial cells. J Nutr Biochem. 2022;101:108924.

    Article  CAS  Google Scholar 

  53. Fu L, Zhang L, Liu L, Yang H, Zhou P, Song F, et al. Effect of heat stress on bovine mammary cellular metabolites and gene transcription related to amino acid metabolism, amino acid transportation and mammalian target of Rapamycin (mTOR) signaling. Animals (Basel). 2021;11(11): 3153. https://0-doi-org.brum.beds.ac.uk/10.3390/ani11113153.

  54. Chen F, Bao H, Xie H, Tian G, Jiang T. Heat shock protein expression and autophagy after incomplete thermal ablation and their correlation. Int J Hyperth. 2019;36(1):95–103.

    Article  CAS  Google Scholar 

  55. Jang KH, Heras CR, Lee G. m(6)A in the signal transduction network. Mol Cell. 2022;45(7):435–43.

    Article  CAS  Google Scholar 

  56. Chen S, Zhou Y, Chen Y, Gu J. Fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–90.

    Article  Google Scholar 

  57. Li H, Durbin R. Fast and accurate short read alignment with burrows-wheeler transform. Bioinformatics. 2009;25(14):1754–60.

    Article  CAS  Google Scholar 

  58. Meng J, Lu Z, Liu H, Zhang L, Zhang S, Chen Y, et al. A protocol for RNA methylation differential analysis with MeRIP-Seq data and exomePeak R/bioconductor package. Methods. 2014;69(3):274–81.

    Article  CAS  Google Scholar 

  59. Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2022. https://0-doi-org.brum.beds.ac.uk/10.1093/nar/gkac963.

Download references

Acknowledgements

Not applicable.

Funding

This study was supported by the Science and Technology Project of Hebei Education Department (ZD2021091), Doctoral Foundation of Langfang Normal University (XBQ202026), and Youth Foundation of Education Committee of Hebei Province (QN2020166).

Author information

Authors and Affiliations

  1. Hebei Key Laboratory of Animal Diversity, Langfang Key Laboratory of Cell Engineering and Application, College of Life Sciences, Langfang Normal University, Langfang, People’s Republic of China

    Ying Qi, Yiming Zhang, Jing Zhang, Jing Wang & Qiuling Li

Authors
  1. Ying Qi
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Yiming Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Jing Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Jing Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Qiuling Li
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

QLL and YQ conceived and designed the study. JZ and JW performed experiments. YMZ and YQ performed bioinformatics analyses. YQ wrote and QLL revised the manuscript. The author(s) read and approved the final manuscript.

Corresponding author

Correspondence to Qiuling Li.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

All the authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

m6A peaks in H and NH groups.

Additional file 2.

m6A peaks related genes in H and NH groups.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Qi, Y., Zhang, Y., Zhang, J. et al. The alteration of N6-methyladenosine (m6A) modification at the transcriptome-wide level in response of heat stress in bovine mammary epithelial cells. BMC Genomics 23, 829 (2022). https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-022-09067-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-022-09067-6

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us