Coordinated transcriptional and post-transcriptional epigenetic regulation during skeletal muscle development and growth in pigs
Journal of Animal Science and Biotechnology volume 13, Article number: 146 (2022)
N6-methyladenosine (m6A) and DNA 5-methylcytosine (5mC) methylation plays crucial roles in diverse biological processes, including skeletal muscle development and growth. Recent studies unveiled a potential link between these two systems, implicating the potential mechanism of coordinated transcriptional and post-transcriptional regulation in porcine prenatal myogenesis and postnatal skeletal muscle growth.
Immunofluorescence and co-IP assays were carried out between the 5mC writers and m6A writers to investigate the molecular basis underneath. Large-scale in-house transcriptomic data were compiled for applying weighted correlation network analysis (WGCNA) to identify the co-expression patterns of m6A and 5mC regulators and their potential role in pig myogenesis. Whole-genome bisulfite sequencing (WGBS) and methylated RNA immunoprecipitation sequencing (MeRIP-seq) were performed on the skeletal muscle samples from Landrace pigs at four postnatal growth stages (days 30, 60, 120 and 180).
Significantly correlated expression between 5mC writers and m6A writers and co-occurrence of 5mC and m6A modification were revealed from public datasets of C2C12 myoblasts. The protein-protein interactions between the DNA methylase and the m6A methylase were observed in mouse myoblast cells. Further, by analyzing transcriptome data comprising 81 pig skeletal muscle samples across 27 developmental stages, we identified a 5mC/m6A epigenetic module eigengene and decoded its potential functions in pre- or post-transcriptional regulation in postnatal skeletal muscle development and growth of pigs. Following integrative multi-omics analyses on the WGBS methylome data and MeRIP-seq data for both m6A and gene expression profiles revealed a genome/transcriptome-wide correlated dynamics and co-occurrence of 5mC and m6A modifications as a consequence of 5mC/m6A crosstalk in the postnatal myogenesis progress of pigs. Last, we identified a group of myogenesis-related genes collaboratively regulated by both 5mC and m6A modifications in postnatal skeletal muscle growth in pigs.
Our study discloses a potential epigenetic mechanism in skeletal muscle development and provides a novel direction for animal breeding and drug development of related human muscle-related diseases.
The skeletal muscle is the largest organ of mammals and is crucial to the metabolism and energy balance of the whole body . Disruption of muscle fiber development may lead to human muscular diseases, e.g., muscle atrophy [2, 3] and a decrease in livestock meat quality and meat production . Pigs represent a major source of meat production worldwide and an ideal animal model in myogenesis studies [5, 6]. Deciphering the molecular mechanisms that regulate skeletal muscle development and growth is significant for both animal husbandry and biomedicine.
The spatiotemporal expression of gene functions in skeletal muscle development and growth is tightly regulated by dynamic epigenetic modifications [7, 8]. DNA 5-methylcytosine (5mC) and RNA N6-methyladenosine (m6A) are two typical epigenetic modifications that play essential roles in transcription and post-transcription regulation, respectively [9,10,11]. Recently, many studies have explored the epigenetic mechanisms underlying skeletal muscle growth and development in mammals, including chromatin landscape [12, 13], microRNA [14, 15], circRNA [16, 17], DNA and histone modifications [14, 18]. The discovery of the RNA methylation system in post-transcriptional regulation of myoblast differentiation and proliferation provides another layer of information for understanding skeletal muscle development and growth [19, 20]. Interestingly, one previous study reported significant positive correlations in genetic variation and gene expression of 5mC and m6A regulators, which include a group of proteins responsible for 5mC and m6A modification addition (writers), identification (readers) and removal (erasers) from DNA or RNA sequence, based on analyzing the gene expression profiles of 11,080 samples from 33 kinds of human cancers in the TCGA database . Another study also indicated for quantitative correlation between 5mC and m6A modifications in the genome of tomatoes . These studies implicated the potential mechanism of coordinated transcriptional and post-transcriptional regulation in various biological processes, which provides new possibilities for deciphering the mechanism of skeletal muscle development and growth.
Previously, we have characterized DNA methylome dynamics of 27 developmental stages in pig skeletal muscle development and revealed the involvement of the DNA methylation/SP1/IGF2BP3 axis in skeletal myogenesis . In addition, we also profiled the RNA methylation landscape during prenatal skeletal muscle development, using the same collection of samples . In the present study, we hypothesized that the transcriptional (5mC) and post-transcriptional (m6A) regulation of hub genes of porcine skeletal muscle development and growth are tightly coordinated via protein-protein interaction of 5mC and m6A methylases. To test this, we first demonstrated the genome/transcriptome-wide association between 5mC and m6A modification as well as their regulators’ expression, following with co-IP and immunofluorescence experiments to confirm the interaction between 5mC and m6A methylases, in the C2C12 myoblasts. Next, we evaluated the key 5mC/m6A regulators that played an important role in regulating prenatal and postnatal skeletal muscle development and identified a group of essential target genes whose transcription and translation are coordinately regulated by these regulators. Our study uncovered a novel epigenetic mechanism in skeletal muscle development and provided novel directions for investigating animal breeding, muscle biology, and related human diseases.
Materials and methods
For immunofluorescence experiments, C2C12 cells were cultured in DMEM supplemented with 10% FBS (Gibco, California, USA) in a humidified incubator with 5% CO2 at 37 °C. For co-IP assays, the cells were transferred to DMEM containing 2% horse serum (Gibco, California, USA) and induced to differentiate for 4 d.
C2C12 cells were fixed with 4% paraformaldehyde for 15 min at room temperature. After permeabilization with 0.25% Triton X-100 for 10 min, the cells were blocked by using 1% bovine serum albumin for 30 min. Then, the cells were incubated overnight with Anti-METTL3 antibody (67733–1-Ig; Protentech, Wuhan, China), and Anti-DNMT3A antibody (ab228691; Abcam, Cambridge, England),or Anti-DNMT3B antibody (ab227883; Abcam, Cambridge, England) diluted in PBST containing 1% bovine serum albumin at 4 °C. The cells were incubated with Alexa Fluor 647 donkey anti-mouse IgG (H + L), (A-31571; Invitrogen, MA, USA) or Alexa Fluor 488 donkey anti-rabbit IgG (H + L) (A-21206; Invitrogen, MA, USA) for 1 h at room temperature. Nuclei were stained with DAPI for 1 min. The cells observed microscopically using an A1HD25 confocal microscope (Nikon Instruments, Tokyo, Japan).
The differentiated C2C12 cells lyzed according to the instructions of the Co-Immunoprecipitation Kit (BersinBio, Guangzhou, China). Anti-Mettl3 antibody (ab195352), Anti-METTL14 antibody (ab252562), Anti-WTAP antibody (ab195380), Anti-DNMT1 antibody (ab19905), Anti-DNMT3A antibody (ab228691) and Anti-DNMT3B antibody (ab227883) were used for IP (Abcam, Cambridge, England). Nonspecific IgG antibodies (ab24584; Abcam, Cambridge, England) precipitated the complex as a control group. Total proteins from C2C12 cells were extracted and incubated overnight with IgG or the Co-IP antibodies (5 μg) at 4 °C, while the extracted proteins without any antibodies were used as Input control. Next, 30 μL of Protein A/G PLUS-Agarose (sc-2027; Santa Cruz Biotechnology, California, USA) was added and incubated for 2 h at 4 °C to form an immune complex. Following centrifugation for 4 min at 3000 r/min at 4 °C, the Protein A/G Plus-Agarose beads were washed 4 times with 1 mL of lysate and boiled in the appropriate protein loading buffer for 5 min. Centrifugation at 3000 r/min again, the supernatants were collected to a new tube for Western blot analysis.
The total proteins from C2C12 cells were lysed in RIPA lysis buffer supplemented with a protease inhibitor cocktail (4693124001; Roche, Mannheim, Germany). The membranes were blocked with 5% BSA for 1 h at room temperature and subsequently probed with primary antibodies overnight at 4 °C. the following dilutions were used for each antibody: METTL3, METTL14, WTAP, DNMT1, DNMT3A and DNMT3B (all 1:1000; Abcam, Cambridge, England). The membranes were then washed with PBS-Tween and incubated for 50 min with horseradish peroxidase-conjugated secondary antibodies (SA00001–1; Proteintech, Wuhan, China). Protein bands were detected after treatment with the SuperSignal west Femto agent (34094; Thermo Scientific, Waltham, USA).
Tissue sample collection
The skeletal muscle (longissimus dorsi, LD) samples were collected from Landrace pigs at four developmental stages, including postnatal days 30, 60, 120, and 180 (abbreviated as LD30, LD60, LD120, and LD180), as we previously reported . The experimental pigs were allowed access to feed and water ad libitum and were housed under identical conditions before slaughtering. After copulation with the boar, the sows and piglets were sacrificed at a commercial slaughter house at the selected stages. At each stage, skeletal muscles from two unrelated pigs were harvested as biological replicates. All samples were stored immediately in liquid nitrogen until further use. All animal procedures were performed according to protocols of the Chinese Academy of Agricultural Sciences and the Institutional Animal Care and Use Committee.
Total RNA was immediately isolated using Trizol reagent according to the manufacturer’s instructions (AM8740; Invitrogen, Carlsbad, USA). All samples were treated with DNase (M0303L; NEB, Ipswich, USA) to avoid DNA contamination. Profiling of m6A sites for each mRNA was performed as described in a previous study . For isolation of mRNA, total RNA was subjected to two rounds of purification using oligo (dT)-coupled magnetic beads according to the manufacturer’s instructions (Invitrogen, Carlsbad, USA). mRNA samples were fragmented into ~ 100-nucleotide fragments by incubating in fragmentation buffer (Invitrogen, Carlsbad, USA) at 90 °C for 1 min. The fragmented mRNA (approximately 50 ng) was used to construct the input library. The remaining mRNA fragments were incubated for 4 h at 4 °C with 3 mg of anti-m6A polyclonal antibody (202003; Synaptic Systems, Goettingen, Germany) combined with protein A beads (Invitrogen, Carlsbad, USA) at room temperature for 1 h in immunoprecipitation purification (IPP) buffer (150 mmol/L NaCl, 0.1% NP-40, 10 mmol/L Tris-HCl). Bound mRNA was eluted from the beads with 0.5 mg/mL of m6A in IPP buffer. Eluted mRNA was precipitated using an ethanol-sodium acetate solution and glycogen (R0551; Thermo Scientific, Waltham, USA) overnight at − 80 °C. mRNA was resuspended in water and used for library generation. Input mRNA and IP mRNA sequencing libraries were generated using a Vazyme mRNA-seq Kit for Illumina (NR601; Vazyme, Nanjing, China). Sequencing was performed using an Illumina HiSeq X Ten system (Illumina, San Diego, USA).
MeRIP sequencing data analysis
For each sample, paired-end reads were used for bioinformatics analysis. Raw data quality control was performed using FastQC software (version 0.11.3). Sequencing data were trimmed using Trimmomatic (version 0.36)  with default parameters to remove the adapter and low-quality data. These high-quality reads were mapped against the ensemble pig genome (Sscrofa11.1) using hisat2 software (version 2.0.1) . m6A-modified RNA regions (m6A peaks) and differentially m6A methylated peaks were identified using exomePeak software (version 2.16.0) [27, 28]. m6A peaks identified in two biological replicates were merged and regarded as highly enriched and selected for further analysis (P-value < 0.05).
Identification of hub 5mC/m6A regulators based on the topology of the co-expression networks
Methods for hub 5mC/m6A regulators and epigenetic module eigengenes (EME) identification were described in Chen et al.’s paper . To identify hub 5mC and m6A regulators for pig skeletal muscle, we introduced the concept of “module” from the weighted gene co-expression network analysis (WGCNA) algorithm and treated the 41 5mC and m6A regulators as a module . The overall expression level of the module was summarized as the module eigengene by the “moduleEigengenes” function in the R package WGCNA. We further calculated the module membership (i.e., module eigengene-based intramodular connectivity) as the correlation between the expression value of a given 5mC/m6A regulator and the module eigengene. Hub 5mC/m6A regulators were then defined as those that achieved a module membership greater than 0.7. The summary expression level of the identified hub 5mC/m6A regulators was again calculated as EME for pig skeletal muscle.
EME correlated co-expression analysis
Genes with the top 25% variant RNA expression level among all the 81 pig skeletal muscle samples were selected for co-expression analysis. Then we identified the co-expression modules by the WGCNA algorithm with default parameters. Correlation analysis was conducted between the EME vector and the identified co-expression modules. Modules with r > 0.8 and P-value < 0.05 were regarded as EME correlated co-expression modules for further analysis.
5mC/m6A interaction effects analysis
To find out 5mC/m6A interaction effects in the 80 up-regulated genes in the postnatal stages, we recruited the general linear model (GLM) using the postnatal days, 5mC methylation level, enrich score of m6A peaks and 5mC/m6A interaction item as explanatory variables:
Here, S denotes the postnatal days, X is the enrich score of an m6A peak, Y is the 5mC methylation level of the m6A peak, X × Y is the interaction item, and ε is the random variable. Additionally, the enrich score of a specified m6A peak is calculated as:
Here, IP is the MeRIP sample and Input is the corresponding control, A is the read counts mapped to a specified peak region in a sample, and ∑A is the total mapped read counts in a sample. Finally, regions or genes with a P-value < 0.05 for the interaction item were selected as of 5mC/m6A modification interaction effect regions or genes. All the fitted GLM were constructed by the lm function in R.
Protein-protein interaction analysis
In addition, we identified the protein-protein interactions among the co-expression genes in the EME negative correlated module based on the STRING database .
Function enrichment analysis
Gene list enrichment analysis was performed using the online tool Enrichr . Enrichment results from MSigDB hallmark pathway, GO biological process, MGI mammalian phenotype level 4, and ARCHS4 TFs co-expression with a P-value < 0.05 were considered nominally significant.
Public data access and analysis
The RNA-seq data of C2C12 myoblasts were downloaded from GEO datasets (GSE131842). Expression levels of DNA 5mC writers (DNMT3A, DNMT3B, and DNMT1) and RNA m6A writers (METTL3, METTL14, and WTAP) were extracted and followed with Pearson’s correlation analysis. The RRBS and MeRIP-seq data of C2C12 cells were also downloaded from GEO datasets with accession numbers GSE131842 and GSE154720, respectively. The CpG methylation matrix and the m6A peak tables were directly used in the further co-enrichment analysis. For 5mC/m6A co-occurrence analysis, we split the upstream 5 kb, the m6A peak region, and the downstream 5 kb sequence into 20, 10, and 20 bins with the same length, respectively. Then the bins covered with no less than 20 CpG sites were selected for average methylation level calculation and visualization. In total, 81 RNA-seq data from 27 prenatal and postnatal developmental stages (including embryonic days 33, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105 and postnatal days 0, 9, 20, 30, 40, 60, 80, 100, 120, 140, 160 and 180, each stage with three replicates) and 8 WGBS data from the four postnatal developmental stages (postnatal days 30, 60, 120 and 180, each stage with two replicates) were obtained from a previously in-house study on pig skeletal muscle development . The whole design of this study was clarified in Additional file 1: Fig. S1.
Protein-protein interactions between 5mC and m6A writers
To investigate the coordinated relations between 5mC and m6A writers, we checked the patterns of expression correlation between DNA 5mC writers (DNMT3A, DNMT3B, and DNMT1) and RNA m6A writers (METTL3, METTL14, and WTAP) in mouse myoblast cells (C2C12) after differentiation for 0, 12, 24, 48 and 60 h (3 biological replicates for each time point). Positive correlations were unveiled between the mRNA level of 5mC and m6A writers (Fig. 1A). In addition, we collected genome-wide 5mC and m6A data of a C2C12 muscle cell line (based on RRBS and MeRIP-seq data from myoblasts that had been in differentiation media for 5 d and 3 d, respectively) and found a significantly higher methylation level in the m6A peaks than in the up- or down-stream regions (Fig. 1B), which indicated co-occurrence of 5mC modification in the genome and m6A modification in the corresponding transcriptome.
We hypothesize that there are protein-protein interactions among the writers of these two types of epigenetic modifications in skeletal muscle. We performed immunofluorescence assays and showed that METTL3 co-localizes with DNMT3A and DNMT3B in the nucleus of C2C12 cells (Fig. 2A), though most of the METTL3 localizes in the cytoplasm. Further, we performed co-immunoprecipitation (co-IP) assays on the C2C12 cells differentiated for 4 d using the main DNA methylases (DNMT1/DNMT3A/DNMT3B) as bait proteins (see Materials and methods). We purified the preys and confirmed the known m6A methylases of (METTL3/METTL14/WTAP) by Western blot analyses (Fig. 2B). To minimize the biases from the experiments, we reversed the co-IP experiments by using RNA methylases as bait proteins and testing the DNA methylases via Western blot. Thereby, the relative intensity of the interaction was calculated as the average values of the two experiments (Additional file 1: Fig. S2). Previous studies have demonstrated the close interactions between the three 5mC methylases [32, 33] and between the three m6A methylases , respectively. We confirmed such close interactions in the present study and further revealed significant interactions between the 5mC and m6A methylases. Specifically, the relative intensities of the interactions between 5mC and m6A methylases are comparable to the interactions among the three individual m6A methylases.
Significantly correlated 5mC and m6A regulators in myogenesis regulation during prenatal and postnatal muscle development of pigs
To investigate the potential 5mC/m6A crosstalk in the skeletal muscle, we first obtained gene expression datasets from a previously published in-house study containing porcine skeletal muscle samples from 27 developmental stages . We checked the patterns of expression correlation between 21 5mC regulators and 20 m6A regulators, including epigenetic writers, readers and erasers that were curated in a previous study . We observed highly correlated expression patterns for the same class of regulators and even higher correlations between the expression of 5mC and m6A regulators in pig skeletal muscle (Fig. 3A), which is in accordance with previous findings in human cancer .
We further analyzed the co-expression gene modules from transcriptome data comprising 81 pig skeletal muscle samples across 27 developmental stages by applying a weighted gene co-expression network analysis (WGCNA) (Additional file 1: Fig. S3). Based on the co-expression gene network, 13 5mC regulators and 16 m6A regulators were identified as hub genes among the 41 m6A and 5mC regulators (see Materials and methods). These hub genes displayed dynamic transcriptional changes across the 27 developmental stages, i.e., mostly down-regulated in postnatal stages compared to prenatal stages (Fig. 3B). We then combined these hub genes to construct an epigenetic module eigengene (EME) to further explore the potential target genes under coordinated pre- and post-transcriptional regulation. We found three co-expression modules displayed a significant correlation with EME, among which the Lightyellow module shows a significantly positive correlation (Fig. 4A, r = 0.98, P-value = 3e-61, consisting of 2414 genes), while the Midnightblue module shows a significantly negative correlation (Fig. 4A, r = − 0.9, P-value = 3e-30, comprised of 80 genes). These results reveal that a group of genes might be under coordinated pre and post-transcriptional regulation by 5mC/m6A crosstalk. The expression heatmap indicated that these genes are highly divergent between prenatal and postnatal development and may be essential to embryonic or postnatal skeletal muscle development (Additional file 1: Fig. S4A and Fig. 4B).
Next, we examined the functional relevance of the EME correlated genes to myogenesis. Enrichment analysis reveals that the Ligthyellow module genes are involved in general RNA processes like splicing, translation and metabolism (Additional file 1: Fig. S4B). Interestingly, the Midnightblue module genes are involved in MSigDB hallmark pathways like myogenesis, hypoxia and glycolysis (Fig. 4C). GO analysis showed they participate in the biological process like glycolytic process, gluconeogenesis and carbohydrate catabolic process (Fig. 4D). They are also relevant to MGI mammalian phenotypes like increased skeletal muscle glycogen level, reticulocytosis, and skeletal muscle fiber degeneration (Fig. 4E). Co-expression analysis of transcript factors in the Midnightblue module revealed crucial known TF genes involved in human muscle development, including LBX1, MYF6, SIX1, and MYOD1 (Fig. 4F). Moreover, we found genes in the Midnightblue module frequently interacted in the protein-protein interaction networks (Additional file 1: Fig. S5).
DNA methylation and RNA m6A dynamics during the postnatal skeletal muscle development and growth in pig
To further address the 5mC and m6A dynamics across the muscle development and growth after birth, we obtained genome-wide 5mC and transcriptome-wide m6A data for the same set of skeletal muscle tissues from Landrace pigs at four postnatal time points (including postnatal days 30, 60, 120, and 180). The DNA methylome data was generated in the previous study . In the present study, MeRIP-seq was conducted, obtaining an average of 36.89 million paired-end reads (150 bp × 2) for each sample, reaching an average mapping rate of 83.92% (range 73.47 to 91.28%) (Additional file 2: Table S1).
Across the four growth stages, m6A modifications were enriched in CDS and 3’UTR regions, especially in the TSS site (Additional file 1: Fig. S6). We then called the m6A consensus peaks and the m6A contained genes in the four growth stages. Intriguingly, both the m6A peak number (8969–12, 221 peaks) and the m6A modified gene number (5324–7003 genes) showed a slight increase in postnatal 60 d and then a persistently decrease onto the postnatal 180 d (Fig. 5A). Corresponding longitudinal gene expression analysis to the prominent m6A writers (METTL3, METTL4 and WTAP) and erasers (ALKBH5 and FTO) reveals that the global m6A dynamic may be mainly driven by the cooperated regulation of METTL3 and ALKBH5 (Fig. 5B). To understand the m6A dynamics between the adjacent stages, we conducted a pairwise comparison to identify differentially methylated m6A peaks and genes. GO analysis reveals that the differentially m6A modification genes of LD30-LD60 and LD60-LD120 are enriched in typical skeletal muscle structures like sarcomere, myofibril, I band and contractile fiber, whereas the differentially m6A modification genes of D120-D180 are enriched in NADH and ATP related energy metabolic progress (Fig. 5F). Similarly, a re-analysis of the WGBS data of the corresponding samples revealed developmental stage-dependent DNA 5mC dynamics (Additional file 1: Fig. S7), as previously reported.
Identification of key candidate genes under coordinated regulation of 5mC/m6A modifications
By integrative analysis of the genome-wide 5mC and m6A sequencing data, we found a significant correlation between global DNA methylation level and m6A peak number during postnatal skeletal muscle development and growth (Fig. 6A, r = 0.954, P-value =0.046). In addition, we also observed significantly higher DNA methylation levels inside the m6A peaks than in the upstream and the downstream regions (Fig. 6B), which implicated an enrichment of DNA 5mC methylation in the corresponding RNA m6A modification sites in pig skeletal muscle. We then introduced a general linear model (GLM) to fit the 5mC/m6A interaction effects on mRNA expression for the Midnightblue module genes in the postnatal stages, under the hypothesis that the observed mRNA levels are affected by both 5mC and m6A modifications (see Materials and methods). We identified 15 regions in 9 genes with significant interaction items in the fitted model (Additional file 2: Table S2, 5mC/m6A interaction item P-value < 0.05). Intriguingly, five of these genes (FBXO40, TMOD4, PGM1, CLIP1 and ACTN3) are involved in myogenesis and three of them (PHKB, PFKFB1 and ENO3) are involved in the glycogen metabolic process (Additional file 2: Table S2).
Nucleotide methylation carries essential information for gene regulation, especially DNA 5-cytosine methylation (5mC) and RNA N6-adenosine methylation (m6A). Several previous studies indicated the essential roles of 5mC [23, 35] and m6A regulation [35, 36] during skeletal muscle development and growth, respectively. Our previous in-house study has revealed an overall decrease in DNA methylation level from the embryo to the adult pig skeletal muscle, which was highly correlated with the downregulated expression of DNMT1 . We also described the m6A methylation landscape during prenatal skeletal muscle development and revealed essential roles in the regulation of myogenesis in another recent study . This study focused on the crosstalk between these two critical epigenetic modifications in the postnatal development and growth of skeletal muscle. We validated the hypothesis that transcriptional (5mC) and posttranscriptional (m6A) regulation of hub genes is coordinated via protein-protein interaction of 5mC and m6A methylases, and disclosed a cooperative regulation of 5mC and m6A modifications on spatiotemporal gene expression in myogenesis.
The potential association of 5mC and m6A regulators was previously examined across 33 human cancer types, revealing their co-occurrence of genetic alterations and co-expression . However, no data on 5mC and m6A modifications were available in that study. By taking advantage of a large-scale study on porcine skeletal muscle across 27 developmental stages, we performed a comprehensive evaluation of 5mC and m6A crosstalk in the present study. Initially, transcriptome data analysis indicated most 5mC and m6A regulators show a positively correlated expression in the skeletal muscle samples, leading to the establishment of 5mC/m6A EME by combining hub 5mC/m6A regulators. With this, we found the EME negatively correlated genes mainly function in myogenesis, hypoxia and glycolysis, which is in line with the physiological and structural changes of the muscle fiber in the postnatal growth stages. According to ATPase activity of the muscle fibers, the muscle fibers are categorized into four fiber types: MyHC I (slow-oxidative), MyHC IIa (fast-oxidative), MyHC IIx (intermediary to MyHC IIa and MyHC IIb), and MyHC IIb (fast-glycolytic) . In pigs, oxidative metabolism represents the principal energy source during fetal life. At birth, glycolytic metabolism dramatically increases during the first postnatal weeks . After that, the diameter of muscle fibers increases, and the muscle fibers lengthen rapidly . Our results indicate that a down-regulated EME, which reflects both the pre- and post-transcriptional epigenetic status, may correlate to the enhanced gene co-expression in the postnatal muscle ontogenesis.
Moreover, immunofluorescence colocalization and co-IP experiments confirmed the interaction between 5mC and m6A methylases (Fig. 2). We further established an m6A epitranscriptome atlas of skeletal muscle and depicted the dynamics during postnatal skeletal muscle development for the first time in pigs. More m6A-contained genes and more variable changes of m6A modification in m6A-contained genes were identified in postnatal stages than prenatal stages . These results indicated the crucial role of m6A modifications in postnatal skeletal muscle development, which is in line with previous documents [36, 39]. The m6A peak number and the DNA methylation level of pig skeletal muscle are changed simultaneously across different development stages. By integrating the 5mC and m6A sequencing data, co-localization of genome-wide 5mC and m6A modification was observed in porcine skeletal muscle. We further identified nine genes from the EME negative correlated module, which showed significant 5mC/m6A interaction effects on their mRNA transcription during postnatal skeletal muscle development and are crucial in myogenesis or glycogen metabolic process in muscle. Among them, ACTN3 (actinin alpha 3) is a structural component of the sarcomeric Z line and is related to sarcopenia and physical fitness inactive older women [40, 41]. PHKB (phosphorylase kinase regulatory subunit beta) is involved in glycogen metabolism and associated with glycogen storage disease IXb, a disease also known as phosphorylase kinase deficiency in liver and muscle .
To our best knowledge, this is the first study suggesting protein-protein interactions between 5mC and m6A methylases, either directly or indirectly. We speculated that proteins like the histone deacetylase HDAC2 and the E3 ubiquitin protein ligase UHRF1 might act as adaptors to guide the interactions, as previous studies indicated they both could interact with the 5mC and m6A writers separately (Additional file 1 Fig. S8) [43,44,45]. Based on current data, we reasoned that 5mC methylase and m6A methylase are likely to add 5mC and m6A modifications simultaneously, thus enabling the synergistic transcriptional and post-transcriptional regulation of gene expression. Considering m6A is an essential mechanism for mRNA metabolism and protein translation, our results suggest that simultaneous nuclear events of transcriptional and translational regulation are coordinated for skeletal muscle cells. Despite this, the 5mC/m6A crosstalk needs further validation in different tissues and cell types from other species. Especially, although we revealed a number of key myogenesis-related genes that may be important targets of this synergistic regulation, further investigations to clarify the spatiotemporal mechanism of regulation and its developmental effects need to be conducted.
In summary, our results lay a foundation for epigenetic regulation in skeletal muscle development and pave a new road for research on either basic muscle biology or animal breeding and muscle-related diseases. Since the development of demethylating drugs is currently a research hot spot, especially for RNA m6A inhibitors, the genes we supposed under tightly epigenetic regulation can serve as candidate targets for the treatment of muscle-related diseases.
Availability of data and materials
All the sequencing data of MeRIP-seq have been deposited to Gene Expression Omnibus (GEO) at the National Center for Biotechnology Information (NCBI) under accession number GSE141943.
Epigenetic module eigengenes
General linear model
Methylated RNA immunoprecipitation sequencing
Reduced representation bisulfite sequencing
Whole-genome bisulfite sequencing
Weighted correlation network analysis
Cabello-Verrugio C, Rivera JC, Garcia D. Skeletal muscle wasting: new role of nonclassical renin–angiotensin system. Curr Opin Clin Nutr Metab Care. 2017;20(3):158–63. https://doi.org/10.1097/MCO.0000000000000361.
Hua Y, Zhang Y, Ren J. IGF-1 deficiency resists cardiac hypertrophy and myocardial contractile dysfunction: role of microRNA-1 and microRNA-133a. J Cell Mol Med. 2012;16(1):83–95. https://doi.org/10.1111/j.1582-4934.2011.01307.x.
Kivelä R, Salmela I, Nguyen YH, Petrova TV, Koistinen HA, Wiener Z, et al. The transcription factor Prox1 is essential for satellite cell differentiation and muscle fibre-type regulation. Nature Comm. 2016;7:13124. https://doi.org/10.1038/ncomms13124.
Lefaucheur L. A second look into fibre typing–relation to meat quality. Meat Sci. 2010;84(2):257–70. https://doi.org/10.1016/j.meatsci.2009.05.004.
Schook L, Beattie C, Beever J, Donovan S, Jamison R, Zuckermann F, et al. Swine in biomedical research: creating the building blocks of animal models. Anim Biotechnol. 2005;16(2):183–90. https://doi.org/10.1080/10495390500265034.
Tang Z, Li Y, Wan P, Li X, Zhao S, Liu B, et al. LongSAGE analysis of skeletal muscle at three prenatal stages in Tongcheng and landrace pigs. Genome Biol. 2007;8:R115. https://doi.org/10.1186/gb-2007-8-6-r115.
Watanabe A, Yamada Y, Yamanaka S. Epigenetic regulation in pluripotent stem cells: a key to breaking the epigenetic barrier. Philos Trans R Soc B Biol Sci. 2013;368(1609):20120292. https://doi.org/10.1098/rstb.2012.0292.
Cardoso-Moreira M, Halbert J, Valloton D, Velten B, Chen C, Shao Y, et al. Gene expression across mammalian organ development. Nature. 2019;571(7766):505–9. https://doi.org/10.1038/s41586-019-1338-5.
Schübeler D. Function and information content of DNA methylation. Nature. 2015;517(7534):321–6. https://doi.org/10.1038/nature14192.
Ru W, Zhang X, Yue B, Qi A, Shen X, Huang Y, et al. Insight into m6A methylation from occurrence to functions. Open Biol. 2020;10(9):200091. https://doi.org/10.1098/rsob.200091.
Chen K, Zhao BS, He C. Nucleic acid modifications in regulation of gene expression. Cell Chem Biol. 2016;23(1):74–85. https://doi.org/10.1016/j.chembiol.2015.11.007.
Yuan R, Zhang J, Wang Y, Zhu X, Hu S, Zeng J, et al. Reorganization of chromatin architecture during prenatal development of porcine skeletal muscle. DNA Res. 2021;28(2):dsab003. https://doi.org/10.1093/dnares/dsab003.
Hernández-Hernández O, Ávila-Avilés RD, Hernández-Hernández JM. Chromatin landscape during skeletal muscle differentiation. Front Genet. 2020;11:578712. https://doi.org/10.3389/fgene.2020.578712.
Li X, Fu L, Cheng H, Zhao S. Advances on microRNA in regulating mammalian skeletal muscle development. Hereditas. 2017;39(11):1046–53. https://doi.org/10.16288/j.yczz.17-112.
Zhang Y, Yao Y, Wang Z, Lu D, Zhang Y, Adetula AA, et al. MiR-743a-5p regulates differentiation of myoblast by targeting Mob1b in skeletal muscle development and regeneration. Genes Dis. 2020;9(4):1038–48. https://doi.org/10.1016/j.gendis.2020.11.018.
Yan J, Yang Y, Fan X, Tang Y, Tang Z. Sp1-mediated circRNA circHipk2 regulates myogenesis by targeting ribosomal protein Rpl7. Genes. 2021;12(5):696. https://doi.org/10.3390/genes12050696.
Yan J, Yang Y, Fan X, Liang G, Wang Z, Li J, et al. circRNAome profiling reveals circFgfr2 regulates myogenesis and muscle regeneration via a feedback loop. J Cachexia Sarcopenia Muscle. 2022;13(1):696–712. https://doi.org/10.1002/jcsm.12859.
Blum R, Dynlacht BD. The role of MyoD1 and histone modifications in the activation of muscle enhancers. Epigenetics. 2013;8(8):778–84. https://doi.org/10.4161/epi.25441.
Kudou K, Komatsu T, Nogami J, Maehara K, Harada A, Saeki H, et al. The requirement of Mettl3-promoted MyoD mRNA maintenance in proliferative myoblasts for skeletal muscle differentiation. Open Biol. 2017;7(9):170119. https://doi.org/10.1098/rsob.170119.
Zhang X, Yao Y, Han J, Yang Y, Chen Y, Tang Z, et al. Longitudinal epitranscriptome profiling reveals the crucial role of N6-methyladenosine methylation in porcine prenatal skeletal muscle development. J Genet Genom. 2020;47(8):466–76. https://doi.org/10.1016/j.jgg.2020.07.003.
Chen Y, Shen J, Chen D, Wu C, Guo R, Zhang P, et al. Identification of cross-talk between m 6 a and 5mC regulators associated with onco-immunogenic features and prognosis across 33 cancer types. J Hematol Oncol. 2020;13:22. https://doi.org/10.1186/s13045-020-00854-w.
Zhou L, Tian S, Qin G. RNA methylomes reveal the m6A-mediated regulation of DNA demethylase gene SlDML2 in tomato fruit ripening. Genome Biol. 2019;20(1):156. https://doi.org/10.1186/s13059-019-1771-7.
Yang Y, Fan X, Yan J, Chen M, Zhu M, Tang Y, et al. A comprehensive epigenome atlas reveals DNA methylation regulating skeletal muscle development. Nucleic Acids Res. 2021;49(3):1313–29. https://doi.org/10.1093/nar/gkaa1203.
Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar L, Osenberg S, et al. Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature. 2012;485(7397):201–6. https://doi.org/10.1038/nature11112.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20. https://doi.org/10.1093/bioinformatics/btu170.
Kim D, Paggi JM, Park C, Bennett C, Salzberg SL. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat Biotechnol. 2019;37(8):907–15. https://doi.org/10.1038/s41587-019-0201-4.
Meng J, Cui X, Rao MK, Chen Y, Huang Y. Exome-based analysis for RNA epigenome sequencing data. Bioinformatics. 2013;29(12):1565–7. https://doi.org/10.1093/bioinformatics/btt171.
Meng J, Lu Z, Liu H, Zhang L, Zhang S, Chen Y, et al. A protocol for RNA methylation differential analysis with MeRIP-Seq data and exomePeak R/Bioconductor package. Methods. 2014;69(3):274–81. https://doi.org/10.1016/j.ymeth.2014.06.008.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559. https://doi.org/10.1186/1471-2105-9-559.
Szklarczyk D, Gable AL, Nastou KC, Lyon D, Kirsch R, Pyysalo S, et al. The STRING database in 2021: customizable protein–protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. 2021;49(D1):D605–D12. https://doi.org/10.1093/nar/gkaa1074.
Xie Z, Bailey A, Kuleshov MV, Clarke DJB, Evangelista JE, Jenkins SL, et al. Gene set knowledge discovery with enrichr. Curr Protoc. 2021;1(3):e90. https://doi.org/10.1002/cpz1.90.
Kim G, Ni J, Kelesoglu N, Roberts RJ, Pradhan S. Co-operation and communication between the human maintenance and de novo DNA (cytosine-5) methyltransferases. EMBO J. 2002;21(15):4183–95. https://doi.org/10.1093/emboj/cdf401.
Rigbolt KTG, Prokhorova TA, Akimov V, Henningsen J, Johansen PT, Kratchmarova I, et al. System-wide temporal characterization of the proteome and phosphoproteome of human embryonic stem cell differentiation. Sci Signal. 2011;4(164):rs3. https://doi.org/10.1126/scisignal.2001570.
Liu J, Yue Y, Han D, Wang X, Fu Y, Zhang L, et al. A METTL3–METTL14 complex mediates mammalian nuclear RNA N6-adenosine methylation. Nat Chem Biol. 2014;10(2):93–5. https://doi.org/10.1038/nchembio.1432.
Pan Y, Chen L, Cheng J, Zhu X, Wu P, Bao L, et al. Genome-wide DNA methylation profiles provide insight into epigenetic regulation of red and white muscle development in Chinese perch Siniperca chuatsi. Comp Biochem Physiol B Biochem Mol Biol. 2021;256:110647. https://doi.org/10.1016/j.cbpb.2021.110647.
Li J, Pei Y, Zhou R, Tang Z, Yang Y. Regulation of RNA N6-methyladenosine modification and its emerging roles in skeletal muscle development. Int J Biol Sci. 2021;17(7):1682. https://doi.org/10.7150/ijbs.56251.
Wang H, Pu J, Chen D, Tian G, Mao X, Yu J, et al. Effects of dietary amylose and amylopectin ratio on growth performance, meat quality, postmortem glycolysis and muscle fibre type transformation of finishing pigs. Arch Anim Nutr. 2019;73(3):194–207. https://doi.org/10.1080/1745039X.2019.1583518.
Picard B, Lefaucheur L, Berri C, Duclos MJ. Muscle fibre ontogenesis in farm animal species. Reprod Nutr Dev. 2002;42(5):415–31. https://doi.org/10.1051/rnd:2002035.
Gheller BJ, Blum JE, Fong EHH, Malysheva OV, Cosgrove BD, Thalacker-Mercer AE. A defined N6-methyladenosine (m6A) profile conferred by METTL3 regulates muscle stem cell/myoblast state transitions. Cell Death Discov. 2020;6:95. https://doi.org/10.1038/s41420-020-00328-5.
Romero-Blanco C, Artiga Gonzalez MJ, Gómez-Cabello A, Vila-Maldonado S, Casajús JA, Ara I, et al. ACTN3 R577X polymorphism related to sarcopenia and physical fitness in active older women. Climacteric. 2021;24:89–94. https://doi.org/10.1080/13697137.2020.1776248.
Kiuchi Y, Makizako H, Nakai Y, Taniguchi Y, Tomioka K, Sato N, et al. Associations of alpha-actinin-3 genotype with thigh muscle volume and physical performance in older adults with sarcopenia or pre-sarcopenia. Exp Gerontol. 2021;154:111525. https://doi.org/10.1016/j.exger.2021.111525.
Beyzaei Z, Ezgu F, Geramizadeh B, Alborzi A, Shojazadeh A. Novel mutations in the PHKB gene in an iranian girl with severe liver involvement and glycogen storage disease type IX: a case report and review of literature. BMC Pediatr. 2021;21:175. https://doi.org/10.1186/s12887-021-02648-6.
Yue Y, Liu J, Cui X, Cao J, Luo G, Zhang Z, et al. VIRMA mediates preferential m6A mRNA methylation in 3′ UTR and near stop codon and associates with alternative polyadenylation. Cell Discov. 2018;4:10. https://doi.org/10.1038/s41421-018-0019-0.
Schwartz S, Mumbach MR, Jovanovic M, Wang T, Maciag K, Bushkin GG, et al. Perturbation of m6A writers reveals two distinct classes of mRNA methylation at internal and 5′ sites. Cell Rep. 2014;8(1):284–96. https://doi.org/10.1016/j.celrep.2014.05.048.
Ignatova VV, Jansen PWTC, Baltissen MP, Vermeulen M, Schneider R. The interactome of a family of potential methyltransferases in HeLa cells. Sci Rep. 2019;9:6584. https://doi.org/10.1038/s41598-019-43010-2.
FG was supported by the Agricultural Science and Technology Innovation Program and The Elite Young Scientists Program of CAAS. ZT was supported by the National Natural Science Foundation of China (31830090), the Basic and Applied Basic Research Foundation of Guangdong province (2019B1515120059), and the Shenzhen Dapeng New District Special Fund for Industry Development (KY20180114), and the Agricultural Science and Technology Innovation Program (CAAS-ZDRW202006).
Ethics approval and consent to participate
All animal procedures were performed according to protocols of the Chinese Academy of Agricultural Sciences and the Institutional Animal Care and Use Committee.
Consent for publication
The authors declare that they have no competing interests.
The whole design of this study. Fig. S2. Relative interaction intensities between the 5mC and m6A writer protein pairs based on the western blot images from co-IP assays. Fig. S3. Identification of WGCNA co-expression modules based on the top 25% expressing variant genes from the 27 prenatal and postnatal skeletal muscle development stages. Fig. S4. Potential regulation role of m6A/5mC EME on prenatal myogenesis of pigs. (A) EME positive correlated Lightyellow module and its gene expression heatmap. (B) GO process enrichment analysis of Midnightblue module genes. Fig. S5. Protein-Protein Interaction network among the Midnightblue module genes. Fig. S6. Distribution of m6A peaks across the whole transcript. Fig. S7. Dynamics of the DNA methylation during postnatal skeletal muscle development in pigs. (A) mRNA expression profiles of 5mC writers (DNMT1/DNMT3A/DNMT3B) and erasers (TET1/TET2/TET3). (B) Numbers of 5mC hypermethylated or hypomethylated genes from inter-stage comparisons (C) Top 5 GO terms enriched by the inter-stage differentially 5mC modification genes. Fig. S8. Possible indirect interactions between 5mC writers and m6A writers. (A) Protein-protein interaction network shows the proteins interacted with 5mC writers (DNMT1/3A/3B) in human beings. (B) Experiment evidence found in the published documents indicated the linker protein (UHRF1 and HDAC1) could interact with the m6A MTC proteins.
Summary of MeRIP-seq data in this study. Table S2. Genes expressed under significant 5mC and m6A interaction effects in postnatal skeletal muscle development. Table S3–S6. RNA expression data were used in this study. Table S7–S10. Consensus m6A peak identification of the four postnatal stages by exomePeak.
About this article
Cite this article
Zhang, D., Wu, S., Zhang, X. et al. Coordinated transcriptional and post-transcriptional epigenetic regulation during skeletal muscle development and growth in pigs. J Animal Sci Biotechnol 13, 146 (2022). https://doi.org/10.1186/s40104-022-00791-3
- DNA methylation
- Epigenetic modification
- Epigenomic analysis
- 5mC regulators
- m6A methylation
- m6A regulators