Regulatory role of N6-Methyladenosine on skeletal muscle development in Hu sheep

N6-Methyladenosine (m6A) RNA modification plays an essential role in many biological processes. To investigate the regulatory role of m6A on the skeletal muscle development in Hu sheep, this study took newborn Hu sheep (b_B Group) and six-month-old Hu sheep (s_B Group) as the objects. MeRIP-seq and RNA-Seq analysis techniques were used to detect differentially methylated genes (DMGs) and differentially expressed genes (DEGs) in the longissimus dorsi muscle of Hu sheep at different months of age. Then, conjoint analysis was further employed to screen for key genes involved in skeletal muscle development that are modified by m6A and expressed by mRNA. According to the results of the MeRIP-seq analysis, there were 285 m6A differentially methylated peaks (DMPs) in total between b_B Group and s_B Group, with 192 significant upregulated peaks and 93 significant downregulated peaks. GO and KEGG analysis revealed that DMGs are mainly enriched in actin-binding, cellular transport, and metabolic pathways. According to the results of the RNA-seq analysis, there were 4,349 DEGs in total between b_B Group and s_B Group, with 2010 upregulated genes and 2,339 downregulated genes. DEGs are found to be mainly enriched in the regulation of actin cytoskeleton tissue, AMPK and FoxO signaling pathways, etc. The conjoint analysis demonstrated that 283 genes were both modified by m6A and expressed by mRNA. Among them, three genes relevant to muscle growth (RGMB, MAPK8IP3, and RSPO3) were selected as candidates for quantitative validation, and the results were in line with the sequencing results. The results mentioned above all suggest that m6A plays a certain role in the skeletal muscle development in Hu sheep.

1 Introduction N6-Methyladenosine (m6A) refers to the methylation modification that occurs on the sixth nitrogen atom of adenine (Wang et al., 2022).m6A methylation is the most common chemical modification of mRNA in eukaryotes, mediated by methyltransferase complexes centered around methyltransferaselike3 (METTL3) and methyltransferaselike14 (METTL14) (Meyer et al., 2012;Liu et al., 2014;Sorci et al., 2018).RNA methylation modification is a regulatory mechanism that controls gene expression in eukaryotic cells (Huang et al., 2020).As a reversible epigenetic modification, m6A is found not only in messenger RNAs but also in non-coding RNAs, which has an impact on the formation of modified RNA molecules and plays an important role in almost all popular biological processes (Huang et al., 2020).
It was found that m6A modification has potential regulatory effects on the growth and development of animal muscle.Deng et al. (2021) performed MeRIP-seq analysis on the skeletal muscle of Haimen goats and resolved the characteristics of m6A modification during muscle development.Zhang et al. (2020) found that m6A methylation and IGF2BP1 play important roles in the regulation of prenatal myogenesis in pig embryos.Ma et al. (2022) found that the m6A gene is mainly involved in regulating longissimus dorsi muscle (LD) differentiation and development in yaks.Chen et al. (2022) found that m6A modification has an impact on ducks' muscle differentiation by regulating gene expression.Xu et al. (2021) uncovered a transcriptome-wide m6A modification pattern that affects embryonic breast muscle development in Dingan goose.As a valuable native sheep breed in China, Hu sheep are characterized by rapid growth, high reproductive rate, tender, and juicy meat (Zhao et al., 2024).There has been no previous research on the regulatory role of m6A in skeletal muscle development in Hu sheep, this study should be conducted, in order to study the role of m6A modification in Hu sheep and provide a theoretical regulatory mechanism for the growth and development of Hu sheep.

Sample collection
Six healthy Hu sheep were selected for this study, three in the newborn (0 day of age) group (b_B Group) and three in the sixmonth-old (180 days of age) group (s_B Group), all of which were female individuals.The Hu sheep used in the experiment were supplied by Huzhou Yihui Eco-Agriculture Co., Ltd.(Zhejiang, China), whose feeding criteria were environmentally friendly.The test sheep's longissimus dorsi muscles were extracted after slaughter and promptly kept in liquid nitrogen before being stored at −80 °C for future studies.All experimental techniques followed the rules established by the Experimental Animal Management Committee of the Zhejiang Academy of Agricultural Sciences.

RNA extraction and fragmentation
Total RNA was extracted and purified with TRIzol reagent (Invitrogen, United States).10 ug of total RNA was mixed with RNA Fragmentation Reagents (Invitrogen, United States) for 10 min of reaction at 70 °C in a Thermomixer to break the RNA into fragments with a size of approximately 100 nt.The fragmented RNA was precipitated using the ethanol technique.

m6A enrichment
The magnetic beads containing protein A and protein G (Invitrogen, United States) were first cleaned with IP buffer (150 mM NaCl, 10 mM Tris-HCl, pH = 7.5) and then treated for 2 h at 4 °C with 5 ug of m6A antibody (Millipore).Second, the magnetic beads were washed twice and resuspended using IP buffer for flipping at 4 °C for 4 h, with fragmented RNA added.Next, incubate the magnetic beads for 1 hour at 4 °C with m6A competitive eluent after washing the magnetic beads three times with IP buffer.The supernatant containing the eluted m6A RNA was transferred to a new test tube for purification using a reagent of phenol, chloroform, and isoamyl alcohol with a ratio of 125:24:1.

Library preparation
Reverse transcription and library preparation on IP and Input samples were carried out using SMARTer ® Stranded Total RNA-seq Kit v2-Pico Input Mammalian User Manual (Takara, JPN).To get the final library, AMPure XP bead [SpeedBead lagnetic Carboxylate lodified Particles (GE, United States)] was employed to choose fragment sizes.After successfully passing the library detection, library pooling was conducted in accordance with the intended sequencing data volume and effective concentration requirements.The Illumina Nova platform was utilized to sequence the libraries, with a strategy of PE150.

Quantitative real-time PCR
Three genes (RGMB, MAPK8IP3, RSPO3) were screened via the conjoint analysis based on differential m6A-associated genes and DEGs for qRT-PCR analysis.Reverse transcription was conducted on the total RNA using Polestar first cDNA Synthesis Kit (Tiosbio, China).Real-time quantitative PCR was performed using SYBR ® Green Realtime PCR Master Mix-Plus (TOYOBO, Japan).The primers were designed using DNAMAN8 with the GAPDH gene of sheep as the internal reference gene (Vorachek et al., 2013).The primer sequences are shown in Supplementary Table S3.Three biological replicates were performed for each sample, and the relative expression was calculated by the 2 −△Ct method (Alcaraz-López et al., 2020).

Data statistics and audit
This study preprocessed the raw data, removing redundant information such as joint sequences and low-quality bases by using Trimmomatic.As shown in Supplementary Table S1, the filtered reads of the six IP libraries were 102,111,130,626.The filtered reads of the Input libraries were 117,937,062-149,706,464.

Data analysis
The filtered data were compared to the reference genome using hisat2 (Kim et al., 2015).It was found that the proportion of read pairs (Uniq_Rate) of all samples that were accurately matched to a position in the reference genome was more than 54.40%, and the results are shown in Supplementary Table S2.
The distribution of reference genome comparison regions showed that the percentage contents of sequenced sequences localized to exon regions were the highest as shown in Figure 1.

Identification of m6A modification sites and analysis of DMGs
The length and location information of the peaks on the genome were determined using the MeRIP-seq data.Reads were found to be abundant near the transcriptome start site (TSS) of the gene.The distribution of reads in the combinable region near the TSS was shown in the form of heatmap plot (as shown in Figure 2A).There were 285 m6A DMPs in total between b_B Group and s_B Group, with 192 being significant upregulated and 93 being significant downregulated (as shown in Figure 2B).The peaks in b_B Group and s_B Group were primarily enriched in the 3′UTR (as shown in Figure 2C).exomePeak analysis found 3,164 and 3,440 specific peaks in b_B Group and s_B Group, respectively, with 11,918 similar peaks between the two groups (as shown in Figure 2D).The peak distributions of the two groups were similar (as shown in Figures 2E, F).The analysis findings revealed that differential peaks were primarily distributed in the 3′UTR (as shown in Figure 2G).Figure 2H demonstrates that differential peaks were found on each chromosome.Table 1 shows the top 20 m6A peaks, with 15 upregulated and 5 downregulated.

Motif analysis
OMER (Hansen et al., 2016) was utilized to identify the motifs in the m6A-modified regions.The motifs were ranked according to the value of P, with the smaller P ranking higher.The results of motif analysis are shown in Figure 3. RRACH, a common motif structure in RNA modification, was present in samples from both groups (where R = A or G; H = A, C or U).

Enrichment analysis of DMGs
GO and KEGG enrichment analyses were performed on DMGs to analyze their potential functions in the skeletal muscle development of Hu sheep.The enriched GO terms for the DMGs mainly include protein modification process, regulation of macromolecule metabolic process, actin-binding, actin cytoskeleton organization, and structural constituent of muscle.KEGG pathway enrichment analysis showed that DMGs were enriched to the cellular transport, metabolism, autophagy, processing of genetic information, ubiquitinmediated protective lysis, glycan biosynthesis and metabolism, as shown in Figure 4.

Analysis of DEGs
The gene expression amount and density were demonstrated in Figures 5A, B. There were 4,349 DEGs detected in total between b_B Group and s_B Group, with 2010 upregulated genes and 2,339 downregulated genes (as shown in Figure 5D).The distribution of DEGs was demonstrated by MA plots, volcano plots, and heat maps, respectively (as shown in Figures 5C-E).Table 2 shows the top 20 DEGs, with 10 upregulated and 10 downregulated.
GO and KEGG enrichment analyses were conducted to further investigate the functions of DEGs.The enriched GO terms for the DEGs mainly include regulation of muscle cell differentiation, actin cytoskeleton, skeletal muscle tissue growth, and carbohydrate metabolism processes (as shown in Figures 6A, B).KEGG pathway enrichment analysis demonstrated that DEGs showed a significant enrichment to the signaling pathways related to muscle development, including the cAMP signaling pathway, AMPK signaling pathway, FoxO signaling pathway, and JAK-STAT signaling pathway, as illustrated in Figures 6C, D. AMPK can be involved in the metabolic regulation of skeletal muscle via activating its downstream target proteins.It plays an important role in regulating skeletal muscle development through its effects on the cellular anabolism and catabolism processes (Thomson, 2018).FoxO signaling pathway is an important pathway involved in skeletal muscle atrophy (García-Prat et al., 2020).It can be suggested that these DEGs may play a crucial role in the skeletal muscle development of Hu sheep.

Conjoint analysis of differential m6A modification and DEGs
The overlapping relationship between differential m6A modification-associated genes and DEGs was analyzed based on the two sequencing results (as shown in Figure 7).It was found that there were 283 genes with m6A methylation modification and significant differential expression, among which there were 11 genes with "m6A_up" and "mRNA_up, 20 genes with m6A_down" and "mRNA_up, 3 genes with m6A_down", and "mRNA_down genes, and 54 genes with m6A_up" and "mRNA_down."

Validation of DEGs
With the purpose of examining the regulation of skeletal muscle development, three candidate genes, RGMB, MAPK8IP3, and RSPO3, were selected for qRT-PCR validation from relevant genes that changed at both m6A and mRNA levels.These three genes were downregulated in m6A modification but upregulated in mRNA expression (Table 3).
The qRT-PCR results revealed that the mRNA expression of the three genes in the b_B Group was significantly higher than that of the s_B Group (Figure 8).The changing trends of the three genes were consistent with the RNA-seq results.Motifs of m6A-modified regions in b_B Group and s_B Group (A) The first six motifs of b_B Group; (B) The first six motifs of s_B Group.
Knowledge of skeletal muscle development is essential to unravel the molecular mechanisms of skeletal muscle formation and disease (Buckingham, 2001;Al-Qusairi and Laporte, 2011).In this study, 285 differentially methylated peaks were found in the skeletal muscles of Hu sheep with different months of age, and the conjoint analysis identified 283 genes with m6A methylation modification and significant differential expression.Dou et al. (2023) found 613 differentially methylated peaks between the QA group (3 days of age) and QN (270 days of age) group of Queshan black pigs when investigating their longissimus dorsi muscle.They also found 88 genes with m6A methylation modification and significant differential expression through conjoint analysis.The results show that IGF1R, CCND, MYOD1, FOS, PHKB, BIN1, and FUT2 were involved in the pathways related to muscle development, may play a regulatory role in the process of muscle growth and development.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.Chen et al. (2022) found 355 differentially methylated peaks in the skeletal muscles of Mountain ducks between the E13 group (Embryo Day 13) and E19 group (Embryo Day 19).They also found 84 genes with m6A methylation modification and significant differential expression through conjoint analysis.They speculate that m6A modification may affect duck muscle development by modulating PDK4, MYLK2 and FBXO40 expression.Xu et al. (2021) found 418 differentially methylated peaks in the skeletal muscles of Dingan goose between the E21 group (Embryo Day 21) and E30 group (Embryo Day 30).Among the m6A-miRNA-genes, they found 10 genes (PDK3, PITPNA, DSTN, BACE, GATM, ITM2A, SOD2, IGFBP4, GAA and TUBB6) are related to breast muscle development which were tightly associated with breast muscle development through affecting the m6A modification levels of their target gene.Wang et al. (2024) found 6,476 differentially methylated peaks in the representatives of leaner Xinghua chicken (XH) and hypertrophic White Recessive Rock chicken (WRR) broilers.They also found 167 genes with m6A methylation modification and significant differential expression through conjoint analysis.As a demethylase of m6A, the highly expression of ALKBH5 in the muscle tissue   of poultry and differential expression between XH and WRR chickens suggest that ALKBH5 may play a crucial role in regulating muscle development.In this study, three candidate genes, RGMB, MAPK8IP3, and RSPO3, were selected for qRT-PCR validation from relevant genes that changed at both m6A and mRNA levels.Results revealed that the mRNA expression of the three genes in the b_B Group was significantly higher than that of the s_B Group, and the changing trends of the three genes were consistent with the RNA-seq results.These results suggested that that RGMB, MAPK8IP3, and RSPO3 play important roles in the skeletal muscle development of Hu sheep.RGMB is the member b of the repulsive guidance molecule (RGM) family (Sartori and Sandri, 2015).Bonemorphogenetic protein (BMP) is a secreted signaling factor belonging to the transforming growth factor β  qRT-PCR results of three differential genes in b_B Group and s_ B Group.
Frontiers in Genetics frontiersin.org10 (TGFβ) superfamily, which plays an important role in bone and cartilage formation (Sartori and Sandri, 2015).RGMB regulates a variety of physiological processes mainly through neogenin-Rho and BMP signaling pathways, and it can promote the activation of the signaling pathways (Samad et al., 2005).MAPK8IP3 (mitogen-activated protein kinase 8 interacting protein 3), also known as JIP3 (JNK-interacting protein3), is a JNK (c-Jun N-terminal kinase) scaffold protein (Akechi et al., 2001;Iwasawa et al., 2019).JIP3 generates a complex with JNK and its upstream signaling kinases that promotes the transduction of the JNK signaling pathway (Kelkar et al., 2000;Matsuura et al., 2002;Song and Lee, 2007).Knockdown of JIP3 significantly reduced IR-induced JNK signaling activation and apoptosis (Xu et al., 2010).R-spondin 3 (RSPO3) is a secreted protein belonging to the RSPO family (RSPO1-4) with a common structure of thrombospondin type 1 domain and N-terminal cysteine-rich region.The RSPO family is expressed in normal human placental, lung and muscle tissues (Kim et al., 2006;Kim et al., 2008), which regulates cell proliferation and differentiation by activating the Wnt signaling pathway and is critical for the growth of bone, muscle, blood vessels, and other tissues (Gu et al., 2020).In sheep, there is no report about regulatory role of N6-Methyladenosine in skeletal muscle Development.This study found that m6A plays a certain role in the skeletal muscle development in Hu sheep and the genes RSPO3, RGMB, MAPK8IP33 play important roles in the skeletal muscle development and differentiation.Subsequently, the three candidate genes RSPO3, RGMB, and MAPK8IP33 can be knocked out, knocked down, overexpressed, etc.Finally, Western blot, qPCR, RNA pull down and other detection and verification are carried out to study role and mechanism.

Conclusion
The findings of this study indicate that the genes RSPO3, RGMB, and MAPK8IP33 may play important roles in the skeletal muscle development and differentiation of Hu sheep.It is inferred that m6A modification has a critical impact on the skeletal muscle development in Hu sheep.
FIGURE 1Comparison region distribution of reference genome.

FIGURE 2
FIGURE 2 Analysis of m6A modification in the longissimus dorsi muscle of Hu sheep (A) Heat map of the enrichment degree of reads near TSS; (B) Visualized volcano plot of differentially m6A-modified regions; (C) Distribution density of differential m6A modifications across gene elements; (D) Veen plots of peak-annotated genomes; (E,F) m6A peak distribution of b_B Group and s_B Group; (G) Differentially m6A peak distribution of b_B Group and s_B Group; (H) Distribution of differential peaks on chromosomes.

FIGURE 5
FIGURE 5 Identification of DEGs (A) Gene expression box plot; (B) Gene expression density plot for DEGs; (C) MA plot for DEGs; (D) Volcano plot for DEGs; (E) Heat map for DEGs.
GO and KEGG enrichment analysis of DEGs function (A,B) DEGs GO enrichment analysis; (C,D) DEGs KEGG enrichment analysis.

FIGURE 7 Correlation
FIGURE 7Correlation Analysis of Differential m6A Modification and DEGs.

TABLE 3 m6A
-modified genes related to skeletal muscle development.