The Comparison Transcriptome Between High and Low Growth Ratio of Meat-Type Rabbit Revealed the Primary Pathway Related to Muscle Development

1. Introduction

In 2021, 570 million meat rabbits provided 900,000 tons of rabbit meat in the world,and the rabbit (Oryctolagus cuniculus) is regarded as an ideal meat-producing animal in theory, because of its short life cycle, short gestation period, large number of litters, and high feed conversion ability [1]. It contains high protein, essential amino acids and a variety of trace elements. In addition, its sensitization and low cholesterol levels have also been demonstrated [1,2]. The thigh muscle and the longissimus dorsi muscle both are the focuses of molecular breeding. With the development and application of molecular biology and bioinformatics, it is possible to reveal the basis of meat rabbit production performance at the molecular level [3]. Therefore, finding muscle growth and development of meat rabbit genes has become particularly significant.

In recent years, RNA sequencing (RNA-Seq) has become a popular choice for gene expression studies [4,5,6,7]. Moreover, sequencing technology is continuously improving, and the costs have decreased. In particular, RNA-Seq has been used to simultaneously measure multitudinous of gene expression levels and to discover differentially expressed genes between two or more groups of a variety of applications [4,8,9]. Therefore, in an increasing number of studies, RNA sequencing has been used to screen and identify key regulatory genes in biological processes such as growth, development, metabolism, and diseases [10,11,12,13,14].

Sichuan linen (Ma) rabbits are an important local domestic rabbit resource with strong adaptability and excellent meat quality, suitable for breeding in hilly areas; Checkered Giant (Ju) rabbits have huge body size and better meat performance [15]. Existing studies have found that differential expression of some genes leads to changes in a subset of skeletal muscle characteristics (e.g., fiber type, oxidative capacity, and insulin binding). The expression of myosin heavy chain (MyHC) gene isoforms (I, IIa, IIx, IIb) [16] and actin (ACTN) [17] and troponin C (TNNC) [18] in skeletal muscle has been reported to be an important feature of a muscle fiber structure and function [19]. These genes such as myostatin (MSTN) [20], muscle regulation factors (MRFs), paired box (PAX), myogenic transcription factor 2C (MEF2C), and insulin-like growth factors (IGFs) [21] form a complex regulatory network involved in regulating physiological activities such as skeletal muscle growth, myofiber type conversion, and metabolism of intramuscular fat deposition [22]. Muscle transcriptional helper studies in a variety of animal models have preliminarily elucidated the transcriptional regulatory network for physiological processes such as muscle fiber type, muscle proliferation and hypertrophy [23]. The objective of the study is to explore the differences in transcriptomic profiles of thigh and longissimus dorsi muscle from Ma and Ju rabbit, thus revealing the hub-gene and primary signaling pathway caused the differences in muscle development between Ju and Ma rabbit.

2. Materials and Methods

2.1. Animals and Muscle Tissue Samples

Ma (n = 9) and Ju (n = 9) were feed with complete diet (Table 1). We collected thirty-six thigh muscle and longissimus dorsi muscle samples of two rabbit breeds (28-day-old, 56-day-old, and 84-day-old, male, same environmental condition and feeding) from the meat rabbit research base of Sichuan Academy of Animal Husbandry Science. All samples were rinsed three times with PBS. Then, the samples were divided into two groups. One group was immediately frozen in liquid nitrogen until RNA extraction was needed, and the other was sectioned into muscle blocks and stored in muscle-specific fixation solution (SAKURA, RAS,USA) to generate paraffin sections.

2.2. Paraffin Section Preparation and Analysis

Muscle tissue samples were collected from the thigh muscle and the longissimus dorsi muscle of rabbit at 28, 56, and 84 days of age. The tissue samples were fixed in 10% formalin in PBS at room temperature for 3 days and then dehydrated, embedded, and sectioned (thickness=4㎛) for routine H&E staining [24]. Images were collected using Olympus BX-53 under 10× eyepiece and 20× objective lens. There were 3 biological replicates in each stage. For each slice, the area of muscle fibers in 3 non-overlapping and non-lost visual fields were measured through image J. The ANOVA in IBM SPSS Statistics 27.0.1 was used for analyzing the average cross-sectional area （CSA）of each stage .

2.3. RNA Sample Detection and Library Construction

RNA integrity was assessed by using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Total RNA was used as an input material for the RNA sample preparations. Briefly, mRNA was purified from total RNA using poly-Toligo-attached magnetic beads. Fragmentation was carried out by using divalent cations under elevated temperature in First Strand Synthesis Reaction Buffer(5X). First strand cDNA was synthesized by using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H-). Second strand cDNA synthesis was subsequently performed by using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3’ ends of DNA fragments, Adaptor with hairpin loop structure were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 370~420 bp in length, the library fragments were purified with AMPure XP system. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. At last, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system.

2.4. Clustering and Quality Control of RNA-Seq Libraries

The clustering of the index-coded samples was performed on a cBot Cluster Generation System by using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) and according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina Novaseq platform and 150 bp paired-end reads were generated.

Raw data (raw reads) of the fastq format were first processed through fastp software. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low-quality reads from raw data. At the same time, Q20, Q30, and GC content the clean data were calculated. All the downstream analyses were based on the clean data with high quality.

2.5. Analysis of Sequencing Data

Reference genome and gene model annotation files were directly downloaded from genome website. Index of the reference genome was built by using Hisat2 v2.0.5 and paired-end clean reads were aligned to the reference genome by using Hisat2 v2.0.5. FPKM, expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced. FeatureCounts v1.5.0-p3 was used to count the reads mapped to each gene. Then, FPKM of each gene was calculated based on the length of the gene and reads mapped to this gene. Differential expression analysis of two conditions/groups (two biological replicates per condition) was performed by using the DESeq2 R package (1.20.0). The resulting P-values were adjusted by using the Benjamini and Hochberg’s approach for controlling the false discovery rate. Genes with an adjusted P-value <=0.05 found by DESeq2 were assigned as differentially expressed. Gene Ontology (GO) enrichment analysis of differentially expressed genes was implemented by the clusterProfiler R package, in which gene length bias was corrected. GO terms with the corrected P-value less than 0.05 were considered as significantly enriched by differential expressed genes. We used clusterProfiler R package to test the statistical enrichment of differential expression genes in KEGG pathways. We used the online Search Tool for the Retrieval of Interacting Genes (STRING database, Version 10.5; http://string-db.org/) to build a protein-protein interaction (PPI) network based on broad function of genome-wide data [25]. After carefully calculating, the threshold >0.150 was applied to build the PPI network, and the cytoscape software (version 3.7.0; http://cytoscape.org/) was used to visualize and analyze the biological networks. The nodes of the networks represented “hub genes” that highly interacted with other genes or proteins.

3. Results

3.1. The Changes in the Area of Muscle Fibers of Rabbits’ Skeletal Muscle at Different Stages

The body weight of Ju rabbit is significantly higher than that of Ma rabbit after 3-week age (Figure 1A). To better observe the changes in the characteristics of rabbits’ thigh muscle and longissimus dorsi muscle, we generated paraffin-sections of the two muscle tissue sites in three developmental stages (Figure 1B and C). We found that as they got older, their cross-sectional areas of muscle fibers increase significantly. On 56-day-old stage, the CAS differed significantly between the same muscle tissue of Ma and Ju (Figure 1D and E) (Table S1), but the area of thigh and longissimus dorsi muscles has no difference (Figure 1F). In the nest experiment, we explored the transcriptome profiles of thigh and longissimus dorsi muscle from two breed at 56-day-old.

3.2. Summarization of RNA-Seq

Twelve libraries of 56-day-old rabbits’ thigh and longissimus dorsi muscle were constructed. All libraries obtained raw data with an average of 40.25 to 47.07 million. Under strict filtering conditions, libraries contained 6.51G effective clean bases, accounting for approximately 98.01 % of the raw reads. Furthermore, the Q30 = 94.58 % of all libraries (Table 2).

Firstly, we analyzed the exons distribution of mRNAs. All candidate mRNAs were found to be mainly spliced by exons (77.33%), introns (7.93 %), and intergenic regions (14.74 %) (Figure 2A) (Table S3). After we calculated the expression values (FPKM) of all genes in each sample, and showed the distribution of gene expression levels in different samples by box plots (Figure 2B). The expression levels of co-expression genes in all samples were analyzed using Pearson’s correlation (Figure 2C), and the average correlation coefficient of duplicate samples within the same tissue and breed (r=0.950, p<0.01) was higher than that within breed (r=0.942, p<0.01). Average correlation coefficient within the breed (r=0.945, p<0.01) and within the same muscle (r=0.937, p<0.01) was higher than that between different breed and muscle (r=0.923, p<0.01). The cluster-heatmap of DEGs showed the expressive pattern of all DEGs was divided into 4 group: Ju_Leg, Ju_Long, Ma_Leg and Ma_Long (Figure 2D). The results suggested RNA-seq high reproducibility and reliability.

3.3. Differential Expression Analysis of mRNAs in Rabbit Skeletal Muscle

We subsequently performed differential expression analysis on genes in the thigh muscle and longissimus dorsi muscle of two rabbit breeds at 56 days post-natal. Venn diagram analysis revealed 11,492 and 11,728 co-expressed genes in the thigh and longissimus dorsi muscles, respectively, across both breeds (Figure 3A and B). Through comparative analysis between the Ma and Ju breeds using established thresholds (|log₂FC| > 0 and padj < 0.05; see Table S4), we identified 1,983 differentially expressed genes (DEGs) in thigh muscle tissue (785 upregulated and 1,198 downregulated) and 1,955 DEGs in longissimus dorsi muscle tissue (734 upregulated and 1,221 downregulated)(Figure 3C).Using established thresholds (|log₂FC| > 1 and padj < 0.05; see Table S5), we identified 284 EGs in thigh muscle tissue (90 upregulated and 194 downregulated) and 305 DEGs in longissimus dorsi muscle tissue (106 upregulated and 199 downregulated)(Figure 3D).

3.4. GO Enrichment and KEGG Pathway of DEGs in Thigh and Longissimus Dorsi Muscle

Firstly, we analyzed the DEGs which were screened using low thresholds (|log₂FC| > 0 and padj < 0.05) between Ma and Ju, and explore the potential functions of DEGs.

In thigh muscle, Gene Ontology (GO) enrichment analysis revealed that these DEGs were significantly associated with 706 biological process terms, 41 cellular component terms, and 62 molecular function terms (Table S5). KEGG pathway analysis indicated DEGs were significantly associated with 64 pathways (p<0.05). These GO enrichment and pathways be related to development of muscles (Table 3).

In longissimus dorsi muscles, Gene Ontology (GO) enrichment analysis revealed that these DEGs were significantly associated with 500 biological process terms, 52 cellular component terms, and 49 molecular function terms (Table S5). KEGG pathway analysis indicated DEGs were significantly associated with 64 pathways (P<0.05). These GO enrichment and pathways be related to development of muscles (Table 4).

In thigh and longissimus dorsi muscle, DEGs were enriched to three common GO term (muscle cell migration, smooth muscle cell migration and myofibril) and two common KEGG pathway (HIF-1 signaling pathway, and Glycolysis/Gluconeogenesis). The result suggested HIF-1 signaling pathway and Glycolysis/Gluconeogenesis is related to muscle development between different rabbit breeds.

3.5. GO Enrichment, KEGG Pathway and PPI of Top-50 DEGs

Secondly, we screened the DEGs using high thresholds (|log₂FC| > 1 and padj < 0.05) between Ma and Ju, and further explore the potential functions of DEGs by GO enrichment, KEGG pathway and PPI analysis.

In thigh muscle, These DEGs were significantly enriched GO term (Figure 4A) related to muscle development (muscle cell migration(GO:0014812), positive regulation of cell proliferation(GO:0008284), smooth muscle cell migration(GO:0014909)) and molecule binding (regulation of cell adhesion(GO:0030155), cell adhesion molecule binding(GO:0050839), fibroblast growth factor binding(GO:0017134), integrin binding(GO:0005178), extracellular matrix binding(GO:0050840), fibronectin binding(GO:0001968)), and significantly enriched KEGG pathway (Figure 4B) related to MicroRNAs in cancer(ocu05206), PI3K-Akt signaling pathway(ocu04151), HIF-1 signaling pathway(ocu04066) and ECM-receptor interaction(ocu04512). PPI result (Figure 4C) showed two group, PRL is one key hub-gene in group 1, and NF-κB, BCL, MYC and GAPDH are important hub-gene in Group 2.

In longissimus dorsi muscle, These DEGs were significantly enriched GO term (Figure 4D) related to muscle development of biological process (muscle organ morphogenesis (GO:0048644), muscle structure development (GO:0061061), muscle tissue morphogenesis (GO:0055008),muscle organ development (GO:0007517), muscle tissue development (GO:0060537)), and fiber (contractile fiber(GO:0043292),myofibril(GO:0030016),contractile fiber part(GO:0044449),stress fiber(GO:0001725)) and filament bundle (contractile actin filament bundle(GO:0097517),actin filament bundle(GO:0032432)) of cellular component.These DEGs were significantly enriched KEGG pathway (Figure 4E) related to dilated cardiomyopathy (ocu05414), hypertrophic cardiomyopathy (ocu05410), motor proteins(ocu04814),HIF-1 signaling pathway (ocu04066), adrenergic signaling in cardiomyocytes (ocu04261), cardiac muscle contraction (ocu04260), glycolysis/gluconeogenesis (ocu00010), arrhythmogenic right ventricular cardiomyopathy (ocu05412), ribosome(ocu03010), proteoglycans in cancer (ocu05205). PPI result (Figure 4F) showed three group, PRL is one key hub-gene in group 1, and GAPDH and MYL/MYH are important hub-gene in Group 2 and 3, respectively. These results indicate that the DEGs among breed exhibits variability between thigh and longissimus dorsi muscle.

3.6. The Analysis of Co-Expressive and Tissue-Specific Expressive DEGs in Thigh and Longissimus Dorsi Muscle

Finally, we analyzed co-expressive and tissue-specific expressive DEGs in thigh and longissimus dorsi muscle. The 124 co-expressive DEGs (Figure 5A) were significantly enriched GO term (Figure 5D) related to Golgi to vacuole transport (GO:0006896), AP-3 adaptor complex (GO:0030123), AP-type membrane coat adaptor complex (GO:0030119), membrane coat (GO:0030117), coated membrane(GO:0048475), NADP binding (GO:0050661), oxidoreductase activity, acting on the aldehyde or oxo group of donors, NAD or NADP as acceptor (GO:0016620), structural constituent of ribosome (GO:0003735), oxidoreductase activity, acting on the aldehyde or oxo group of donors (GO:0016903), NAD binding (GO:0051287), and GAPDH is a important hub-gene in PPI network of these co-expressive DEGs (Figure 5E). The result suggested that the difference in energy metabolism related to glucose might be the reason for the difference in muscle development between Ma and Ju rabbit.

The longissimus dorsi muscle-specific expresive DEGs were significantly enriched GO term (Figure 5B) related to purine and nucleoside t biosynthetic process (ATP biosynthetic process (GO:0006754), purine ribonucleoside triphosphate biosynthetic process (GO:0009206), purine nucleoside triphosphate biosynthetic process (GO:0009145), ribonucleoside triphosphate biosynthetic process (GO:0009201), nucleoside triphosphate biosynthetic process (GO:0009142), purine nucleoside monophosphate biosynthetic process (GO:0009127), purine ribonucleoside monophosphate biosynthetic process (GO:0009168), ribonucleoside monophosphate biosynthetic process (GO:0009156), nucleoside monophosphate biosynthetic process (GO:0009124), purine ribonucleotide biosynthetic process (GO:0009152), purine nucleotide biosynthetic process (GO:0006164), purine-containing compound biosynthetic process (GO:0072522), ribonucleotide biosynthetic process (GO:0009260), ribose phosphate biosynthetic process (GO:0046390)), and proton transport related molecules (energy coupled proton transport, down electrochemical gradient (GO:0015985), ATP synthesis coupled proton transport (GO:0015986), mitochondrial proton-transporting ATP synthase complex, coupling factor (GO:0000276), mitochondrial proton-transporting ATP synthase complex (GO:0005753), proton-transporting ATP synthase complex, coupling factor (GO:0045263), proton-transporting ATP synthase complex (GO:0045259), proton-transporting two-sector ATPase complex, proton-transporting domain (GO:0033177), proton-transporting two-sector ATPase complex (GO:0016469), structural constituent of ribosome (GO:0003735)). PPI of these DEGs cluster into three independent network (Figure 5C), and SOX, MYH/MYL and PRL are hub-genes of networks.

The thigh muscle-specific expressive DEGs were significantly enriched GO term (Figure 5F) related to structural constituent of ribosome (GO:0003735), chemokine activity (GO:0008009), chemokine receptor binding (GO:0042379), cytokine activity (GO:0005125), G protein-coupled receptor binding (GO:0001664), Hsp70 protein binding (GO:0030544), oxidoreductase activity, acting on the aldehyde or oxo group of donors, NAD or NADP as acceptor (GO:0016620), oxidoreductase activity, acting on the aldehyde or oxo group of donors (GO:0016903), NAD binding (GO:0051287), cytokine receptor binding (GO:0005126), and unfolded protein binding (GO:0051082). These DEGs form a regulatory network consisting of NF-κB, BCL, CXCL10 and PRL hub genes.

4. Discussion

Skeletal muscle development is a dynamic developmental process. It gradually differentiates from mesoderm-derived cells into myoblasts, merges into myotubes, and matures to form muscle fibers [26]. Skeletal muscle morphology has been extensively measured in pigs [27,28], chickens [29] and sheep [30], but little research has been done on rabbits. In one study, hind leg muscles of fetus, child and adult rabbits were collected for paraffin section, founding extremely significant differences in the diameter, number, area and density of skeletal muscle fibers of the fetus, child, adult rabbit hind legs [31]. In addition, it was found that the leg muscles of fetuses’ rabbits were not fully developed at 2 weeks [32]. There were only a few primary muscle cells. When in childhood and adulthood, secondary muscle fibers have been completely formed, and the diameter of muscle fibers has increased significantly. Exploration of skeletal muscle structure reveals changes in leg muscle fibers during fetal, childhood, and adult stages, providing a theoretical basis for better understanding of rabbit skeletal muscle development.

With the development of high-throughput sequencing, RNA sequencing technology has been applied to systematically study animal growth and development regulatory genes, and massive amounts of transcriptional data have been produced, especially for mammals. RNA-Seq was used to study the transcriptome of the breast muscle of Jinghai yellow chickens at different growth phase and found that RAC2 might regulate cell proliferation by regulating the PAKs/MAPK8 pathway [33]. Additionally, RNA-Seq was used for tissues transcriptome analysis such as skeletal muscles of crossbred beef steers, found that variation among DEG identified by cohort suggests that environment and breed play large roles in the expression of genes associated with feed efficiency in the muscle of beef cattle [34].

To uncover genes that regulate muscle tissue development, in this study, we used RNA-Seq to analyze and identify differentially expressed genes in the muscle tissue transcriptome between Sichuan linen rabbits and Checkered Giant rabbits. Hence, we collected the thigh muscle and longissimus dorsi muscle, from Sichuan linen rabbit and Checkered Giant rabbit, of 28, 56, and 84 days for RNA-seq to explore mRNAs that might have potential regulatory effects on skeletal muscle development. The sample correlation and HCA verified the reliability of the samples. We identified 1983 differentially expressed genes, including 785 upregulated and 1198 downregulated genes, in the thigh muscle of Sichuan linen rabbits compared with that of Checkered Giant rabbits. GO enrichment analysis of these genes showed significant enrichment in related to muscle cell migration, negative regulation of cell death, smooth muscle cell migration, onsarcomere, Myofibril, actin filament bundle, organic acid binding, calmodulin binding, and NADP binding. KEGG pathway annotation revealed that these genes were enriched in Glycolysis/Gluconeogenesis, HIF-1 signaling pathway, PI3K-Akt signaling pathway, FoxO signaling pathway, and MAPK signaling pathway in the thigh muscle of the two breeds of rabbit.

We found that these up-regulated genes GAPDH, PGM1, IGF1, and SGK1 were significantly higher and the down-regulated gene ACSS1 was lower in Ma Leg compared with Ju Leg. These results may suggest that they play critical roles in muscle development of rabbit. Studies have shown that these differentially expressed genes play important roles in different biological processes such as cell proliferation, apoptosis, and migration. GAPDH were expressed in breast muscles, with the prolongation of ages in breast muscle, the number of DETs decreased progressively [35]. PGM1 in regulating the growth and development of myoblasts and specific IMP deposition in Jingyuan chickens [36]. IGF1 have been demonstrated to play essential roles in protein synthesis and fish muscle growth [37]. That serum- and glucocorticoid-inducible kinase 1 (SGK1) regulates muscle mass maintenance via downregulation of proteolysis and autophagy as well as increased protein synthesis during hibernation [38]. Acss1 may constitute a regulatory effect in mitochondrial dynamics and muscle wasting, and was negatively related to the weight of gastrocnemius [39]. So we speculate these genes involved in rabbit muscle development.

We identified 1955 differentially expressed genes, including 734 upregulated and 1221 downregulated genes, in the longissimus dorsi muscle of Sichuan linen rabbits compared with that of Checkered Giant rabbits. GO enrichment analysis showed that these genes were significantly enriched in muscle cell migration, inflammatory response, actin cytoskeleton, contractile fiber, myofibril, sarcomere calcium ion binding, fibronectin binding, and actin filament binding. KEGG pathway analysis results revealed that HIF-1 signaling pathway, Glycolysis/Gluconeogenesis, regulation of actin cytoskeleton, and PPAR signaling pathway were significantly enriched. We found that these up-regulated genes GAPDH, ALDOA, ACTN1, and MYH10 were significantly higher and the down-regulated gene ACSS1 was lower in Ma Long compared with Ju Long. ALDOA genes related to glycogen degradation, glucose transport, and glycolysis [40]. That proper regulation of ACTN1 levels is essential for normal muscle and neuromuscular junction development [41]. Bisulfite sequencing was utilized to analyze and compare the muscle DNA methylation profiles of Larimichthys crocea inhabiting different environments, finding muscle-related genes such as myh10, myf5 and myf6 etc [42].

Through analyzing the GO enrichment and KEGG pathways of related up-regulated and down-regulated Top-50 DEGs.In this part we also found some GO enrichment and KEGG pathway involved in rabbit muscle development,like muscle cell cellular homeostasis,muscle structure development,sarcolemma,glycoprotein complex,fibroblast growth factor bindingand growth factor binding, Circadian rhythm and PPAR signaling pathway,etc.These results surmise that genes in these pathways may play a key role in muscle development.

At present, although the function of these DEGs in rabbit muscle development has not been reported, according to our analysis results in this study, we hypothesize that these DEGs may play an important role in determining the differences in muscle development between Ma and Ju rabbits. The function of these candidate genes should be investigated in future studies.

5. Conclusions

In this study, RNA-Seq was used to analyze the transcriptome of the thigh muscle and the longissimus dorsi muscle between Sichuan linen rabbits and Checkered Giant rabbits, and our data showed that the differentially expressed genes in muscle tissues of the two rabbit breeds were mainly enriched in some important pathways such as Glycolysis/Gluconeogenesis and HIF-1 signaling pathway. These DEGs may play an important role in muscle development. These results provide an important theoretical basis for future clarifying the regulation mechanism of muscle development and molecular breeding of Ma rabbit.