Xinran Yang
Jianfang Wang
Xinhao Ma
Jiawei Du
Chugang Mei
Linsen Zan- 1College of Animal Science and Technology, Northwest A&F University, Yangling, China
- 2National Beef Cattle Improvement Center, Northwest A&F University, Yangling, China
N6-methyladenosine (m6A) is the most prevalent methylation modification of eukaryotic mRNA, and it plays an important role in regulating gene expression. Previous studies have found that m6A methylation plays a role in mammalian skeletal muscle development. However, the effect of m6A on bovine skeletal myogenesis are still unclear. Here, we selected proliferating myoblasts (GM) and differentiated myotubes (on the 4th day of differentiation, DM) for m6A-seq and RNA-seq to explore the m6A methylation modification pattern during bovine skeletal myogenesis. m6A-seq analysis revealed that m6A methylation was an abundant modification of the mRNA in bovine myoblasts and myotubes. We scanned 5,691–8,094 m6A-modified transcripts, including 1,437 differentially methylated genes (DMGs). GO and KEGG analyses revealed that DMGs were primarily involved in transcriptional regulation and RNA metabolism, as well as insulin resistance and metabolic pathways related to muscle development. The combined analysis further identified 268 genes that had significant changes at both m6A and mRNA levels, suggesting that m6A modification may regulate myoblast differentiation by mediating the expression of these genes. Furthermore, we experimentally confirmed four genes related to myogenesis, including MYOZ2, TWIST1, KLF5 and MYOD1, with differential changes in both m6A and mRNA levels during bovine myoblast differentiation, indicating that they can be potential candidate targets for m6A regulation of skeletal myogenesis. Our results may provide new insight into molecular genetics and breeding of beef cattle, and provide a reference for investigating the mechanism of m6A regulating skeletal muscle development.
Introduction
Skeletal muscle is an important factor that regulates livestock muscle quality and maintains metabolic homeostasis (Picard et al., 2010). The growth and development of skeletal muscle are extremely complex biological processes, which successively include directional differentiation of progenitor cells, myoblast proliferation, differentiation and fusion of myocytes, and, finally, formation of multinucleated muscle fibers with contractile function (Braun and Gautel, 2011). Besides a series of specific transcription factors, epigenetic modifications such as DNA methylation and histone methylation also play an important role in skeletal myogenesis (McKinnell et al., 2008; Yang et al., 2021). Nonetheless, molecular selection breeding in beef cattle mostly focused on the exploration of some key genes, and rarely improved the breeding process from the perspective of RNA.
More than 150 kinds of chemical RNA modifications have been identified, and RNA methylation accounts for more than 60% of all modifications (Cantara et al., 2011). Among these types of modification, N6-methyladenosine (m6A) is considered the most prevalent internal mRNA modification in eukaryotes (Wei et al., 1975, 1976; Schibler et al., 1977; Zhong et al., 2008; Jia et al., 2013; Schwartz et al., 2013; Luo et al., 2014; Meyer and Jaffrey, 2014; Haussmann et al., 2016; Shen et al., 2016). m6A is a dynamic and reversible posttranscriptional methylation modification (Fu et al., 2014; Meyer and Jaffrey, 2014; Zhao et al., 2017), that is catalyzed by m6A writer proteins (methyltransferase complexes composed of methyltransferase-like 3 (METTL3), METTL14 and Wilms tumour 1-associated protein (WTAP)) and is demethylated by m6A eraser proteins [fat mass and obesity-associated protein (FTO) and AlkB homolog H5 (ALKBH5)] (Jia et al., 2011; Zheng et al., 2013; Liu et al., 2014; Ping et al., 2014; Shi et al., 2019). The m6A modification is functionally interpreted by m6A “reader” proteins, such as the widely studied YTH-domain family proteins (Zhang et al., 2010; Wang et al., 2014). m6A modification plays a variety of roles in mRNA metabolism, including mRNA translation efficiency, stability, splicing, and nuclear export (Dominissini et al., 2012; Fu et al., 2014; Wang et al., 2014; Zhao et al., 2014; Alarcon et al., 2015; Meyer et al., 2015; Wang et al., 2015; Shi et al., 2019). Accumulating evidence suggests that m6A affects different developmental and biological processes, such as multiple cancer processes, mESC differentiation, antitumor immunity, embryonic and postembryonic development, cell rhythms, cell fate determination, and adipogenesis (Fustin et al., 2013; Batista et al., 2014; Zhao et al., 2014; Lin et al., 2016; Bertero et al., 2018; Yang et al., 2018; Han et al., 2019; Song et al., 2019; Yao et al., 2019; Wang et al., 2020).
Recent advances in the field have shown that m6A modification plays an important role in processes related to muscle growth and development, such as myocyte stem maintenance, myocyte proliferation, cell differentiation and myocardial function (Kudou et al., 2017; Wang et al., 2017; Dorn et al., 2019; Mathiyalagan et al., 2019; Gheller et al., 2020; Lin et al., 2020; Zhang et al., 2020). In primary mouse myoblasts and C2C12 cells, m6A is essential for skeletal muscle differentiation (Kudou et al., 2017; Li et al., 2021) and regulates the transitions of muscle stem cells/myoblasts (Gheller et al., 2020). In farm animals, using m6A-specific methylated RNA immunoprecipitation coupled with next-generation sequencing (m6A-seq), Tao revealed the m6A modification maps in porcine adipose and muscle tissue (Tao et al., 2017), Cheng found the m6A modification pattern in chicken fat deposition (Cheng et al., 2021), and Xu uncovered the m6A regulation in goose embryonic muscle development (Xu et al., 2021). These studies indicated that m6A modification plays a potentially important role in adipogenesis and myogenesis of animals. However, the studies on m6A in livestock and poultry mainly focus on the tissue development level, while the transcriptome-wide m6A methylome maps in cellular level have been seldom reported. Moreover, m6A modification in bovine skeletal muscle development and myogenic differentiation has not been reported.
The present study was aimed to uncover the m6A modification pattern in bovine skeletal myoblast differentiation and explore the potential function of m6A modification during myoblast differentiation. We used bovine skeletal myoblasts as a research target to investigate the abundance, function and mechanism of m6A modification in the process of myogenic differentiation. Thus, we performed m6A-seq and RNA-seq in pre-differentiation (GM, myoblasts) and post-differentiation (DM, myotubes) cells. Our results indicate that m6A modifications were highly enriched in mRNA, especially in the 3′UTR and CDS regions, and were likely to participate in the regulation of myogenic differentiation. Finally, we identified, screened and verified four skeletal muscle development-related genes (MYOZ2, TWIST1, KLF5 and MYOD1), which showed significant differences in both m6A methylation and mRNA expression. Our study first revealed the mRNA m6A modification map during bovine skeletal myogenic differentiation in vitro, which could contribute to further understand the roles of m6A in bovine skeletal muscle development.
Materials and Methods
Ethical Statement
The animal experiments of this study were conducted in the light of the protocol of the Experimental Animal Management Committee of Northwest A&F University (Protocol NWAFAC1120), and in accordance with the Regulations on Administration of Animals Used as Subjects of Experiments issued by the State Council of China in 2017.
Isolation and Culture of Bovine Myoblasts
Primary myoblasts were isolated from the longissimus dorsi of a newborn Qinchuan beef cattle according to the method previously reported in our laboratory (Wang et al., 2018). The isolated myoblasts were cultured to 80% confluence in growth medium, and then, myogenic differentiation was induced with differentiation medium. The culture conditions were a humidified incubator (Thermo Fisher Scientific, MA, USA) containing 5% carbon dioxide at 37°C. The myoblast growth medium was composed of Dulbecco’s modified Eagle’s medium: Nutrient Mixture F-12 (DMEM/F12; HyClone, USA), 20% fetal bovine serum (FBS; GIBCO, USA) and 1% penicillin/streptomycin. The myoblast differentiation medium consisted of DMEM/F12 containing 2% horse serum (HS; GIBCO, USA) and 1% penicillin/streptomycin. The medium was changed every 2 days.
Immunofluorescence
Cultured myoblasts and myotubes were washed briefly with PBS and fixed with PBS containing 4% paraformaldehyde for 20 min at room temperature, and then permeabilized with PBS containing 0.2% Triton X-100 (Solarbio, Beijing, China) for 10 min. The cells were subsequently washed with PBS 3 times. The cells were blocked with 10% donkey serum, 1% BSA and 0.3 M glycine in PBS at room temperature for 30 min. The primary antibodies were diluted to different concentrations in blocking buffer according to the protocols, and then, the cells were incubated overnight at 4°C. After washing three times with PBS, the cells were incubated with fluorescent dye-conjugated secondary antibodies diluted in blocking buffer for 1 h at room temperature, and this step was performed in the dark. The cells were washed 3 times with PBS, stained with 0.1% DAPI (Sigma-Aldrich, USA) for 15 min and then visualized under a fluorescence microscope (Olympus IX71, Japan). The primary antibodies used were anti-PAX7 (1:200, ab187339, Abcam) and anti-MyoD1 (1:200, ab16148, Abcam). The secondary antibodies used were Alexa Fluor 555-conjugated donkey anti-rabbit IgG (1:1,000, ab150074, Abcam) and Alexa Fluor 488-conjugated goat anti-mouse IgG (1:1,000, ab150113, Abcam).
RNA Isolation and Fragmentation
Proliferating myoblasts (named GM; 80% confluence, cultured in GM) and differentiated myotubes (named DM; cultured in DM for 4 days) were harvested, and total RNA was extracted using RNAiso reagent (TaKaRa, Dalian, China). Then the RNA samples were sent to LC-BIO Bio-tech ltd. (Hangzhou, China) for RNA sequencing and m6A sequencing. The quality and quantity of the total RNA were analyzed by Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, CA, USA) with RIN >7.0. Poly (A) mRNA was isolated from total RNA over 200 ug using poly-T oligo attached magnetic beads (Invitrogen, Massachusetts, USA). Following purification, the poly(A) mRNA fractions are fragmented into about 100-nt-long oligonucleotides by divalent cations at high temperatures.
m6A Immunoprecipitation, Library Construction and Sequencing
The cleaved RNA fragments were incubated with m6A-specific antibody (202003, Synaptic Systems, Germany) in IP buffer (50 mM Tris-HCl, 750 mM NaCl and 0.5% Igepal CA-630) supplemented with BSA (0.5 μg/μl) at 4°C for 2 h. Then the mixture was incubated with protein-A beads and eluted with elution buffer (1 × IP buffer and 6.7 mM m6A). The eluted RNA was precipitated with 75% ethanol. According to the strand-specific library prepared by dUTP method, the eluted fragments containing m6A (IP) and the untreated input control fragments were converted to final cDNA library. The average insert size for the paired terminal libraries was ∼100 ± 50 bp. And then we performed the paired-end 2 × 150 bp sequencing on an Illumina Novaseq™ 6000 platform at the LC-BIO Bio-tech ltd. (Hangzhou, China) following the vendor’s recommended protocol.
RNA-Seq Data Analysis
RNA-seq alignment. Trimmomatic (Bolger et al., 2014) were used to remove the contained adaptor contamination and low-quality bases, Then, we used the fastQC software to verify the sequence quality of each sample. We mapped valid reads to the reference genome of Bos taurus (ARS-UCD1.2) published on Ensembl website using HISAT2 (Kim et al., 2015). Then StringTie (Pertea et al., 2015) was accessed to quantify the expression level of all genes and transcripts by calculating FPKM [total exon fragments/mapped reads (millions) × exon length (kB)]. And reads were counted by featureCounts software (Liao et al., 2014).
Differential analysis. Differential expression analysis was performed using edgeR (https://bioconductor.org/packages/edgeR), a Bioconductor package in R, and marked with significant parameter. The differentially expressed genes (DEGs) were selected with Fold Change (FC) > 1.5 or FC < −1.5 and p value <0.05. The R package ggplot2 was used to generate differential volcano maps.
GO and KEGG analysis. Gene Ontology (GO) classification and enrichment analysis was performed by g:Profiler online tool (Raudvere et al., 2019). Kyoto Encyclopedia of Genes and Genomes (KEGG) was analyzed by KOBAS online website (Bu et al., 2021). Bonferroni and Hochberg corrected p value <0.05 in individual genes was considered to be statistically significant.
m6A-Seq Data Analysis
m6A-seq alignment and peak calling. After obtaining raw sequence data, we first conduct quality control for raw reads by fastQC and trim low-quality reads and adapter sequences by Trimmomatic (Bolger et al., 2014). Trimmed reads were called clean reads, and were aligned to the Bos taurus reference genome (ARS-UCD1.2) using HISAT2 software (Kim et al., 2015). The Uniquely aligned sequences were extracted by Sambamba (Tarasov et al., 2015) and only the uniquely mapped and non-duplicated alignments were further analyzed. MACS2 (Zhang et al., 2008) was used to identify the m6A-modification peaks of each sample with the default parameters, which identifies m6A peaks with bed or bam format that can be adapted for visualization on the Integrative Genomics Viewer (IGV) software (http://www.igv.org/). De novo and known motif were found using MEME (Bailey et al., 2009) and HOMER (Heinz et al., 2010), and Perl scripts in house were used to locate the motif with respect to peak summit. Meanwhile, the input RNA sequencing (RNA-seq) data were used as the background when calling peaks.
Overall analysis of m6A data. Deeptools (Ramirez et al., 2016) was used to analyze the correlation between samples and reads enrichment signal. ChIPseeker (Yu et al., 2015) were used to annotate the peaks. Bedtools (Quinlan and Hall, 2010) was used to count the peak number of each bin, and the counts were employed to plot the patterns by R. Next, Guitar (Cui et al., 2016) was used to examine the distribution pattern of the m6A peaks throughout different regions of the transcripts, the mRNA transcripts were divided into five non-overlapping segments: the 5′UTR, start codon (100 nucleotides centered on the start codon), CDS, stop codon (100 nucleotides centered on the stop codon), and 3′UTR. Each area was separated into 20 bins. Circos analysis was performed using the OmicStudio tools at https://www.omicstudio.cn/tool/.
Analysis of differential peaks between GM and DM. Correlation Heatmap was generated by DiffBind (Stark and Brown, 2011). All m6A peaks in growth and differentiation conditions for myoblasts were identified by the exomePeak (Meng et al., 2014). For differential analysis, we retrieved all peaks with >1.5-fold differences for downstream analysis. The analysis of GO and KEGG was the same as RNA-seq. Finally, we analyzed the correlation between m6A abundance and mRNA level during myoblast differentiation using R packages. The scatter plots of correlation were generated using ggplot2 (v3.3.5), the Pearson correlation coefficient was calculated using ggpubr (v0.4.0), and the fitting equation was added to the plots by ggpmisc (v0.4.2).
cDNA Synthesis and Real-Time Quantitative PCR
The PrimeScript RT reagent kit (TaKaRa) was used to synthesize cDNA. The residual genomic DNA was removed at 42°C for 2 min, and then, the reverse transcription reaction was conducted at 37°C for 15 min and then at 85°C for 5 s. Real-time quantitative PCR (RT-qPCR) was performed using the TB Green Premix Ex Taq II Kit (TaKaRa) and a CFX Connect Real-Time PCR Detection System (BIO-RAD, CA, USA). Bovine GAPDH was used as the internal reference to standardize the data. Each sample analyzed by RT-qPCR was subjected to at least three biological repeats. Relative mRNA expression was calculated and analyzed using the 2−ΔΔCt method (Livak and Schmittgen, 2001). All of the primers used in the RT-qPCR are listed in Supplementary Table S1.
m6A-IP-qPCR
m6A immunoprecipitation assays were performed as previously described (Dominissini et al., 2012). In brief, 48 h after transfection, RNA from the cells was chemically digested into 200-nt fragments, and more than 200 μg of total RNA was subjected to immunoprecipitation with affinity-purified m6A-specific antibody (202003, Synaptic Systems, Germany). The RNA fragments that bound to m6A were separated by TRIzol reagent. Following ethanol precipitation, the input RNA and eluted m6A RNA were reverse transcribed by random hexamers, and then, the enriched sequences were detected by RT-qPCR. The ΔΔCt between the 10% input and the immunoprecipitated RNA was determined, and the relative enrichment was calculated as 2−ΔΔCt. The primers used to amplify the m6A peak region are listed in Supplementary Table S1.
Statistical Analysis
All data were presented with the means ± standard deviation (SD) of at least three biological repeat samples. Student’s t-test (between two groups) or ANOVA (among multiple groups) were used to compare significance of the mean values. Differences were considered to be very significant or significant at p < 0.01 or p < 0.05, respectively. The results were analyzed using GraphPad Prism 7.00 (CA, USA) software and images were generated.
Results
Identification of Bovine Skeletal Myoblasts
To verify whether the cells isolated from bovine longissimus dorsi muscles could undergo myogenic differentiation, we seeded the isolated cells in culture dishes, grew them to 80–90% confluence, and passaged them to a 6-well plate for further culture. After 48 h in the growth medium, immunofluorescence showed that PAX7 and MYOD1 were simultaneously expressed, while the expression of MYOD1 was relatively low (Figure 1A). Therefore, we preliminarily identified the isolated cells as myoblasts. The myoblasts were cultured in growth medium to 90% confluence, and myogenic differentiation was induced with differentiation medium. Microscopic observation revealed that the myotubes formed by myoblast fusion gradually increased and became longer as the number of days of differentiation increased (Figure 1B). We evaluated the differentiation status of the myoblasts by detecting the mRNA levels of MYOD1, MYOG, MYH3 (myosin heavy chain 3), MYMK (myomaker, myoblast fusion factor), MRF4 (myogenic regulatory factor 4) and CKM (creatine kinase), which are widely recognized marker genes of differentiated myoblasts and fused myotubes (Bentzinger et al., 2012; Millay et al., 2013). Notably, the levels of MYOG, MYH3, MYMK, MRF4 and CKM gradually increased during myogenic differentiation, and the trends in the expression of these pivotal genes were consistent with the differentiation stage (Figure 1C). Alternatively, the levels of MYOD1 peaked on D2 (Figure 1C), which was consistent with previous studies showed that MYOD1 plays a vital role in the proliferation and early differentiation of myoblasts (Weintraub et al., 1991; Weintraub, 1993). These results indicated that the isolated bovine skeletal myoblasts in this study were capable of myogenic differentiation and could serve as a model for our follow-up study.
FIGURE 1. Identification of bovine skeletal myoblasts. (A) Identification of myoblasts based on PAX7 and MYOD1 expression in growth medium (GM) (scale bar: 20 μm). (B) Observation of myotube formation on days 0, 1, 2, 3 and 4 of culture in differentiation medium (DM) after the induced myogenic differentiation (scale bar: 100 μm). (C) Relative mRNA expression of specific myogenic genes (MYOD1, MYOG, MYH3, MYMK, MRF4 and CKM) during bovine skeletal myoblast differentiation. The results were normalized to the GAPDH levels and were presented as the means ± SD. Different capital letters indicate very significant differences (p < 0.01), different lowercase letters indicate significant differences (p < 0.05), and the same letters indicate no significant differences (p > 0.05).
Overview of Samples and N6-Methyladenosine Methylation of mRNA in Bovine Sketetal Myoblasts and Myotubes
To investigate the role of m6A in myoblast differentiation, mRNA was extracted from pre-differentiation (GM, myoblasts, D0) and post-differentiation (DM, myotubes, D4) cells for m6A-seq and RNA-seq. Pearson correlation analysis showed that there was a strong correlation between the three biological repeat samples in GM and DM groups, respectively (Figure 2A). RNA-seq and m6A-seq produced 46785592–66324558 clean reads in Input and IP groups, of which more than 95% were mapped to reference genome of Bos taurus (Table 1). After eliminating low-quality reads, more than 92% of unique mapped reads were obtained in each group of clean reads (Table 1). These data show that high-throughput sequencing in this study has been carried out successfully. Furthermore, we identified 23287–23912 total transcripts in GM and DM groups (Table 2, Supplementary Table S2). Then the methylated mRNA was mapped to the transcriptome. 6,042–8,094 m6A-modified transcripts were identified in GM groups and 5,691–7,361 m6A-modified transcripts were found in DM groups. There were 9,246–13596 and 9,128–11861 m6A peaks in the two groups, respectively (Table 2, Supplementary Table S3). Therefore, we found that there were ∼1.60 m6A peaks per m6A transcript in bovine skeletal myoblasts and myotubes, which were similar to those detected in human HepG2 (∼1.7 peaks per gene) (Dominissini et al., 2012), pig muscle (∼1.7 peaks per m6A transcript) (Jiang et al., 2019), pig adipose tissue (∼1.3 peaks per gene) (Tao et al., 2017) and chicken fat (∼1.5 peaks per m6A transcript) (Cheng et al., 2021). The results of the three biological replicates were similar, which also showed the accuracy of m6A-seq results and the similarity between duplicate samples. Summary plots and heatmaps were generated by deepTools using normalized read coverages from m6A-seq in bovine skeletal myoblasts (Figure 2B). The Summary plots on top of the heatmap (Figure 2B, top panel) indicate that the enrichment intensity between and near TSS (transcription start site) and TES (transcription end site) in the genome region of the immunoprecipitation (IP) group is higher than that in the Input group, and peaks signal is most strongly enriched around the start and end of genes which is also visible in the heatmaps. Clustering of significantly enriched sequences found m6A consensus motif RRACH was highly enriched in GM and DM (Figure 2C), which was consistent with the general pattern of mammals.
FIGURE 2. Overview of samples and m6A methylation of mRNAs in bovine skeletal myoblasts and myotubes. (A) Pearson correlation analysis showing the correlation between the groups, represented as values in modules. (B) Distribution of peaks on gene functional elements in all groups. (C) Sequence motifs identified within m6A peaks. (D) Metagene profiles of m6A peaks distribution across the 5′ UTR, CDS, and 3′ UTR in GM and DM. (E) Percentage of m6A peak in 5′ UTR, CDS, and 3′ UTR of mRNA. (F) Percentage of m6A-methylated genes with different number of m6A peaks. (G) Circos plot showing the global distribution of m6A, in cow genome (left) and chromosome 3 (right) of GM and DM. (H) Cumulative curves of m6A abundance in GM and DM.
TABLE 1. Summary of sequencing data and read-alignment statistics from m6A-seq in bovine skeletal myoblasts and myotubes.
TABLE 2. Number of m6A peaks detected in bovine skeletal myoblasts and myotubes.
To investigate the modification position of m6A in transcripts, we then studied the metagene profiles of the m6A peak in the entire transcriptome of GM and DM (Figure 2D). Three distinct m6A peaks were observed in the start codon, stop codon and 3′UTR (untranslated region). Simultaneously, the peak near the stop codon is significantly higher than the other two peaks. Interestingly, the m6A density of DM in 3′UTR is higher than GM, while the density in CDS is lower than GM. Further, to evaluate m6A enrichment systematically in mRNAs, we calculated the enrichment proportion of m6A peak in 5′UTR, CDS (coding sequence) and 3′UTR. More than 50% of m6A peaks are in the CDS region, nearly 40% are in the 3′UTR, and just under 10% are in the 5′UTR (Figure 2E). Moreover, the distribution pattern of mRNA m6A modification was highly similar in myoblasts and myotubes, and coincided with the typical m6A peak distribution in mammals. To further determine the distribution of m6A modifications in the transcriptome, we analyzed the number of m6A peaks contained in each gene. There was only one m6A peak in more than 60% of the genes, and ∼23 and 9% of the genes with 2 and 3 peaks, respectively. Only about 2–4% of the genes contain more than 4 peaks (Figure 2F). Furthermore, to investigate the abundance of m6A modification at the chromosome level and the difference between GM and DM, the circle diagram of m6A peaks was produced using OmicStudio online tools (https://www.omicstudio.cn/tool) (Figure 2G). The circular visualization of peaks on chromosomes suggested that m6A methylation is extensively modified in myoblasts and myotubes, but there is an inconspicuous difference in m6A abundance between GM and DM. Cumulative curves also showed that there was no visible difference in m6A abundance between GM and DM (Figure 2H).
Analysis of Differentially Methylated Peaks (DMPs) Between Bovine Skeletal Myoblasts (GM) and Myotubes (DM)
Correlation heatmap showed that the immunoprecipitated samples were obviously divided into two treatment groups, GM_IP and DM_IP, with weak intra-group differences and significant inter-group differences visible (Figure 3A). To investigate the difference in the abundance of methylated peaks between the two groups, firstly, we identified 7,140 common peaks in both DM and GM (Supplementary Table S4). Then, 1,559 differentially methylated peaks (DMPs) were further screened, representing the annotated 1,437 differentially methylated genes (DMGs). The expression of 949 peaks augmented (correspond to 881 genes with up-regulated m6A abundance) and 610 peaks diminished (correspond to 571 genes with down-regulated m6A abundance) (Figure 3B, Supplementary Table S4). The analysis of DMPs enrichment sites showed that 44.39% of these DMPs were enriched in 3′UTR and about 41% in the exon region, while 14.33% of m6A modification occurs in 5′UTR (Figure 3C).
FIGURE 3. Analysis of differentially methylated peaks (DMPs) between bovine skeletal myoblasts and myotubes. (A) Correlation heatmap showing the correlation between GM_IP and DM_IP. (B) MA plot showing DMPs between DM_IP and GM_IP. (C) Pie chart showing the fraction of m6A peaks in 4 non-overlapping transcript segments (5′ UTRs, 1st exon, other exons, and 3′ UTRs). (D) Gene ontology (GO) enrichment analysis of molecular function for differentially methylated genes (DMGs). (E) GO analysis of biological process for DMGs. (F) KEGG analysis of DMGs.
Furthermore, to analyze the potential function of m6A modified genes in myoblast differentiation, we performed GO [terms of molecular function (MF) and biological process (BP)] and KEGG pathway enrichment analysis for DMGs. As shown in Figures 3D,E, DMPs were mainly clustered in GO-MFs related to transcriptional regulation, including DNA binding, protein binding, transcription regulator activity and cis-regulatory region sequence (Supplementary Table S5). Likewise, they were also enriched in the GO-BPs of transcriptional regulation. Moreover, they were mainly involved in nucleic acid metabolic-related and biosynthetic-related Biological Processes. KEGG analysis indicated that DMGs enriched in disease-related signal pathways, but also participated in some signal pathways related to myogenesis, including ubiquitin mediated proteolysis, insulin resistance, cell cycle and metabolic pathways (Figure 3F, Supplementary Table S6). These results suggested that DMPs (representing DMGs) might play a role in gene transcription regulation and cell metabolism during bovine myoblast differentiation.
Analysis of Differentially Expressed Genes (DEGs) Between Bovine Skeletal Myoblasts and Myotubes
m6A abundance has been reported to affect mRNA levels (Wang et al., 2014; Shen et al., 2016; Zhao et al., 2017; Zhou et al., 2019). To evaluate whether there is a potential correlation between m6A mRNA methylation and gene transcript levels during myoblast differentiation, RNA-seq analyses were simultaneously performed in all samples. The Volcano Plot of DEGs data were shown in Figure 4A. 2257 DEGs were identified using RNA-seq, of which 986 genes were up-regulated and 1,271 genes were down-regulated (p < 0.05, FC > 1.5) (Figure 4A, Supplementary Table S7). To determine the clustering pattern of genes under different experimental conditions, we took the top 50 up-regulated and top 50 down-regulated genes with the lowest p-value for differential gene cluster analysis (heat map) on Omicstudio tools at https://www.omicstudio.cn/tool. The hierarchical clustering of DEGs data were displayed in Figure 4B. We found that MEF2C, MYL1, MYOM2, ACTN2, MYH8 and other genes related to myoblast differentiation and myotube formation were significantly up-regulated, while the expression of genes inhibiting muscle development was significantly decreased, including ID1 (Jen et al., 1992), IGFBP5 (Salih et al., 2004), SOX8 (Schmidt et al., 2003) and PRC1 (Bracken et al., 2006).
FIGURE 4. Analysis of differentially expressed genes (DEGs) between bovine skeletal myoblasts and myotubes. (A) The Volcano plot of DEGs. (B) Hierarchical clustering heatmap of DEGs. (C) GO analysis of biological process for up-regulated DEGs. (D) KEGG analysis of up-regulated DEGs. (E) GO analysis of biological process for down-regulated DEGs. (F) KEGG analysis of down-regulated DEGs.
To further reveal the role of DEGs in myoblast differentiation, we performed GO and KEGG analysis on up-regulated and down-regulated DEGs, respectively. The results indicated, on the one hand, that the up-regulated genes are mainly enriched in biological processes closely related to myogenesis, such as muscle contract, muscle tissue development and muscle structure development (Figure 4C, Supplementary Table S8). The major signaling pathways involved include cardiomyopathy, insulin secret, cardiac muscle contract, calcium signaling pathway, Wnt signaling pathway, PI3K-Akt signaling pathway and cardiology pathway, which directly regulated muscle development (Figure 4D, Supplementary Table S9). On the other hand, the down-regulated genes were mainly clustered in cell cycle-related and cell division-related biological processes. Similarly, they were enriched in related signaling pathways of Cell cycle and Cell division, including DNA replication, cell cycle, and Rap1 signaling Pathway (Figures 4E,F, Supplementary Tables S8, 9). Meanwhile, the TNF signaling pathway and PI3K-Akt signaling pathway involved in the regulation of myogenic differentiation were enriched. These data verified that RNA-seq results were consistent with the biological process of the samples (GM and DM).
Integrated Analysis of m6A-Seq and RNA-Seq Data
To demonstrate the potential regulation of m6A modification on gene expression in bovine myoblast differentiation, we analyzed genes with significant changes in both mRNA and m6A levels. As shown in Figure 5A, the Venn diagram of DEGs and DMPs found that a total of 268 genes had significant changes at both levels, accounting for 18.4% of 1,438 DMGs and 11.9% of 2257 DEGs. This result implied that m6A modification may regulate the expression of these genes during the course of myogenic differentiation. Additionally, the overlapping results of DEGs and DMPs showed that there were 67 common genes in both “m6A_up” and “mRNA_up” (means hyper-up), 83 common genes in both “m6A_up” and “mRNA_down” (means hyper-down), 65 common genes in both “m6A_down” and “mRNA_up” (means hypo-up), and 46 common genes in both “m6A_down” and “mRNA_down” (means hypo-down) (Figure 5A, Supplementary Table S10). Intriguingly, m6A modification abundances of 15 genes exhibited both up-regulation and down-regulation, indicating that multiple peaks with m6A methylation modification within these genes had different significant changes in myoblast differentiation.
FIGURE 5. Conjoint analysis of m6A-seq and RNA-seq data. (A) Venn diagram showing the overlapping genes between DMPs and DEGs. (B) GO analysis of molecular function for common genes between DMGs and DEGs. (C) GO analysis of biological process for common genes between DMGs and DEGs. (D) KEGG analysis of common genes between DMGs and DEGs. (E) Heatmap of 15 genes mRNA expression in Table 3. (F) Heatmap of 15 genes mRNA expression in Table 4.
Likewise, GO and KEGG analyses were performed to predict the role of the genes with significant changes in both m6A and mRNA levels during myoblast differentiation. GO-MF analysis showed that these genes were predominantly concentrated in binding-related pathways, including protein binding, enzyme binding, DNA binding, ion binding (Figure 5B, Supplementary Table S11). This result was similar to our results in DMPs (Figure 3D), suggesting that m6A modification may participate in the regulation of transcripts. Meanwhile, the results of the GO-BP analysis revealed that the genes were enriched in signal transduction, cellular process, biological process and developmental process (Figure 5C, Supplementary Table S11). KEGG analysis indicated that these genes were involved in signaling pathways related to muscle development such as cardiomyopathy and insulin resistance, besides being enriched in disease and cancer-related pathways (Figure 5D, Supplementary Table S12). As shown in Tables 3, 4, in view of the results of differential integration analysis, we listed top 15 genes with the highest up-regulated and down-regulated of m6A level (with significant mRNA expression level (p < 0.05, FC > 1.5 and FPKM >0.5) in all three biological replications of GM and DM), and Heatmaps showed differences in their mRNA expression (Figures 5E,F). Strikingly, we found that m6A-modified peaks of 24 of the 30 genes were enriched in 3′UTR and exons, which supported the results of the m6A-enriched region in this study and suggested that m6A may influence bovine myoblast differentiation by mediating the expression of these genes.
TABLE 3. Top 15 up-regulated m6A methylated genes between DMGs and DEGs.
TABLE 4. Top 15 down-regulated m6A methylated genes between DMGs and DEGs.
Finally, to determine whether m6A modification could regulate gene expression in bovine skeletal myoblasts, we analyzed the correlation between mRNA level and m6A peaks abundance in bovine myoblasts (GM) and myotubes (DM). Our results revealed that mRNA expression was strongly positively correlated with the abundance of m6A peaks in each of GM and DM (Figures 6A,B). However, there was no striking correlation between the change of m6A level and the differential expression of mRNA during myoblast differentiation (Figure 6C). These data suggested that the higher the expression level of the genes, the higher the abundance of m6A methylation modification. However, there was no significant positive or negative correlation between the changes of mRNA m6A methylation and differential expression during bovine skeletal myoblast differentiation.
FIGURE 6. Correlation analysis of m6A peak enrichment and mRNA abundance in bovine myoblasts and myotubes. (A) Scatter plots of m6A peak enrichment and mRNA abundance in three biological replicates of GM showing the same tendency. (B) Scatter plots of m6A peak enrichment and mRNA abundance in three biological replicates of DM showing the same tendency. (C) Scatter plots showing a weak correlation between the fold changes in m6A methylation and mRNA expression levels for genes.
Validation of m6A and mRNA Levels of 4 Specific Myogenesis-Related Genes in Bovine Myoblast Differentiation
In the m6A-seq and RNA-seq analysis, we found that transcripts of several well-known myogenesis-related genes, including MYOZ2, TWIST1, KLF5 and MYOD1, exhibited differential both m6A and mRNA levels in bovine myoblast differentiation. MYOZ2 is a muscle-specific protein in the Z-band of the sarcomere, which plays an important role in maintaining muscle fiber structure and myotube formation (Takada et al., 2001). TWIST1 is a transcription factor that inhibits skeletal myocyte differentiation (Spicer et al., 1996), while KLF5 and MYOD1 are key transcription factors that promote myogenesis (Hayashi et al., 2016).
We generated the mRNA m6A peaks of these four genes by using IGV software (Figure 7A). In all three biological replicates of GM and DM, m6A-modified peaks of MYOZ2 and TWIST1 were observed to be enriched near the start codon in mRNA, while m6A-modified peaks of KLF5 and MYOD1 were enriched within CDS and 3′UTR, respectively (Figure 7A). We performed m6A-IP-qPCR to verify the results of m6A-seq and indicated that the mRNAs of MYOZ2 and TWIST1 displayed higher levels of m6A enrichment in DM compared with GM, and the mRNAs of KLF5 and MYOD1 exhibited lower levels of m6A enrichment in DM (Figure 7B). Additionally, RNA-seq and RT-qPCR revealed and verified the transcript levels of MYOZ2 and KLF5 were enhanced significantly in DM compared with GM, and TWIST1 and MYOD1 were diminished significantly in DM (Figure 7C). These results confirmed the accuracy of the m6A-seq and RNA-seq data. What’s more, these data provided possible mechanisms for m6A methylation to regulate myoblast differentiation, and these four genes may also be considered as downstream candidate targets for m6A modification during myoblast differentiation.
FIGURE 7. Changes in m6A and mRNA levels in transcripts of 4 myogenesis-related genes during bovine myoblast differentiation. (A) IGV tracks displaying the distribution of m6A peaks in MYOZ2, TWIST1, KLF5 and MYOD1 transcripts in all six groups. (B) Validations of the m6A enrichment by m6A-IP-qPCR. (C) mRNA expression levels of MYOZ2, TWIST1, KLF5 and MYOD1 were revealed and verified by RNA-seq and RT-qPCR, respectively. Data were presented as the means ± SD. *p < 0.05, ***p < 0.001, **p < 0.01; using Student’s t test.
Discussion
The discovery of the first m6A demethylase FTO (Jia et al., 2011) and the development of m6A-specific high-throughput sequencing technology (Dominissini et al., 2012; Meyer et al., 2012) provided theoretical basis and technical support for the study of m6A in plants and animals, respectively. Related studies in farm animals have also been uncovered in recent years. However, as far as we know, the research on m6A modification in cattle has not been reported, and the enrichment pattern and potential role of m6A methylation in beef cattle muscle development were still unclear. In the present study, we identified transcriptome-wide N6-methyladenosine profiling of bovine skeletal myoblast differentiation by performing m6A-seq and RNA-seq. Bioinformatics analysis suggested the potential role of m6A modification in regulating myoblast differentiation. Further experiments verified the accuracy of the sequencing results, and screened some myogenesis-related genes with obvious m6A methylation modification.
First of all, we verified that the bovine skeletal muscle myoblasts isolated in this study could carry out normal myoblast differentiation by observing myotube formation and using RT-qPCR. The results of RNA-seq later indicated that DEGs were mainly involved in the process of muscle growth and development, which also proved the reliability of sequencing. Our results showed that m6A methylation occurred in 24.3–33.9% transcripts during bovine skeletal myoblast differentiation, and these transcripts with m6A modification contained more than one m6A peak on average. The results suggest that m6A is widely modified in bovine myoblasts and may play an important role in skeletal myogenesis. Similar to other mammals, m6A peaks in myoblasts and myotubes were also mainly enriched in motif RRACH. This result was consistent with the findings in pigs (Tao et al., 2017; Jiang et al., 2019), chickens (Cheng et al., 2021) and geese (Xu et al., 2021), but different from the results in plants (Wei et al., 2018; Zhou et al., 2019), suggesting that m6A modification in mammals may be similar to birds, but different from plants. Exploring the sequence and location of motifs could provide a reference for subsequent molecular mechanism studies. The motifs are catalyzed by m6A methylases and recognized and bound by m6A reader proteins. Many studies have shown that m6A modification is reduced after specific mutation of m6A motif (Guo et al., 2020; Ye et al., 2020), which would provide support for subsequent investigations into the potential role of m6A modification. In our study, four genes were selected for m6A-IP-qPCR validation, and the primers used were also designed in view of m6A motif positions.
The distribution of m6A peaks in bovine myoblasts is similar to that in humans, mice and pigs, mainly abundant near the stop codon, CDS and 3′UTR (Dominissini et al., 2012; Tao et al., 2017; Gheller et al., 2020; Zhang et al., 2020), but is inconsistent with the result that m6A peaks were mainly distributed near the start codon in chicken fat and goose muscle (Cheng et al., 2021; Xu et al., 2021). The results demonstrated a difference in the distribution of m6A methylation in mammals and birds. Previous studies have shown that m6A modification in different regions of mRNA may have different mechanisms. m6A within the 5′UTR regulated cap-independent translation in stress response, and CDS and 3′UTR m6A was recognized by YTHDF2 or IGF2BPs to degrade and stabilize the target mRNAs, respectively (Shi et al., 2019). Besides, YTHDF1 was more inclined to bind to the m6A site of 3′UTR to promote translation (Shi et al., 2019). Our results revealed that among the TOP15 genes with the highest up-regulated or down-regulated m6A level, a total of 24 genes had m6A peaks distribution in the 3′UTR and CDS regions, and m6A modification occurred in 5′UTR of the only six genes. As shown in Table 4, CKM mRNA was upregulated in the bovine myoblast differentiation, while the mRNA m6A abundance was diminished. It is reasonable to speculate that YTHDF1 or IGF2BPs may recognize the m6A site in 3′UTR of CKM mRNA and promote its mRNA translation or stabilization. CKM is a terminal differentiation gene for myogenesis (Millay et al., 2013), suggesting that m6A may influence myoblast differentiation by mediating the m6A level of CKM. Subsequent experiments will be conducted to verify the molecular mechanism.
GO analysis revealed that DMGs were commonly enriched in DNA binding, protein binding and transcription-related terms of Molecular Function, and Biological Process enrichment analysis demonstrated that DMGs were mainly associated with transcriptional regulation, nucleic acid metabolism and RNA metabolism. These results are consistent with previous results of m6A modification in C2C12 myoblast differentiation (Gheller et al., 2020) and porcine skeletal muscle development (Tao et al., 2017; Zhang et al., 2020), suggesting that m6A modification may be involved in gene transcription regulation and RNA metabolism during myoblast differentiation. KEGG analysis in this study showed that DMGs were involved in insulin resistance and metabolic pathways, among which STAT3, JAK2, IGF2, CKM, PIK3CA and other genes that regulate skeletal muscle development were scanned, implying a new potential mechanism of myogenesis-related genes regulating the bovine myoblast differentiation. The previous study has found that the m6A level of both C2C12 and primary mouse myoblasts on the third day of differentiation was significantly lower than that in the proliferation phase (Gheller et al., 2020). In our study, there was no significant difference in the m6A abundance of GM and DM, which may be related to the interference of some genes with low abundance on the sequencing results. Subsequent experiments, such as LC-MS/MS or dot blot, are required to determine the changes of m6A level during bovine myoblast differentiation.
The associated analyses of m6A-seq and RNA-seq revealed that there was a strong positive correlation between mRNA m6A abundance and expression level in each GM and DM group. The higher the gene expression level, the higher the m6A abundance. The result is contrary to the idea that m6A and mRNA expression in chicken adipose tissue was negatively correlated (Cheng et al., 2021). We speculate that it may be attributed to differences in species and cell types. Notably, integrated analysis for differential expression showed that there was no marked correlation between fold changes of mRNA expression and fold change of m6A level during bovine myoblast differentiation. Inconsistent with our findings, changes in m6A methylation were negatively correlated with gene expression during goose embryonic muscle development (Xu et al., 2021). We analyzed the m6A-seq data of chicken fat deposition and found that there was a positive correlation between differential m6A abundance and differential gene expression level (Cheng et al., 2021). In addition, the two modules with positive or negative correlation with body traits were screened by WGCNA during porcine embryonic muscle development, and the combined analysis displayed that there was either a significant positive or negative correlation between m6A methylation changes and gene expression multiples in the individual module (Zhang et al., 2020). These results suggest that the association between mRNA m6A methylation abundance and gene expression level may be conservative in species, but also related to the differences between tissue level in vivo and cellular level in vitro.
Finally, the m6A level and mRNA relative expression of four well-known myogenesis-related genes were verified by m6A-IP-qPCR and RT-qPCR. The differentially expressed m6A peaks in TWIST1 mRNA were enriched near 5 ′UTR, and the m6A level was up-regulated while the mRNA level was down-regulated, suggesting that m6A may inhibit TWIST1 mRNA expression via promoting m6A methylation of TWIST1 mRNA to promote myoblast differentiation due to TWIST1 inhibited skeletal muscle differentiation (Spicer et al., 1996). Previous studies showed that no m6A modification was found in Myod1 mRNA in mESC (Batista et al., 2014), whereas m6A methylation of Myod1 in C2C12 was significantly enriched in 5′UTR, and siMETTL3 led to a restraint of myoblast differentiation by reducing MYOD1 mRNA expression (Kudou et al., 2017). In contrast, our study found that m6A modification of MYOD1 mRNA was enriched in 3′UTR and its mRNA expression was decreased, which may have different or similar regulatory mechanisms in bovine myoblast differentiation. In view of previous reports and our data, it is reasonable to speculate that the role of m6A modification in skeletal muscle differentiation may not be a pure promoting or inhibiting mechanism, but rather may be complex. Our ongoing research project will be to investigate the molecular mechanism of m6A modification regulating skeletal myogenesis in vitro using bovine myoblasts as a model.
Conclusion
In conclusion, we first revealed the mRNA m6A modification map in bovine skeletal myoblasts and myotubes. We found that m6A methylation may regulate myogenesis via mediating the gene expression. Further, four candidate target genes were identified and screened, including MYOZ2, TWIST1, KLF5 and MYOD1. These comprehensive analyses open a new perspective for the genetic improvement and molecular breeding of beef cattle, and provide a theoretical basis for studying the functional and molecular mechanism of m6A methylation in regulating skeletal muscle development and myogenesis.
Data Availability Statement
The datasets supporting our findings are included in the article and Supplementary Material. The m6A-seq and RNA-seq data generated by this study have been deposited to GEO database under accession number GEO: GSE173477.
Ethics Statement
The animal study was reviewed and approved by the animal experiments of this study were conducted in the light of the protocol of the Experimental Animal Management Committee of Northwest A&F University (Protocol NWAFAC1120), and in accordance with the Regulations on Administration of Animals Used as Subjects of Experiments issued by the State Council of China in 2017. Written informed consent was obtained from the owners for the participation of their animals in this study.
Author Contributions
LZ and XY designed the experiment. XY performed experiments, analyzed the data and wrote the manuscript. JW carried out most of the bioinformatic analysis. XM and JD conducted partial cellular experiments. CM contributed to discussions and draft writing. All authors read and approved the final manuscript.
Funding
This work was supported by the National Key Research and Development Program of China (2018YFD0501700), National Natural Science Foundation of China (31972994), Key Research and Development Program of Ningxia Province (2019BEF02004), National Beef and Yak Industrial Technology System (CARS-37), Agricultural Science and Technology Innovation and Transformation Project of Shaanxi Province (NYKJ-2018-LY09).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
We are grateful to our colleagues of Zan’s laboratory who provided expertise and contributed to this research.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.785380/full#supplementary-material
References
Alarcón, C. R., Goodarzi, H., Lee, H., Liu, X., Tavazoie, S., and Tavazoie, S. F. (2015). HNRNPA2B1 Is a Mediator of m6A-Dependent Nuclear RNA Processing Events. Cell 162 (6), 1299–1308. doi:10.1016/j.cell.2015.08.011
Bailey, T. L., Boden, M., Buske, F. A., Frith, M., Grant, C. E., Clementi, L., et al. (2009). MEME SUITE: Tools for Motif Discovery and Searching. Nucleic Acids Res. 37 (Web Server issue), W202–W208. doi:10.1093/nar/gkp335
Batista, P. J., Molinie, B., Wang, J., Qu, K., Zhang, J., Li, L., et al. (2014). m6A RNA Modification Controls Cell Fate Transition in Mammalian Embryonic Stem Cells. Cell Stem Cell 15 (6), 707–719. doi:10.1016/j.stem.2014.09.019
Bentzinger, C. F., Wang, Y. X., and Rudnicki, M. A. (2012). Building Muscle: Molecular Regulation of Myogenesis. Cold Spring Harbor Perspect. Biol. 4 (2), a008342. doi:10.1101/cshperspect.a008342
Bertero, A., Brown, S., Madrigal, P., Osnato, A., Ortmann, D., Yiangou, L., et al. (2018). The SMAD2/3 Interactome Reveals that TGFβ Controls m6A mRNA Methylation in Pluripotency. Nature 555 (7695), 256–259. doi:10.1038/nature25784
Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a Flexible Trimmer for Illumina Sequence Data. Bioinformatics 30 (15), 2114–2120. doi:10.1093/bioinformatics/btu170
Bracken, A. P., Dietrich, N., Pasini, D., Hansen, K. H., and Helin, K. (2006). Genome-Wide Mapping of Polycomb Target Genes Unravels Their Roles in Cell Fate Transitions. Genes Dev. 20 (9), 1123–1136. doi:10.1101/gad.381706
Braun, T., and Gautel, M. (2011). Transcriptional Mechanisms Regulating Skeletal Muscle Differentiation, Growth and Homeostasis. Nat. Rev. Mol. Cel. Biol. 12 (6), 349–361. doi:10.1038/nrm3118
Bu, D., Luo, H., Huo, P., Wang, Z., Zhang, S., He, Z., et al. (2021). KOBAS-i: Intelligent Prioritization and Exploratory Visualization of Biological Functions for Gene Enrichment Analysis. Nucleic Acids Res. 49 (W1), W317–W325. doi:10.1093/nar/gkab447
Cantara, W. A., Crain, P. F., Rozenski, J., McCloskey, J. A., Harris, K. A., Zhang, X., et al. (2011). The RNA Modification Database, RNAMDB: 2011 Update. Nucleic Acids Res. 39 (Database issue), D195–D201. doi:10.1093/nar/gkq1028
Cheng, B., Leng, L., Li, Z., Wang, W., Jing, Y., Li, Y., et al. (2021). Profiling of RNA N6-Methyladenosine Methylation Reveals the Critical Role of m6A in Chicken Adipose Deposition. Front. Cel Dev. Biol. 9, 590468. doi:10.3389/fcell.2021.590468
Cui, X., Wei, Z., Zhang, L., Liu, H., Sun, L., Zhang, S.-W., et al. (2016). Guitar: An R/Bioconductor Package for Gene Annotation Guided Transcriptomic Analysis of RNA-Related Genomic Features. Biomed. Res. Int. 2016, 1–8. doi:10.1155/2016/8367534
Dominissini, D., Moshitch-Moshkovitz, S., Schwartz, S., Salmon-Divon, M., Ungar, L., Osenberg, S., et al. (2012). Topology of the Human and Mouse m6A RNA Methylomes Revealed by m6A-Seq. Nature 485 (7397), 201–206. doi:10.1038/nature11112
Dorn, L. E., Lasman, L., Chen, J., Xu, X., Hund, T. J., Medvedovic, M., et al. (2019). The N(6) -Methyladenosine mRNA Methylase METTL3 Controls Cardiac Homeostasis and Hypertrophy. Circulation 139 (4), 533–545. doi:10.1161/CIRCULATIONAHA.118.036146
Fu, Y., Dominissini, D., Rechavi, G., and He, C. (2014). Gene Expression Regulation Mediated Through Reversible m6A RNA Methylation. Nat. Rev. Genet. 15 (5), 293–306. doi:10.1038/nrg3724
Fustin, J.-M., Doi, M., Yamaguchi, Y., Hida, H., Nishimura, S., Yoshida, M., et al. (2013). RNA-Methylation-Dependent RNA Processing Controls the Speed of the Circadian Clock. Cell 155 (4), 793–806. doi:10.1016/j.cell.2013.10.026
Gheller, B. J., Blum, J. E., Fong, E. H. H., Malysheva, O. V., Cosgrove, B. D., and Thalacker-Mercer, A. E. (2020). A Defined N6-Methyladenosine (m6A) Profile Conferred by METTL3 Regulates Muscle Stem Cell/myoblast State Transitions. Cell Death Discov. 6, 95. doi:10.1038/s41420-020-00328-5
Guo, X., Li, K., Jiang, W., Hu, Y., Xiao, W., Huang, Y., et al. (2020). RNA Demethylase ALKBH5 Prevents Pancreatic Cancer Progression by Posttranscriptional Activation of PER1 in an m6A-YTHDF2-Dependent Manner. Mol. Cancer 19 (1), 91. doi:10.1186/s12943-020-01158-w
Han, D., Liu, J., Chen, C., Dong, L., Liu, Y., Chang, R., et al. (2019). Anti-tumour Immunity Controlled through mRNA m6A Methylation and YTHDF1 in Dendritic Cells. Nature 566 (7743), 270–274. doi:10.1038/s41586-019-0916-x
Haussmann, I. U., Bodi, Z., Sanchez-Moran, E., Mongan, N. P., Archer, N., Fray, R. G., et al. (2016). m6A Potentiates Sxl Alternative Pre-mRNA Splicing for Robust Drosophila Sex Determination. Nature 540 (7632), 301–304. doi:10.1038/nature20577
Hayashi, S., Manabe, I., Suzuki, Y., Relaix, F., and Oishi, Y. (2016). Klf5 Regulates Muscle Differentiation by Directly Targeting Muscle-Specific Genes in Cooperation with MyoD in Mice. Elife 5, e17462. doi:10.7554/eLife.17462
Heinz, S., Benner, C., Spann, N., Bertolino, E., Lin, Y. C., Laslo, P., et al. (2010). Simple Combinations of Lineage-Determining Transcription Factors Prime Cis-Regulatory Elements Required for Macrophage and B Cell Identities. Mol. Cel. 38 (4), 576–589. doi:10.1016/j.molcel.2010.05.004
Jen, Y., Weintraub, H., and Benezra, R. (1992). Overexpression of Id Protein Inhibits the Muscle Differentiation Program: In Vivo Association of Id with E2A Proteins. Genes Dev. 6 (8), 1466–1479. doi:10.1101/gad.6.8.1466
Jia, G., Fu, Y., Zhao, X., Dai, Q., Zheng, G., Yang, Y., et al. (2011). N6-methyladenosine in Nuclear RNA Is a Major Substrate of the Obesity-Associated FTO. Nat. Chem. Biol. 7 (12), 885–887. doi:10.1038/nchembio.687
Jia, G., Fu, Y., and He, C. (2013). Reversible RNA Adenosine Methylation in Biological Regulation. Trends Genet. 29 (2), 108–115. doi:10.1016/j.tig.2012.11.003
Jiang, Q., Sun, B., Liu, Q., Cai, M., Wu, R., Wang, F., et al. (2019). MTCH2 Promotes Adipogenesis in Intramuscular Preadipocytes via an M 6 A‐YTHDF1‐dependent Mechanism. FASEB j. 33 (2), 2971–2981. doi:10.1096/fj.201801393RRR
Kim, D., Langmead, B., and Salzberg, S. L. (2015). HISAT: a Fast Spliced Aligner with Low Memory Requirements. Nat. Methods 12 (4), 357–360. doi:10.1038/nmeth.3317
Kudou, K., Komatsu, T., Nogami, J., Maehara, K., Harada, A., Saeki, H., et al. (2017). The Requirement of Mettl3-Promoted MyoD mRNA Maintenance in Proliferative Myoblasts for Skeletal Muscle Differentiation. Open Biol. 7 (9), 170119. doi:10.1098/rsob.170119
Li, J., Pei, Y., Zhou, R., Tang, Z., and Yang, Y. (2021). Regulation of RNA N6-Methyladenosine Modification and its Emerging Roles in Skeletal Muscle Development. Int. J. Biol. Sci. 17 (7), 1682–1692. doi:10.7150/ijbs.56251
Liao, Y., Smyth, G. K., and Shi, W. (2014). featureCounts: an Efficient General Purpose Program for Assigning Sequence Reads to Genomic Features. Bioinformatics 30 (7), 923–930. doi:10.1093/bioinformatics/btt656
Lin, S., Choe, J., Du, P., Triboulet, R., and Gregory, R. I. (2016). The M(6)A Methyltransferase METTL3 Promotes Translation in Human Cancer Cells. Mol. Cel. 62 (3), 335–345. doi:10.1016/j.molcel.2016.03.021
Lin, J., Zhu, Q., Huang, J., Cai, R., and Kuang, Y. (2020). Hypoxia Promotes Vascular Smooth Muscle Cell (VSMC) Differentiation of Adipose-Derived Stem Cell (ADSC) by Regulating Mettl3 and Paracrine Factors. Stem Cell Int. 2020, 1–11. doi:10.1155/2020/2830565
Liu, J., Yue, Y., Han, D., Wang, X., Fu, Y., Zhang, L., et al. (2014). A METTL3-METTL14 Complex Mediates Mammalian Nuclear RNA N6-Adenosine Methylation. Nat. Chem. Biol. 10 (2), 93–95. doi:10.1038/nchembio.1432
Livak, K. J., and Schmittgen, T. D. (2001). Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2−ΔΔCT Method. Methods 25 (4), 402–408. doi:10.1006/meth.2001.1262
Luo, G.-Z., MacQueen, A., Zheng, G., Duan, H., Dore, L. C., Lu, Z., et al. (2014). Unique Features of the m6A Methylome in Arabidopsis thaliana. Nat. Commun. 5, 5630. doi:10.1038/ncomms6630
Mathiyalagan, P., Adamiak, M., Mayourian, J., Sassi, Y., Liang, Y., Agarwal, N., et al. (2019). FTO-dependent N 6 -Methyladenosine Regulates Cardiac Function during Remodeling and Repair. Circulation 139 (4), 518–532. doi:10.1161/CIRCULATIONAHA.118.033794
McKinnell, I. W., Ishibashi, J., Le Grand, F., Punch, V. G. J., Addicks, G. C., Greenblatt, J. F., et al. (2008). Pax7 Activates Myogenic Genes by Recruitment of a Histone Methyltransferase Complex. Nat. Cel. Biol. 10 (1), 77–84. doi:10.1038/ncb1671
Meng, J., Lu, Z., Liu, H., Zhang, L., Zhang, S., Chen, Y., et al. (2014). A Protocol for RNA Methylation Differential Analysis with MeRIP-Seq Data and exomePeak R/Bioconductor Package. Methods 69 (3), 274–281. doi:10.1016/j.ymeth.2014.06.008
Meyer, K. D., and Jaffrey, S. R. (2014). The Dynamic Epitranscriptome: N6-Methyladenosine and Gene Expression Control. Nat. Rev. Mol. Cel. Biol. 15 (5), 313–326. doi:10.1038/nrm3785
Meyer, K. D., Saletore, Y., Zumbo, P., Elemento, O., Mason, C. E., and Jaffrey, S. R. (2012). Comprehensive Analysis of mRNA Methylation Reveals Enrichment in 3′ UTRs and Near Stop Codons. Cell 149 (7), 1635–1646. doi:10.1016/j.cell.2012.05.003
Meyer, K. D., Patil, D. P., Zhou, J., Zinoviev, A., Skabkin, M. A., Elemento, O., et al. (2015). 5′ UTR m6A Promotes Cap-Independent Translation. Cell 163 (4), 999–1010. doi:10.1016/j.cell.2015.10.012
Millay, D. P., O’Rourke, J. R., Sutherland, L. B., Bezprozvannaya, S., Shelton, J. M., Bassel-Duby, R., et al. (2013). Myomaker Is a Membrane Activator of Myoblast Fusion and Muscle Formation. Nature 499 (7458), 301–305. doi:10.1038/nature12343
Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T.-C., Mendell, J. T., and Salzberg, S. L. (2015). StringTie Enables Improved Reconstruction of a Transcriptome from RNA-Seq Reads. Nat. Biotechnol. 33 (3), 290–295. doi:10.1038/nbt.3122
Picard, B., Berri, C., Lefaucheur, L., Molette, C., Sayd, T., and Terlouw, C. (2010). Skeletal Muscle Proteomics in Livestock Production. Brief. Funct. Genomics 9 (3), 259–278. doi:10.1093/bfgp/elq005
Ping, X.-L., Sun, B.-F., Wang, L., Xiao, W., Yang, X., Wang, W.-J., et al. (2014). Mammalian WTAP Is a Regulatory Subunit of the RNA N6-Methyladenosine Methyltransferase. Cell Res. 24 (2), 177–189. doi:10.1038/cr.2014.3
Quinlan, A. R., and Hall, I. M. (2010). BEDTools: a Flexible Suite of Utilities for Comparing Genomic Features. Bioinformatics 26 (6), 841–842. doi:10.1093/bioinformatics/btq033
Ramírez, F., Ryan, D. P., Grüning, B., Bhardwaj, V., Kilpert, F., Richter, A. S., et al. (2016). deepTools2: a Next Generation Web Server for Deep-Sequencing Data Analysis. Nucleic Acids Res. 44 (W1), W160–W165. doi:10.1093/nar/gkw257
Raudvere, U., Kolberg, L., Kuzmin, I., Arak, T., Adler, P., Peterson, H., et al. (2019). g:Profiler: a Web Server for Functional Enrichment Analysis and Conversions of Gene Lists (2019 Update). Nucleic Acids Res. 47 (W1), W191–W198. doi:10.1093/nar/gkz369
Salih, D. A. M., Tripathi, G., Holding, C., Szestak, T. A. M., Gonzalez, M. I., Carter, E. J., et al. (2004). Insulin-like Growth Factor-Binding Protein 5 (Igfbp5) Compromises Survival, Growth, Muscle Development, and Fertility in Mice. Proc. Natl. Acad. Sci. 101 (12), 4314–4319. doi:10.1073/pnas.0400230101
Schibler, U., Kelley, D. E., and Perry, R. P. (1977). Comparison of Methylated Sequences in Messenger RNA and Heterogeneous Nuclear RNA from Mouse L Cells. J. Mol. Biol. 115 (4), 695–714. doi:10.1016/0022-2836(77)90110-3
Schmidt, K., Glaser, G., Wernig, A., Wegner, M., and Rosorius, O. (2003). Sox8 Is a Specific Marker for Muscle Satellite Cells and Inhibits Myogenesis. J. Biol. Chem. 278 (32), 29769–29775. doi:10.1074/jbc.M301539200
Schwartz, S., Agarwala, S. D., Mumbach, M. R., Jovanovic, M., Mertins, P., Shishkin, A., et al. (2013). High-resolution Mapping Reveals a Conserved, Widespread, Dynamic mRNA Methylation Program in Yeast Meiosis. Cell 155 (6), 1409–1421. doi:10.1016/j.cell.2013.10.047
Shen, L., Liang, Z., Gu, X., Chen, Y., Teo, Z. W. N., Hou, X., et al. (2016). N 6 -Methyladenosine RNA Modification Regulates Shoot Stem Cell Fate in Arabidopsis. Dev. Cel. 38 (2), 186–200. doi:10.1016/j.devcel.2016.06.008
Shi, H., Wei, J., and He, C. (2019). Where, when, and How: Context-dependent Functions of RNA Methylation Writers, Readers, and Erasers. Mol. Cel. 74 (4), 640–650. doi:10.1016/j.molcel.2019.04.025
Song, T., Yang, Y., Wei, H., Xie, X., Lu, J., Zeng, Q., et al. (2019). Zfp217 Mediates m6A mRNA Methylation to Orchestrate Transcriptional and post-transcriptional Regulation to Promote Adipogenic Differentiation. Nucleic Acids Res. 47 (12), 6130–6144. doi:10.1093/nar/gkz312
Spicer, D. B., Rhee, J., Cheung, W. L., and Lassar, A. B. (1996). Inhibition of Myogenic bHLH and MEF2 Transcription Factors by the bHLH Protein Twist. Science 272 (5267), 1476–1480. doi:10.1126/science.272.5267.1476
Stark, R., and Brown, G. (2011). DiffBind: Differential Binding Analysis of ChIP-Seq Peak Data. [Online]. Available at http://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf [Accessed June 18, 2021].
Takada, F., Vander Woude, D. L., Tong, H. Q., Thompson, T. G., Watkins, S. C., Kunkel, L. M., et al. (2001). Myozenin: An Alpha -actinin- and Gamma -Filamin-Binding Protein of Skeletal Muscle Z Lines. Proc. Natl. Acad. Sci. U S A. 98 (4), 1595–1600. doi:10.1073/pnas.041609698
Tao, X., Chen, J., Jiang, Y., Wei, Y., Chen, Y., Xu, H., et al. (2017). Transcriptome-wide N 6 -methyladenosine Methylome Profiling of Porcine Muscle and Adipose Tissues Reveals a Potential Mechanism for Transcriptional Regulation and Differential Methylation Pattern. BMC Genomics 18 (1), 336. doi:10.1186/s12864-017-3719-1
Tarasov, A., Vilella, A. J., Cuppen, E., Nijman, I. J., and Prins, P. (2015). Sambamba: Fast Processing of NGS Alignment Formats. Bioinformatics 31 (12), 2032–2034. doi:10.1093/bioinformatics/btv098
Wang, X., Lu, Z., Gomez, A., Hon, G. C., Yue, Y., Han, D., et al. (2014). N6-Methyladenosine-Dependent Regulation of Messenger RNA Stability. Nature 505 (7481), 117–120. doi:10.1038/nature12730
Wang, X., Zhao, B. S., Roundtree, I. A., Lu, Z., Han, D., Ma, H., et al. (2015). N6-methyladenosine Modulates Messenger RNA Translation Efficiency. Cell 161 (6), 1388–1399. doi:10.1016/j.cell.2015.05.014
Wang, X., Huang, N., Yang, M., Wei, D., Tai, H., Han, X., et al. (2017). FTO Is Required for Myogenesis by Positively Regulating mTOR-PGC-1α Pathway-Mediated Mitochondria Biogenesis. Cell Death Dis. 8 (3), e2702. doi:10.1038/cddis.2017.122
Wang, Y.-N., Yang, W.-C., Li, P.-W., Wang, H.-B., Zhang, Y.-Y., and Zan, L.-S. (2018). Myocyte Enhancer Factor 2A Promotes Proliferation and its Inhibition Attenuates Myogenic Differentiation via Myozenin 2 in Bovine Skeletal Muscle Myoblast. PLoS One 13 (4), e0196255. doi:10.1371/journal.pone.0196255
Wang, Y., Gao, M., Zhu, F., Li, X., Yang, Y., Yan, Q., et al. (2020). METTL3 Is Essential for Postnatal Development of Brown Adipose Tissue and Energy Expenditure in Mice. Nat. Commun. 11 (1), 1648. doi:10.1038/s41467-020-15488-2
Wei, C.-M., Gershowitz, A., and Moss, B. (1975). Methylated Nucleotides Block 5′ Terminus of HeLa Cell Messenger RNA. Cell 4 (4), 379–386. doi:10.1016/0092-8674(75)90158-0
Wei, C. M., Gershowitz, A., and Moss, B. (1976). 5'-Terminal and Internal Methylated Nucleotide Sequences in HeLa Cell mRNA. Biochemistry 15 (2), 397–401. doi:10.1021/bi00647a024
Wei, L.-H., Song, P., Wang, Y., Lu, Z., Tang, Q., Yu, Q., et al. (2018). The m6A Reader ECT2 Controls Trichome Morphology by Affecting mRNA Stability in Arabidopsis. Plant Cel. 30 (5), 968–985. doi:10.1105/tpc.17.00934
Weintraub, H., Davis, R., Tapscott, S., Thayer, M., Krause, M., Benezra, R., et al. (1991). The myoD Gene Family: Nodal point during Specification of the Muscle Cell Lineage. Science 251 (4995), 761–766. doi:10.1126/science.1846704
Weintraub, H. (1993). The MyoD Family and Myogenesis: Redundancy, Networks, and Thresholds. Cell 75 (7), 1241–1244. doi:10.1016/0092-8674(93)90610-3
Xu, T., Xu, Z., Lu, L., Zeng, T., Gu, L., Huang, Y., et al. (2021). Transcriptome-wide Study Revealed m6A Regulation of Embryonic Muscle Development in Dingan Goose (Anser cygnoides Orientalis). BMC Genomics 22 (1), 270. doi:10.1186/s12864-021-07556-8
Yang, D., Qiao, J., Wang, G., Lan, Y., Li, G., Guo, X., et al. (2018). N 6-Methyladenosine Modification of lincRNA 1281 Is Critically Required for mESC Differentiation Potential. Nucleic Acids Res. 46 (8), 3906–3920. doi:10.1093/nar/gky130
Yang, Y., Fan, X., Yan, J., Chen, M., Zhu, M., Tang, Y., et al. (2021). A Comprehensive Epigenome Atlas Reveals DNA Methylation Regulating Skeletal Muscle Development. Nucleic Acids Res. 49 (3), 1313–1329. doi:10.1093/nar/gkaa1203
Yao, Y., Bi, Z., Wu, R., Zhao, Y., Liu, Y., Liu, Q., et al. (2019). METTL3 Inhibits BMSC Adipogenic Differentiation by Targeting the JAK1/STAT5/C/EBPβ Pathway via an M 6 A‐YTHDF2-dependent Manner. FASEB J. 33 (6), 7529–7544. doi:10.1096/fj.201802644R
Ye, J., Wang, Z., Chen, X., Jiang, X., Dong, Z., Hu, S., et al. (2020). YTHDF1-enhanced Iron Metabolism Depends on TFRC m6A Methylation. Theranostics 10 (26), 12072–12089. doi:10.7150/thno.51231
Yu, G., Wang, L.-G., and He, Q.-Y. (2015). ChIPseeker: an R/Bioconductor Package for ChIP Peak Annotation, Comparison and Visualization. Bioinformatics 31 (14), 2382–2383. doi:10.1093/bioinformatics/btv145
Zhang, Y., Liu, T., Meyer, C. A., Eeckhoute, J., Johnson, D. S., Bernstein, B. E., et al. (2008). Model-based Analysis of ChIP-Seq (MACS). Genome Biol. 9 (9), R137. doi:10.1186/gb-2008-9-9-r137
Zhang, Z., Theler, D., Kaminska, K. H., Hiller, M., de la Grange, P., Pudimat, R., et al. (2010). The YTH Domain Is a Novel RNA Binding Domain. J. Biol. Chem. 285 (19), 14701–14710. doi:10.1074/jbc.M110.104711
Zhang, X., Yao, Y., Han, J., Yang, Y., Chen, Y., Tang, Z., et al. (2020). Longitudinal Epitranscriptome Profiling Reveals the Crucial Role of N6-Methyladenosine Methylation in Porcine Prenatal Skeletal Muscle Development. J. Genet. Genomics 47 (8), 466–476. doi:10.1016/j.jgg.2020.07.003
Zhao, X., Yang, Y., Sun, B.-F., Shi, Y., Yang, X., Xiao, W., et al. (2014). FTO-dependent Demethylation of N6-Methyladenosine Regulates mRNA Splicing and Is Required for Adipogenesis. Cel Res. 24 (12), 1403–1419. doi:10.1038/cr.2014.151
Zhao, B. S., Roundtree, I. A., and He, C. (2017). Post-transcriptional Gene Regulation by mRNA Modifications. Nat. Rev. Mol. Cel. Biol. 18 (1), 31–42. doi:10.1038/nrm.2016.132
Zheng, G., Dahl, J. A., Niu, Y., Fedorcsak, P., Huang, C.-M., Li, C. J., et al. (2013). ALKBH5 Is a Mammalian RNA Demethylase that Impacts RNA Metabolism and Mouse Fertility. Mol. Cel. 49 (1), 18–29. doi:10.1016/j.molcel.2012.10.015
Zhong, S., Li, H., Bodi, Z., Button, J., Vespa, L., Herzog, M., et al. (2008). MTA Is an Arabidopsis Messenger RNA Adenosine Methylase and Interacts with a Homolog of a Sex-Specific Splicing Factor. Plant Cel. 20 (5), 1278–1288. doi:10.1105/tpc.108.058883
Keywords: N6-methyladenosine, cattle, myoblast differentiation, m6A-seq, differentially methylated genes
Citation: Yang X, Wang J, Ma X, Du J, Mei C and Zan L (2021) Transcriptome-wide N6-Methyladenosine Methylome Profiling Reveals m6A Regulation of Skeletal Myoblast Differentiation in Cattle (Bos taurus). Front. Cell Dev. Biol. 9:785380. doi: 10.3389/fcell.2021.785380
Received: 29 September 2021; Accepted: 12 November 2021;
Published: 06 December 2021.
Edited by:
Alaguraj Veluchamy, St. Jude Children’s Research Hospital, United StatesReviewed by:
Jiangbo Wei, University of Chicago, United StatesGeorgia Dermentzaki, Columbia University, United States
Copyright © 2021 Yang, Wang, Ma, Du, Mei and Zan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Linsen Zan, zanlinsen@163.com