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

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.

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.

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.