Transcriptome-wide analysis of RNA m6A methylation regulation of muscle development in Queshan Black pigs
BMC Genomics volume 24, Article number: 239 (2023)
N6-methyladenosine (m6A) refers to the methylation modification of N6 position of RNA adenine, a dynamic reversible RNA epigenetic modification that plays an important regulatory role in a variety of life processes. In this study, we used MeRIP-Seq and RNA-Seq of the longissimus dorsi (LD) muscle of adult (QA) and newborn (QN) Queshan Black pigs to screen key genes with m6A modification involved in muscle growth by bioinformatics analysis.
A total of 23,445 and 25,465 m6A peaks were found in the whole genomes of QA and QN, respectively. Among them, 613 methylation peaks were significantly different (DMPs) and 579 genes were defined as differentially methylated genes (DMGs). Compared with the QN group, there were 1,874 significantly differentially expressed genes (DEGs) in QA group, including 620 up-regulated and 1,254 down-regulated genes. In order to investigate the relationship between m6A and mRNA expression in the muscle of Queshan Black pigs at different periods, a combined analysis of MeRIP-Seq and RNA-Seq showed that 88 genes were significantly different at both levels. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes results showed that DEGs and DMGs were mainly involved in skeletal muscle tissue development, FoxO signaling pathway, MAPK signaling pathway, insulin signaling pathway, PI3K–Akt signaling pathway, and Wnt signaling pathway. Four DEGs (IGF1R, CCND2, MYOD1 and FOS) and four DMGs (CCND2, PHKB, BIN1 and FUT2), which are closely related to skeletal muscle development, were selected as candidate genes for verification, and the results were consistent with the sequencing results, which indicated the reliability of the sequencing results.
These results lay the foundation for understanding the specific regulatory mechanisms of growth in Queshan Black pigs, and provide theoretical references for further research on the role of m6A in muscle development and breed optimization selection.
With the advancement of epigenetic studies, chemical RNA modifications have received increasing attention from researchers. To date, more than 150 types of RNA post-transcriptional modifications have been identified in all organisms . The common modifications are N6-methyladenosine (m6A), N1-methyladenosine (m1A), 5-methylcytosine (m5C), and pseudouridine (ψ) . In the 1970s, the presence of m6A modification was first identified in hepatocellular carcinoma cells by Desrosiers et al. . m6A refers to the methylation modification of N6 position of RNA adenine, which occurs in the highly conserved consensus sequence RRACH (R = G or A; H = A, C, or U)  and is most enriched near the 3′ untranslated regions (UTRs), followed by the coding sequence (CDS) region and 5′ UTRs . It is one of the most abundant and highly conserved forms of post-transcriptional modifications in eukaryotes and is widely present in mRNAs, tRNAs, rRNAs, mRNAs, and lncRNAs ; and plays a role in post-transcriptional regulation, affecting mRNA splicing, expression, translation, and stability [7, 8]. In addition, it is involved in a variety of biological processes, such as energy metabolism and adipogenesis .
To date, studies on m6A methylation modifications have focused on humans, mice, and other model animals; however, little is known about its specific mechanisms of action. Similar to DNA and histone methylation classes, m6A methylation is dynamically reversible in mammals  and relies on the co-regulation of multiple proteins, including m6A methyltransferases (writers), demethylases (erasers), and m6A methylated reading proteins (readers) [11, 12]. m6A methyltransferases, such as METTL3, METTL14, and WTAP are the core proteins of methyltransferases that form the m6A methyltransferase complex and play a catalytic role . m6A-methylated methyl is derived from S-adenosylmethionine (SAM); METTL3 binds SAM and transfers the methyl to a specific site in the mRNA . METTL3 exerts catalytic activity, whereas METTL14 does not have catalytic activity but can facilitate the binding of METTL3 to substrates. WTAP does not have methylation activity but can interact with the METTL3–METTL14 complex and participate in the m6A modification process together. Only two types of demethylases have been identified, namely, FTO and AlkB homolog 5 (ALKBH5). FTO is the earliest identified demethylase. ALKBH5 can remove methyl directly from m6A-modified adenine without going through the oxidation process and can effectively remove m6A modification from mRNA . m6A-binding proteins mostly carry conserved YTH structural domains, such as the YTH homologous structural domain protein family (YTHDF1, YTHDF2, and YTHDF3) and YTH structural domain proteins (YTHDC1, YTHDC2). YTHDF1 is able to bind to m6A modification sites to improve mRNA translation efficiency, and YTHDF2 promotes mRNA degradation and dynamically regulates mRNA abundance . YTHDF3 can further enhance mRNA translation or degradation by binding to YTHDF1 or YTHDF2 . YTHDC1 is an intranuclear binding protein that is capable of participating in mRNA processing and nuclear localization. YTHDC2 possesses a specific helix and protein repeat structural domain where binding to the m6A modification site is achieved, which in turn promotes mRNA translation . In addition to proteins containing the YTH structural domain, an increasing number of binding proteins play important functions in the m6A modification process.
Queshan Black pig is a local excellent pig breed with high fertility, adaptability, excellent meat quality, and stable genetic performance in Henan Province, China. Our laboratory has conducted long-term research around this breed to understand its specific physiological characteristics, growth, and development mechanism [18, 19]. Based on the indispensable role played m6A methylation modifications in regulating gene expression and participating in various biological processes, we hypothesized that m6A modifications are involved in the muscle growth and development of Queshan Black pigs. Understanding the molecular mechanisms of muscle growth is essential for maintaining meat yield and quality. Therefore, this study sequenced QN and QA LD samples by MeRIP-Seq and RNA-Seq to identify differential methylation peaks (DMPs) and differentially expressed genes (DEGs), so as to further explore their function and mechanism. This study may provide a theoretical basis for further research on the specific regulatory mechanisms of growth of Queshan Black pigs and the optimal selection of this breed.
Materials and methods
All of the experiments involving animals were carried out in accordance with the guidelines for the care and use of experimental animals established by the Ministry of Science and Technology of the People’s Republic of China (Approval Number DWLL20211193). The animal study was reviewed and approved by the Ethics Committee of Henan Agricultural University. In addition, all experiments were conducted in accordance with the relevant approved guidelines and regulations during slaughter, sampling, and sample conservation.
Animals and tissue collection
Three QA and three QN pigs, all male, were selected in this study. The Queshan Black pigs used in this experiment were all from Queshan Black pig breeding farm in Henan Province. Each group of pigs were fed under the same conditions. The piggery type is completely enclosed piggery, the breeding environment is suitable, the pigs are healthy, no genetic diseases. Piglets were slaughtered at 3 days of age with a weight of 2.2 − 2.5 kg. The body length of QN1 − QN3 were 34, 35 and 35 cm, respectively. Adult pigs were slaughtered at 270 days of age and weighing 100 − 102 kg. The body length of QA1 − QA3 were 103, 111 and 103 cm, respectively. The longissimus dorsal (LD) muscle between the 6th and 7th ribs was collected and immediately stored in liquid nitrogen. Afterward, the muscle samples were placed at − 80 °C for subsequent experiments.
RNA isolation, library preparation and sequencing
Total RNA was isolated and purified using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's procedure. The total RNA quality and quantity were determined by Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, CA, USA) with RIN number > 7.0. Poly(A) RNA was specifically captured from 50 μg total RNA using Oligo-dT magnetic beads and fragmented using Magnesium RNA Fragmentation Module (NEB, cat.e6150, USA). The cleaved RNA fragments were incubated for 2 h at 4 °C with m6A-specific antibody (No. 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 mg/ ml). The eluted RNA was precipitated with 75% ethanol. According to the chain-specific library prepared by the dUTP method, the eluted m6A fragment (IP) and the unprocessed input control fragment were converted into the final cDNA library. The average insert size of the paired-end libraries was 100 ± 50 bp. Lastly, we performed paired-end sequencing (PE150) on an Illumina Novaseq™ 6000 platform (LC-Bio Technology Co., Ltd., Hangzhou, China) following the vendor's recommended protocol.
Bioinformatics analysis process
Fastp software (v0.19.4, https://github.com/OpenGene/fastp) was used to remove the reads with adaptor contamination, low-quality bases, and undetermined bases with default parameter . The sequence quality of the IP and input samples were verified using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and RseQC (http://rseqc.sourceforge.net/) [21, 22]. Then, HISAT2 (v2.0.4, http://daehwankimlab.github.io/hisat2) was used to map reads to the Sus scrofa (version 11.1) . Peak calling and diff peak analysis were performed by R package exomePeak2 (v1.5.0, https://bioconductor.org/packages/release/bioc/html/exomePeak2.html). The threshold settings of differential peak and differential expression were |log2FC|≥ 1 and p-adjust < 0.05. and peaks were annotated by intersection with gene architecture using R package ANNOVAR (http://www.openbioinformatics.org/annovar/) [24, 25]. MEME (v5.3.3, http://meme-suite.org) and HOMER (v4.10, http://homer.ucsd.edu/homer/motif) were used for de novo and known motif finding, followed by motif localization with respect to the peak summit [26, 27]. StringTie (v2.1.2, https://ccb.jhu.edu/software/stringtie) was used to determine the expression levels of all transcripts and genes from input libraries by calculating the FPKM (total exon fragments / mapped reads [millions] × exon length [kB]) . The differentially expressed transcripts and DEGs were selected with |log2FC|≥ 1 and p-adjust < 0.05 by R package edgeR (v4.1, https://bioconductor.org/packages/edgeR) . The gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichments of DEGs were performed by KOBAS (v3.0) online analyses. P < 0.05 was considered statistically significant. All KEGG pathways shown in the manuscript can be found on the official KEGG website [30-32]. Protein − protein interaction (PPI) analysis was performed for DEGs and DMGs. The selected genes were imported into STRING (v11.5) for online analysis, and corresponding data were imported into Cytoscape (v.3.9.1) for visualization. Cytoscape (v.3.9.1) was also used to visualize the network of pathways and genes. Gene set enrichment analysis (GSEA) was performed using the OmicStudio tools at https://www.omicstudio.cn/tool. The pathway with |NES|> 1, NOM p-val < 0.05, FDR q-val < 0.25 was considered statistically significant.
According to the manufacturer's instructions, total RNA was extracted from the LD muscle tissue, the RNA fragments were converted into ~ 200nt fragments using the riboMeRIP™ m6A Transcriptome Profiling Kit (10 Assay) (RN: R11096.6, RiboBio, Guangzhou, China), and anti-m6A magnetic beads were prepared. Part of the RNA samples were taken as Input, and the rest were used as IP group for immunoprecipitation. Protein A/G magnetic beads were added into IP group and incubated at 4℃. Then the MagenTM Hipure Serum/plasma miRNA kit (R4317-03, Magen, Guangzhou,China) was used for elution and RNA recovery, and the obtained RNA was used for subsequent reverse transcription and q-PCR verification.
Quantitative Real-Time PCR
We verified four DEGs and four DMGs to study the status of m6A in LD tissues of Queshan Black pigs at different periods. RNA extracted from muscle was reverse-transcribed into cDNA using the Evo M-MLV RT Kit with gDNA Clean for qPCR kit (AG11705, Accurate Biotechnology (Hunan) Co.,Ltd, Changsha, China). q-PCR was performed using the SYBR Green Premix Pro Taq HS qPCR Kit (AG11701, Accurate Biotechnology (Hunan) Co.,Ltd, Changsha, China) using the CFX96 real-time PCR detection system (Thermo Fisher Scientific, USA) according to the instructions. The GAPDH was used as the internal reference gene to normalize the gene expression levels. The relative expression of genes was calculated using the 2−ΔΔCt method. MeRIP-qPCR does not require internal reference genes, and the IP/input ratio is calculated by 2−ΔCt (ΔCt = CtIP − Ctinput) and the ratio of IP RNA template and input RNA template to initial RNA when reverse transcription is introduced. Detailed primers are shown in supplementary Table 1 (Table S1).
The data are expressed as mean ± standard deviation (n = 3). The student’s t-test and one-way analysis of variance (ANOVA) was performed using GraphPad Prism 8 software to determine the significance of differences between comparison groups. The difference between the means was considered statistically significant when p-value ≤ 0.05.
Sequence statistics, quality control and reference genome alignment
After using fastp to filter out unqualified sequences from the raw data, the clean data was used for MeRIP-Seq and RNA-Seq analysis. We obtained two sets of QA and QN muscle sample data readings with three biological replicates per set. For the valid data of each group of samples, base proportion with mass value ≥ 20 (sequencing error rate less than 0.01), base proportion with mass value ≥ 30 (sequencing error rate less than 0.001) and GC content proportion are shown in Table S2.
The IP samples from the longissimus dorsi (LD) muscles of QA and QN pigs are denoted as QA_IP and QN_IP, respectively, in the m6A-Seq library in Figure S1. The LD muscle samples from QA and QN pigs are denoted as QA_input and QN_input in the RNA-Seq library, respectively. Each set of samples was repeated three times. The unique mapped reads are shown in Table S3. According to the regional distribution information of the reference genome, which can be defined compared with exon, intron, and spacer regions, the percentage contents of sequenced sequences localized to exon regions should be the highest under normal conditions. The analysis results are shown in Figure S1.
Identification of m6A modification sites and analysis of differential methylation peaks
Based on MeRIP-seq (IP) and RNA-seq (Input) sequencing data, the position information and length information of peaks on the genome were obtained. Reads near TSS were abundant at the transcriptome start of genes, and the peak distribution is shown as a heat map in Fig. 1A. Next, with p-adjust < 0.05 and |log2FC|≥ 1 for threshold selection differential methylation peaks (DMPs) and differentially methylated genes (DMGs). A total of 613 DMPs were screened and 579 DMGs were annotated in the QA group compared with the QN group (Table S4). A total of 176 peaks showed increased expression (corresponding to 167 genes upregulated in m6A abundance), and 437 peaks showed decreased expression (corresponding to 418 genes downregulated in m6A abundance) as shown in Figs. 1B and C. The m6A peaks in QA and QN were enriched in the CDS near the stop codon (Fig. 1D). Through exomePeak2 analysis, 8,825 and 10,845 peaks were specifically expressed in the QA and QN groups, respectively, and 14,620 peaks were common between the two groups (Fig. 1E). The transcription products were divided into four regions, namely, 5′ UTR, 3′ UTR, exonic region, and intronic region. The distribution of m6A peaks in QA and QN was similar (Figs. 1F and G). The analysis of DMPs enrichment sites showed that 53.51% of the DMPs were enriched in the 3′ UTR, about 30.83% were in the exonic region, and 15.17% of the m6A modifications occurred in the 5′ UTR (Fig. 1H). By analyzing the distribution of m6A peaks for each mRNA or gene, we found that most of the mRNAs or genes had one m6A peak (mRNAs with upregulated peaks: 362/384, mRNAs with downregulated peaks: 714/758, genes with upregulated peaks: 360/383,; genes with downregulated peaks: 693/747. Figs. 1I and J). All variable m6A peaks were mapped to human chromosomes. The presence of m6A peaks was found in all chromosomes, especially chr1, chr3, and chr6 (Fig. 1K). In Table 1, the differential m6A peaks were concentrated in the 3′ UTR. The Table 1 showed the top 20 differential m6A peaks, where log2FC < 0 represents hypomethylation, and log2FC ≥ 0 represents hypermethylation.
RNA methylation and demethylation are initiated by the combined action of multiple binding proteins to the motifs of methylation sites. A motif is a biologically important nucleic acid sequence pattern that is highly conserved. We performed motif prediction for the two groups of samples and ranked them according to p-value. The smaller the p-value, the higher the ranking (Fig. 2). A common motif structure in RNA modification is RRACH (where R = A or G; H = A, C, or U).
Enrichment analysis of differentially methylated genes
GO and KEGG pathways of DMGs were performed to analyse the potential function of m6A-modified genes in the skeletal muscle growth and development of Queshan Black pigs. The enriched GO terms for the DMGs mainly included protein phosphorylation, regulation of glucose metabolic process, positive regulation of actin filament polymerization and regulation of cell cycle (Fig. 3A, Table S5). KEGG pathway enrichment analysis showed that DMGs were enriched to the PI3K–Akt signaling pathway, Wnt signaling pathway, p53 signaling pathway, Thyroid hormone signaling pathway, ECM-receptor interaction and other pathways related to myogenesis and development (Fig. 3B, Table S6). PPI analysis was performed for DMGs in GO terms and KEGG pathways shown in Figs. 3A and B. As shown in Fig. 3C, larger nodes indicate more connections. The network diagram of partial pathways and DMGs is shown in Fig. 3D. The results indicate that DMGs have a potential role in regulating gene expression and biological metabolism during the skeletal muscle growth and development of Queshan Black pigs.
Analysis of differentially expressed genes
RNA-Seq analysis was performed in all Input samples, and the overall distribution of DEGs could be understood through the volcano map in Fig. 4A. In this study, FPKM was used to measure the abundance of gene expression in different samples. As shown in Fig. 4B, this study found that among the 19,023 genes identified in the two groups of samples, 1,874 genes were differentially expressed in the QA group compared with the QN group, including 620 upregulated genes, 1254 downregulated genes (|log2FC|≥ 1 and p-adjust < 0.05, Table S7), and 17,149 genes without remarkable difference. As shown in Fig. 4C, the clustering patterns of genes between samples in two periods. Gene expression and expression density plots are shown in Figs. 4D and E, respectively. The top 20 DEGs are shown in Table 2.
GO and KEGG analyses were performed to further reveal the functions of DEGs. GO enrichment analysis of DEGs showed that including muscle organ development, fatty acid metabolic process, fat cell differentiation, muscle contraction, skeletal muscle tissue development and myotube differentiation, were significantly enriched (Fig. 5A, Table S8). KEGG analysis showed that DEGs were significantly enriched in such as PPAR signaling pathway, calcium signaling pathway, AMPK signaling pathway, Fatty acid degradation, Fatty acid metabolism, FoxO signaling pathway, MAPK signaling pathway, mTOR signaling pathway, Wnt signaling pathway, and other pathways related to fat deposition and muscle development regulation (Fig. 5B, Table S9). As shown in Fig. 5C, PPI analysis was performed for DEGs in GO terms and KEGG pathways shown in Figs. 5A and B. The network diagram of partial pathways and DEGs is shown in Fig. 5D.
A conjoint analysis of MeRIP-Seq and RNA-Seq data
In order to explore the potential regulatory effect of m6A modification on gene expression during skeletal muscle growth and development, the data of m6A-Seq and RNA-Seq were jointly analyzed in this study to further screen genes with substantial changes at the mRNA and m6A levels. The results revealed a negative correlation between methylation peak and gene expression level (Fig. 6A). This study found that 11,491 genes were modified by m6A in the QA group and 12,384 genes were modified by m6A in the QN group, among which 1,874 genes were significantly differentially expressed. Based on the results, 88 genes were screened out with remarkable alterations at both levels (Fig. 6B, Table S10). This result implies that m6A modifications may affect the expression of these genes during muscle growth and development. The overlapping results of DEGs and DMGs are shown in Fig. 6C, including 22 genes in common between "m6A_up" and "mRNA_up" (hyper-up) and 11 genes in common between "m6A_up" and "mRNA_down" (hyper-down). 29 genes (hypo-up) were common in "m6A_down" and "mRNA_up" and 27 genes (hypo-down) were common in "m6A_down" and "mRNA_down".
GO analysis showed that these genes were enriched in terms, such as long-chain fatty acid transport, thyroid hormone transport, response to muscle activity, adipose tissue development and regulation of I-kappaB kinase/NF-kappaB signaling (Fig. 6D, Table S11). KEGG pathway enrichment analysis showed that these genes were enriched in Calcium signaling pathway, FoxO signaling pathway, insulin signaling pathway, and PI3K–Akt signaling pathway and Thyroid hormone signaling pathway, which are remarkably associated with muscle development (Fig. 6E, Table S12). The network diagram of partial pathways and codifferential genes is shown in Fig. 6F. The differentially methylated sites in QA and QN showed altered intensity around the corresponding m6A peaks, according to Integrative Genomics Viewer (IGV) software (Fig. 6G).
Gene set enrichment analysis
The results of gene set enrichment analysis (GSEA) showed that it was consistent with the above KEGG enrichment results, skeletal muscle cell differentiation, Insulin signaling pathway and FoxO signaling pathway were also highly enriched (Fig. 7A-C). In addition, Fatty acid-related Fatty acid biosynthesis, Fatty acid degradation and fatty acid metabolism were also significantly enriched (Fig. 7D-F).
Validation by MeRIP-qPCR and qRT-PCR
As shown in Fig. 8A, qRT-PCR was used in this study to detect the expression levels of methylation-related enzymes in QA and QN. The expressions of methylated transferase METTL3, METTL14, demethylated enzyme FTO, ALKBH5 and methylated reading protein YTHDF2, YTHDF3 in QN were significantly higher than those in QA group. In order to verify the accuracy of the sequencing results, four DMGs and four DEGs that may have potential regulatory relationships with muscle growth and development were screened in this study (Table 3), among which CCND2 showed significant differences at both levels. The methylation level of DMGs and the gene expression level of DEGs were detected by MeRIP-qPCR and qRT-PCR respectively (Fig. 8B-D). The expression trends obtained by the experiment were consistent with the sequencing results. The consistency of the results confirmed the existence of m6A modification in the muscle of Queshan Black pig, which verified the reliability of the sequencing data.
With the rapid development of science and technology in recent years, RNA methylation research has become a cutting-edge research hotspot in the field of epigenetics based on DNA and protein modification research. RNA m6A methylation is a highly conserved epigenomic modification that is widely found in various eukaryotes, such as yeast, plants, drosophila, and mammals, but has been studied rarely in livestock production animals, such as pigs. RNA m6A modifications can play a variety of biological functions; it is involved in disease development, stem cell differentiation, mRNA metabolism, animal energy metabolism, fat deposition, muscle growth and development [33-37]. In 2017, Tao et al. constructed the first transcriptome-wide m6A methylation modification map of porcine muscle tissue with the help of methylation RNA IP and MeRIP-Seq techniques, revealing that m6A methylation modification sites are distributed in the CDS region, stop codon, and 3′ UTR region in the porcine transcriptome . It has been previously reported that m6A peaks are enriched near the stop codon in human transcripts, and many methylation sites are conserved in mouse and human transcripts . The methylation modification sites in yak are concentrated in the stop codon, CDS, TSS, 3′ UTR, and, to some extent, in the 5′ UTR . The results of the present study are consistent with this finding, as methylation modifications were located near the stop codon and in the 3′ UTR. Thus, these results suggest that the overall distribution of m6A modification sites is similar in the mammalian transcriptome, which again demonstrates that m6A modifications are conserved among species. In addition, our data showed that a large number of m6A methylation modifications existed in muscle tissue during the growth of the Queshan Black pig, which may have important effects on muscle fiber type, muscle cell maturation, muscle structural changes and further play an important regulatory role in the muscle growth and development of pigs by regulating gene expression levels.
RNA m6A modifications may play a key regulatory role in the differentiation and development of animal cells; for example, genes consistently modified by m6A are associated with myoblast growth and differentiation during three different developmental stages in yak . In the transcriptome of embryonic stem cells, m6A-modified genes are involved in the regulation of embryonic stem cell pluripotency . Therefore, in Queshan Black pig, two sequencing libraries, namely, m6A-Seq (IP) and RNA-Seq (Input), and the results were analyzed bioinformatically. By MeRIP-Seq, we detected a large number of m6A methylation peaks in the transcriptome of the muscle tissues of Queshan Black pigs, with 613 DMPs were detected. In this study, through GO, KEGG and GSEA, it is speculated that DEGs and DMGs have potential important functions in the regulation of skeletal muscle development and are involved in many important pathways, such as AMPK signaling pathway, FoxO signaling pathway, PI3K-Akt signaling pathway and Wnt signaling pathway. Among them, the AMPK signaling pathway is an important signaling pathway of energy metabolism that is heavily involved in skeletal muscle metabolic processes; controls skeletal muscle growth and development by regulating many downstream targets; is related to myocyte energy metabolism, protein synthesis, and catabolism; and plays an important role in regulating muscle mass and regeneration . The MAPK signaling pathway is a class of transcription factors that promote skeletal muscle development, and genes on this pathway are continuously expressed from myogenic cell development to myotube formation. The Wnt signaling pathway is an important pathway that regulates skeletal muscle development and differentiation. It is mainly reflected in myogenic regulatory factors (MRFs) gene family myogenic regulatory factor 5(Myf5) and myogenic determination gene (MyoD) in the early stage of myogenesis . The PI3K–Akt signaling pathway is involved in the regulation of important biological processes, such as muscle growth and development, metabolic regulation and homeostasis maintenance. The aberrant expression of genes related to the PI3K–Akt signaling pathway or the abnormal phosphorylation of proteins can be detected in diseases, such as antimyotrophic protein-deficient muscle and Duchenne muscular dystrophy . Members of the Fox family, especially the FOXO1 gene, are closely associated with muscle growth and development, as they play a key role in myogenic cell fusion and myofiber type conversion . PPAR is a member of the nuclear receptor transcription factor family involved in growth and development, metabolism, inflammation, and many cellular processes in different organs; three PPAR isoforms, namely, PPARα, -β/δ, and -γ, are expressed in skeletal muscle in different degrees, but all have a remarkable impact on muscle homeostasis, directly or indirectly .
Skeletal muscle development is a complex biological process that is regulated by multiple mechanisms. The study of the regulatory role of myogenic regulators and the regulation of skeletal muscle development through epigenetic modifications has provided a preliminary understanding of the regulatory network of skeletal muscle development . Based on previous studies on m6A methylation , m6A is hypothesized to be involved in the regulation of skeletal muscle development process. m6A-related modifying enzymes are able to regulate muscle development by regulating the expression of genes related to muscle cell differentiation. METTL3 is a key regulatory factor of skeletal muscle differentiation, and participates in the process of skeletal muscle differentiation of myogenic progenitor cells by regulating the mRNA expression level of myogenic transcription factors such as MYOD in myoblasts . METTL14 is involved in the regulation of skeletal muscle differentiation, and METTL14 knockdown leads to the decreased expression of MYOD and MYOG in C2C12 cells, inhibiting cell differentiation and promoting cell proliferation . FTO positively regulates the mTOR–PGC-1α pathway by affecting mTOR activity through m6A demethylase activity, participating in skeletal muscle differentiation . In this study, m6A methylation was analyzed in LD muscle tissue of Chinese native breed Queshan Black pig at both newborn and adult stages. The experimental results showed that the expressions of METTL3, METTL14, FTO, ALKBH5, YTHDF2 and YTHDF3 in QN were higher than those in QA, indicating that methylation may play a regulatory role in the early growth stage, thus mediating the growth and development process of muscle. Subsequently, four DEGs (IGF1R, CCND2, MYOD1 and FOS) and four DMGs (CCND2, PHKB, BIN1 and FUT2) were selected as candidate genes through a series of analysis, and were verified by experiments.
From the above analysis, it was found that some DEGs with methylation modification and DMGs were involved in the pathways related to muscle development, may play a regulatory role in the process of muscle growth and development. Some DEGs modified by m6A methylation were found to be involved in pathways related to muscle development, and they play a regulatory role in muscle growth and development. MYOD is a member of the MRF family; myogenic differentiation 1 (MYOD1) promotes muscle growth and development, enhances muscle cell metabolism, and plays an important role in lean muscle mass improvement and meat quality enhancement . FOS is enriched in the MAPK signaling pathway; is involved in cell growth and differentiation and muscle production. It can regulate skeletal muscle cell proliferation, differentiation, and transformation . Insulin-like growth factor 1 receptor (IGF1R) is an important component of the insulin-like growth factor (IGF) system that is widely expressed in animal muscles and is a key gene affecting growth and development . IGF1R can regulate the MYOD activation, promote muscle differentiation through Akt signaling, and induce MYOG expression to stimulate the terminal differentiation of myogenic cells . As a cyclin, Cyclin D2 (CCND2) is involved in a number of pathways related to muscle growth and development. It has been reported that CCND2 can significantly enhance the myogenic differentiation of muscle progenitor cells and is an effective regulator of muscle fiber generation . Phosphorylase b kinase (PHKB) affects the disorder of skeletal muscle glycogen metabolism and is involved in the process of glycogen decomposition, which leads to the occurrence of related myopathy . Lack of Fucosyltransferase 2 (FUT2) increases energy expenditure and heat production in brown adipose tissue . Bridging Integrator 1 (BIN1) is a key player in muscle development and its muscle-specific isoforms are required for skeletal muscle development and function at birth and muscle regeneration in adulthood .
From the above analysis, it can be seen that the genes screened in this study participate in several pathways related to muscle growth and development, and all exist in the LD muscle tissue of Queshan Black pigs. Meanwhile, the regulation level of m6A is negatively correlated with the transcription level in muscle, indicating that m6A modification not only participates in the process of muscle growth and development, but also may regulate gene expression level. These results laid a foundation for further exploration of the role of m6A modification in muscle growth and development.
This study was the first to discover the transcriptome-wide m6A methylation modification pattern affecting skeletal muscle development in Queshan Black pigs. The m6A map revealed the distribution characteristics of m6A modification in the transcriptome of Queshan Black pigs. A total of 1,874 DEGs and 613 DMPs were identified, including 176 up-regulated and 437 down-regulated peaks. Through bioinformatics analysis, four DEGs (IGF1R, CCND2, MYOD1 and FOS) and four DMGs (CCND2, PHKB, BIN1 and FUT2), which are closely related to skeletal muscle development, were selected as candidate genes for verification. The results of this study can lay a foundation for further determining the potential effect of m6A RNA modification on the regulation of muscle growth of Queshan Black pig, and provide a theoretical reference for the optimization of this breed.
Availability of data and materials
All raw data of high-throughput sequencing have been deposited to the National Genomics Data Center (NGDC, https://bigd.big.ac.cn) with the dataset accession number CRA009120.
Boccaletto P, Machnicka MA, Purta E, Piatkowski P, Baginski B, Wirecki TK, de Crécy-Lagard V, Ross R, Limbach PA, Kotter A, et al. MODOMICS: a database of RNA modification pathways. 2017 update. Nucleic Acids Res. 2018;46(D1):D303-d307.
Gilbert WV, Bell TA, Schaening C. Messenger RNA modifications: form, distribution, and function. Science. 2016;352(6292):1408–12.
Desrosiers R, Friderici K, Rottman F. Identification of methylated nucleosides in messenger RNA from Novikoff hepatoma cells. Proc Natl Acad Sci U S A. 1974;71(10):3971–5.
Zhang H, Shi X, Huang T, Zhao X, Chen W, Gu N, Zhang R. Dynamic landscape and evolution of m6A methylation in human. Nucleic Acids Res. 2020;48(11):6251–64.
Yue Y, Liu J, He C. RNA N6-methyladenosine methylation in post-transcriptional gene expression regulation. Genes Dev. 2015;29(13):1343–55.
Luo GZ, MacQueen A, Zheng G, Duan H, Dore LC, Lu Z, Liu J, Chen K, Jia G, Bergelson J, et al. Unique features of the m6A methylome in Arabidopsis thaliana. Nat Commun. 2014;5:5630.
Xiao W, Adhikari S, Dahal U, Chen YS, Hao YJ, Sun BF, Sun HY, Li A, Ping XL, Lai WY, et al. Nuclear m(6)A Reader YTHDC1 Regulates mRNA Splicing. Mol Cell. 2016;61(4):507–19.
Wang X, Zhao Boxuan S, Roundtree Ian A, Lu Z, Han D, Ma H, Weng X, Chen K, Shi H, He C. N6-methyladenosine modulates messenger RNA translation efficiency. Cell. 2015;161(6):1388–99.
Wang X, Zhu L, Chen J, Wang Y. mRNA m6A methylation downregulates adipogenesis in porcine adipocytes. Biochem Biophys Res Commun. 2015;459(2):201–7.
Wang Y, Zheng Y, Guo D, Zhang X, Guo S, Hui T, Yue C, Sun J, Guo S, Bai Z, et al. m6A methylation analysis of differentially expressed genes in skin tissues of coarse and fine type liaoning cashmere goats. Front Genet. 2019;10:1318.
Zhang C, Fu J, Zhou Y. A review in research progress concerning m6A methylation and immunoregulation. Front Immunol. 2019;10:922.
Meyer KD, Jaffrey SR. The dynamic epitranscriptome: N6-methyladenosine and gene expression control. Nat Rev Mol Cell Biol. 2014;15(5):313–26.
Liu J, Yue Y, Han D, Wang X, Fu Y, Zhang L, Jia G, Yu M, Lu Z, Deng X, et al. A METTL3-METTL14 complex mediates mammalian nuclear RNA N6-adenosine methylation. Nat Chem Biol. 2014;10(2):93–5.
Bokar JA, Shambaugh ME, Polayes D, Matera AG, Rottman FM. Purification and cDNA cloning of the AdoMet-binding subunit of the human mRNA (N6-adenosine)-methyltransferase. RNA. 1997;3(11):1233–47.
Zheng G, Dahl JA, Niu Y, Fedorcsak P, Huang CM, Li CJ, Vågbø CB, Shi Y, Wang WL, Song SH, et al. ALKBH5 is a mammalian RNA demethylase that impacts RNA metabolism and mouse fertility. Mol Cell. 2013;49(1):18–29.
Shi H, Wang X, Lu Z, Zhao BS, Ma H, Hsu PJ, Liu C, He C. YTHDF3 facilitates translation and decay of N(6)-methyladenosine-modified RNA. Cell Res. 2017;27(3):315–28.
Shi H, Zhang X, Weng YL, Lu Z, Liu Y, Lu Z, Li J, Hao P, Zhang Y, Zhang F, et al. m(6)A facilitates hippocampus-dependent learning and memory through YTHDF1. Nature. 2018;563(7730):249–53.
Li X, Qiao R, Ye J, Wang M, Zhang C, Lv G, Wang K, Li X, Han X. Integrated miRNA and mRNA transcriptomes of spleen profiles between Yorkshire and Queshan black pigs. Gene. 2019;688:204–14.
Qi K, Liu Y, Li C, Li X, Li X, Wang K, Qiao R, Han X. Construction of circRNA-related ceRNA networks in longissimus dorsi muscle of Queshan Black and Large White pigs. Mol Genet Genomics. 2022;297(1):101–12.
Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–90.
de Sena Brandine G, Smith AD. Falco: high-speed FastQC emulation for quality control of sequencing data. F1000Res. 1874;2019:8.
Wang L, Wang S, Li W. RSeQC: quality control of RNA-seq experiments. Bioinformatics. 2012;28(16):2184–5.
Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.
Meng J, Lu Z, Liu H, Zhang L, Zhang S, Chen Y, Rao MK, Huang Y. A protocol for RNA methylation differential analysis with MeRIP-Seq data and exomePeak R/Bioconductor package. Methods. 2014;69(3):274–81.
Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res. 2010;38(16): e164.
Bailey TL, Boden M, Buske FA, Frith M, Grant CE, Clementi L, Ren J, Li WW, Noble WS. MEME SUITE: tools for motif discovery and searching. Nucleic Acids Res. 2009;37(Web Server issue):W202-208.
Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, Cheng JX, Murre C, Singh H, Glass CK. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. 2010;38(4):576–89.
Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51.
Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2023;51(D1):D587-d592.
Zhang S, Zhao BS, Zhou A, Lin K, Zheng S, Lu Z, Chen Y, Sulman EP, Xie K, Bögler O, et al. m(6)A demethylase ALKBH5 maintains tumorigenicity of glioblastoma stem-like cells by sustaining FOXM1 expression and cell proliferation program. Cancer Cell. 2017;31(4):591-606.e596.
Boissel S, Reish O, Proulx K, Kawagoe-Takaki H, Sedgwick B, Yeo GS, Meyre D, Golzio C, Molinari F, Kadhom N, et al. Loss-of-function mutation in the dioxygenase-encoding FTO gene causes severe growth retardation and multiple malformations. Am J Hum Genet. 2009;85(1):106–11.
Tao X, Chen J, Jiang Y, Wei Y, Chen Y, Xu H, Zhu L, Tang G, Li M, Jiang A, et al. 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. 2017;18(1):336.
Batista PJ, Molinie B, Wang J, Qu K, Zhang J, Li L, Bouley DM, Lujan E, Haddad B, Daneshvar K, et al. m(6)A RNA modification controls cell fate transition in mammalian embryonic stem cells. Cell Stem Cell. 2014;15(6):707–19.
Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR. Comprehensive analysis of mRNA methylation reveals enrichment in 3’ UTRs and near stop codons. Cell. 2012;149(7):1635–46.
Ma X, La Y, Bao P, Chu M, Guo X, Wu X, Pei J, Ding X, Liang C, Yan P. Regulatory role of N6-methyladenosine in longissimus dorsi development in Yak. Front Vet Sci. 2022;9: 757115.
Thomson DM. The role of AMPK in the regulation of skeletal muscle size, hypertrophy, and regeneration. Int J Mol Sci. 2018;19(10):3125.
Abu-Elmagd M, Robson L, Sweetman D, Hadley J, Francis-West P, Münsterberg A. Wnt/Lef1 signaling acts via Pitx2 to regulate somite myogenesis. Dev Biol. 2010;337(2):211–9.
Alexander MS, Casar JC, Motohashi N, Myers JA, Eisenberg I, Gonzalez RT, Estrella EA, Kang PB, Kawahara G, Kunkel LM. Regulation of DMD pathology by an ankyrin-encoded miRNA. Skelet Muscle. 2011;1:27.
Sin TK, Yung BY, Siu PM. Modulation of SIRT1-Foxo1 signaling axis by resveratrol: implications in skeletal muscle aging and insulin resistance. Cell Physiol Biochem. 2015;35(2):541–52.
Manickam R, Duszka K, Wahli W. PPARs and microbiota in skeletal muscle health and wasting. Int J Mol Sci. 2020;21(21):8056.
Bharathy N, Ling BM, Taneja R. Epigenetic regulation of skeletal muscle development and differentiation. Subcell Biochem. 2013;61:139–50.
Kudou K, Komatsu T, Nogami J, Maehara K, Harada A, Saeki H, Oki E, Maehara Y, Ohkawa Y. The requirement of Mettl3-promoted MyoD mRNA maintenance in proliferative myoblasts for skeletal muscle differentiation. Open Biol. 2017;7(9):170119.
Zhang X, Yao Y, Han J, Yang Y, Chen Y, Tang Z, Gao F. Longitudinal epitranscriptome profiling reveals the crucial role of N(6)-methyladenosine methylation in porcine prenatal skeletal muscle development. J Genet Genomics. 2020;47(8):466–76.
Wang X, Huang N, Yang M, Wei D, Tai H, Han X, Gong H, Zhou J, Qin J, Wei X, et al. FTO is required for myogenesis by positively regulating mTOR-PGC-1α pathway-mediated mitochondria biogenesis. Cell Death Dis. 2017;8(3): e2702.
Lee EA, Kim JM, Lim KS, Ryu YC, Jeon WM, Hong KC. Effects of variation in porcine MYOD1 gene on muscle fiber characteristics, lean meat production, and meat quality traits. Meat Sci. 2012;92(1):36–43.
Almada AE, Horwitz N, Price FD, Gonzalez AE, Ko M, Bolukbasi OV, Messemer KA, Chen S, Sinha M, Rubin LL, et al. FOS licenses early events in stem cell activation driving skeletal muscle regeneration. Cell Rep. 2021;34(4): 108656.
O’Neill BT, Lauritzen HP, Hirshman MF, Smyth G, Goodyear LJ, Kahn CR. Differential Role of Insulin/IGF-1 Receptor Signaling in Muscle Growth and Glucose Homeostasis. Cell Rep. 2015;11(8):1220–35.
Sun L, Liu L, Yang XJ, Wu Z. Akt binds prohibitin 2 and relieves its repression of MyoD and muscle differentiation. J Cell Sci. 2004;117(Pt 14):3021–9.
Khanjyan MV, Yang J, Kayali R, Caldwell T, Bertoni C. A high-content, high-throughput siRNA screen identifies cyclin D2 as a potent regulator of muscle progenitor cell fusion and a target to enhance muscle regeneration. Hum Mol Genet. 2013;22(16):3283–95.
Tarnopolsky MA. Myopathies related to glycogen metabolism disorders. Neurotherapeutics. 2018;15(4):915–27.
Zhou R, Llorente C, Cao J, Zaramela LS, Zeng S, Gao B, Li SZ, Welch RD, Huang FQ, Qi LW, et al. Intestinal α1-2-fucosylation contributes to obesity and steatohepatitis in mice. Cell Mol Gastroenterol Hepatol. 2021;12(1):293–320.
Prokic I, Cowling BS, Kutchukian C, Kretz C, Tasfaout H, Gache V, Hergueux J, Wendling O, Ferry A, Toussaint A, et al. Differential physiological roles for BIN1 isoforms in skeletal muscle development, function and regeneration. Dis Model Mech. 2020;13(11):dmm044354.
This research was funded by the 14th Five-Year National Key R&D Program, grant number 2021YFD1301202; Agricultural Breeds Research Project of Henan Province, grant number 2022020101.
Ethics approval and consent to participate
All of the experiments involving animals were carried out in accordance with the guidelines for the care and use of experimental animals established by the Ministry of Science and Technology of the People’s Republic of China (Approval Number DWLL20211193). The animal study was reviewed and approved by the Ethics Committee of Henan Agricultural University. All experiments were performed in accordance with relevant guidelines and followed the ARRIVE guidelines for the reporting of animal experiments.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1. The primers used for the validation.
Table S2. Summary of reads quality control.
Table S3. Summary of reads mapped to the Sus scrofa.
Table S4. Differentially methylated peaks (DMPs) in QA vs. QN.
Table S5. All GO terms for the DMGs.
Table S6. All KEGG pathways for the DMGs.
Table S7. Differentially expressed genes (DEGs) in QA vs QN.
Table S8. All GO terms for the DMEGs.
Table S9. All KEGG pathways for the DEGs.
Table S10. Intersection of DMGs and DEGs.
Table S11. All GO terms for the intersection of DMGs and DEGs.
Table S12. All KEGG pathways for the intersection of DMGs and DEGs.
Figure S1. Refer to the genome to compare the regional distribution.
About this article
Cite this article
Dou, Y., Wei, Y., Zhang, Z. et al. Transcriptome-wide analysis of RNA m6A methylation regulation of muscle development in Queshan Black pigs. BMC Genomics 24, 239 (2023). https://doi.org/10.1186/s12864-023-09346-w