Profiling Analysis of N6-Methyladenosine mRNA Methylation Reveals Differential m6A Patterns during the Embryonic Skeletal Muscle Development of Ducks

Simple Summary Recent studies show that N6-methyladenosine (m6A) modification, the most common RNA chemical modification, influences the modification, processing, transport, and translation of RNA. N6-methyladenosine is an epigenetic modification that influences skeletal myogenesis and skeletal muscle development. However, the N6-methyladenosine modification profile and its function during poultry muscle development is unclear, and there is only one report about m6A modification in ducks, which focuses on duck hepatitis A virus infection. Here, we compared the m6A modification profiles between E13 (embryonic day 13) and E19 (embryonic day 19) in duck breast muscle differentiation using MeRIP-seq, and evaluated the expression profile of the methyl transferase METTL14 and its cofactors during breast muscle development. This is the first study of N6-methyladenosine modification patterns in duck muscle tissue. The current study not only elucidates the regulation mechanisms of duck skeletal muscle development, but also lays the groundwork for studying the role of RNA modification in poultry muscle development. Abstract N6-Methyladenosine is a reversible epigenetic modification that influences muscle development. However, the m6A modification profile during poultry skeletal muscle development is poorly understood. Here, we utilized m6A-specific methylated RNA immunoprecipitation sequencing to identify m6A sites during two stages of breast muscle development in ducks: embryonic days 13 (E13) and E19. MeRIP-seq detected 19,024 and 18,081 m6A peaks in the E13 and E19 groups, respectively. Similarly to m6A distribution in mammalian transcripts, our results revealed GGACU as the main m6A motif in duck breast muscle; they also revealed that m6A peaks are mainly enriched near the stop codons. In addition, motif sequence analysis and gene expression analysis demonstrated that m6A modification in duck embryo skeletal muscles may be mediated by the methyltransferase-like 14. GO and KEGG analysis showed that m6A peaks containing genes at E19 were mainly enriched in muscle-differentiation- and muscle-growth-related pathways, whereas m6A peaks containing genes in E13 were mainly enriched in embryonic development and cell proliferation pathways. Combined analysis of MeRIP-seq and RNA-seq showed that the mRNA expression may be affected by m6A modification. Moreover, qRT-PCR analysis of the expression of METTL14 and its cofactors (WTAP, ZC3H13, RBM15 and VIRMA) during duck embryonic skeletal muscle development in breast and leg muscle samples revealed a significant downward trend as the developmental age progressed. Our results demonstrated that m6A mRNA methylation modifications control muscle development in ducks. This is the first study of m6A modification patterns in duck muscle tissue development, and it lays the foundation for the study of the effects of RNA modification on poultry skeletal muscle development.


Introduction
Poultry meat contributes to human nutrition by providing high-quality protein [1]. Poultry muscle development is a complex process which is influenced by many factors, such as genetics, nutrition, and disease [2][3][4]. However, genetics plays the main role [5,6]. Myogenesis occurs in the following stages: the differentiation of somites into myoblasts, the migration of myoblasts to the limb buds and the expression of specific myogenic transcription factors, and the differentiation of myoblasts into myotubes under the influence of these myogenic transcription factors [7]. The transition point from the proliferation to the fusion of duck embryonic breast muscle development occurs at embryonic day 19 (E19), and this coincides with the highest expression level of myogenic markers, such as myogenin, MRF4, and myostatin [8].
Numerous studies have investigated the regulatory role of m6A during muscle development [18]. Although m6A plays an important role in livestock growth and development [18][19][20][21][22], its role in poultry muscle development is still unknown. The purpose of the current study was to profile the mRNA m6A modification at different stages in duck breast muscle development. Here, undifferentiated (E13) and differentiated (E19) duck breast muscles were collected and subjected to MeRIP-seq and mRNA-seq analysis to profile m6A modification sites during duck breast muscle development. The expression levels of METTL14, WTAP, ZC3H13, RBM15, and VIRMA were determined at different stages during breast and leg muscles development. The data generated from this study can be utilized in future studies of m6A function in poultry skeletal muscle development.
from the incubator. The interval for egg turning was set at 2 h during the incubation, and the egg-turning was halted after all eggs were placed in the hatcher tray. Tissue collection was conducted as described in our previous work [23]. Breast muscles and leg muscles from 15 embryos were collected every 3 days from embryonic day 10 (E10) to 1 day post hatch (E10, E13, E16, E19, E22, and 1 day post hatch). The embryos were carefully taken out and their developmental stages were confirmed by morphological observation. Next, breast muscles and leg muscles were sampled, snap-frozen in liquid nitrogen, and stored at −80 • C. Total DNA from duck embryo breast muscles was extracted using a Trelief TM Animal Genomic DNA Kit (Tsingke, Beijing, China) following the manufacturer's instructions. The sex of the embryos was determined by PCR amplification analysis using sex-specific primers of chromodomain helicase DNA-binding protein 1 (CHD1) gene [24]. A total of 6, 9, 5, 6, and 7 female embryos and 9 female ducklings were identified from E10, E13, E16, E19, E22 and 1 day post hatch (P1), respectively. Five females from each period were selected and used in this study.

RNA Extraction and m6A-Specific Methylated RNA Immunoprecipitation (MeRIP)
Total RNA was extracted from the breast muscle and leg muscle tissues of E10, E13, E16, E19, E22, and P1 ducks using TRIzol reagent (Invitrogen, CA, USA), following the manufacturer's instructions. Three breast muscle samples of E13 and E19 were chosen for MeRIP-seq. To remove any DNA contamination from the samples, DNase I was used after RNA extraction. RNA quality was determined on Nanodrop 1000 (NanoDrop, Wilmington, DE, USA). RNA integrity was confirmed by 1.5% agarose gel electrophoresis and quantified by Qubit3.0 using a QubitTM RNA Broad Range assay kit (Life Technologies, Carlsbad, CA, USA). In total, 50 µg of total RNA was used for polyadenylated RNA enrichment using VAHTS mRNA capture beads (Vazyme, Nanjing, China). Then, 20 mM ZnCl 2 was added to the mRNA and incubated at 95 • C for 10 min to generate 100-200 nt RNA fragments. Next, 10% of the RNA fragments were saved as "input" and the rest were used for m6A immunoprecipitation (IP) using a specific anti-m6A antibody (Synaptic Systems, Göttingen, Germany). The IP experiment was carried out according to the protocol described in previous research [25].

Library Constructions and Sequencing
For MeRIP-seq, the stranded RNA sequencing library was constructed using a KC-DigitalTM Stranded mRNA Library Prep kit from Illumina ® (Seqhealth, Wuhan, China), following the manufacturer's instructions. The kit eliminates duplication bias in the PCR and sequencing steps by using a unique molecular identifier (UMI) of 8 random bases to label pre-amplified cDNA. The library products were enriched, quantified, and sequenced on a Novaseq 6000 sequencer (Illumina, San Diego, CA, USA) with pair-end 150.
To analyze mRNA expression, in addition to the RNA-seq with input, total RNAs were used to construct a library of 200-500 nt RNA fragments to enhance the accuracy of gene expression results. Then, 2 µg of RNA were used for stranded RNA sequencing library preparation, using a protocol similar to the one used for MeRIP-seq, except for 200-500 nt library products. All sequencing data in this study can be found in the NCBI Sequence Read Archive (SRA) (https://submit.ncbi.nlm.nih.gov/, accessed on 22 September 2022) with the accession numbers PRJNA725663 (MeRIP-seq data) and PRJNA726590 (RNA-seq data).

MeRIP-Seq and RNA-Seq Data Analysis
For MeRIP-seq: first, raw sequencing data were filtered using Trimmomatic [26]. Lowquality reads were then discarded and reads contaminated with adaptor sequences were trimmed. Clean reads were further treated to minimize duplication. In brief, clean reads were first clustered according to the UMI sequence, whereby reads with the same UMI sequence were grouped into the same cluster. Reads in the same cluster were compared using pairwise alignment and reads with sequence identities >95% were extracted to a new sub-cluster. After all sub-clusters were generated, multiple sequence alignment was used to obtain one consensus sequence for each sub-cluster.
The de-duplicated consensus sequences were used for m6A peak analysis. They were mapped to NCBI's duck (Anas platyrhynchos) reference genome (https://www.ncbi. nlm.nih.gov/assembly/GCF_003850225.1/, accessed on 22 September 2022) using STAR software (Version 2.5.3a) with default parameters. ExomePeak version 3.8 was used for peak calling. The m6A peaks were annotated using bedtools version 2.25. DeepTools version 2.4.1 was used for peak distribution analysis. Next, peaks present in both groups were classified as intersection peaks, and group-unique peaks were classified as specific peaks. Differential m6A peaks were identified in the intersection peaks using Fisher's test with p < 0.05 and |log2FC| > 1 as the cutoff threshold. Sequence motifs enriched in m6A peak regions were identified using Homer version 4.10. Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses for annotated genes were carried out using KOBAS version: 2.1.1 with a corrected p value cutoff of 0.05 to judge statistically significant enrichment. Protein-protein interaction analysis was carried out on STRING version 11.0 and adjusted on Cytoscape version 3.7.2).
For RNA-seq analysis, raw data and reads mapping were processed same as for MeRIP-seq. Reads mapped to exons were summarized by featureCounts (Subread-1.5.1; Bioconductor) and reads per kilobase per million mapped reads (RPKM) were calculated. Differentially expressed genes (DEGs) in E19 vs. E13 were identified using edge R package (version 3.12.1) using p < 0.05 and |log2FC| > 1 as cutoffs. Enrichment and proteinprotein interaction analyses were conducted using the procedure used for MeRIP-seq. For Pearson correlation analysis of mRNA expression and m6A abundance, RPKM and m6A enrichment were log2 transformed and m6A enrichment was calculated by dividing IP RPKM by input RPKM.

qRT-PCR Assay
Total RNAs were reverse transcribed using the Monad MonScript™ All-in-One Kit with DNase (Biopro, Shanghai, China). The expression level was determined using 2 × T5 Fast qPCR Mix (TsingKe, Beijing, China) in a final volume of 20 µL. GAPDH served as the internal control. The qRT-PCR reaction was carried out on an ABI QuantStudio 5 system (Thermo Fisher, Waltham, MA, USA) under the following protocol: 95 • C for 3 min; 40 cycles of 95 • C for 10 s, Tm for 30 s min, and fluorescence collection at 65-95 • C. Each sample was analyzed in triplicate and the relative mRNA expression level was calculated using the 2 −∆∆Ct method. The qRT-PCR primers' sequences are included in Table S1.
The qPCR results are presented as the mean ± standard error of the mean (S.E.M.). For multiple comparison analysis, groups were compared with a one-way ANOVA test followed by a Duncan test using SPSS 26.0 (https://www.ibm.com/support/pages/ downloading-ibm-spss-statistics-26, accessed on 22 September 2022). The different letters between the two groups represent significant differences (p < 0.05). For two-group comparison, the results were subjected to statistical analysis using the two-tailed student's t-test. The level of significance was presented as * (p < 0.05), ** (p < 0.01), and *** (p < 0.001).

Transcriptome-Wide Detection of m6A Peaks in Duck Embryonic Breast Muscle
MeRIP-seq analysis of six embryonic breast muscle tissues generated 91.7 Gb in data (1,426,191,010 raw reads). After data quality control, 672,172,477 clean reads were aligned to the duck reference genome. This analysis revealed a unique mapping ratio of 83.46 to 88.22% (Table 1). RNA-seq analysis of the same samples used in MeRIP-seq produced 38.6 Gb of data and a total of 429,253,704 reads were mapped to the reference genome with a unique mapping range of 82.78-88.06% (Table 2). Furthermore, the correlation coefficient value of MeRIP-seq samples of the same group was high (about 0.97) ( Figure S1). MeRIP-seq identified 19,024 m6A peaks in E13 breast muscle tissues and 18,081 m6A peaks in E19 group ( Figure 1A), with 17,741 peaks common to the two groups, 1283 peaks unique to E13, and 340 peaks unique to E19 ( Figure 1B). Using homer to scan motifs between peaks, we identified GGACU as the most likely motif sequence in duck embryo breast muscle tissue ( Figure 1C). A count of each gene's peak number found that >73% of the mRNAs had 1-2 peaks, and about 13% of the genes had ≥4 m6A peaks in both groups ( Figures 1D and S2). Counting the number of m6A peak-associated genes in each group, we detected 8915 m6A peak genes in E13 tissue and 8546 m6A peak genes in E19 samples. Comparison of the two sets of genes identified 7703 intersection genes (Table S2), while 1212 and 843 genes were specific to E13 and E19, respectively ( Figure 1E, Tables S3 and S4). Information about the functional regions of genes, where all peaks were located, can be determined by peak location. In embryonic duck breast muscle, m6A peaks were mainly distributed in the anterior and posterior regions of the stop codon (the CDS and 3 UTR regions close to the stop codon, Figure 1F). According to annotation and statistics, the distribution of peaks in different regions were introns beyond the coding gene (PR intron, about 37%), CDS region (about 29%), and 3 UTR region (about 20%, Figure 1G).

GO and KEEG Pathway Analysis of Differentially m6A-Methylated mRNAs in E13 and E19
To explore the important biological roles and pathways m6A actively participates in during the embryonic development of duck breast muscle, we subjected differentially m6A-methylated mRNAs in E13 and E19 to GO and KEGG enrichment analysis. GO biological process (BP) analysis of all the differentially expressed genes revealed that the m6A-containing genes are significantly enriched in RNA splicing, cell metabolism, and cell component organization (Figure 2A, Table S5). In addition, the KEGG pathway revealed m6A-containing genes to be significantly involved in TGF-beta, VEGF, and cell proliferation and embryonic development pathways ( Figure 2B, Table S6). In E19, the KEGG enrichment analysis found that m6A-containing genes were significantly enriched in MAPK and TNF signaling ( Figure 2C, Table S7), while in E13, the m6A-containing genes were mainly enriched in signaling pathways regulating stem cell pluripotency, cAMP signaling, and signal pathways involved in cell pluripotency and proliferation ( Figure 2D, Table S8).

GO and KEEG Pathway Analysis of Differentially m6A-Methylated mRNAs in E13 and E19
To explore the important biological roles and pathways m6A actively participates in during the embryonic development of duck breast muscle, we subjected differentially m6A-methylated mRNAs in E13 and E19 to GO and KEGG enrichment analysis. GO biological process (BP) analysis of all the differentially expressed genes revealed that the m6A-containing genes are significantly enriched in RNA splicing, cell metabolism, and cell component organization (Figure 2A, Table S5). In addition, the KEGG pathway revealed m6A-containing genes to be significantly involved in TGF-beta, VEGF, and cell proliferation and embryonic development pathways ( Figure 2B, Table S6). In E19, the KEGG enrichment analysis found that m6A-containing genes were significantly enriched in MAPK and TNF signaling ( Figure 2C, Table S7), while in E13, the m6A-containing genes were mainly enriched in signaling pathways regulating stem cell pluripotency, cAMP signaling, and signal pathways involved in cell pluripotency and proliferation (Figure 2D, Table S8).  To investigate the dynamic changes in m6A modification at different stages (E19 vs. E13) of embryonic duck skeletal muscle development, the abundance of m6A modification at intersection peaks between E13 and E19 were compared. Thresholds of p value < 0.05 and |log2FC| > 1 were used to identify differentially methylated peaks (DMPs). Differentially methylated genes (DMGs), which are genes associated with DMPs, were submitted to the enrichment analysis. Based on the threshold, 355 m6A peaks were differentially abundant in the 17,741 intersection peaks in E13 and E19 samples ( Figure 3A, Table S9). The DMPs had 230 peaks that were significantly hypo-methylated and 125 peaks that were noticeably hyper-methylated ( Figure 3B). DMPs annotation identified 331 DMGs (Table S10). KEGG pathway enrichment analysis of the DMGs revealed that the hypermethylated peak genes were significantly enriched in 21 pathways, including the AMPK signaling pathway and the neuroactive ligand-receptor interaction pathway ( Figure 3C, Table S11). The hypo-methylated peak genes were significantly enriched in the pathways related to cell proliferation and growth, including cAMP and GnRH signaling ( Figure 3D, Table S12). Protein-protein interaction network analysis of DMGs on STRING identified SRY-Box Transcription Factor 2 (SOX2), Distal-Less Homeobox 5 (DLX5), and F-Box Protein 40 (FBXO40), which are involved in myogenesis and muscle growth, as central factors in the network ( Figure 3E).

Identification of Differentially Expressed Genes (DEGs) in E19 vs. E13 Samples
To investigate gene expression changes in E13 vs. E19 samples, total RNA from the samples subjected to MeRIP-seq analysis was used to construct the RNA-seq libraries and DEGs were identified using a threshold of p value < 0.05 and |log2FC| > 1. This analysis revealed 2880 DEGs in E19 vs. E13. Of these, 1381 were upregulated and 1499 were downregulated ( Figure 4A,B, Table S13). Heatmap visualization of the DEGs revealed good repeatability across the three samples in each group and separate clustering of upregulated and downregulated genes ( Figure 4C). Next, the DEGs were subjected to enrichment and protein-protein interaction network analysis. GO biological process enrichment analysis showed that the DEGs were significantly enriched in muscle organ development, embryonic skeletal system morphogenesis, and negative regulation of the canonical Wnt signaling pathway ( Figure 4D, Table S14). KEGG enrichment analysis of the DEGs revealed that they were significantly enriched in signaling pathways, including the calcium, Wnt, TGF-beta, and melanogenesis pathways ( Figure 4E, Table S15). Interestingly, melanogenesis was also enriched in the hypo-methylated peak genes (Table S15). Subsequently, the top 300 most significant DEGs (sorted by p value from small to large) were submitted to protein-protein interaction analysis. DEGs were mainly clustered into two parts based on up-and downregulation, and proteins including myozenin (MYOZ) family members and myomesin 2 (MYOM2), which are associated with skeletal muscle differentiation and muscle fiber fusion, were at the core of the upregulated cluster ( Figure 4F). Bone morphogenetic protein (BMP) and fibroblast growth factor receptor (FGFR) families, which regulate the cell cycle and cell proliferation, were at the center of the downregulated cluster ( Figure 4F).

Conjoint Analysis of MeRIP-Seq and mRNA-Seq in E19 and E13
To investigate the relationship between m6A modification and mRNA expression, m6A peak regions in E19 and E13 samples were cross-compared with the DEGs from RNA-seq data. This analysis identified 129 upregulated and 57 downregulated m6A peakcontaining DEGs in E19, and 45 upregulated and 304 downregulated m6A peak-containing DEGs in E13 ( Figure 5A, Table S16). In the E19 m6A peak-containing genes, pyruvate dehydrogenase kinase 4 (PDK4) and myosin light chain kinase 2 (MYLK2) are related to muscle differentiation and growth. Animals 2022, 12, x 9 of 19  nodes stand for the hyper-methylated peak genes; the pink nodes indicate the hypo-methylated peak genes.  development ( Figure 6D), whereas its expression in breast muscle samples showed a downward tendency ( Figure 6C). Moreover, we performed qRT-PCR validation for the RNA-seq and the results were consistent with the RNA-seq data (Figure 7). These data showed that METTL14 and the RNA methyltransferase cofactors' expression were drastically downregulated from pre-differentiation (before E13) to post-differentiation (after E19) in duck skeletal muscles. We speculate that METTL14-mediated m6A modification plays a crucial role in embryonic skeletal muscle development. The m6A abundances in IRX5 mRNA for the E19 and E13 groups. The m6A peak, which was significantly higher in E19 than in E13, is shown in the black rectangle (p < 0.05 and |log2FC| > 1). (D) A plot of m6A peak abundance and mRNA expression in E19 and E13.
Next, dynamically changed peak genes and their expression levels were plotted, and the genes with |log2FC| > 1 in both MeRIP-seq and RNA-seq analyses were screened out ( Figure 5B). This analysis revealed that 44 hyper-methylated genes were upregulated (hyper-up genes) and 49 were downregulated (hyper-down genes). For hypo-methylated genes, 86 were up-regulated (hypo-up genes) and 237 were down-regulated (hypo-down genes). These genes were then filtered in the MeRIP-seq and RNA-seq data using p < 0.05 as the threshold, and 15 hyper-up, 9 hyper-down, 2 hypo-up, and 58 hypo-down genes with significant differences were obtained (Tables 3-6 and S17). Notably, some genes associated with embryonic development and muscle growth were significantly different in terms of both expression and m6A modification data. The expression and m6A abundance of FBXO40 in E19 samples were significantly higher than in E13 samples. The m6A abundance of iroquois homeobox 5 (IRX5) was significantly elevated, but mRNA levels were significantly reduced ( Figure 5C). Pearson analysis of the correlation between m6A modification and expression regulation, using m6A abundance and mRNA expression data, revealed a negative correlation between m6A abundance and the mRNA expression level (E13, R = −0.43680, p < 0.01; E19, R = −0.45243, p < 0.01), indicating that m6A modification affected gene expression ( Figure 5D).

Expression Patterns of the m6A Modification Regulators during Duck Embryonic Skeletal Muscle Development
Motif sequence analysis showed that m6A modification in duck embryo skeletal muscles may be mediated by METTL14. qRT-PCR analysis of METTL14 expression during the embryonic skeletal muscle development of duck breast and leg muscles from E10, E13, E16, E19, E22, and 1 day post hatch (P1, 5 samples in each stage) revealed significant changes from E10 to E22 in breast muscles, and marked downregulation from E10 to E19 in leg muscles ( Figure 6A,B). In addition, we evaluated the expression levels of METTL14 cofactors, including WTAP, ZC3H13, RBM15, and VIRMA during duck embryonic skeletal muscle development (as above). ZC3H13, RBM15, and VIRMA showed similar downward expression patterns to METTL14 during embryonic skeletal muscle development ( Figure 6E-J). Meanwhile, the expression levels of WTAP peaked on E13 during the leg muscle development ( Figure 6D), whereas its expression in breast muscle samples showed a downward tendency ( Figure 6C). Moreover, we performed qRT-PCR validation for the RNA-seq and the results were consistent with the RNA-seq data ( Figure 7). These data showed that METTL14 and the RNA methyltransferase cofactors' expression were drastically downregulated from pre-differentiation (before E13) to post-differentiation (after E19) in duck skeletal muscles. We speculate that METTL14-mediated m6A modification plays a crucial role in embryonic skeletal muscle development. seq data in E19 and E13. (C) The m6A abundances in IRX5 mRNA for the E19 and E13 groups. The m6A peak, which was significantly higher in E19 than in E13, is shown in the black rectangle (p < 0.05 and |log2FC| > 1). (D) A plot of m6A peak abundance and mRNA expression in E19 and E13.   presented as the mean ± S.E.M. Different letters between two groups represent differences that are significant (p < 0.05).

Discussion
Muscle tissue drives body movement, and animal muscle tissue is a source of high quality protein [27,28]. Muscle development is regulated by many genetic factors, such as DNA methylation genes and non-coding RNAs [29,30]. Epigenetic RNA modification, especially that of m6A, influences various biological processes [9]. Although m6A has been identified in various animals, including pigs, chickens, geese, sheep, cattle, and fish [18][19][20][21][22]25,31,32], it has not been described in ducks, and its role in poultry muscle development is unknown. Here, we selected two important stages of duck embryo skeletal muscle differentiation for MeRIP-seq and found differences in m6A modification patterns before and after skeletal muscle differentiation. Together, our data suggest that m6A methylation plays an important role in duck skeletal muscle development by regulating gene expression.
We found that the m6A modification pattern changed after skeletal muscle differentiation, with a reduction of 1123 m6A peaks and 369 peak genes. More than 73% of the transcripts had one or two m6A peaks and about 13% had ≥ 4 peaks, which is higher than the rates seen in humans, pigs, and chickens [22,33,34]. Motif analysis identified the classical GGACU sequence as the one modified by m6A in duck embryonic skeletal muscle, and METTL3 and METTL14 as potential methyltransferases [11]. However, duck METTL3 has not been annotated and it is possible that METTL14 mediates m6A modification in duck embryo skeletal muscles. The m6A distribution peaks in duck breast muscle mainly occur in the region before and after the stop codon. This region is proposed to possess transport and translocation features, as well as protein translation regulatory elements, which affect RNA stability, transport, and translocation signal transduction [35]. m6A modification in this region may affect the above regulatory functions by changing RNA conformation. The m6A peak distribution was similar to that observed in pig longissimus dorsi muscles [36], but is different from that seen in chicken abdominal fat tissue [22]. In chicken fat tissue, m6A peaks are mainly distributed before and after the start codon, in the CDS region, and before and after the stop codon. m6A peak distribution patterns may be tissue specific.
Previous studies have shown that different stages of the same tissue have different m6A modification patterns [19,37]. In the E13 and E19 stages of breast muscle tissues, 78.9% of the genes were modified by m6A in both stages, and modification of 21.1% of the

Discussion
Muscle tissue drives body movement, and animal muscle tissue is a source of high quality protein [27,28]. Muscle development is regulated by many genetic factors, such as DNA methylation genes and non-coding RNAs [29,30]. Epigenetic RNA modification, especially that of m6A, influences various biological processes [9]. Although m6A has been identified in various animals, including pigs, chickens, geese, sheep, cattle, and fish [18][19][20][21][22]25,31,32], it has not been described in ducks, and its role in poultry muscle development is unknown. Here, we selected two important stages of duck embryo skeletal muscle differentiation for MeRIP-seq and found differences in m6A modification patterns before and after skeletal muscle differentiation. Together, our data suggest that m6A methylation plays an important role in duck skeletal muscle development by regulating gene expression.
We found that the m6A modification pattern changed after skeletal muscle differentiation, with a reduction of 1123 m6A peaks and 369 peak genes. More than 73% of the transcripts had one or two m6A peaks and about 13% had ≥4 peaks, which is higher than the rates seen in humans, pigs, and chickens [22,33,34]. Motif analysis identified the classical GGACU sequence as the one modified by m6A in duck embryonic skeletal muscle, and METTL3 and METTL14 as potential methyltransferases [11]. However, duck METTL3 has not been annotated and it is possible that METTL14 mediates m6A modification in duck embryo skeletal muscles. The m6A distribution peaks in duck breast muscle mainly occur in the region before and after the stop codon. This region is proposed to possess transport and translocation features, as well as protein translation regulatory elements, which affect RNA stability, transport, and translocation signal transduction [35]. m6A modification in this region may affect the above regulatory functions by changing RNA conformation. The m6A peak distribution was similar to that observed in pig longissimus dorsi muscles [36], but is different from that seen in chicken abdominal fat tissue [22]. In chicken fat tissue, m6A peaks are mainly distributed before and after the start codon, in the CDS region, and before and after the stop codon. m6A peak distribution patterns may be tissue specific.
Previous studies have shown that different stages of the same tissue have different m6A modification patterns [19,37]. In the E13 and E19 stages of breast muscle tissues, 78.9% of the genes were modified by m6A in both stages, and modification of 21.1% of the genes was stage specific. We speculate that group-specific genes may play important roles in skeletal muscle development. E19-specific m6A genes were significantly enriched in the MAPK signaling pathway. MAPKs are part of a highly conserved network, and the kinases of MAPKs, Mitogen-Activated Protein Kinase 1 (MAPK1), and MAPK3 induce slow fiber-type switching in skeletal muscles [38]. E13 m6A peak-containing genes were mainly enriched in embryo-development-and cell-proliferation-related signaling pathways, including pathways that regulate stem cell pluripotency and cAMP signaling [39,40]. Intersection m6A-containing genes between E19 and E13 were involved in RNA splicing and embryo development. Modifications of the RNA processing genes may affect gene expression, globally affecting mRNA splicing and protein synthesis [41].
Here, 355 differentially methylated peaks (DMPs) and 331 differentially methylated genes (DMGs) were identified. Enrichment analysis revealed that the hyper-DMGs were significantly enriched in AMPK signaling. AMPK induces mitochondrial biogenesis and skeletal muscle glucose uptake during training adaptation, and AMPK signaling mediates muscle-fiber-type transformation [42,43]. Protein-protein interaction network analysis indicated that the DMGs, including FBXO40, SOX2, and DLX5, are involved in skeletal muscle development. Interestingly, the FBXO40, SOX2, and DLX5 genes were also differentially expressed, and the differential abundance of m6A is consistent with this expression pattern (Table S13). FBXO40, a member of the F-box protein family, ubiquitinates Insulin Receptor Substrate 1 (IRS1) for degradation, blocks IGF1 signaling, and causes muscle atrophy [44,45]. MeRIP-seq and RNA-seq analyses revealed that, among the m6A peak-associated genes and the hyper-methylated peak genes in E19, many upregulated genes, including PDK4, MYLK2, and FBXO40, are involved in muscle differentiation and growth. Forkhead box O1 (FOXO1) and peroxisome proliferator-activated receptor (PPAR) α/β regulate the PDK4 response to insulin, elevate free fatty acids or hunger-mediated signal transduction, and mediate muscle fiber atrophy [46,47]. MYLK2 is mainly expressed in striated muscle and mediates the switch of fast and slow muscle fibers in skeletal muscles [48]. We speculate that m6A modification may affect duck muscle development by modulating PDK4, MYLK2 and FBXO40 expression.
The methyltransferases METTL3 and METTL14 catalyze RNA m6A modification using S-adenosylmethionine as the methylation donor [10,11]. METTL14 knockout inhibits the differentiation of C2C12 myoblasts and promotes their proliferation [49]. Motif analysis showed that m6A modification in duck embryo skeletal muscle may be mediated by METTL14. The METTL14 expression level was high before duck embryonic skeletal muscle differentiation (E10 to E16) and decreased significantly between E16 to E19, and remained low at E19 to P1, indicating that METTL14 influences duck muscle differentiation. We speculate that METTL14 mediates m6A modification in duck embryonic skeletal muscles and regulates muscle development. The underlying mechanism of METTL14 regulation of skeletal muscle development and the m6A genes induced by METTL14 need to be further explored.

Conclusions
In summary, we profiled the m6A modification pattern in duck embryonic skeletal muscle, presented the difference between the pre-and post-differentiation m6A patterns in breast muscle, and proposed that m6A modification influences muscle differentiation by regulating gene expression. The current study not only explored the profiling effects of RNA methylation on duck skeletal muscle development, but also laid the groundwork for studying RNA modification in poultry muscle development.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ani12192593/s1, Figure S1: The correlation coefficient of MeRIPseq samples; Figure S2: The proportion of 1-2 m6A peaks in each group; Table S1: Information of the qPCR primers; Table S2: the list of intersection peak genes in E13 and E19; Tables S3 and S4: the group-specific peak genes in E13 and E19, respectively; Table S5: the GO biological process analysis with all intersection genes in E13 and E19; Table S6: the KEGG analysis with all intersection genes in E13 and E19; Tables S7 and S8: the KEGG analysis with group-specific peak genes in E13 and E19, respectively; Tables S9 and S10: the differentially methylated peaks (DMPs) and the DMPs associated genes in E19 vs. E13; Tables S11 and S12: the KEGG analysis with the hyper-methylated peak genes and hypo-methylated peak genes in E19 vs. E13, respectively; Table S13: the differentially expressed genes (DEGs) in E19 vs. E13; Table S14: the GO biological process analysis with all DEGs in E19 vs. E13; Table S15: the KEGG analysis with all DEGs in E19 vs. E13; Table S16: the conjoint analysis of group-specific peak genes and DEGs in E19 vs. E13; Table S17: List of hypo-down methylated genes that were significantly different between E19 and E13.