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.

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.

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

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 m⁶A 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.

The cleaved RNA fragments were incubated with m⁶A-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 m⁶A). The eluted RNA was precipitated with 75% ethanol. According to the strand-specific library prepared by dUTP method, the eluted fragments containing m⁶A (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 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.

m⁶A-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 m⁶A-modification peaks of each sample with the default parameters, which identifies m⁶A 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 m⁶A 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 m⁶A 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 m⁶A 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 m⁶A 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).

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.

m⁶A 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 m⁶A-specific antibody (202003, Synaptic Systems, Germany). The RNA fragments that bound to m⁶A were separated by TRIzol reagent. Following ethanol precipitation, the input RNA and eluted m⁶A 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 m⁶A peak region are listed in Supplementary Table S1.

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.

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.

To investigate the role of m⁶A in myoblast differentiation, mRNA was extracted from pre-differentiation (GM, myoblasts, D0) and post-differentiation (DM, myotubes, D4) cells for m⁶A-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 m⁶A-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 m⁶A-modified transcripts were identified in GM groups and 5,691–7,361 m⁶A-modified transcripts were found in DM groups. There were 9,246–13596 and 9,128–11861 m⁶A peaks in the two groups, respectively (Table 2, Supplementary Table S3). Therefore, we found that there were ∼1.60 m⁶A peaks per m⁶A 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 m⁶A transcript) (Jiang et al., 2019), pig adipose tissue (∼1.3 peaks per gene) (Tao et al., 2017) and chicken fat (∼1.5 peaks per m⁶A transcript) (Cheng et al., 2021). The results of the three biological replicates were similar, which also showed the accuracy of m⁶A-seq results and the similarity between duplicate samples. Summary plots and heatmaps were generated by deepTools using normalized read coverages from m⁶A-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 m⁶A consensus motif RRACH was highly enriched in GM and DM (Figure 2C), which was consistent with the general pattern of mammals.

To investigate the modification position of m⁶A in transcripts, we then studied the metagene profiles of the m⁶A peak in the entire transcriptome of GM and DM (Figure 2D). Three distinct m⁶A 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 m⁶A density of DM in 3′UTR is higher than GM, while the density in CDS is lower than GM. Further, to evaluate m⁶A enrichment systematically in mRNAs, we calculated the enrichment proportion of m⁶A peak in 5′UTR, CDS (coding sequence) and 3′UTR. More than 50% of m⁶A 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 m⁶A modification was highly similar in myoblasts and myotubes, and coincided with the typical m⁶A peak distribution in mammals. To further determine the distribution of m⁶A modifications in the transcriptome, we analyzed the number of m⁶A peaks contained in each gene. There was only one m⁶A 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 m⁶A modification at the chromosome level and the difference between GM and DM, the circle diagram of m⁶A peaks was produced using OmicStudio online tools (https://www.omicstudio.cn/tool) (Figure 2G). The circular visualization of peaks on chromosomes suggested that m⁶A methylation is extensively modified in myoblasts and myotubes, but there is an inconspicuous difference in m⁶A abundance between GM and DM. Cumulative curves also showed that there was no visible difference in m6A abundance between GM and DM (Figure 2H).

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 m⁶A abundance) and 610 peaks diminished (correspond to 571 genes with down-regulated m⁶A 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 m⁶A modification occurs in 5′UTR (Figure 3C).

Furthermore, to analyze the potential function of m⁶A 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.

m⁶A 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 m⁶A 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).

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

To demonstrate the potential regulation of m⁶A modification on gene expression in bovine myoblast differentiation, we analyzed genes with significant changes in both mRNA and m⁶A 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 m⁶A 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 “m⁶A_up” and “mRNA_up” (means hyper-up), 83 common genes in both “m⁶A_up” and “mRNA_down” (means hyper-down), 65 common genes in both “m⁶A_down” and “mRNA_up” (means hypo-up), and 46 common genes in both “m⁶A_down” and “mRNA_down” (means hypo-down) (Figure 5A, Supplementary Table S10). Intriguingly, m⁶A modification abundances of 15 genes exhibited both up-regulation and down-regulation, indicating that multiple peaks with m⁶A methylation modification within these genes had different significant changes in myoblast differentiation.

Likewise, GO and KEGG analyses were performed to predict the role of the genes with significant changes in both m⁶A 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 m⁶A 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 m⁶A 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 m⁶A-modified peaks of 24 of the 30 genes were enriched in 3′UTR and exons, which supported the results of the m⁶A-enriched region in this study and suggested that m⁶A may influence bovine myoblast differentiation by mediating the expression of these genes.

Finally, to determine whether m⁶A modification could regulate gene expression in bovine skeletal myoblasts, we analyzed the correlation between mRNA level and m⁶A peaks abundance in bovine myoblasts (GM) and myotubes (DM). Our results revealed that mRNA expression was strongly positively correlated with the abundance of m⁶A peaks in each of GM and DM (Figures 6A,B). However, there was no striking correlation between the change of m⁶A 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 m⁶A methylation modification. However, there was no significant positive or negative correlation between the changes of mRNA m⁶A methylation and differential expression during bovine skeletal myoblast differentiation.

In the m⁶A-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 m⁶A 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 m⁶A peaks of these four genes by using IGV software (Figure 7A). In all three biological replicates of GM and DM, m⁶A-modified peaks of MYOZ2 and TWIST1 were observed to be enriched near the start codon in mRNA, while m⁶A-modified peaks of KLF5 and MYOD1 were enriched within CDS and 3′UTR, respectively (Figure 7A). We performed m⁶A-IP-qPCR to verify the results of m⁶A-seq and indicated that the mRNAs of MYOZ2 and TWIST1 displayed higher levels of m⁶A enrichment in DM compared with GM, and the mRNAs of KLF5 and MYOD1 exhibited lower levels of m⁶A 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 m⁶A-seq and RNA-seq data. What’s more, these data provided possible mechanisms for m⁶A methylation to regulate myoblast differentiation, and these four genes may also be considered as downstream candidate targets for m⁶A modification during myoblast differentiation.