Skip to main content
Download PDF

Comprehensive characterization of lncRNA N6-methyladenosine modification dynamics throughout bovine skeletal muscle development

Journal of Animal Science and Biotechnology volume 16, Article number: 36 (2025) Cite this article

Abstract

Background

N6-methyladenosine (m6A) methylation is a key epigenetic modification that can modulate gene expression and strongly affect mammalian developmental processes. However, the genome-wide methylation of long non-coding RNAs (lncRNAs) and its implications for the development of skeletal muscle remain poorly understood. Bovine skeletal muscle samples from five developmental stages were analyzed in this study to establish lncRNA methylome and transcriptomic maps.

Results

Globally, 59.67% of lncRNAs in skeletal muscle with m6A modifications, and this percentage decreased progressively during development. lncRNA expression levels were positively associated with the number of m6A peaks, with lncRNAs possessing 3 or more peaks showing significantly higher expression levels than those with 1 or 2 peaks. Specific lncRNAs involved in skeletal muscle development were identified through two analytical approaches. The first approach employed weighted gene co-expression network analysis (WGCNA) of transcriptomic data to identify correlations between annotated lncRNAs and growth-related traits, resulting in 21 candidate hub lncRNAs. The intersection of these 21 hub lncRNAs with 151 differentially methylated lncRNAs (DM-lncRNAs) identified 10 shared candidate lncRNAs. The second approach integrated MeRIP-seq and RNA-seq data to identify 36 lncRNAs that were both differentially m6A modified and differentially expressed (dme-lncRNAs). GO and KEGG enrichment analyses of cis-target genes associated with these dme-lncRNAs identified eight candidate lncRNAs. Combining the results from the two approaches identified 16 key m6A-modified lncRNAs likely involved in skeletal muscle development.

Conclusions

These findings highlight the regulatory and functional significance of dynamic lncRNA methylation in skeletal muscle development.

Background

In addition to being the most abundant tissue type in the mammalian body, skeletal muscle growth and development are the primary determinants of meat yield and meat quality for livestock [1, 2]. Studies of skeletal muscle development are thus an extremely important topic in the animal genetics and breeding field. Skeletal muscle growth involves the coordinated regulation of processes such as myoblast proliferation and differentiation across various stages of development [3]. Skeletal muscle of the trunk and limbs is of mesodermal origin, developing from dorsal portions of somites arising from mesoderm segmentation [3, 4]. This process involves the regulatory control of many transcription factors, signaling molecules, long non-coding RNAs (lncRNAs), and N6-methyladenosine (m6A) methylation residues. The influence of novel regulatory factors, including m6A methylation, on this process, remains partially understood, requiring further research to improve livestock meat quality through genetic improvements.

lncRNAs, a large class of RNAs over 200 nucleotides long, generally lack protein-coding potential but often exhibit poly-A tailing and splicing. They are abundant in the cytosol and nuclear compartments and have complex spatial structures and a range of functions [5, 6]. Specific lncRNAs regulate various biological processes, including development and disease progression. Altered lncRNA expression is associated with muscle-related diseases and plays a significant role in regulating muscle tissue development [7]. lncRNAs influence biological processes by competitively binding molecules, altering target protein stability, regulating translation, acting as molecular sponges, activating transcription or modifying chromosomes [7]. The competing endogenous RNA (ceRNA) mechanism is thought to be the primary pathway by which lncRNAs influence muscle development through the sequestration of miRNAs with sequence complementarity, thereby altering target mRNA expression to shape developmental processes [8,9,10,11,12]. IGF2-AS [13] and lncMGPF [14] are established ceRNAs that respectively bind to miR-503 and miR-135a-5p, emphasizing their importance in the context of post-transcriptional regulation. However, many lncRNAs related to the development of muscle tissues in domestic animals remain to be characterized.

m6A modifications involve the methylation of the amino group at the 6th position of adenine residues in RNA molecules. This reversible epigenetic modification is regulated in a coordinated manner by methyltransferases (writers) and demethylases (erasers) in the nucleus, thereby controlling the overall levels of m6A modification and the specific sites that are modified. The m6A reader proteins can then recognize these modifications to control the splicing, transport, translation and degradation of the modified RNAs [15,16,17,18]. m6A modification is catalyzed by methyltransferase complexes, primarily METTL3, METTL14 and WTAP. FTO and ALKBH5 are m6A demethylases capable of removing the m6A modifications from RNA molecules that have been methylated [19]. The m6A reader proteins are generally classified into three groups, including (1) YTH domain-containing readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2 and YTHDF3) capable of recognizing transcripts that have undergone m6A modification; (2) the 43S translation initiation complex protein eIF3 (eukaryotic initiation factor 3), which modulates translational initiation by binding to the 5′ untranslated region (UTR) of mRNAs that have undergone m6A modification [15]; and (3) factors capable of stabilizing mRNAs that have undergone m6A modification, including insulin-like growth factor binding proteins 1–3 (IGF2BP1–3) and Prrc2a [20, 21]. A growing body of evidence has attested to the important role that m6A modifications play in the coordination of the development of skeletal muscle. In pigs, for instance, whole-transcriptome m6A mapping of the prenatal skeletal muscle development process highlighted dynamic shifts in the m6A methylome over time, with most impacted genes being closely associated with pathways relevant to skeletal muscle development [22]. In goats, knocking down FTO is associated with an increase in GADD45B m6A levels and a reduction in GADD45B mRNA stability during muscle development, ultimately suppressing myogenic differentiation [23]. One recent report highlighted the effects of m6A modification during bovine muscle cell development, and METTL3 and METTL14 silencing were shown to improve bovine myoblast proliferation while inhibiting the differentiation of these cells [24].

Bohai black cattle, previously referred to as Wudi black cattle, are one of the three most prominent breeds of black cattle globally and are members of the yellow cattle family. In this study, the dynamic profiling of the lncRNA methylome and transcriptome was performed in bovine skeletal muscle samples collected from five stages of development ranging from the newborn stage to 30 months of age. The potential biological importance of lncRNA methylation during development of skeletal muscle was assessed. These findings provide insights into the molecular processes governing ruminants skeletal muscle development, laying a theoretical foundation for future studies on the specific mechanisms that regulate skeletal muscle development and advancing the selection and breeding of Bohai black cattle.

Materials and methods

Tissue samples

Samples of skeletal muscle (longissimus dorsi muscle, LDM) were collected from Bohai black cattle at 5 stages of development, including 0, 6, 12, 20 and 30 months postnatally (M0, M6, M12, M20 and M30). At each of these stages, samples were collected from three bulls as biological replicates and stored in liquid nitrogen.

MeRIP-seq and RNA-seq library construction

Total RNA was extracted from 15 samples using TRIzol (Cat. 15596-026, Invitrogen, Carlsbad, CA, USA) as directed. A NanoDrop spectrophotometer (NanoDrop ND-1000, NanoDrop Technologies, Wilmington, DE, USA) was used to quantify RNA and check purity. A bioanalyzer (Agilent 2100 Bioanalyzer, Agilent Technologies Inc., Santa Clara, CA, USA) was then used to establish the quality of the isolated RNA, retaining samples with a RIN < 7.0 and confirming these results through denaturing agarose gel electrophoresis. More than 25 μg of total RNA from a particular sample was used to deplete rRNA with an Epicentre Ribo-Zero Gold Kit (Cat. RZH1046, Epicentre, Madison, WI, USA). Purified rRNA-depleted RNA was then fragmented with a Magnesium RNA Fragmentation Module (Cat. E6150, New England Biolabs, Ipswich, MA, USA) at 86 °C for 7 min and incubated with an m6A-specific antibody (Cat. 202003, Synaptic Systems, Göttingen, Lower Saxony, GER) in IP buffer (50 mmol/L Tris-HCl, 750 mmol/L NaCl and 0.5% Igepal CA-630) for immunoprecipitation. The precipitated RNA was used to prepare cDNA with SuperScript™ II Reverse Transcriptase (Cat. 1896649, Invitrogen, Carlsbad, CA, USA), which was then used for the preparation of U-labeled second-stranded DNA using E. coli DNA polymerase I (Cat. M0209, New England Biolabs, Ipswich, MA, USA), RNase H (Cat. M0297, New England Biolabs, Ipswich, MA, USA) and a dUTP Solution (Cat. R0133, Thermo Fisher Scientific, Waltham, MA, USA). The blunt ends of each strand then had an A-base added to prepare for index adapter ligation, with each adapter harboring a T-base overhand to allow for ligation of the A-tailed DNA fragments. These fragments were ligated to single- or dual-index adapters, after which AMPureXP beads were used for size selection. U-labeled second-stranded DNA was then treated with the heat-labile UDG enzyme (Cat. M0280, New England Biolabs, Ipswich, MA, USA), and ligated products were subjected to PCR amplification (95 °C for 3 min; 8 cycles of 98 °C for 15 s, 60 °C for 15 s, and 72 °C for 30 s; 72 °C for 5 min). The final cDNA library had an average insert size of 300 ± 50 bp. At last, 2 × 150 bp paired-end sequencing (PE150) of these samples was performed using an Illumina Novaseq™ 6000 (LC-Bio Technology Co., Ltd., Hangzhou, China).

Bioinformatics analyses

Adaptor-containing, low-quality and undetermined reads were removed using Fastp (https://github.com/OpenGene/fastp) with default settings [25]. FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and RseQC (http://rseqc.sourceforge.net/) [26, 27] were used for the sequence validation of IP and input samples, while the read mapping to the Bos taurus reference genome (ARS-UCD1.2, Ensemblv107) was performed with HISAT2 (http://daehwankimlab.github.io/hisat2) [28]. The R exomePeak package (http://bioconductor.jp/packages/3.17/bioc/html/exomePeak2.html) was used for peak calling and differential peak analyses, with peak annotation being achieved based on the overlap with gene architecture using the R ANNOVAR package (https://annovar.openbioinformatics.org/en/latest/) [29, 30]. De novo and known motifs were identified with HOMER (http://homer.ucsd.edu/homer/motif), followed by motif localization relative to the peak summit [31]. Transcript and gene expression analyses were conducted using StringTie (https://ccb.jhu.edu/software/stringtie), calculating FPKM values (total exon fragments/mapped reads (millions) × exon length (kb)) to quantify expression levels [32]. lncRNAs considered expressed such that they were retained for subsequent analysis were those with an average FPKM ≥ 0.1 in three samples. The threshold for differential transcript expression was |log2 (fold change)| ≥ 0.585 and P < 0.05, as determined with the R edgeR package (https://bioconductor.org/packages/release/bioc/html/edgeR.html) [33].

lncRNA identification

Unqualified sequences were filtered out Cutadapt (https://cutadapt.readthedocs.io/en/stable/) [34], after which HISAT2 (https://daehwankimlab.github.io/hisat2/) [28, 35, 36] was used to remove portions of reads with low quality and to map the remaining reads to the reference genome. StringTie (https://ccb.jhu.edu/software/stringtie) [32, 36, 37] was then used to establish the transcripts, with gffcompare (http://ccb.jhu.edu/software/stringtie/gffcompare.shtml) [38, 39] being used to identify novel transcripts. For this approach, transcripts overlapping with known mRNAs and lncRNAs or transcripts < 200 bp long were screened. Then, CPC0.9-r2 (http://cpc2.cbi.pku.edu.cn) and CNCI2.0 (https://github.com/www-bioinfo-org/CNCI#install-cnci) were used to predict new transcripts with coding potential using the default analytical parameters, retaining those transcripts with a CPC score < 0.5 and a CNCI score < 0 as putative novel lncRNAs. The remaining transcripts with class codes (I, j, o, u, x) were regarded as novel lncRNAs. The known and novel lncRNA datasets were then combined into a final lncRNA dataset for further analysis.

Weighted gene co-expression network analysis (WGCNA)

An R WGCNA package was used to perform a WGCNA using the provided tutorials. Euclidean distances calculated using gene expression data and integrated with growth and development-related parameters were used to cluster the 15 cattle LDM samples in this study. Network topology analyses ensured a scale-free topology, using a soft-thresholding power of 5. A dynamic tree-cutting algorithm with parameters minModuleSize at 30 and mergeCutHeight at 0.25 identified four modules. The eigengene (defined as the first component expression of genes in that module) was determined, and correlations between these eigengenes and cattle body weight, withers height, hip height, body length, chest circumference, abdominal circumference and cannon bone circumference were assessed. Genes exhibiting a high degree of connectivity within the established modules were regarded as hub genes.

Functional enrichment analyses

The potential functional roles of the analyzed lncRNAs were explored by selecting expressed mRNAs within 100 kb as putative cis-target genes [40]. Significantly differentially expressed lncRNAs and mRNAs were identified based on a threshold of |log2FC|> 1 and P < 0.05. mRNAs within 100 kb upstream or downstream of lncRNAs were identified as potential target genes. Finally, these mRNAs were used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses using OmicStudio (https://www.omicstudio.cn/tool), with terms having P < 0.05 considered statistically significant.

MeRIP-qPCR

A NEB/EpiMark® N6-Methyladenosine Enrichment Kit (Cat. E1610S, New England Biolabs, Ipswich, MA, USA) was used for MeRIP analyses establishing individual transcript m6A modification status. Briefly, 150 μg of total RNA from pretreated samples was randomly fragmented to yield < 100 nucleotide fragments, which were then immunoprecipitated using magnetic beads that had been coated with 10 μg of anti-m6A (Cat. 202003, Synaptic Systems, Göttingen, GER). The precipitated m6A-modified fragments of RNA were then eluted with N6-methyladenosine 5′-monophosphate sodium salt (6.7 mmol/L), after which MeRIP-qPCR analyses were performed with appropriate primers developed with MeRIP-Seq data and the motif-dependent m6A site predictor SRAMP (http://www.cuilab.cn/sramp) (Table S3). Relative m6A enrichment was normalized against the input RNA as follows: %input = 1/10 × 2Ct [IP] – Ct [input].

qPCR

TRIzol (Cat. 15596-026, Invitrogen, Carlsbad, CA, USA) was used to extract total RNA from PSCs and tissue samples, and a PrimeScript™ RT reagent Kit with gDNA Eraser (Cat. RR047A; Takara, Shiga, Japan) was used to prepare cDNA that was subsequently used for qPCR analyses on an ABI 2720 Thermal Cycler (Thermo Fisher Scientific, Waltham, MA, USA) with PowerUp SYBR Green Master Mix (Cat. A25742, Thermo Fisher Scientific, Waltham, MA, USA). Relative expression was established with the Ct (2−ΔΔCT) method. Primers for selected genes (Table S4) were designed using the NCBI Primer BLAST software. Primer pairs were selected to minimize non-specific amplification, with high ΔG values to avoid self or pair dimers and hairpin formations. Primer pairs were designed to amplify products spanning exon-exon junctions to avoid annotated variants from public SNP databases (https://www.ncbi.nlm.nih.gov/snp/). After the primer design, the predicted product was BLAST-searched against the bovine database to ensure the specificity of the primers.

Statistical analyses

qPCR data were analyzed using SAS (version 9.2; SAS Institute, Cary, NC, USA). GraphPad Prism (version 6.0; GraphPad Software, San Diego, CA, USA) was used to prepare figures. Data are presented as mean ± standard error of the mean (SEM), with significance set at P < 0.05.

Results

Characterization of the dynamic RNA methylome and transcriptome during skeletal muscle development

In an effort to systematically examine the regulatory and functional importance of RNA methylation in the context of skeletal muscle development, the dynamic changes in the RNA methylome and transcriptome were characterized via MeRIP-seq and RNA-seq in skeletal muscle tissue samples over five stages of development (Fig. 1A). MeRIP-seq analyses yielded 1,121,349,968 raw reads covering 160.43 Gb of sequence (Table S1), while RNA-seq analyses yielded 1,266,898,132 raw reads covering 181.69 Gb of sequence (Table S2). Clean MeRIP-seq and RNA-seq reads had average mapping rates of 95.64% (94.03%–96.55%) and 96.46% (95.84%–96.88%) to the ARS-UCD1.2 reference genome (Ensemblv107), respectively.

Fig. 1
figure 1

RNA methylome and transcriptomic profiles during skeletal muscle development. A The experimental workflow for profiling the RNA methylome and transcriptome across five stages of skeletal muscle development in cattle (illustration by Figdraw, ID: ISIIO07141). B Methylation levels of lncRNAs at each developmental stage. C Distribution map of lncRNA expression. D Principal component analysis (PCA) of the samples. E Pearson correlation analysis of RNA-seq data for each sample pair, accompanied by hierarchical clustering. Darker colors represent stronger correlations between samples. *P < 0.05

Full size image

Dynamic changes in lncRNA m6A methylation were assessed over the skeletal muscle developmental process using methylomic sequencing data, revealing the lowest and highest levels of lncRNA methylation in the M6 and M20 samples, respectively (Fig. 1B). Transcriptomic data showed high intersample correlation, with no clear differentiation among groups based on Pearson's correlation coefficients, principal component analyses and density curves (Fig. 1C–E).

Developing skeletal muscle exhibits dynamic changes in lncRNA methylation

Dynamic changes in the lncRNA methylome were examined throughout skeletal muscle development. We detected 419, 416, 417, 380 and 402 lncRNAs expressed in the M0, M6, M12, M20 and M30 groups respectively (Fig. 2A), of which 255, 244, 256, 244 and 214 were respectively m6A-modified (Fig. 2B), for respective m6A modification proportions of 60.86%, 58.65%, 61.39%, 64.21% and 53.23% (average: 59.67%, Fig. 2C). These dynamic changes in lncRNA m6A modification over the developmental process also coincided with a general reduction in the overall level of RNA methylation, declining from 60.86% at the newborn (M0) stage to 53.23% at 30 months of age (M30).

Fig. 2
figure 2

Global lncRNA methylation dynamics during skeletal muscle development. A Number of lncRNAs expressed at each developmental stage. B Number of m6A-modified lncRNAs across stages. C Proportion of lncRNAs with m6A modifications at each stage. D Distribution of peak counts in m6A-modified lncRNAs. E Correlations between lncRNA expression levels and m6A peak counts. F–G Gene Ontology (GO) (F) and KEGG pathway (G) enrichment analyses of lncRNAs with single m6A peaks. H–I RNA-seq-based expression analysis of RNA methyltransferase and demethylase genes (H) and their corresponding heatmaps (I). *P < 0.05, ***P < 0.01

Full size image

Numbers of m6A peaks for each lncRNA were analyzed, revealing that most had 1–2 such peaks (Fig. 2D). Correlation analyses showed that lncRNAs with three or more m6A peaks had higher expression levels than those with one or two peaks (Fig. 2E). The relationship between lncRNA function and m6A peak numbers was also assessed. As lncRNAs cannot code for proteins under normal circumstances, they function primarily through their effects on gene targets. Cis-target genes of these lncRNAs were thus leveraged in an effort to understand their potential functions through GO and KEGG enrichment analyses of these genes. These analyses showed no significant differences in enriched pathways or functions based on m6A peak numbers. These cis-target genes were enriched in pathways such as transcriptional regulation, DNA binding, RNA polymerase II transcriptional function, cAMP, calcium signaling, Apelin, cGMP-PKG and TNF signaling (Fig. 2F and G, Fig. S1).

Changes in RNA methylation status are under the dynamic control of RNA methyltransferases and demethylases [41,42,43,44]. RNA-seq data were thus used to assess the expression of these genes. The most highly expressed methyltransferase and demethylase genes identified in this study were VIRMA and ALKBH5, respectively (Fig. 2H). VIRMA levels were gradually downregulated with the progression of skeletal muscle development, whereas ALKBH5 levels were upregulated. These trends coincided with decreased muscle cell proliferation as indicated by the reduced expression of the proliferation biomarker MKI67 throughout development (Fig. 2I). These results indicated that the machinery responsible for maintaining methylation and generating de novo methylation were both downregulated with the development of skeletal muscle.

Dynamic changes in differentially methylated peaks and lncRNAs are related to skeletal muscle development

Initial analysis identified 316, 317, 323, 320 and 271 m6A peaks across the five stages of development (M0, M6, M12, M20 and M30) that were respectively associated with 255, 244, 256, 244 and 214 methylated lncRNAs (Fig. 3A). To fully explore the changes in m6A methylation status over the process of skeletal muscle development, pairwise comparisons were employed to analyze m6A peaks and methylated lncRNAs between consecutive stages. For the M0–M6, M6–M12, M12–M20 and M20–M30 comparisons, 41, 54, 75 and 67 differentially methylated peaks (DMPs), respectively, were identified that were associated with 40, 53, 69 and 63 differentially methylated lncRNAs (DM-lncRNAs) (Fig. 3B). A total of 225 DM-lncRNAs were identified, with 151 retained for further analysis after duplicate removal. The four control groups had 21, 39, 64 and 16 upregulated DMPs (20, 38, 60 and 15 DM-lncRNAs) and 20, 15, 11 and 51 downregulated DMPs (20, 15, 9 and 48 DM-lncRNAs) (Fig. 3C and D). Only a small number of common lncRNAs were detected among these DM-lncRNAs, suggesting that the development of skeletal muscle is characterized by very dynamic shifts in m6A levels (Fig. 3E). Potential DM-lncRNA functions in this context were next evaluated through GO and KEGG enrichment analyses of the 224 cis-target genes associated with these DM-lncRNAs. In GO analyses, these cis-target genes were primarily enriched in the regulation of DNA-templated transcription and DNA binding biological processes (Fig. 3F). KEGG pathway analyses indicated that these cis-target genes were primarily enriched in metabolic and human disease-related pathways, in addition to the PI3K-Akt, cAMP and calcium signaling pathways relevant to skeletal muscle development (Fig. 3G). Together, these analyses revealed many DM-lncRNAs expressed in the skeletal muscle in different stages of development, emphasizing the roles that these lncRNAs may play in this developmental process.

Fig. 3
figure 3

Differential methylation of peaks and lncRNAs in skeletal muscle. A Number of m6A peaks and methylated lncRNAs at specific developmental stages. B Count of DMPs and DM-lncRNAs. C–D Count of hypomethylated and hypermethylated DMPs (C) and DM-lncRNAs (D). E Venn diagram showing overlaps among DM-lncRNAs. F–G GO (F) and KEGG (G) enrichment analyses of cis-target genes associated with DM-lncRNAs

Full size image

Skeletal muscle development is characterized by dynamic changes in the lncRNA transcriptome

Transcriptomic analysis identified 24,865 lncRNAs, classified as intronic 45.29% (11,261), intergenic 44.64% (11,099), sense-overlapping 3.62% (900), antisense 3.31% (823) and bidirectional 2.41% (599) (Fig. 4A). These included 23,435 novel lncRNAs (Fig. 4B). To better characterize these novel lncRNAs, comparisons with gene structure and expression were performed among novel lncRNAs, annotated lncRNAs and mRNAs. Novel and annotated lncRNAs showed similar features, including shorter transcript lengths, fewer exons and shorter ORFs than mRNAs. Most of these novel and annotated lncRNAs had expression levels from 0–0.5, whereas most mRNAs had expression values in the 0–0.25 range (Fig. 4C–F).

Fig. 4
figure 4

Characterization of lncRNAs identified during skeletal muscle development. A Proportion of various lncRNA types identified in this study. B Classification of identified lncRNAs. C–E Distribution of exon numbers (C), transcript lengths (D) and open reading frame (ORF) lengths (E) in lncRNAs and mRNAs. F Expression level distribution (log10(FPKM + 1)) of lncRNAs and mRNAs

Full size image

Transcriptomic data were used in a WGCNA to identify correlations between annotated lncRNAs and traits relevant to cattle development. A network with scale-free topology was achieved at β = 5, with a scale independence value of 0.85 and lower levels of mean connectivity (Fig. 5A). lncRNAs showing similar expression dynamics were grouped into modules via hierarchical clustering, with a height threshold of 0.25, merging highly similar modules until ultimately obtaining a final set of four modules (Fig. 5B). Correlations between these modules and cattle traits were then assessed, revealing that the turquoise model was the most strongly negatively correlated with all of these traits (r = −0.97 to −0.91, P = 9e-09 to 1e-06) (Fig. 5C). This suggests that genes in the turquoise module may regulate these critical growth-related traits. Hub lncRNAs within the turquoise module were then used for functional enrichment analysis, and correlations between module membership and all growth and development-related traits were established. A preliminary assessment led to the selection of 21 hub lncRNAs from this module (Fig. 5D, Fig. S2). GO enrichment analyses of cis-target genes associated with these 21 lncRNAs were then performed, revealing their enrichment in terms including leukotriene receptor activity, troponin complex and negative regulation of miRNA transcription (Fig. 5E). They were also enriched in the KEGG Wnt, cGMP-PKG, cAMP and calcium signaling pathways (Fig. 5F). These 21 hub lncRNAs may thus serve as important mediators of the dynamic regulation necessary for appropriate skeletal muscle development.

Fig. 5
figure 5

Weighted Gene Co-expression Network Analysis (WGCNA) of annotated lncRNAs in skeletal muscle development. A Scale-free network coefficient (R2) analysis for soft thresholding (β) and corresponding mean connectivity. The red line indicates R2 = 0.85, with β = 5. B Cluster dendrogram based on topological overlap dissimilarity, identifying four modules. C Heatmap showing correlations between module eigengenes and traits associated with growth and development. Each cell includes the correlation coefficient and P value. D Scatter plot of module eigengenes, highlighting the turquoise module, where hub lncRNAs were associated with body weight (module membership = 0.8, gene significance = 0.2). E–F GO (E) and KEGG (F) enrichment analyses of hub lncRNAs

Full size image

Establishment of key skeletal muscle development-related lncRNAs

Two methods were next leveraged for the identification of key lncRNAs during the development of skeletal muscle. In the initial approach (Method 1), the intersection of the 151 DM-lncRNAs and 21 hub lncRNAs identified above yielded a list of 10 shared genes established as candidate lncRNAs (Fig. 6A).

Fig. 6
figure 6

Selection of key lncRNAs associated with skeletal muscle development. A Venn diagram illustrating the overlap between DM-lncRNAs and hub lncRNAs. B A four-quadrant diagram showing 45 dme-lncRNAs identified in M0–M6, M6–M12, M12–M20 and M20–M30 comparisons. C–D GO (C) and KEGG (D) enrichment analyses for cis-target genes of dme-lncRNAs. E Overview of lncRNAs and their associated cis-target genes. F Venn diagram showing overlap between Methods 1 and 2

Full size image

In the second approach (Method 2), a conjoint analysis of the MeRIP-seq and RNA-seq datasets was performed. Pairwise comparisons identified 5, 10, 27 and 3 lncRNAs that were both differentially m6A modified and differentially expressed (dme-lncRNAs). Of the 45 total dme-lncRNAs, 12 upregulated ones were significantly methylated, with 9 hyper-methylated (hyper-up) and 3 hypo-methylated (hypo-up). Meanwhile, 33 down-regulated lncRNAs were significantly methylated, with 27 hyper-methylated (hyper-down) and 6 hypo-methylated (hypo-down) (Fig. 6B). Following duplicate lncRNAs removal, 36 dme-lncRNAs were retained as candidate lncRNAs associated with 78 putative cis-target genes. GO enrichment analyses of these genes indicated that they were enriched in the DNA-templated transcription and DNA binding terms (Fig. 6C). Furthermore, three cis-target genes associated with two dme-lncRNAs (MSTRG.27754 and MSTRG.18394) were enriched in pathways related to skeletal muscle development. Cis-target genes for MSTRG.18394 were identified as the HOX family members HOXD3, HOXD9 and HOXD10, all of which control key developmental processes [45,46,47]. KEGG pathway analyses of the cis-target genes of these lncRNAs also indicated that they are enriched in muscle development-related pathways including the Hedgehog, Wnt and cGMP-PKG signaling pathways (Fig. 6D). In total, 10 cis-target genes associated with 8 dme-lncRNAs (MSTRG.27754, MSTRG.8738, MSTRG.23985, MSTRG.18394, MSTRG.7208, MSTRG.16924, MSTRG.28502 and ENSBTAG00000052793) were analyzed. This approach yielded 8 dme-lncRNAs that may help shape skeletal muscle development (Fig. 6E). In conclusion, 16 total candidate lncRNAs were identified by combining results derived from these two analytical approaches (Fig. 6F).

MeRIP-qPCR and qPCR validation of study results

Three dme-lncRNAs (MSTRG.8738, MSTRG.18394 and ENSBTAG00000052793) were selected for MeRIP-qPCR and qPCR to validate their methylation and expression levels (Fig. 7A and B). The results showed that these three dme-lncRNAs had m6A enrichment across five developmental stages. The m6A methylation levels of MSTRG.8738 and ENSBTAG00000052793 were significantly upregulated, and expression levels were significantly downregulated in M12 samples compared to M6 samples. The m6A methylation level of MSTRG.18394 was significantly upregulated, and the expression level showed an upward trend in M12 samples compared to M6 samples.

Fig. 7
figure 7

Validation of findings using MeRIP-qPCR and qPCR. A–B Validation of three dme-lncRNAs using MeRIP-qPCR (A) and qPCR (B). C Validation of five cis-target genes using qPCR. Relative expression levels were normalized to GAPDH. Data are presented as mean ± SD for three independent biological replicates. *P < 0.05, **P < 0.01, ***P < 0.001

Full size image

The expression levels of the cis-target genes linked to these three dme-lncRNAs were also analyzed. Compared with M6, the SFMBT2 expression level in M12 samples showed an upward trend. However, HOXD10 expression was significantly downregulated, while KAT7 expression was significantly upregulated (Fig. 7C). These results confirm the m6A modification status of specific lncRNAs in Bohai black bovine skeletal muscle.

Discussion

This study presents dynamic genome-wide transcriptomic and methylomic maps for five stages of skeletal muscle development. This is the first reported effort to systemically profile the lncRNA m6A methylome in bovine skeletal muscle to the best of our knowledge. While it is important to note that the multi-omics analyses in this study relied on the use of bulk tissue instead of single muscle fibers such that the results may have been influenced by differences in fiber time composition and the populations of non-muscle cells present over the course of development, such forms of bias have been suggested to have minimal confounding effects on study conclusions [48]. These datasets are thus comprehensive tools for efforts to understand the roles that lncRNA methylation plays in the development of skeletal muscle.

In addition to being ubiquitously present among mRNA transcripts, m6A methylation is closely tied to the control of gene expression [49]. lncRNA m6A methylation has been a significant focus of research. Here, MeRIP-seq identified several m6A peaks in lncRNAs from bovine skeletal muscle, with results validated by MeRIP-qPCR analyses. The m6A methylation of mRNAs can impact their transport, splicing, translation and stability [15]. Similarly, the m6A methylation of lncRNAs can affect expression levels thereof. THAP7-AS1, for instance, can be transcriptionally induced by SP1 and m6A modified by METTL3, whereupon it can support oncogenesis by promoting NLS interactions with importin α1 and enhancing nuclear CUL4B protein entry, repressing miR-22-3p and miR-320a transcription [50]. The stability of the lncRNA DIAPH1-AS1 can also be improved by m6A modification mediated by WTAP through a pathway dependent on IGF2BP2, whereupon it functions as a molecular adaptor capable of promoting the formation of the MTDH-LASP1 complex and upregulating LASP1 to facilitate the growth and metastasis of nasopharyngeal carcinoma [51].

Approximately 59.67% of lncRNAs showed m6A modifications, with levels declining during development, likely due to reduced VIRMA expression and increased ALKBH5 expression. The m6A modification of Traf6 was predicted to regulate translation, with no effect observed for transcripts with 1–2 m6A peaks, while transcripts with 3 peaks activated the YTHDF1 regulatory mechanism [52]. In this analysis, lncRNAs with three or more m6A peaks showed significantly higher expression levels than those with one or two peaks. The mechanistic basis for this observation is uncertain and warrants further research.

RNA-seq analyses of LDM samples from Bohai black cattle at five developmental stages (0, 6, 12, 20 and 30 months) characterized transcriptomic dynamics with high confidence [53]. Novel lncRNAs identified herein were consistent with previous reports regarding length, exon numbers and ORF length [54,55,56,57,58]. Most of the novel and annotated lncRNAs had expression values in the 0–0.5 range. However, most mRNAs showed expression values in the 0–0.25 range, differing from previous findings [54,55,56,57,58].

Two approaches were subsequently employed to identify specific lncRNAs that shape the development of skeletal muscle. Initially, transcriptomic data were analyzed using WGCNA to identify correlations between annotated lncRNAs and developmental traits, resulting in a list of 21 putative hub lncRNAs. The intersection of these 21 hub lncRNAs with 151 DM-lncRNAs identified 10 overlapping candidate lncRNAs. The second approach integrated MeRIP-seq and RNA-seq datasets to identify 36 dme-lncRNAs. Cis-target gene analyses have provided valuable insights into lncRNA functions [56]. Accordingly, 78 cis-target genes associated with these dme-lncRNAs were herein identified and used to conduct additional functional enrichment analyses. This strategy ultimately identified 8 putative dme-lncRNAs that may influence skeletal muscle development. These two methods identified 16 lncRNAs as final candidates, representing potential regulators of skeletal muscle formation.

MeRIP-qPCR and qPCR were finally used to determine the methylation and gene expression levels of three dme-lncRNAs (MSTRG.8738, MSTRG.18394 and ENSBTAG00000052793) and their cis-target genes. In comparison to M6, the methylation level of MSTRG.8738 in M12 was significantly upregulated, and the expression level was significantly downregulated, while the expression of its cis-target gene SFMBT2 demonstrated an increasing trend. These results suggest that increased m6A methylation inhibited MSTRG.8738 expression while elevating SFMBT2 expression, contributing to skeletal muscle development. Nevertheless, the exact mechanism by which MSTRG.8738, with m6A modification and its cis-target gene SFMBT2, regulates myoblast proliferation and differentiation remains unclear and warrants further investigation. Compared with M6, MSTRG.18394 methylation level in M12 was significantly upregulated, and the expression level showed an upward trend, whereas the expression level of its cis-target gene HOXD10 was significantly downregulated. These results indicate that increased m6A methylation enhanced MSTRG.18394 expression while inhibiting HOXD10 expression, supporting skeletal muscle growth and development. The 39 HOX genes in mammals are organized into four clusters labeled A through D. These genes can be categorized into 13 paralogous groups (1–13) according to their sequence similarities and positions within the clusters [59, 60]. For the proper growth and skeletal design of tetrapod limbs, HOX genes are needed, especially the HOXA and HOXD clusters, which are vital for the development of both forelimbs and hindlimbs [61, 62]. In the limb development of mice, Shh expression is primarily driven by the genes HOXA9, 10, 11 and HOXD9, 10, 11 [63]. Studies indicate that Shh is vital for the early induction of the myogenic determination genes Myf5 and MyoD in epaxial somite cells, which lead to the formation of deep back muscle progenitors [64]. Compared with M6, ENSBTAG00000052793 methylation level in M12 was significantly upregulated, and the expression level was downregulated, while the expression level of its cis-target gene KAT7 was significantly upregulated. These findings suggest that increased m6A methylation inhibited ENSBTAG00000052793 expression and increased KAT7 expression, facilitating skeletal muscle growth and development. Previous studies have shown that lncRNA ADAMTS9-AS is competitively bound to miR-185-5p to upregulate KAT7 and thus inhibit cardiomyocyte hypertrophy [65]. Another investigation found that circFoxo3 alleviated myocardial ischemia/reperfusion injury by reducing autophagy, achieved by inhibiting HMGB1 through the suppression of KAT7 in myocardial infarction [66]. The three genes identified in this study are potential candidates for regulating skeletal muscle development. Further investigation is required to shed light on how these genes regulate and influence skeletal muscle development.

Conclusions

In conclusion, this study analyzed the expression and m6A methylation profiles of lncRNAs linked to skeletal muscle development, ultimately revealing 16 m6A-modified lncRNAs that may play a key regulatory role during this process. These findings may offer evidence for future studies to clarify the mechanistic functions of these m6A-modified lncRNAs, providing an opportunity for more comprehensive analyses of the epigenetic modification of RNA during skeletal muscle development.

Data availability

The datasets generated and analyzed during the current study are not publicly available. They can be made available from the corresponding author upon reasonable request.

Abbreviations

ALKBH5:

AlkB homolog 5

ceRNA:

Competing endogenous RNA

dme-lncRNAs:

Differentially m6A modified and differentially expressed lncRNAs

DM-lncRNAs:

Differentially methylated lncRNAs

DMP:

Differentially methylated peaks

FTO:

Fat mass and obesity related

GO:

Gene Ontology

KEGG :

Kyoto Encyclopedia of Genes and Genomes

IGF2BP1–3 :

Insulin-like growth factor binding proteins 1–3

LDM :

Longissimus dorsi muscle

lncRNA :

Long non-coding RNA

MeRIP-qPCR :

Methylated RNA immunoprecipitation-real-time quantitative PCR

MeRIP-seq :

Methylated RNA immunoprecipitation sequencing

MKI67 :

Marker of proliferation Ki-67

METTL3 :

Methyltransferase-like 3

METTL14 :

Methyltransferase-like 14

m6A :

N6-methyladenosine

Prrc2a :

Proline-rich coiled-coil2A

qPCR :

Real-time quantitative PCR

RNA-seq :

RNA sequencing

UTR :

Untranslated region

VIRMA :

Vir like N6-methyladenosine methyltransferase associated protein

WGCNA :

Weighted gene co-expression network analysis

WTAP :

Wilms tumor 1-associating protein

YTHDC1–2 :

YTH domain-containing 1–2

YTHDF1–3 :

YTH-domain family 1–3

References

  1. Buckingham M. Myogenic progenitor cells and skeletal myogenesis in vertebrates. Curr Opin Genet Dev. 2006;16(5):525–32. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.gde.2006.08.008.

    Article  PubMed  CAS  Google Scholar 

  2. Cagnazzo M, te Pas MFW, Priem J, de Wit AA, Pool MH, Davoli R, et al. Comparison of prenatal muscle tissue expression profiles of two pig breeds differing in muscle characteristics. J Anim Sci. 2006;84(1):1–10. https://doiorg.publicaciones.saludcastillayleon.es/10.2527/2006.8411.

    Article  PubMed  CAS  Google Scholar 

  3. Buckingham M, Rigby PW. Gene regulatory networks and transcriptional mechanisms that control myogenesis. Dev Cell. 2014;28(3):225–38. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.devcel.2013.12.020.

    Article  PubMed  CAS  Google Scholar 

  4. Endo T. Molecular mechanisms of skeletal muscle development, regeneration, and osteogenic conversion. Bone. 2015;80:2–13. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.bone.2015.02.028.

    Article  PubMed  CAS  Google Scholar 

  5. Chen LL, Carmichael GG. Decoding the function of nuclear long non-coding RNAs. Curr Opin Cell Biol. 2010;22(3):357–64. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.ceb.2010.03.003.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  6. Xie SJ, Diao LT, Cai N, Zhang LT, Xiang S, Jia CC, et al. mascRNA and its parent lncRNA MALAT1 promote proliferation and metastasis of hepatocellular carcinoma cells by activating ERK/MAPK signaling pathway. Cell Death Discov. 2021;7:110. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41420-021-00497-x.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  7. Wang S, Jin J, Xu Z, Zuo B. Functions and regulatory mechanisms of lncRNAs in skeletal myogenesis, muscle disease and meat production. Cells. 2019;8(9):1107. https://doiorg.publicaciones.saludcastillayleon.es/10.3390/cells8091107.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  8. Song C, Yang Z, Jiang R, Cheng J, Yue B, Wang J, et al. lncRNA IGF2 AS regulates bovine myogenesis through different pathways. Mol Ther Nucleic Acids. 2020;21:874–84. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.omtn.2020.07.002.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  9. Cheng X, Li L, Shi G, Chen L, Fang C, Li M, et al. MEG3 promotes differentiation of porcine satellite cells by sponging miR-423-5p to relieve inhibiting effect on SRF. Cells. 2020;9(2):449. https://doiorg.publicaciones.saludcastillayleon.es/10.3390/cells9020449.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  10. Zou C, Li L, Cheng X, Li C, Fu Y, Fang C, et al. Identification and functional analysis of long intergenic non-coding RNAs underlying intramuscular fat content in pigs. Front Genet. 2018;9:102. https://doiorg.publicaciones.saludcastillayleon.es/10.3389/fgene.2018.00102.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  11. Wu T, Wang S, Wang L, Zhang W, Chen W, Lv X, et al. Long noncoding RNA (lncRNA) CTTN-IT1 elevates skeletal muscle satellite cell proliferation and differentiation by acting as ceRNA for YAP1 through absorbing miR-29a in Hu sheep. Front Genet. 2020;11:843. https://doiorg.publicaciones.saludcastillayleon.es/10.3389/fgene.2020.00843.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  12. Li Z, Cai B, Abdalla BA, Zhu X, Zheng M, Han P, et al. LncIRS1 controls muscle atrophy via sponging miR-15 family to activate IGF1-PI3K/AKT pathway. J Cachexia Sarcopenia Muscle. 2019;10(2):391–410. https://doiorg.publicaciones.saludcastillayleon.es/10.1002/jcsm.12374.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Huang J, Chen YX, Zhang B. IGF2-AS affects the prognosis and metastasis of gastric adenocarcinoma via acting as a ceRNA of miR-503 to regulate SHOX2. Gastric Cancer. 2020;23(1):23–38. https://doiorg.publicaciones.saludcastillayleon.es/10.1007/s10120-019-00976-2.

    Article  PubMed  CAS  Google Scholar 

  14. Lv W, Jin J, Xu Z, Luo H, Guo Y, Wang X, et al. lncMGPF is a novel positive regulator of muscle growth and regeneration. J Cachexia Sarcopenia Muscle. 2020;11(6):1723–46. https://doiorg.publicaciones.saludcastillayleon.es/10.1002/jcsm.12623.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Meyer KD, Jaffrey SR. The dynamic epitranscriptome: N6-methyladenosine and gene expression control. Nat Rev Mol Cell Biol. 2014;15(5):313–26. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/nrm3785.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  16. Lence T, Paolantoni C, Worpenberg L, Roignant JY. Mechanistic insights into m6A RNA enzymes. Biochim Biophys Acta Gene Regul Mech. 2019;1862(3):222–9. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.bbagrm.2018.10.014.

    Article  PubMed  CAS  Google Scholar 

  17. Gerken T, Girard CA, Tung YC, Webby CJ, Saudek V, Hewitson KS, et al. The obesity-associated FTO gene encodes a 2-oxoglutarate-dependent nucleic acid demethylase. Science. 2007;318(5855):1469–72. https://doiorg.publicaciones.saludcastillayleon.es/10.1126/science.1151710.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  18. Falnes P, Johansen RF, Seeberg E. AlkB-mediated oxidative demethylation reverses DNA damage in Escherichia coli. Nature. 2002;419(6903):178–82. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/nature01048.

    Article  PubMed  CAS  Google Scholar 

  19. Yue Y, Liu J, He C. RNA N6-methyladenosine methylation in post-transcriptional gene expression regulation. Genes Dev. 2015;29(13):1343–55. https://doiorg.publicaciones.saludcastillayleon.es/10.1101/gad.262766.115.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  20. Huang H, Weng H, Sun W, Qin X, Shi H, Wu H, et al. Recognition of RNA N6-methyladenosine by IGF2BP proteins enhances mRNA stability and translation. Nat Cell Biol. 2018;20(3):285–95. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41556-018-0045-z.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  21. Wu R, Li A, Sun B, Sun JG, Zhang J, Zhang T, et al. A novel m6A reader Prrc2a controls oligodendroglial specification and myelination. Cell Res. 2019;29(1):23–41. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41422-018-0113-8.

    Article  PubMed  CAS  Google Scholar 

  22. 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 Genomics. 2020;47(8):466–76. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.jgg.2020.07.003.

    Article  PubMed  CAS  Google Scholar 

  23. Deng K, Fan Y, Liang Y, Cai Y, Zhang G, Deng M, et al. FTO-mediated demethylation of GADD45B promotes myogenesis through the activation of p38 MAPK pathway. Mol Ther Nucleic Acids. 2021;26:34–48. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.omtn.2021.06.013.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  24. Yang X, Mei C, Ma X, Du J, Wang J, Zan L. m6A methylases regulate myoblast proliferation, apoptosis and differentiation. Animals (Basel). 2022;12(6):773. https://doiorg.publicaciones.saludcastillayleon.es/10.3390/ani12060773.

    Article  PubMed  Google Scholar 

  25. Chen S, Zhou Y, Chen Y, Gu J. Fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–90. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/bioinformatics/bty560.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  26. de Sena Brandine G, Smith AD. Falco: high-speed FastQC emulation for quality control of sequencing data. F1000Res. 2019;8:1874. https://doiorg.publicaciones.saludcastillayleon.es/10.12688/f1000research.21142.2.

  27. Wang L, Wang S, Li W. RSeQC: quality control of RNA-seq experiments. Bioinformatics. 2012;28(16):2184–5. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/bioinformatics/bts356.

    Article  PubMed  CAS  Google Scholar 

  28. Kim D, Langmead B, Salzberg SL. HISAT: A fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/nmeth.3317.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  29. 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://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.ymeth.2014.06.008.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  30. Wang K, Li M, Hakonarson H. ANNOVAR: Functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res. 2010;38(16):e164. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/nar/gkq603.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  31. Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. 2010;38(4):576–89. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.molcel.2010.05.004.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  32. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/nbt.3122.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  33. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/bioinformatics/btp616.

    Article  PubMed  CAS  Google Scholar 

  34. Kechin A, Boyarskikh U, Kel A, Filipenko M. cutPrimers: A new tool for accurate cutting of primers from reads of targeted next generation sequencing. J Comput Biol. 2017;24(11):1138–43. https://doiorg.publicaciones.saludcastillayleon.es/10.1089/cmb.2017.0096.

    Article  PubMed  CAS  Google Scholar 

  35. 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://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41587-019-0201-4.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  36. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. 2016;11(9):1650–67. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/nprot.2016.095.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  37. Kovaka S, Zimin AV, Pertea GM, Razaghi R, Salzberg SL, Pertea M. Transcriptome assembly from long-read RNA-seq alignments with StringTie2. Genome Biol. 2019;20:278. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s13059-019-1910-1.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  38. Kang YJ, Yang DC, Kong L, Hou M, Meng YQ, Wei L, et al. CPC2: A fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res. 2017;45(W1):W12–6. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/nar/gkx428.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  39. Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, et al. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;41(17):e166. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/nar/gkt646.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  40. Bao Z, Yang Z, Huang Z, Zhou Y, Cui Q, Dong. LncRNADisease 2.0: an updated database of long non-coding RNA-associated diseases. Nucleic Acids Res. 2019;47(D1):D1034–7. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/nar/gky905.

  41. Chao M, Wang M, Han H, Liu Y, Sun X, Tian T, et al. Profiling of m6A methylation in porcine intramuscular adipocytes and unravelling PHKG1 represses porcine intramuscular lipid deposition in an m6A-dependent manner. Int J Biol Macromol. 2024;272(Pt 1):132728. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.ijbiomac.2024.132728.

  42. Wang Y, Li Y, Toth JI, Petroski MD, Zhang Z, Zhao JC. N6-methyladenosine modification destabilizes developmental regulators in embryonic stem cells. Nat Cell Biol. 2014;16(2):191–8. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/ncb2902.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  43. Jia G, Fu Y, Zhao X, Dai Q, Zheng G, Yang Y, et al. N6-methyladenosine in nuclear RNA is a major substrate of the obesity-associated FTO. Nat Chem Biol. 2011;7(12):885–7. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/nchembio.687.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  44. Mauer J, Luo X, Blanjoie A, Jiao X, Grozhik AV, Patil DP, et al. Reversible methylation of m6Am in the 5’ cap controls mRNA stability. Nature. 2017;541(7637):371–5. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/nature21022.

    Article  PubMed  CAS  Google Scholar 

  45. Wang M, Xie H, Shrestha S, Sredni S, Morgan GA, Pachman LM. Methylation alterations of WT1 and homeobox genes in inflamed muscle biopsy samples from patients with untreated juvenile dermatomyositis suggest self-renewal capacity. Arthritis Rheum. 2012;64(10):3478–85. https://doiorg.publicaciones.saludcastillayleon.es/10.1002/art.34573.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  46. de Wilde J, Hulshof MF, Boekschoten MV, de Groot P, Smit E, Mariman EC. The embryonic genes Dkk3, Hoxd8, Hoxd9 and Tbx1 identify muscle types in a diet-independent and fiber-type unrelated way. BMC Genomics. 2010;11:176. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/1471-2164-11-176.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  47. Turner DC, Gorski PP, Maasar MF, Seaborne RA, Baumert P, Brown AD, et al. DNA methylation across the genome in aged human skeletal muscle tissue and muscle-derived cells: the role of HOX genes and physical activity. Sci Rep. 2020;10(1):15360. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41598-020-72730-z.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  48. Terry EE, Zhang X, Hoffmann C, Hughes LD, Lewis SA, Li J, et al. Transcriptional profiling reveals extraordinary diversity among skeletal muscle tissues. Elife. 2018;7:e34613. https://doiorg.publicaciones.saludcastillayleon.es/10.7554/eLife.34613.

    Article  PubMed  PubMed Central  Google Scholar 

  49. 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://doiorg.publicaciones.saludcastillayleon.es/10.1038/nature11112.

    Article  PubMed  CAS  Google Scholar 

  50. Liu HT, Zou YX, Zhu WJ, Sen L, Zhang GH, Ma RR, et al. lncRNA THAP7-AS1, transcriptionally activated by SP1 and post-transcriptionally stabilized by METTL3-mediated m6A modification, exerts oncogenic properties by improving CUL4B entry into the nucleus. Cell Death Differ. 2022;29(3):627–41. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41418-021-00879-9.

    Article  PubMed  CAS  Google Scholar 

  51. Li ZX, Zheng ZQ, Yang PY, Lin L, Zhou GQ, Lv JW, et al. WTAP-mediated m6A modification of lncRNA DIAPH1-AS1 enhances its stability to facilitate nasopharyngeal carcinoma growth and metastasis. Cell Death Differ. 2022;29(6):1137–51. https://doiorg.publicaciones.saludcastillayleon.es/10.1038/s41418-021-00905-w.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  52. Zong X, Xiao X, Shen B, Jiang Q, Wang H, Lu Z, et al. The N6-methyladenosine RNA-binding protein YTHDF1 modulates the translation of TRAF6 to mediate the intestinal immune response. Nucleic Acids Res. 2021;49(10):5537–52. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/nar/gkab343.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  53. Liu R, Han M, Liu X, Yu K, Bai X, Dong Y. Genome-wide identification and characterization of long non-coding RNAs in longissimus dorsi skeletal muscle of shandong black cattle and luxi cattle. Front Genet. 2022;13:849399. https://doiorg.publicaciones.saludcastillayleon.es/10.3389/fgene.2022.849399.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  54. Sun J, Xie M, Huang Z, Li H, Chen T, Sun R, et al. Integrated analysis of non-coding RNA and mRNA expression profiles of 2 pig breeds differing in muscle traits. J Anim Sci. 2017;95(3):1092–103. https://doiorg.publicaciones.saludcastillayleon.es/10.2527/jas.2016.0867.

    Article  PubMed  CAS  Google Scholar 

  55. Zhan S, Dong Y, Zhao W, Guo J, Zhong T, Wang L, et al. Genome-wide identification and characterization of long non-coding RNAs in developmental skeletal muscle of fetal goat. BMC Genomics. 2016;17:666. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-016-3009-3.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  56. Hong L, Hu Q, Zang X, Xie Y, Zhou C, Zou X, et al. Analysis and screening of reproductive long non-coding RNAs through genome-wide analyses of goat endometrium during the pre-attachment phase. Front Genet. 2020;11:568017. https://doiorg.publicaciones.saludcastillayleon.es/10.3389/fgene.2020.568017.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  57. Shen L, Gan M, Tang Q, Tang G, Jiang Y, Li M, et al. Comprehensive analysis of lncRNAs and circRNAs reveals the metabolic specialization in oxidative and glycolytic skeletal muscles. Int J Mol Sci. 2019;20(12):2855. https://doiorg.publicaciones.saludcastillayleon.es/10.3390/ijms20122855.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  58. Wang S, Tan B, Xiao L, Zhao X, Zeng J, Hong L, et al. Comprehensive analysis of long noncoding RNA modified by m6A methylation in oxidative and glycolytic skeletal muscles. Int J Mol Sci. 2022;23(9):4600. https://doiorg.publicaciones.saludcastillayleon.es/10.3390/ijms23094600.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  59. Mallo M, Alonso CR. The regulation of Hox gene expression during animal development. Development. 2013;140(19):3951–63. https://doiorg.publicaciones.saludcastillayleon.es/10.1242/dev.068346.

    Article  PubMed  CAS  Google Scholar 

  60. Alexander T, Nolte C, Krumlauf R. Hox genes and segmentation of the hindbrain and axial skeleton. Annu Rev Cell Dev Bi. 2009;25:431–56. https://doiorg.publicaciones.saludcastillayleon.es/10.1146/annurev.cellbio.042308.113423.

    Article  CAS  Google Scholar 

  61. Wellik DM, Capecchi MR. Hox10 and Hox11 genes are required to globally pattern the mammalian skeleton. Science. 2003;301(5631):363–7. https://doiorg.publicaciones.saludcastillayleon.es/10.1126/science.1085672.

    Article  PubMed  CAS  Google Scholar 

  62. Zakany J, Duboule D. The role of Hox genes during vertebrate limb development. Curr Opin Genet Dev. 2007;17(4):359–66. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/j.gde.2007.05.011.

    Article  PubMed  CAS  Google Scholar 

  63. Raines AM, Magella B, Adam M, Potter SS. Key pathways regulated by Hox A9,10,11/HoxD9,10,11 during limb development. BMC Dev Biol. 2015;15:28. https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12861-015-0078-5.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  64. Borycki A-G, Brunk B, Tajbakhsh S, Buckingham M, Chiang C, Emerson CP Jr. Sonic hedgehog controls epaxial muscle determination through Myf5 activation. Development. 1999;126(18):4053–63. https://doiorg.publicaciones.saludcastillayleon.es/10.1242/dev.126.18.4053.

    Article  PubMed  CAS  Google Scholar 

  65. Song B, Li W, Xu X, Dang H, Dong R. The lncRNA ADAMTS9-AS1/miR-185-5p/KAT7 ceRNA network inhibits cardiomyocyte hypertrophy in hypertrophic obstructive cardiomyopathy. Biomed Res. 2023;44(3):105–15. https://doiorg.publicaciones.saludcastillayleon.es/10.2220/biomedres.44.105.

    Article  PubMed  CAS  Google Scholar 

  66. Sun G, Shen JF, Wei XF, Qi GX. Circular RNA Foxo3 relieves myocardial ischemia/reperfusion injury by suppressing autophagy via inhibiting HMGB1 by repressing KAT7 in myocardial infarction. J Inflamm Res. 2021;14:6397–407. https://doiorg.publicaciones.saludcastillayleon.es/10.2147/JIR.S339133.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

Download references

Acknowledgements

The authors would like to thank all the reviewers who participated in the review and MJEditor (https://www.mjeditor.com) for their linguistic assistance during the preparation of this manuscript.

Funding

This work was supported by the National Key R&D Program of China (2023YFD1300103), the Science and Technology Plan Project of Yantai City (2023ZDCX024), the National Natural Science Foundation of China (32372852), the Science Fund for Distinguished Young Scholars of Shaanxi Province (2024JC-JCQN-30), Shaanxi Provincial Innovation Leadership Program in Sciences and Technologies for Young and Middle-aged Scientists (2023SR205), the China Agriculture Research System-beef (CARS-37), the Innovation Team of Cattle Industry in Technological System of Shandong Modern Agriculture Industry (SDAIT-09-03), the Key Research and Development Project in Shandong Province (Competitive Innovation Platform) (2022CXPT010).

Author information

Authors and Affiliations

  1. Key Laboratory of Animal Genetics, Breeding and Reproduction of Shaanxi Province, College of Animal Science and Technology, Northwest A&F University, Yangling, 712100, China

    Cui Mao, Yuta Yang & Xianyong Lan

  2. Institute of Animal Science and Veterinary Medicine, Shandong Academy of Agricultural Sciences, Jinan, 250100, China

    Wei You, Haijian Cheng, Xin Hu & Enliang Song

  3. Key Laboratory of Livestock and Poultry Multi-Omics of MARA, Jinan, 250100, China

    Wei You, Haijian Cheng, Xin Hu & Enliang Song

Authors
  1. Cui Mao
    View author publications

    You can also search for this author inPubMed Google Scholar

  2. Wei You
    View author publications

    You can also search for this author inPubMed Google Scholar

  3. Yuta Yang
    View author publications

    You can also search for this author inPubMed Google Scholar

  4. Haijian Cheng
    View author publications

    You can also search for this author inPubMed Google Scholar

  5. Xin Hu
    View author publications

    You can also search for this author inPubMed Google Scholar

  6. Xianyong Lan
    View author publications

    You can also search for this author inPubMed Google Scholar

  7. Enliang Song
    View author publications

    You can also search for this author inPubMed Google Scholar

Contributions

ES and XL conceived and designed the project. CM, WY, and YY performed the experiments and collected samples. CM analyzed the data and wrote the manuscript. HC and XH revised the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Xianyong Lan or Enliang Song.

Ethics declarations

Ethics approval and consent to participate

The experimental procedures complied with the guidelines of the Animal Experimental Ethics Committee of the Institute of Animal Science and Veterinary Medicine, Shandong Academy of Agricultural Sciences (Jinan, China).

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Supplementary Information

Additional file 1: Table S1. Sequencing statistics summary of the methylome data.

Additional file 2: Table S2. Sequencing statistics summary of the transcriptome data.

Additional file 3: Table S3. Primers used for MeRIP-qPCR experiments.

Additional file 4: Table S4. Primers used for qPCR experiments.

40104_2025_1164_MOESM5_ESM.tif

Additional file 5: Fig. S1. GO and KEGG enrichment analyses of lncRNAs with two or more m6A peaks. A GO enrichment analyses of lncRNAs with two m6A peaks. B KEGG enrichment analyses of lncRNAs with two m6A peaks. C GO enrichment analyses of lncRNAs with 3 + m6A peaks. D KEGG enrichment analyses of lncRNAs with 3 + m6A peaks.

40104_2025_1164_MOESM6_ESM.tif

Additional file 6: Fig. S2. Scatter plots of module eigengenes. A–F The genes of the turquoise module in the upper right were chosen as hub lncRNAs associated with withers height (A), hip height (B), body length (C), chest circumference (D) abdominal circumference (E) and cannon bone circumference (F) (module membership = 0.8, gene significance = 0.2).

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Mao, C., You, W., Yang, Y. et al. Comprehensive characterization of lncRNA N6-methyladenosine modification dynamics throughout bovine skeletal muscle development. J Animal Sci Biotechnol 16, 36 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s40104-025-01164-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s40104-025-01164-2

Keywords

Journal of Animal Science and Biotechnology

ISSN: 2049-1891