Skip to main content
Download PDF

Comprehensive transcriptome of muscle development in Sichuan white rabbit

BMC Genomic Data volume 26, Article number: 32 (2025) Cite this article

Abstract

Background

The Sichuan white rabbit is a unique domestic breed and is famous for its high meat production. Muscle development is a complicated biological process, but the underlying regulatory mechanisms have not been elucidated. Here, we generated comprehensive transcriptome datasets (i.e., mRNAs, miRNAs and lncRNAs) in three developmental stages of Sichuan white rabbits, and aim to systematically explore the regulatory network in myogenesis.

Results

We generated extensive transcriptome datasets (mRNAs, miRNAs and lncRNAs) revealing the myogenic regulatory network at different time points. Our differential expression analysis identified 2,995 DE genes, 1,211 DE-lncRNAs, and 305 DE-miRNAs with distinct expression patterns across developmental stages. In addition, functional enrichment analysis of DE mRNAs and miRNAs indicates their involvement in muscle growth, development, and regeneration, highlighting biological processes and muscle-specific functions. Interaction analysis between DE-lncRNAs and mRNAs uncovered a complex regulatory network, especially between 21 and 27 days of development. These findings contribute to better understanding of the transcriptomic changes during muscle development and have implications for breeding improvement in Sichuan white rabbits.

Conclusions

Our study provides a comprehensive overview of the transcriptomic changes during muscle development in Sichuan white rabbits. The identification and functional annotation of DE genes, miRNAs, and lncRNAs provide valuable insights into the molecular mechanisms underlying this process. These findings pave the way for targeted investigations into the role of non-coding RNAs in muscle biology.

Peer Review reports

Background

Rabbit meat becomes increasingly popular in China because of its remarkable nutritional value and unique flavour. On the one hand, rabbit meat has moderately high energy values, low fat contents, and low cholesterol levels [1] and greatly satisfies modern consumers’ desire for a healthy lifestyle. On the other hand, the delectable texture and taste of rabbit meat make it more prevalent worldwide, especially in China. As the largest rabbit meat producer worldwide, Chinese rabbit meat production has steadily increased from 690,000 tons in 2010 to 849,150 tons in 2016 (http://www.fao.org/faostat/en/#data/QL). China has various domesticated rabbit breeds, while most local breeds achieve a low meat production rate [2]. Sichuan white rabbit, due to its high-yield meat production, strong adaptability and fecundity [2, 3], has become one of the well-known breeds in China. Thus, improving the quality and production of Sichuan white rabbits is critical for breeding improvement.

Muscle is the largest organ of body mass in humans and other animals, and its functions include movement, postural support and thermogenesis [4]. The development of muscle is a complicated biological process that includes distinct embryonic and postnatal phases. Many diffusible signaling molecules, transcription factors and non-coding RNAs [e.g., microRNAs (miRNAs) and long non-coding RNAs (lncRNAs)] that contribute to muscle development have been identified. These regulators serve as direct templates for protein synthesis, which is fundamental to the growth and repair of muscle fibers [5]. In vertebrate embryos, numerous miRNAs are expressed in the developing somites of zebrafish [6], Xenopus [7] and chicks [8]. Muscle-specific mRNAs, such as miR- 1 and miR- 133, which are known to regulate muscle differentiation by targeting specific transcription factors and structural proteins, can regulate muscle maturation [9]. In addition, lncRNAs can act as scaffolds for certain regulatory proteins or be involved in the regulation of other non-coding RNAs [10,11,12]. Previous studies reported that lncRNAs associated with myogenesis include steroid receptor RNA activator (Sra), which co-activates MyoD [13], and LncMyoD, which is itself activated by MyoD, together with linc-MD1, which regulates miR- 133 to further enhance differentiation [14]. In addition, a novel lncRNA lncMGR, which promotes myoblast differentiation and muscle fiber hypertrophy, can recruit cyclin-dependent kinase 9 (CDK9) and sponge miRNAs, such as miR- 2131 - 5p, to regulate the expression of skeletal muscle myosin heavy chain 1 A (MYH1 A) [4].

Current studies are limited to the regulation of either mRNAs or miRNAs related to the proliferation and differentiation of skeletal muscle cells. Thus, the identification of more potential muscle-associated protein-coding genes (i.e., mRNAs) and non-coding RNAs (e.g., miRNAs and lncRNAs) could better dissect the underlying regulatory mechanisms. In this study, we generated comprehensive transcriptome datasets (i.e., mRNAs, miRNAs and lncRNAs) in different developmental stages of Sichuan white rabbit and aimed to systematically explore the regulatory network in myogenesis. Our study reveals an epigenetic-mediated myogenic regulatory mechanism and provides insights into the roles of non-coding RNA in myogenesis.

Methods

Animals and sample collection

The Sichuan white rabbit (SWR) used in this study were raised on the farm of Sichuan Animal Sciences Academy (Chengdu, Sichuan Province) under standard and uniform housing conditions (temperature: 22–26 °C; humidity: 60–70%). All the animals were healthy male rabbits (i.e., siblings) with similar body conditions. The animals were fed twice per day with formula diets containing 1.2% crude protein, 16% crude fiber, 8% crushed ash, 0.6% calcium, 1.2% lysine, 0.4% phosphorus, and 0.6% sodium chloride and had ad libitum access to water. The longissimus dorsi of rabbits is no longer expressed 27 days after delivery. For long non-coding RNA sequencing, longissimus muscle tissues were collected at 21 days, 24 days and 27 days after birth, with three biological replicates for each time point. Hereinafter, we used “21 d”, “24 d” and “27 d” to represent the above time points (Fig. 1A). For the mRNA and miRNA sequencing samples, muscle tissues were collected at 0 day, 1 month and 6 months after delivery to represent the fetal, child and adult stages separately, and were marked as “0 d”, “1 mon” and “6 mon” hereinafter. Due to the limitations of sample collection, we obtained three biological replicates for 0 d, one biological replicate for 1 month and one biological replicate for 6 months (Fig. 1A). Muscle tissues were collected immediately after slaughter. The tissues were cut into small pieces and rinsed with PBS. The collected tissues were stored in liquid nitrogen and then transferred to − 80 °C for subsequent high-throughput sequencing.

Fig. 1
figure 1

Experimental design and overview of the data. A Schematic diagram of the experiment and data generation. B Length distribution of sRNAs. C Classification of miRNAs. Colors represent different types of sRNAs. D Classification of lncRNAs. Colors represent different types of genomic annotations. E The intersection of lncRNAs predicted by the coding potential calculator (CPC), coding-non-coding index (CNCI) and protein families database (PFAM)

Full size image

RNA extraction, library construction and sequencing

Total RNA from muscle was extracted using TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer’s protocol. After purification, the quality was checked via agarose gel electrophoresis and a NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). The RNA was measured with an Agilent 2100 RNA 6000 Nano Kit (Agilent Technologies, Waldbronn, Germany). All the samples had high-quality RNA with an RNA integrity number (RIN) > 6. For mRNAs and lncRNAs, strand-specific sequencing libraries were constructed via the ribosomal RNA (rRNA) removal method following the manufacturer’s instructions using the Illumina Standard RNA Sample Library Preparation Kit (Illumina, San Diego, CA, USA). For miRNAs, RNA molecules ranging from 18–30 nt in size were enriched from total RNA by the polyacrylamide gel electrophoresis (PAGE). The 3′ adaptors were then added, followed by enrichment of RNAs with lengths of 36–44 nt and ligation of 5’ adaptors to the RNAs. RNAs were converted to cDNA, and PCR amplification was performed for library construction. All the above libraries were sequenced on the Illumina NovaSeq 6000 platform, and generated 150-bp paired-end (PE150) reads.

Data pre-processing

All the raw RNA-seq (i.e., mRNA, miRNA and lncRNA) reads were filtered using Trimmomatic (v0.36) software [15]. Specifically, adapters and reads of low quality, in which more than half of the bases had quality < = 20 or more than 10% of the bases were missing (Ns > 10%) were discarded to obtain clean reads. The clean reads were assessed with FastQC (v0.11.9) (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) for quality check.

For the mRNA sequencing reads, the clean reads were mapped to the rabbit reference genome OryCun2.0 (GCA_000003625.1) by HISAT2 software (v2.1.0) [16] with the following parameters: −n-ceil: L, 0, 0.15; −mp: MX = 6, MN = 2; −np: 1; −rdg: 5, 3; −efg: 5, 3; −score-min: L, 0, − 0.2. Mapped reads were quantified with the FeatureCounts program from the Subread suite (v2.0) [17]. The raw counts of genes were normalized to transcripts per million (TPM) with the in-house script for subsequent analysis.

For miRNA sequencing, clean reads were aligned to the rabbit (OryCun2.0) mature miRNA database in miRBase (v20) [18] and matched with known miRNAs, ribosomal RNA (rRNA), small nuclear RNA (snRNA), small nucleolar RNA (snoRNA), and transfer RNA (tRNA) sequences. miRNA abundance was quantified as counts and normalized to the TPM. Unmapped clean reads were further aligned to the rabbit reference genome OryCun2.0 using bowtie2 (v2.4.4) [19] software. To remove tags derived from protein-coding genes, repeat sequences, ncRNA, rRNA, tRNA, snRNA, and snoRNA, small RNA tags were mapped to RepeatMasker software [20] and the Rfam database. The distribution of the alignments was summarized using the software miREvo [21]. Thek novel miRNAs were predicted by miRDeep2 [22].

For the lncRNA sequencing reads, the clean reads were mapped to the rabbit reference genome via the HISAT2 (v2.1.0) program [16]. Then alignments were transferred to StringTie (v1.3.3) [23] and Cuffcompare (v2.2.1) [24] software for transcript assembly.

Differential expression analysis

For each type of RNA-seq data in this study, we performed differential expression analysis via the R package DESeq2 (v1.32.0) [25] to identify differentially expressed (DE) genes, miRNAs and lncRNAs between time points (mRNA-seq: 0 d vs. 1 mon, 0 d vs. 6 mon, and 1 mon vs. 6 mon; miRNA-seq: 0 d vs. 1 mon, 0 d vs. 6 mon, and 1 mon vs. 6 mon; lncRNA-seq: 21 d vs. 24 d, 21 d vs. 27 d, and 24 d vs. 27 d). In each comparison, the former group (e.g., “0 d” in “0 d vs. 1 mon”) was used as the control group when we mentioned up-regulated and down-regulated expressions/genes. The matrix of normalized TPM was used for DE analysis. The P-values were corrected via the Benjamini–Hochberg method. Corrected P-value < 0.05 and |log2FC| > 1 were set as the thresholds.

Target gene prediction for differentially expressed MiRNAs and LncRNAs

Target gene prediction was performed between three groups of differentially expressed miRNAs and lncRNAs. The prediction of target genes of miRNAs was performed by miRanda [26]. For DE lncRNAs, target genes were predicted by the positional relationship (co-location) and expression correlation (co-expression) of lncRNAs with protein-coding genes. lncRNA–mRNA co-regulated pairs (Pearson’s correlation coefficient > 0.8 and P- value < 0.05) were screened for Gene Ontology (GO) analysis.

Functional enrichment analysis

For DE RNAs (i.e., mRNAs, miRNAs and lncRNAs) from all comparisons, we conducted Gene Ontology (GO) enrichment analysis with the R package Goseq (v3.19) [27]. We built a reference database for rabbits with the Oryctolagus_cuniculus GTF file (Ensembl genome browser 113). Additionally, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was also performed by KOBAS (https://bio.tools/kobas). P-values were calculated by hyper-geometric test. The above biological processes and pathways were considered statistically significant with P-values < 0.05.

Results

Summary of the mRNA, MiRNA and LncRNA sequencing data

A total of 19 RNA sequencing datasets were generated in this study, including 5 mRNA-seq, 5 miRNA-seq and 9 lncRNA-seq datasets. Each sequencing data contained 3 developmental time points (Fig. 1A, Table S1). After the quality control of the raw reads, we obtained 282,575,600 mapped paired-end reads with an average mapping rate of 85.49% (ranging from 82.52 to 90.28%) for the mRNA datasets (Table S2), 38,485,181 mapped paired-end reads with an average mapping rate of 91.72% (ranging from 85.36 to 95.97%) for the miRNA datasets (Table S3) and 794,268,228 mapped paired-end reads with an average mapping rate of 90.51% (ranging from 89.37 to 91.21%) for the lncRNA datasets (Table S4). In addition, the length of all the miRNAs was 18–35 nt, most of which were 21–23 nt (Fig. 1B), and rRNA accounted for less than 1.35% (Fig. 1C). The classification of lncRNAs is shown in Fig. 1D, with 48.7% of the reads were lncRNAs (i.e., antisense lncRNAs and lncRNAs). We identified 17,959 lncRNA transcripts from the intersection of the coding potential calculator (CPC), coding-non-coding index (CNCI), and protein families database (Pfam) (Fig. 1E).

Differential expression analysis of mRNAs, MiRNAs and LncRNAs

To explore the characteristics of different developmental stages, we first verified the reliability of the identified RNAs with principal component analysis (PCA) and hierarchical cluster analysis (HCA) (Fig. 2). Although the sample size is limited, samples in the same group cluster together roughly. Next, we performed differential expression analysis between groups for each sequencing dataset separately. The former group (e.g., “0 d” in “0 d vs. 1 mon”) was used as the control group when we mentioned up-regulated and down-regulated expressions/genes.

Fig. 2
figure 2

Principal component analysis (PCA) (A-C) and hierarchical cluster analysis (HCA) (D-F) of samples in mRNA, miRNA and lncRNA sequencing

Full size image

Among the comparisons, we detected a total of 2,995 differentially expressed (DE) genes from RNA-seq, 1,211 DE-lncRNAs from lncRNA-seq and 305 DE-miRNAs from miRNA-seq (Fig. 3A-C). For mRNA-seq, 44 genes were detected in all three comparisons, whereas 308, 231 and 1,777 DEGs were exclusively differentially expressed in the comparison of “0 d vs. 1 mon”, “1 mon vs. 6 mon” and “0 d vs. 1 mon”, respectively (Fig. 3D). Besides, there was small difference in the number of up and down-regulated DEGs across the three comparisons (Fig. 3A). For miRNA-seq, we detected the greatest number of DEGs in the comparison of “0 d vs. 6 mon”, followed by “0 d vs. 1 mon” and “1 mon vs. 6 mon” (Fig. 3B), which included 93 and 74 DE-miRNAs that were up- and down-regulated in the comparison of “0 d vs. 6 mon”, 39 and 30 DE-miRNAs that were up- and down-regulated in the comparison of “0 d vs. 1 mon”, and 53 and 16 DE-miRNAs that were up- and down-regulated in the comparison of “1 mon vs. 6 mon” (Fig. 3B). The greatest number of unique DE-miRNAs was detected in the comparison of “0 d vs. 6 mon”, followed by “1 mon vs. 6 mon” and “0 d vs. 1 mon” (Fig. 3E). We also observed that the number of up-regulated DE-miRNAs was generally greater than that of down-regulated DE-miRNAs (Fig. 3B). In terms of lncRNA-seq, we found the greatest number of DE-lncRNAs in the comparison of “21 d vs. 27 d”, including 389 up-regulated and 703 down-regulated DE-lncRNAs (Fig. 3C). Additionally, only 8 DE-lncRNAs were detected in all three comparisons, whereas 19,865 and 71 lncRNAs were uniquely differentially expressed in the comparison of “21 d vs. 24 d”, “21 d vs. 27 d” and “24 d vs. 27 d”, respectively (Fig. 3F).

Fig. 3
figure 3

Differentially expressed DE genes, DE-miRNAs and DE-lncRNAs between different developmental stages. A-C The number of up- and down-regulated DEGs (A), DE-miRNAs (B) and DE-lncRNAs (C) in each pairwise comparison. D-F Venn diagram of DEGs (D), DEMs (E) and DELs (F) in each pairwise comparison

Full size image

Functional enrichment of differentially expressed mRNAs, MiRNAs and LncRNAs

To investigate the functions of the differentially expressed genes, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses for comparisons of each sequencing dataset separately (Tables S5-S9). We mainly focused on the comparison of which harbored the greatest number of DE-RNAs. We found the highest number of DE genes in the “1 mon vs. 6 mon” comparison, which included 1,178 up-regulated and 954 down-regulated DE genes (Fig. 4A). We then performed functional enrichment for these genes and found that the functions of the up- and down-regulated DEGs were different (Fig. 4B and C). For example, up-regulated DE genes were associated with the GO terms “action filament-based process” (e.g., CACNA2D1, AKAP9 and RYR2) and “anatomical structure development” (e.g., HOMER1, MYOC and GSK3B) (Fig. 4B, Table S5 and Table S6), whereas down-regulated DEGs were largely associated with macromolecule metabolic processes (e.g., ATF3, JAK2 and IDE) (Fig. 4C).

Fig. 4
figure 4

Differential expression analysis of mRNAs between 1 month and 6 months. A Volcano plots of up-regulated and down-regulated DEGs between 1 month and 6 months. The blue dots represent non-significantly DEGs. B, C The top 10 GO terms enriched from down- and up-regulated DEGs in the “1 mon vs. 6 mon” comparison. The blue bar represents terms of down-regulated genes (B) and the red bar represents terms of up-regulated genes (C)

Full size image

For miRNA-seq, the hierarchical clustering of the DE-miRNAs heatmap showed their expression dynamics at three time points (Fig. 5A). We detected the greatest number of DE-miRNAs in the comparison of “0 d vs. 6 mon” (Fig. 5B). The top 10 GO terms were associated with basic biological processes (Table S7), such as “transport” (e.g., SPTBN2 and TGFB3), “single-organism transport” (e.g., CAMK1 and CLCN3) and “establishment of localization” (e.g., CHMP7 and RTN2). Besides, we also found several muscle-related GO terms (Fig. 5C), such as “positive regulation of growth” (e.g., SMO, TBX2 and WNT3 A) and “striated muscle cell differentiation” (e.g., MYOG, MYPN and EDN1). For lncRNA-seq, we examined the expression of the DE-lncRNAs and observed large changes in the comparison of “21 d vs. 27 d” (Fig. 6A), which obtained the greatest number of DE-lncRNAs (Fig. 6B). Surprisingly, we found that the top 10 GO terms from this comparison were closely related to muscle development (Fig. 6C, Table S8, Table S9), such as “muscle structure development” (e.g., TNNC1, MYF6 and MYL3), “muscle organ development” (e.g., ACTN3, IGF1 and MYLK2) and “skeletal muscle organ development” (e.g., MSTN, MYOG and CXCL9).

Fig. 5
figure 5

Differential expression analysis of miRNAs between 0 day and 6 months. A Heatmap showing the expression of DE-miRNAs in all the miRNA samples. B Volcano plots of up-regulated and down-regulated between 0 day and 6 months. The blue dots represent non-significant DEGs. C Muscle-related GO terms enriched from DE-miRNAs in the “0 d vs. 6 mon” comparison

Full size image
Fig. 6
figure 6

Differential expression analysis of lncRNAs between 21 days and 27 days. A Heatmap showing the expression of DE-lncRNAs in all the lncRNA samples. B Volcano plots of up-regulated and down-regulated between 21 days and 27 days. The blue dots represent non-significant DEGs. C The top 10 GO terms enriched DE-lncRNAs in the “21 d vs. 27 d” comparison

Full size image

Interaction of differentially expressed LncRNAs between key time points

To address how miRNAs and lncRNAs interact with their target genes (mRNAs) to regulate rabbit muscle development, we predicted potential target genes and investigated their functions. The results of the above differential expression analysis, we found the greatest number of DE-lncRNAs in the comparison of “21 d vs. 27 d”, and their function was associated with muscle growth. Previous studies have confirmed that lncRNAs regulate the expression of neighboring protein-coding genes through cis-acting mechanisms [28]. Thus, we performed co-expressed gene analysis between lncRNAs and mRNAs with a correlation coefficient > 0.95 as the threshold and concentrated on this comparison. The directed acyclic graph (DAG) plot for GO terms enriched from the target genes of the DE-lncRNAs showed complex regulatory network (Fig. 7).

Fig. 7
figure 7

Directed acyclic graph (DAG) plot of the enriched GO terms for the target genes of the DE-lncRNAs

Full size image

Discussion

The growth and development of muscle is a dynamic and complex process that largely affects the meat production of livestock. In this study, we systemically investigated the complex regulatory landscape of muscle development by conducting a comprehensive sequencing analysis of mRNAs, miRNAs, and lncRNAs. Myogenesis is a highly coordinated developmental process. Our transcriptome data included three time points that can cover different stages. For mRNAs and miRNAs, we selected 0 day, 1 month and 6 months to represent the fetal, child, and adult periods of rabbits, respectively [29, 30]. The total number of fibers is fixed during the first postnatal month of myogenesis in rabbits [30]. In particular, longissimus dorsi is no longer expressed 27 days after birth [31]. Hence we selected 21 days, 24 days and 27 days to explore lncRNA changes. By performing PCA and HCA analyses, we found that samples in the same group tended to cluster together (Fig. 2) and verified the reliability of the data for downstream analysis.

Differential expression analysis revealed that gene expression changes corresponded to distinct developmental stages. The comparison of mRNA and miRNA datasets revealed a total of 2,995 DE genes and 305 DE-miRNAs, respectively, with notable variations in expression patterns. The lncRNA-seq analysis identified 1,211 DE-lncRNAs, and the most significant changes were observed between 21 d and 27 d. We found gene PAX7 and MEF2 were down-regulated in the comparison of “21 d vs. 27 d” (Table S9) only. PAX7 is expressed in satellite cells during the postnatal development of rabbits and plays an important role in muscle regeneration and repair. Similar functions have also been reported in mice and chickens [32]. MEF2 serves as a transcription factor that can regulate muscle fiber identity and maintenance. The down-regulation of these two genes suggested that the number of fibers stopped increasing in the first postnatal month in rabbits [30]. Moreover, muscle also can change their functional characteristics in response to the physiological stage (i.e., growing, maintaining and senescing) of animals. For instance, the dynamic expression of miRNAs and lncRNAs in 7 different periods in goats was reported to affect skeletal muscle development [33, 34]. In black Muscovy duck, differentially expressed genes like MyoG, FBXO1, MEF2 A, and FoxN2 in leg muscle were enriched in growth-related biological processes [35]. These findings suggest that myogenesis is relatively conserved among species.

Functional enrichment analysis of the DE mRNAs and miRNAs shed light on their biological roles, with a particular emphasis on the “1 mon vs. 6 mon” comparison for mRNAs. This analysis revealed that the up-regulated DE genes were involved in processes such as actin filament-based movement (e.g., CACNA2D1, AKAP9 and RYR2) and anatomical structure development (e.g., HOMER1, MYOC and GSK3B) (Fig. 4B), while down-regulated DE genes were predominantly associated with macromolecule metabolism (e.g., ATF3, JAK2 and IDE) (Fig. 4C). miRNA-seq analysis revealed a dynamic expression pattern, with the “0 d vs. 6 mon” comparison yielding the greatest number of DE-miRNAs (Fig. 3B). The GO terms associated with these miRNAs pointed towards fundamental biological processes and muscle-specific functions (Fig. 5C), such as positive regulation of growth (e.g., SMO, TBX2 and WNT3 A) and striated muscle cell differentiation (e.g., MYOG, MYPN and EDN1). Notably, we detected the greatest number of DE-lncRNAs in the comparison of “21 d vs. 27 d” (Fig. 3C), and the functions of these DE-lncRNAs were closely associated with muscle growth (Fig. 6C). MyoG and Myf5, these two important MRFs were largely enriched in multiple muscle-related GO terms, such as “muscle organ development” and “striated muscle tissue development” (Fig. 6, Table S9) in the comparison of “21 d vs. 27 d”. In addition, interaction analysis between DE-lncRNAs and their target genes unveiled a complex regulatory network. Our co-expression analysis used a stringent threshold of correlation > 0.95, and the subsequent DAG plot of GO terms enriched from the target genes clearly showed the intricate interplay between lncRNAs and mRNAs in modulating muscle development (Fig. 7).

Furthermore, the above mentioned target genes play important roles in muscle development. Myogenic cell specification and differentiation are determined by the master transcription factor MyoD in concert with other myogenic regulatory factors (MRFs) [36]. In particular, MyoG plays a central role in the terminal differentiation of myoblasts into mature muscle fibers [37]. Besides, MyoG is regulated by growth factor signaling pathways, especially the IGF1 pathway. IGF1 can enhance the expression of above-mentioned MRFs by activating the PI3 K/Akt pathway, which works collectively to promote muscle growth and repair [38]. Previous study reported that the interplay between IGF signaling and MRFs can also be modulated by miRNAs like miR- 1 and miR- 133 [39]. These results suggest dramatic lncRNA changes in the expression of myosin heavy chain isoforms in the first postnatal month. Besides, gene MYOC encodes a protein that belongs to the olfactomedin family and is expressed in various tissues, including skeletal muscle. Although its specific function in muscle is not well characterized, it is known to be involved in the development and maintenance of muscle mass, potentially through interactions with other muscle regulatory factors [40]. WNT3 A belongs to the Wnt signaling pathway, this pathway regulates muscle formation and the maintenance of adult tissue homeostasis [41, 42]. Besides, JAK2 is part of the JAK-STAT signaling pathway. This pathway is involved in the regulation of muscle growth and differentiation, with JAK2 being a key mediator of these processes. Mutations in JAK2 have been associated with myeloproliferative disorders, which can have secondary effects on muscle function [43]. The RYR2 gene encodes the type 2 ryanodine receptor, a calcium release channel found in the sarcoplasmic reticulum of muscle cells. It is essential for the regulation of calcium ions during muscle contraction and relaxation. Dysfunction of RYR2 has been implicated in various muscle diseases, including malignant hyperthermia and central core disease, highlighting its importance in muscle function [44]. Overall, these genes and their protein products are integral to the complex processes of muscle growth, development, and maintenance.

The present study has several limitations. First, the small sample size for each time point has potential impact on expression patterns, although unrelated factors were controlled. Previous studies have reported that sufficient biological replicates enable to capture of authentic biological variability and avoid technical artifacts [45, 46]. Proper statistical models or decreased multiple comparisons can also reduce false positive errors. In addition, quantitative polymerase chain reaction (qPCR) is one of the most reliable methods for validating gene expression [47] and has been widely used in transcriptome research. Due to the restriction of sample collection, our study only explored expression patterns during muscle development based on sequencing data. In future work, we will cautiously consider the impact of sample size and add a qPCR experiment to validate the expression of candidate genes. Second, tissues were collected from males only. Thus, sex-based differences in muscle expression cannot be addressed. Future research could benefit from the increased sample size of both male and female rabbits and the use of different layers of multi-omics data (e.g., chromatin accessibility and DNA methylation). It will greatly enhance the understanding of the regulation of rabbit muscle development.

Conclusions

This study provides a comprehensive overview of the transcriptomic changes that occur during muscle development. The identification and functional annotation of DE genes, miRNAs, and lncRNAs provide valuable insights into the molecular mechanisms underlying this process. These findings pave the way for investigations into the role of non-coding RNAs in muscle biology.

Data availability

The datasets generated in this study are available in the Sequence Read Archive (SRA) with the primary accession code PRJNA1173484 of the NCBI database.

Abbreviations

miRNAs:

Micro RNA

lncRNAs:

Long non-coding RNAs

SWR:

Sichuan white rabbit

PBS:

Phosphate buffered saline

TPM:

Transcripts per million

DE:

Differential expressed

GO:

Gene Ontology

CPC:

Coding potential calculator

CNCI:

Coding-non-coding index

Pfam:

Protein families database

KEGG:

Kyoto Encyclopedia of Genes and Genomes

DAG:

Directed acyclic graph

qPCR:

quantitative polymerase chain reaction

References

  1. Dalle A, Szendrő Z. The role of rabbit meat as functional food. Meat Sci. 2011;88:319–31.

    Article  Google Scholar 

  2. Li S, Zeng W, Li R, Hoffman LC, He Z, Sun Q, Li H. Rabbit meat production and processing in China. Meat Sci. 2018;145:320–8.

    Article  PubMed  Google Scholar 

  3. Zhang X, Zhang C, Yang C, Kuang L, Zheng J, Tang L, Lei M, Li C, Ren Y, Guo Z, et al. Circular RNA, microRNA and Protein Profiles of the Longissimus Dorsi of Germany ZIKA and Sichuan White Rabbits. Front Genet. 2021;12:777232.

  4. Guo Y, Geng W, Chen Z, Zhi Y, Zhang K, Li Z, Li G, Kang X, Tian W, Li H, et al. LncRNA LncMGR regulates skeletal muscle development and regeneration by recruiting CDK9 and sponging MiRNAs. Int J Biol Macromol. 2024;266:131049.

    Article  CAS  PubMed  Google Scholar 

  5. Yang Y, Wu J, Liu W, Zhao Y, Chen H. The function and regulation mechanism of Non-Coding RNAs in muscle development. Int J Mol Sci. 2023;24:14534.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Wienholds E, Kloosterman WP, Miska E, Alvarez-Saavedra E, Berezikov E, De Bruijn E, Horvitz HR, Kauppinen S. Plasterk R H A. MicroRNA expression in zebrafish embryonic development. Science. 2005;309:310–1.

    Article  CAS  PubMed  Google Scholar 

  7. Ahmed A, Ward NJ, Moxon S, Lopez-Gomollon S, Viaut C, Tomlinson ML, Patrushev I, Gilchrist MJ, Dalmay T, Dotlic D, et al. A Database of microRNA Expression Patterns in Xenopus Laevis. Plos One. 2015;10(10):e13813.

  8. Darnell DK, Kaur S, Stanislaw S, Konieczka JK, Yatskievych TA, Antin PB. MicroRNA expression during chick embryo development. Dev Dyn. 2006;235:3156–65.

    Article  CAS  PubMed  Google Scholar 

  9. Horak M, Novak J, Bienertova-Vasku J. Muscle-Specific MicroRNAs in skeletal muscle development. Dev Biol. 2016;410:1–13.

    Article  CAS  PubMed  Google Scholar 

  10. Chen LL, Kim VN. Small and long Non-Coding RNAs: past, present, and future. Cell. 2024;187:6451–85.

    Article  CAS  PubMed  Google Scholar 

  11. Cesana M, Cacchiarelli D, Legnini I, Santini T, Sthandier O, Chinappi M, Tramontano A, Bozzoni I. A long noncoding RNA controls muscle differentiation by functioning as a competing endogenous RNA. Cell. 2011;147:358–69.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Wilusz JE, Sunwoo H, Spector DL. Long noncoding RNAs: functional surprises from the RNA world. Genes Dev. 2009;23:1494–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Wang X, Song X, Glass CK, Rosenfeld MG. The long arm of long noncoding RNAs: roles as sensors regulating gene transcriptional programs. Cold Spring Harb Perspect Biol. 2011;3:a003756.

    Article  PubMed  PubMed Central  Google Scholar 

  14. Gong C, Li Z, Ramanujan K, Clay I, Zhang Y, Lemire-Brachat S, Glass DJ. A long Non-Coding RNA, LncMyoD, regulates skeletal muscle differentiation by blocking IMP2-Mediated mRNA translation. Dev Cell. 2015;34:181–91.

    Article  CAS  PubMed  Google Scholar 

  15. Bolger AM, Lohse M, Usadel B. Trimmomatic: A flexible trimmer for illumina sequence data. Bioinformatics. 2014;30:2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Kim D, Langmead B, Salzberg SL. HISAT: A fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Liao Y, Smyth GK, Shi W. The Subread aligner: fast, accurate and scalable read mapping by Seed-and-Vote. Nucleic Acids Res. 2013;41:e108–108.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Kozomara A, Birgaoanu M, Griffiths-Jones S. MiRBase: from MicroRNA sequences to function. Nucleic Acids Res. 2019;47:D155–62.

    Article  CAS  PubMed  Google Scholar 

  19. Langmead B, Salzberg SL. Fast Gapped-Read alignment with bowtie 2. Nat Methods. 2012;9:357–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Tarailo-Graovac M, Chen N. Using RepeatMasker to identify repetitive elements in genomic sequences. Curr Protoc Bioinformatics. 2009;25:4.10.1–14.

  21. Wen M, Shen Y, Shi S, Tang T. MiREvo: an integrative MicroRNA evolutionary analysis platform for next-Generation sequencing experiments. BMC Bioinformatics. 2012;13:140.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Friedländer MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel MicroRNA genes in seven animal clades. Nucleic Acids Res. 2012;40:37–52.

    Article  PubMed  Google Scholar 

  23. 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:290–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, Van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28:511–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Love MI, Huber W, Anders S. Moderated Estimation of fold change and dispersion for RNA-Seq data with DESeq2. Genome Biol. 2014;15:550.

    Article  PubMed  PubMed Central  Google Scholar 

  26. John B, Enright AJ, Aravin A, Tuschl T, Sander C, Marks DS. Human MicroRNA targets. PLoS Biol. 2004;2:e363.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-Seq: accounting for selection Bias. Genome Biol. 2010;11:R14.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Bu Q, Hu Z, Chen F, Zhu R, Deng Y, Shao X, Li Y, Zhao J, Li H, Zhang B, et al. Transcriptome analysis of long Non-Coding RNAs of the nucleus accumbens in Cocaine-Conditioned mice. J Neurochem. 2012;123:790–9.

    Article  CAS  PubMed  Google Scholar 

  29. Jing J, Jiang X, Zhu C, Zheng Q, Ji Q, Yin H, Huang J, Zhu Y, Wang J, Qin S, et al. Dynamic changes of MiRNAs in skeletal muscle development at new Zealand rabbits. BMC Genomics. 2021;22:577.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Picard B, Lefaucheur L, Berri C, Duclos MJ. Muscle fibre ontogenesis in farm animal species. Reprod Nutr Dev. 2002;42:415–31.

    Article  PubMed  Google Scholar 

  31. Gondret F, Lefaucheur L, D’Albis A, Bonneau M. Myosin isoform transitions in four rabbit muscles during postnatal growth. J Muscle Res Cell Motil. 1996;17:657–67.

    Article  CAS  PubMed  Google Scholar 

  32. Relaix F, Montarras D, Zaffran S, Gayraud-Morel B, Rocancourt D, Tajbakhsh S, Mansouri A, Cumano A, Buckingham M. Pax3 and Pax7 have distinct and overlapping functions in adult muscle progenitor cells. J Cell Biol. 2005;172:91–102.

    Article  PubMed  Google Scholar 

  33. Ling Y, Zheng Q, Sui M, Zhu L, Xu L, Zhang Y, Liu Y, Fang F, Chu M, Ma Y, et al. Comprehensive analysis of LncRNA reveals the Temporal-Specific module of goat skeletal muscle development. Int J Mol Sci. 2019;20:3950.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Ling Y, Zheng Q, Jing J, Sui M, Zhu L, Li Y, Zhang Y, Liu Y, Fang F, Zhang X. RNA-Seq Reveals miRNA Role Shifts in Seven Stages of Skeletal Muscles in Goat Fetuses and Kids. Front Genet. 2020;11:684.

  35. Hu Z, Cao J, Liu G, Zhang H, Liu X. Comparative transcriptome profiling of skeletal muscle from black muscovy Duck at different growth stages using RNA-Seq. Genes. 2020;11:1228.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Yu X, Zhang Y, Li T, Ma Z, Jia H, Chen Q, Zhao Y, Zhai L, Zhong R, Li C, et al. Long Non-Coding RNA Linc-RAM enhances myogenic differentiation by interacting with myod. Nat Commun. 2017;8:14016.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Hasty P, Bradley A, Morris JH, Edmondson DG, Venuti JM, Olson EN, Klein WH. Muscle deficiency and neonatal death in mice with a targeted mutation in the Myogenin gene. Nature. 1993;364:501–6.

    Article  CAS  PubMed  Google Scholar 

  38. Tapscott SJ. The circuitry of a master switch: myod and the regulation of skeletal muscle gene transcription. Development. 2005;132:2685–95. https://doiorg.publicaciones.saludcastillayleon.es/10.1242/dev.01874.

    Article  CAS  PubMed  Google Scholar 

  39. Chen JF, Mandel EM, Thomson JM, Wu Q, Callis TE, Hammond SM, Conlon FL, Wang DZ. The role of microRNA-1 and microRNA-133 in skeletal muscle proliferation and differentiation. Nat Genet. 2006;38:228–33.

    Article  CAS  PubMed  Google Scholar 

  40. Xie S, Liu Q, Fu C, Chen Y, Li M, Tian C, Li J, Han M, Li C. Molecular regulation of Porcine skeletal muscle development: insights from research on CDC23 expression and function. Int J Mol Sci. 2024;25:3664.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Maltzahn J, Chang NC, Bentzinger CF, Rudnicki MA. Wnt signaling in myogenesis. Trends Cell Biol. 2012;22:602–9.

    Article  Google Scholar 

  42. Girardi F, Le Grand F. Wnt Signaling in Skeletal Muscle Development and Regeneration. Prog Mol Transl Sci. 2018;153:157–79.

  43. Liu J, Jing X, Gan L, Sun C. The JAK2/STAT3 signal pathway regulates the expression of genes related to skeletal muscle development and energy metabolism in mice and mouse skeletal muscle cells. Biosci Biotechnol Biochem. 2012;10:1866–70.

    Article  Google Scholar 

  44. Sleiman Y, Lacampagne A, Meli AC. Ryanopathies and RyR2 dysfunctions: can we further Decipher them using in vitro human disease models?? Cell Death Dis. 2021;12:1–19.

    Article  Google Scholar 

  45. Khang TF, Lau CY. Getting the most out of RNA-seq data analysis. Peer J. 2015;3:e1360.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Stark R, Grzelak M, Hadfield J. RNA sequencing: the teenage years. Nat Rev Genet. 2019;20:631–56.

    Article  CAS  PubMed  Google Scholar 

  47. Everaert C, Luypaert M, Maag JL, et al. Benchmarking of RNA-sequencing analysis workflows using whole-transcriptome RT-qPCR expression data. Sci Rep. 2017;7:1559.

    Article  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This work was financially supported by grants from Natural Science Foundation of Sichuan Province (2023 NSFSC0171); Sichuan Rabbit Innovation Team of National Modern Agricultural Industry Technology System (sccxtd-2024-26), Sichuan Province science and technology plan project (2021YFYZ0009), Basic Research Funds for Public Welfare Research Institutes in Sichuan Province in 2023-Construction and application of genome selection method for growth and development traits of Rex rabbits (2023 JDKY0033), Supported by the Earmarked Fund for China Agriculture Research System (CARS-43-A-3).

Author information

Authors and Affiliations

  1. Sichuan Academy of Science Academy, Chengdu, 610066, China

    Xiangyu Zhang, Dengping Huang, Shangjun Yang, Min Zhang & Qin Yin

  2. Animal Breeding and Genetics Key Laboratory of Sichuan Province, Chengdu, 610066, China

    Xiangyu Zhang, Dengping Huang, Shangjun Yang, Min Zhang & Qin Yin

  3. Sichuan Academy of Grassland Sciences, Chengdu, 611743, China

    Kai Zhang

Authors
  1. Xiangyu Zhang
    View author publications

    You can also search for this author inPubMed Google Scholar

  2. Kai Zhang
    View author publications

    You can also search for this author inPubMed Google Scholar

  3. Dengping Huang
    View author publications

    You can also search for this author inPubMed Google Scholar

  4. Shangjun Yang
    View author publications

    You can also search for this author inPubMed Google Scholar

  5. Min Zhang
    View author publications

    You can also search for this author inPubMed Google Scholar

  6. Qin Yin
    View author publications

    You can also search for this author inPubMed Google Scholar

Contributions

X.Y.Z., M.Z., and Q.Y. conceived, designed and supervised the study. X.Y.Z. and K.Z. performed bioinformatic analysis. X.Y.Z, K.Z, D.P.H., S.J.Y., Q.Y. and M.Z contributed to the sample collection and resource generation. X.Y.Z. wrote and revised the manuscript. All authors reviewed, edited, and approved the final manuscript.

Corresponding authors

Correspondence to Min Zhang or Qin Yin.

Ethics declarations

Ethics approval and consent to participate

The experiment was approved by the Animal Care and Use Committee of the Sichuan Animal Sciences Academy. The identification number was not required since the commercial animal sampling was approved. The application form for welfare and ethical review in animal experimentation was approved by the Sichuan Animal Sciences Academy.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

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/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhang, X., Zhang, K., Huang, D. et al. Comprehensive transcriptome of muscle development in Sichuan white rabbit. BMC Genom Data 26, 32 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12863-025-01322-5

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12863-025-01322-5

Keywords

BMC Genomic Data

ISSN: 2730-6844

Contact us