Front. Genet. Frontiers in Genetics Front. Genet. 1664-8021 Frontiers Media S.A. 777232 10.3389/fgene.2021.777232 Genetics Original Research Circular RNA, microRNA and Protein Profiles of the Longissimus Dorsi of Germany ZIKA and Sichuan White Rabbits Zhang et al. Omics Data of Rabbit Muscle Zhang Xiangyu 1 2 * Zhang Cuixia 1 2 Yang Chao 1 2 Kuang Liangde 1 2 Zheng Jie 1 2 Tang Li 1 2 Lei Min 1 2 Li Congyan 1 2 Ren Yongjun 1 2 Guo Zhiqiang 1 2 Ji Yang 1 2 Deng Xiaodong 1 Huang Dengping 1 2 Wang Gaofu 3 Xie Xiaohong 1 2 * 1 Sichuan Animal Sciences Academy, Chengdu, China 2 Animal Breeding and Genetics Key Laboratory of Sichuan Province, Chengdu, China 3 Chongqing Academy of Animal Sciences, Chongqing, China

Edited by: Angela Cánovas, University of Guelph, Canada

Reviewed by: Lifan Zhang, Nanjing Agricultural University, China

Emily Louise Clark, University of Edinburgh, United Kingdom

 *Correspondence: Xiangyu Zhang, zhangxiangyu216@sina.com; Xiaohong Xie, 928742851@qq.com

This article was submitted to Livestock Genomics, a section of the journal Frontiers in Genetics

 24 12 2021 2021 12 777232 15 09 2021 29 11 2021 Copyright © 2021 Zhang, Zhang, Yang, Kuang, Zheng, Tang, Lei, Li, Ren, Guo, Ji, Deng, Huang, Wang and Xie. 2021 Zhang, Zhang, Yang, Kuang, Zheng, Tang, Lei, Li, Ren, Guo, Ji, Deng, Huang, Wang and Xie

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Due to the dietetic properties and remarkable nutritive value of rabbit meat, its industry is increasing rapidly. However, the association between circular RNAs, microRNAs, and proteins and muscle fiber type, and meat quality of rabbit is still unknown. Here, using deep sequencing and iTRAQ proteomics technologies we first identified 3159 circRNAs, 356 miRNAs, and 755 proteins in the longissimus dorsi tissues from Sichuan white (SCWrabs) and Germany great line ZIKA rabbits (ZIKArabs). Next, we identified 267 circRNAs, 3 miRNAs, and 29 proteins differentially expressed in the muscle tissues of SCWrabs and ZIKArabs. Interaction network analysis revealed some key regulation relationships between noncoding RNAs and proteins that might be associated with the muscle fiber type and meat quality of rabbit. Further, miRNA isoforms and gene variants identified in SCWrabs and ZIKArabs revealed some pathways and biological processes related to the muscle development. This is the first study of noncoding RNA and protein profiles for the two rabbit breeds. It provides a valuable resource for future studies in rabbits and will improve our understanding of the molecular regulation mechanisms in the muscle development of livestock. More importantly, the output of our study will benefit the researchers and producers in the rabbit breeding program.

 circular RNA miRNA muscle longissimus dorsi rabbit 1 Introduction

As a functional food, rabbit meat can provide properties and be remarkably nutritive (Dalle Zotte and Szendro, 2011). The rabbit industry for meat production is developing rapidly (van Harten and Cardoso, 2010; Paci et al., 2013) and has reached 200 million tons of meat every year, because rabbit meat has characteristics of rich protein, low cholesterol, and low fat. Thus, improving the yield and quality of rabbit meat is important for the breeding program. Sichuan white rabbits (SCWrabs) are a unique Chinese domestic breed which has excellent characteristics, including strong adaptability, fecundity, and resistance to disease, easy to feed, and tender meat (Li S. et al., 2018). ZIKA rabbits (ZIKArabs) are a rabbit breed developed in Germany as a high yielding hybrid for the meat industry and can attend a weight of 3.2 kg in 84 days (Kuang et al., 2014). With our observation, SCWrabs and ZIKArabs have significant differences in the growth rate, adult weight, and meat quality.

Muscle fiber type is important for the muscle growth and meat production of rabbits and the muscle fiber property can be affected by the growth rate (Ruusunen and Puolanne, 2004). It has been reported that many factors, such as age, gender, genetic heritage, muscle function, and production system, can affect the muscle fiber type (Korfage et al., 2009; Park et al., 2009), which can be classified by the isoforms of myosin proteins that are well known for their roles in muscle contraction and intracellular motility (McKoy et al., 1998; Sweeney and Hammers, 2018). In addition, some studies have been demonstrated to identify genes related to the muscle fiber type and meat quality in livestock animals. For example, the skeletal muscle transcriptome study in cattle showed the importance of energy and protein metabolism in determining meat quality after the aging process (Kee et al., 2008); development of longissimus muscle related genes were reported to be different in two different cattle breeds, such as FSTL1, IGFBP5, and FABP5 (Lehnert et al., 2007); MSTN was reported as a regulator of the skeletal muscle development in terms of the lean mass, fat mass, and glucose metabolism (Guo et al., 2009). In rabbit, some genes like FATSO (Xing et al., 2013), MYF5 (Wang et al., 2017a), POU1F1 (Wang et al., 2015), CAST (Wang et al., 2016), and MYPN (Wang et al., 2017b) have been found to be associated with the rabbit meat quality; and Germany great line of ZIKA rabbits (ZIKArabs) were found to be subjected to less growth inhibition from MSTN at slaughter age, which occurred most possibly in skeletal muscles (Kuang et al., 2014). SCWrabs have slower growth rate, lower adult weight, and tender meat, compared to the imported ZIKArabs. However, our knowledge about the molecular basis of muscle fiber type and meat quality of rabbits is still poor, especially for SCWrabs.

Circular RNA (circRNA) and microRNA (miRNA) are two types of regulatory noncoding RNAs. CircRNAs are generated by splicing events occurring during the maturation of the corresponding pre-mRNAs and can be characterized by a closed ring structure without 3 and 5′ ends (Li X. et al., 2018). They have been hypothesized with multiple functions, including regulating the transcription and splicing of their parental genes, acting as miRNA sponges, regulating protein functions through the direct interaction circRNA/protein, and being translated into proteins with a cap-independent mechanism (Greco et al., 2018). It has been shown that circRNAs are abundant and dynamically expressed during muscle development (Greco et al., 2018). In rabbit, a study has identified 9418 circRNAs and seven of them were related to the atherosclerosis (Zhang et al., 2018). miRNAs are small (∼22 nt) noncoding RNAs that can degrade or repress the expression of their target genes (John et al., 2004). Some studies have been demonstrated to identify rabbit miRNAs (Holcomb et al., 2015), however, the rabbit miRNA sequences recorded in miRBase, a database of miRNAs, are much less than other model species like mouse and chicken. More importantly, very little is known about circRNAs and miRNAs in the longissimus dorsi of rabbits.

Transcriptome sequencing has been widely used to study the expression profiles of protein coding and noncoding genes in the muscle development. Some long noncoding RNAs, another type of noncoding RNA, were identified to be differentially expressed in SCWrabs and ZIKArabs and might play important roles in skeletal muscle development (Kuang et al., 2018). Noce and colleagues analyzed the gene expression profile of the ovine skeletal muscle and characterized the genetic variation of transcripts expressed in five Spanish meat sheep breeds (Noce et al., 2018). Sun compared the longissimus muscle transcriptomes of Qianhua Mutton Merino and Small Tail Han sheep (Sun et al., 2016). While iTRAQ (isobaric Tags for Relative and Absolute Quantitation) is a powerful and popular proteomics technology used for protein identification and quantification due to its sensitivity, accuracy, and high throughput (Yang et al., 2016). It has been used to study the protein profiles and functions in multiple tissues from rabbits (Qiu et al., 2011; Zhou et al., 2013; Liu et al., 2015). In the present study, we aimed to use transcriptome sequencing and iTRAQ proteomics to study the difference of circRNAs, miRNAs, and proteins in the longissimus muscle tissues from SCWrabs and ZIKArabs. The output of this study is a valuable resource for future studies and will benefit researchers and producers of rabbits in the breeding program.

 2 Materials and Methods 2.1 Animals

The experiments of this study were conducted based on the guidelines established by the National Institute of Health. The protocols and proposals of this study were approved by The Medical Ethics Committee on Animal Research of the Sichuan Academy of Animal Science. ZIKArabs and SCWrabs were maintained in the experimental area of Sichuan Academy of Animal Science, Sichuan, China. Three male ZIKArabs and three male SCWrabs from the same litter were selected for this project. After they were maintained under the same experimental conditions for 6 months (both reached individual and sexual maturity), the longissimus dorsi tissues were collected and stored in liquid nitrogen for RNA extraction.

 2.2 Total RNA Extraction

The longissimus dorsi tissue (1 g) of each animal was homogenized in liquid nitrogen, and 10 mg of each homogenized sample was used for total RNA extraction using TRIzol reagent, as described (Chen et al., 2016). In brief, the tissue sample was mixed with 1 ml TRIzol reagent, homogenized using a power homogenizer, and centrifuged (12,000 ×g) for 10 min at 4°C. Then, the fatty layer was removed, and the supernatant was transferred into a new tube. After 0.2 ml chloroform (0.2 ml) was added to the tube and an incubation under room temperature for 3 min, the tube was centrifuged (12,000 × g) for 15 min at 4°C and the aqueous phase was transferred to another RNA tube. For RNA precipitation, we added RNase-free glycogen (10 μg) and 100% isopropanol (0.5 ml) to the tube, and incubated the sample at room temperature for 10 min. After a centrifugation (12,000 ×g, 10 min, 4°C), the RNA pellet was washed using 75% ethanol (1 ml) and air-dried. Then, it was suspended in RNase-free water and water-bathed at 60°C for 10 min. The quality and quantity of the total RNA were determined using the Agilent 2100 Bioanalyzer.

 2.3 Strand Transcriptome Sequencing With rRNA Removal and Small RNA Sequencing

Total RNA (8 μg) of each sample was processed to remove the ribosomal RNA using the Epicentre Ribo-zero™ rRNA Removal Kit (Epicentre, USA). Then, the rRNA-depleted RNA sample was used for the library construction using the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina (NEB, USA), as described (Zhang et al., 2018). Small RNA sequencing libraries of the two rabbit breeds were constructed using the Illumina TruSeq Small RNA Sample Preparation Kit v2, as described (Chen et al., 2019b). Both stranded transcriptome sequencing and small RNA sequencing were performed on the HiSeq 2500 platform in NovoGene (Beijing, China).

 2.4 Bioinformatics Analysis of circRNAs and miRNAs

We initially cleaned the raw data of transcriptome sequencing and small RNA sequencing using SOAPnuke (Chen et al., 2018). Then, we used “find_circ” to predict the potential circRNAs in each sample according to the protocol and filtered the circRNAs less than 5 reads in at least two replicates (Memczak et al., 2013). To obtain the gene expression profiles of muscle tissues of SCWrabs and ZIKArabs, the clean reads were aligned to the rabbit genome (oryCun2) using HISAT2 and the gene expression was profiled using StringTie, according to the protocols (Pertea et al., 2016). The raw read counts aligned to all genes were measured using htseq-count, as described (Chen et al., 2019a). For miRNA identification, we first aligned the clean reads to the rabbit genome (oryCun2) using SOAP2 software (Li et al., 2009) and predicted the rabbit miRNAs using mireap (https://github.com/liqb/mireap) using the genome mapping results. Next, we aligned the predicted mature miRNA sequences to the miRBase (v22) to identify miRNAs homolog to other species. Both known and novel miRNA sequences were used as the reference for miRNA screen in each sample. miRNAs identified in at least two replicates were kept for downstream analyses. edgeR (Robinson et al., 2010) was used to identify differentially expressed circRNAs, genes, and miRNAs in the muscle tissues of SCWrabs and ZIKArabs with the following cut-offs: log2 fold change (Log2FC) > 1 or < −1, p-value < 0.05 and FDR < 0.05.

 2.5 Protein Preparation and iTRAQ Analysis

Proteins of each sample were prepared from 1 mg muscle tissue, and these were used for the iTRAQ proteomics analysis. The procedures of this experiment were conducted the same as a previous study (Yang et al., 2016). Samples were labelled with the iTRAQ tags as follows: SCWrabs (113, 114, 115) and ZIKArabs (1s16, 117, 118). After the data acquisition was performed on a Q EXACTIVE (Thermo Fisher Scientific, CA, United States) coupled online to the HPLC, we performed the database search and protein quantification using Proteome Discoverer (Thermo Fisher Scientific, CA, United States). For differential protein expression, we first calculated the average abundance of a protein in each group and calculated the ratio using the average protein abundance. Then, t-test was used to calculate the p-value. Differentially expressed proteins in the two rabbit groups should satisfy ratio > 1.2 and p-value < 0.05.

 2.6 miRNA, circRNA, mRNA, and Protein Interaction Network

The interactions of miRNA-mRNA and miRNA-circRNA in this study were predicted using the miranda (v3.3a) (John et al., 2004). Then, we predicted the enriched motifs in the mature miRNA sequences using the DREME software (Bailey, 2011) and searched the motifs against the ATtRACT (a database of RNA binding proteins and associated motifs) (Giudice et al., 2016). For circRNAs, we identified their RNA binding proteins (RBPs) using the human data in starBaseV3 (Li et al., 2014). After the miRNA-mRNA, miRNA-circRNA, and circRNA_RBPs were prepared, we visualized their interaction in Cytoscape (v3.7.1).

 2.7 Mutation Analysis

The mutations in the muscle tissues of SCWrabs and ZIKArabs were characterized using the Strelka pipeline (v2.9.10) (Kim et al., 2018). Then, we removed the variants that passed the filter of Strelka. Next, variants identified in all three replicates were annotated using the ANNOVAR pipeline (Wang et al., 2010) with the Ensembl gene annotation.

 3 Results 3.1 CircRNAs in the Rabbit Muscle

Stranded transcriptome sequencing produced 1095.66 million reads (182.61 million reads on average) for all the samples and we obtained 167.37 ∼ 194.56 million reads (12.55 ∼ 14.59 G clean bases) for each sample after data cleaning. Then, the clean data was aligned to the rabbit genome (oryCun2) and it showed that 89.14 ∼ 91.36% of the clean reads mapped. Next, we systematically identified 26,529 circRNAs with 1 or more reads supported in the muscle tissues of SCWrabs and ZIKArabs. After lowly expressed circRNAs (>5 read count in at least two replicates) were filtered, we obtained 3159 circRNAs, of which 2331 and 2733 distributed in SCWrabs and ZIKArabs, respectively (Supplementary Table S1). We showed the numbers of circRNAs identified in all chromosomes of the rabbit genome (Figure 1A) and found that chromosomes chr13, chr1, and chr2 were the top three which can produce circRNAs. Notably, we found that 382 (12.09%), 411 (13.01%), and 854 (27.03%) circRNAs were derived from the intergenic, exon, and intron regions, respectively, and that 1512 (47.86%) circRNAs were formed by both exon and intron parts. Figure 1B showed that 64.95% of the rabbit circRNAs were longer than 1400 bp. Then, we found that 1905 circRNAs were expressed in the muscle tissues of both rabbit breeds and that 426 and 828 circRNAs were specifically expressed in SCWrabs and ZIKArabs, respectively (Figure 1C). Further, we found 31 highly expressed (TPM > 1000) circRNAs in the muscle tissues of two rabbit breeds and that 26 were shared. Using edgeR we identified 120 up-regulated and 147 down-regulated circRNAs in the muscle tissues of SCWrabs compared to ZIKArabs (Figure 1D; Supplementary Table S2). Notably, myosin heavy chain genes are a class of host genes that can derive dysregulated circRNAs, such as MYH1, MYH7, and MYH13, and up-regulated circRNAs in SCWrabs were derived from PI4KA and NEB genes. Functional analysis (Figure 1E) of the circRNA host genes showed that many of them were related to the muscle development related biological processes and pathways, such as “ATP binding,” “Tight junction,” “Regulation of actin cytoskeleton” and “myofibril.”

circRNAs identified in the muscle tissue of SCWrabs and ZIKArabs. (A) Distribution of circRNAs identified in the rabbit chromosomes. (B) Length distribution of circRNAs identified in this study. (C) Venn diagram of circRNAs identified in the muscle tissues of SCWrabs and ZIKArabs. (D) Volcano plot showing differentially expressed circRNAs in the muscle tissues of SCWrabs and ZIKArabs. (E) Functional annotation of the host genes for differentially expressed circRNAs.

 3.2 Gene Expression Profile

Using the sequencing data, we profiled the protein coding gene expression in the muscle tissues of SCWrabs and ZIKArabs. Initially, we identified 7846 and 7193 genes in SC and ZIKA rabbit muscle tissues, respectively. SCWrabs and ZIKArabs shared 7004 genes while 842 and 189 genes were only detected in SCWrabs and ZIKArabs, respectively (Figure 2A). Compared to ZIKArabs edgeR identified 153 differentially expressed protein coding genes (DEGs) in the muscle tissues of SCWrabs (Figure 2B; Supplementary Table S3), including 55 up-regulated and 98 down-regulated genes. It is notable that 4 myosin genes (e.g., MYH4, MYH13, MYLK4, MYL6B) were up-regulated while 9 heat shock protein and 14 histone genes were down-regulated between SCWrabs and ZIKArabs (Table 1). Functional analysis showed 9, 3, and 3 DEGs involved in “negative regulation of cell proliferation,” “negative regulation of growth,” and “fat cell differentiation,” respectively.

Gene expression profiles and miRNAs identified in the muscle tissues of SCWrabs and ZIKArabs. (A) Venn diagram of genes identified in the muscle tissues of SCWrabs and ZIKArabs. (B) Volcano plot of differentially expressed genes in the muscle tissues of SCWrabs and ZIKArabs. (C) Distribution of miRNAs identified in the chromosomes of rabbit. (D) Base bias of the mature miRNAs identified in the muscle tissues of SCWrabs and ZIKArabs. (E) Length distribution of miRNAs identified in this study. (F) Comparison of target genes for the differentially expressed miRNAs.

Significant gene groups differentially expressed in the muscle tissues of SCWrabs and ZIKArabs.

GeneID Log2FC p-value FDR Regulation Gene Name Description
ENSOCUG00000005076 1.48 1.07E-04 4.93E-03 UP MYLK4 Myosin light chain kinase family member 4
ENSOCUG00000017368 1.29 7.30E-04 2.29E-02 UP MYL6B Myosin light chain 6B
ENSOCUG00000010207 1.59 3.38E-05 1.88E-03 UP MYH13 Myosin heavy chain 13
ENSOCUG00000029652 2.13 4.62E-08 5.15E-06 UP MYH4 Myosin heavy chain 4
ENSOCUG00000023743 −2.94 2.85E-13 8.37E-11 DOWN Heat shock 70 kda protein 1B
ENSOCUG00000025112 −3.24 1.94E-15 7.53E-13 DOWN Heat shock 70 kda protein 1B
ENSOCUG00000001280 −1.82 2.32E-06 1.76E-04 DOWN HSP90AA1 Heat shock protein 90 alpha family class A member 1
ENSOCUG00000015533 −1.25 1.01E-03 3.00E-02 DOWN HSPA8 Heat shock protein family A (Hsp70) member 8
ENSOCUG00000004833 −1.98 3.45E-07 3.23E-05 DOWN DNAJA1 DNAJ heat shock protein family (Hsp40) member A1
ENSOCUG00000014773 −2.06 1.19E-07 1.24E-05 DOWN HSPH1 Heat shock protein family H (Hsp110) member 1
ENSOCUG00000012690 −1.68 1.26E-05 7.93E-04 DOWN Heat shock protein family B (small) member 1
ENSOCUG00000027940 −1.91 1.01E-06 8.53E-05 DOWN HSPA4L Heat shock protein family A (Hsp70) member 4 like
ENSOCUG00000021817 −4.70 4.68E-27 8.05E-24 DOWN Heat shock 70 kda protein 1B-like
ENSOCUG00000014459 −1.45 2.72E-04 1.06E-02 DOWN Histone H3.1
ENSOCUG00000024086 −2.21 3.07E-07 2.90E-05 DOWN Histone H2A type 1
ENSOCUG00000023129 −1.21 1.79E-03 4.67E-02 DOWN Histone H2A type 1-E
ENSOCUG00000026810 −2.03 2.45E-05 1.41E-03 DOWN Histone H2B type 1
ENSOCUG00000027315 −1.54 1.09E-04 5.00E-03 DOWN Histone H3
ENSOCUG00000027038 −2.70 5.66E-09 7.58E-07 DOWN Histone H2A type 1
ENSOCUG00000025069 −1.70 2.01E-05 1.19E-03 DOWN Histone H2A type 1-E
ENSOCUG00000026226 −1.42 9.86E-04 2.94E-02 DOWN Histone H2A type 1-E
ENSOCUG00000026771 −1.48 3.24E-04 1.21E-02 DOWN HIST1H2AM Histone H2A type 1
ENSOCUG00000025038 −1.96 3.53E-05 1.94E-03 DOWN Histone H3.1
ENSOCUG00000014458 −2.01 1.29E-06 1.06E-04 DOWN Histone H2A type 1-E
ENSOCUG00000025297 −5.09 2.55E-15 9.32E-13 DOWN Histone H2A type 1
ENSOCUG00000027227 −1.46 2.73E-04 1.06E-02 DOWN HIST1H2BB Histone H2B type 1-B
ENSOCUG00000021035 −2.85 3.03E-11 6.41E-09 DOWN Histone H2B type 1
3.3 miRNA Expression Profiles of Rabbit Muscle Tissues

Using small RNA sequencing, we identified 291 (185 known and 106 novel) miRNA precursor sequences in all the chromosomes of rabbit genome and these miRNA precursors can produce 356 mature miRNAs (Figure 2C; Supplementary Table S3). Among them, 68 were found with both mature and passenger miRNAs (Supplementary Table S3). Bias analysis of nucleotide base showed that the first base preference of rabbit miRNAs was uracil (Figure 2D), and the rabbit miRNAs were enriched in 22 nt in length (Figure 2E). As expected, no specific miRNAs were found in SCWrabs or ZIKArabs. The top 10 highly expressed miRNAs remained the same in the muscle tissues of SCWrabs and ZIKArabs, including ocu-miR-1-3p, ocu-miR-206-3p, ocu-miR-133a-3p, ocu-miR-26b-5p, ocu-miR-101-3p, ocu-miR-378-3p, ocu-miR-27-3p, ocu-let-7d-5p, ocu-let-7i-5p, and ocu-miR-378d-5p. We next identified 3 miRNAs down-regulated in muscle tissues of SCWrabs compared to ZIKArabs (Table 2). Target prediction identified 2124, 2068, and 4694 candidate genes for ocu-miR-182-5p, ocu-miR-146b, and ocu-N-miR-10-5p, respectively. We next analyzed the gene regulation of their targets using the RNA-seq data. It was revealed that 7, 5, and 11 up-regulated genes might be regulated by ocu-miR-182-5p, ocu-miR-146b, and ocu-N-miR-10-5p, respectively (Table 2). It is interesting that some target genes were shared by these three miRNAs (Figure 2F). For example, SH3RF2 was regulated by all three miRNAs; ANO2 and ABCA1 were targeted by ocu-miR-182-5p and ocu-miR-146b-5p; DUSP1 and MFSD8 were shared by ocu-miR-146b-5p and ocu-N-miR-10-5p; and MYH13 was the target gene of ocu-miR-182-5p and ocu-N-miR-10-5p. The regulation relationship between these three miRNAs and their target genes requires further experiments to be explored.

Differentially expressed miRNAs in the muscle tissues of SCWrabs and ZIKArabs.

miRNA edgeR values Target number Gene Log2FC FDR Gene_symbol Description
ocu-miR-182-5p Log2FC = -1.40, p = 0.0004 2124 ENSOCUG00000010207 1.59 1.88E-03 MYH13 Myosin heavy chain 13
ENSOCUG00000006775 2.28 1.12E-06 ANO2 Anoctamin 2
ENSOCUG00000010659 1.51 5.00E-03 ADAM23 ADAM metallopeptidase domain 23
ENSOCUG00000007700 1.76 5.48E-04 GDAP1 Ganglioside induced differentiation associated protein 1
ENSOCUG00000003254 1.33 1.74E-02 CDKL2 Cyclin dependent kinase like 2
ENSOCUG00000009279 1.29 2.15E-02 ABCA1 ATP binding cassette subfamily A member 1
ENSOCUG00000015667 1.98 4.83E-05 SH3RF2 SH3 domain containing ring finger 2
ocu-miR-146b-5p Log2FC = -1.23, p = 0.0016 2068 ENSOCUG00000006775 2.28 1.12E-06 ANO2 Anoctamin 2
ENSOCUG00000007999 1.65 1.03E-03 DUSP1 Dual specificity phosphatase 1
ENSOCUG00000009279 1.29 2.15E-02 ABCA1 ATP binding cassette subfamily A member 1
ENSOCUG00000015003 1.22 4.34E-02 MFSD8 Major facilitator superfamily domain containing 8
ENSOCUG00000015667 1.98 4.83E-05 SH3RF2 SH3 domain containing ring finger 2
ocu-N-miR-10-5p Log2FC = -1.01, p = 0.0096 4694 ENSOCUG00000005076 1.48 4.93E-03 MYLK4 Myosin light chain kinase family member 4
ENSOCUG00000017507 1.49 5.69E-03 PAX3 Paired box 3
ENSOCUG00000027746 3.68 3.34E-15
ENSOCUG00000011746 1.43 7.59E-03 NR4A1 Nuclear receptor subfamily 4 group A member 1
ENSOCUG00000007999 1.65 1.03E-03 DUSP1 Dual specificity phosphatase 1
ENSOCUG00000015003 1.22 4.34E-02 MFSD8 Major facilitator superfamily domain containing 8
ENSOCUG00000015667 1.98 4.83E-05 SH3RF2 SH3 domain containing ring finger 2
ENSOCUG00000023102 2.17 3.36E-06 LOC100349173 Uncharacterized protein c15orf52 homolog
ENSOCUG00000010207 1.59 1.88E-03 MYH13 Myosin heavy chain 13
ENSOCUG00000029652 2.13 5.15E-06 MYH4 Myosin heavy chain 4
ENSOCUG00000008198 1.29 2.21E-02 IFFO1 Intermediate filament family orphan 1
3.4 Protein Expression Profile and Differential Expression

Then, we used iTRAQ to study the protein profiles of the muscle tissues of SCWrabs and ZIKArabs. Initially, we obtained 269,250 spectra from the liquid chromatography coupled to mass spectrometry (LC-MS/MS) analysis, of which 13,232 (5221 unique peptides) were aligned to 755 rabbit proteins. Length distribution of peptides showed that they were enriched between 9 and 14 amino acids (Figure 3A). The analysis of protein molecular weight (Figure 3B) revealed that 116, 114, 111, and 100 proteins were in 30 ∼ 40 kDa, 20 ∼ 30 kDa, 10 ∼ 20 kDa, and > 100 kDa, respectively. Like miRNAs, no specific proteins were found in the muscle tissues of SCWrabs and ZIKArabs, probably because many rabbit proteins are unknown currently. Next, we showed the top 20 proteins with the most unique peptides identified in the muscle tissues of SCWrabs and ZIKArabs (Figure 3C). TTN was identified with the most unique peptides (865 peptides), followed by NEB (237 peptides) and MYH4 (135 peptides). Differential expression analysis identified 19 up-regulated and 10 down-regulated proteins in the muscle tissues of SCWrabs compared to ZIKArabs (Figure 3D; Table 3). Among them, 3 myosin proteins (MYH7B, MYH8, and MYH13) and 2 proteasome proteins (PSMB1 and PSMC5) were up-regulated in SCWrabs compared to ZIKArabs. Then, we performed KEGG pathway enrichment analysis and found that 5, 3, 2, and 2 proteins were involved by the pathways of “tight junction,” “hippo signaling pathway,” “proteasome,” and “ribosome,” respectively (Figure 3E).

iTRAQ proteomics identified differentially expressed proteins in the longissimus muscle tissues of SCWrabs and ZIKArabs. (A) Length distribution of peptides detected in the longissimus muscle tissues by iTRAQ. (B) Distribution of protein mass. (C) Top highly expressed proteins in the longissimus muscle tissue of rabbits. (D) Volcano plot of the differentially expressed proteins in the muscle tissues of SCWrabs and ZIKArabs. (E) KEGG pathway enrichment analysis for the differentially expressed proteins.

Differentially expressed proteins in the longissimus muscle tissues of SCWrabs and ZIKArabs.

ProteinID Log2FC p-value Regulation Gene name Description
ENSOCUP00000021805.1 −0.393 0.043 down AAMDC Adipogenesis associated Mth938 domain containing
ENSOCUP00000020925.2 −0.667 0.048 down AKR1B10 Aldo-keto reductase family 1 member B10
ENSOCUP00000006694.2 −0.197 0.050 down CTNNA1 Catenin alpha 1
ENSOCUP00000014440.2 −0.383 0.011 down CYB5R1 NADH-cytochrome b5 reductase 1
ENSOCUP00000010958.3 −0.412 0.015 down HSPD1 Heat shock protein family D (Hsp60) member 1
ENSOCUP00000012421.1 −0.180 0.019 down LOC100348835 Histone H1.4
ENSOCUP00000007190.2 −0.302 0.014 down LOC100354435 Heat shock 70 kDa protein 1B
ENSOCUP00000001099.3 −0.157 0.034 down NNT Nicotinamide nucleotide transhydrogenase
ENSOCUP00000005290.3 −0.235 0.031 down USP36 Ubiquitin specific peptidase 36
ENSOCUP00000011737.2 −0.179 0.020 down
ENSOCUP00000002957.2 0.146 0.043 Up COPS4 COP9 signalosome subunit 4
ENSOCUP00000009625.2 0.302 0.001 Up COPS7A COP9 signalosome subunit 7A
ENSOCUP00000021470.1 0.251 0.037 Up DHRS7B Dehydrogenase/reductase 7B
ENSOCUP00000005141.2 0.223 0.012 Up GPX1 Glutathione peroxidase 1
ENSOCUP00000007030.2 0.243 0.045 Up H2AFV H2A histone family, member V
ENSOCUP00000010917.2 0.311 0.005 Up HSPB1 Heat shock protein family B (small) member 1
ENSOCUP00000021175.1 0.382 0.038 Up LOC108175352 60S ribosomal protein L17
ENSOCUP00000013763.2 0.201 0.009 Up LRRC20 Leucine rich repeat containing 20
ENSOCUP00000023874.2 0.812 0.004 Up MYH13 Myosin heavy chain 13
ENSOCUP00000001968.2 0.341 0.021 Up MYH7B Myosin heavy chain 7B
ENSOCUP00000008830.2 0.321 0.007 Up MYH8 Myosin-8
ENSOCUP00000021074.2 0.171 0.004 Up PPP2R2A Protein phosphatase 2 regulatory subunit Balpha
ENSOCUP00000006889.2 0.215 0.001 Up PSMB1 Proteasome subunit beta 1
ENSOCUP00000003741.2 0.152 0.013 Up PSMC5 Proteasome 26S subunit, Atpase 5
ENSOCUP00000017814.1 0.165 0.045 Up SHMT1 Serine hydroxymethyltransferase 1
ENSOCUP00000010672.2 0.466 0.011 Up TMEM38A Transmembrane protein 38A
ENSOCUP00000022564.1 0.267 0.017 Up YWHAG Tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein gamma
ENSOCUP00000021670.1 0.297 0.003 Up
ENSOCUP00000025416.1 0.299 0.033 Up
3.5 Noncoding RNA-Protein Interaction Network

To build the interaction network of rabbit noncoding RNAs, genes, and proteins, we first predicted the miRNA targets in rabbit circRNAs and genes. It showed that 3152 circRNAs and 20262 genes can be regulated by the rabbit mature miRNAs. Among them, 129 differentially expressed circRNAs (64 up-regulated and 65 down-regulated) and 160 genes (57 up-regulated and 103 down-regulated) were the potential targets of differentially expressed miRNAs in the muscle tissues of SCWrabs compared to ZIKArabs. Using the iTRAQ results we found 39 protein coding genes (16 up-regulated and 23 down-regulated) regulated by the differentially expressed miRNAs. Next, using the DREME software we found two motifs—AGGUAG and GGAGUCA in 23 and 19 mature miRNA sequences, respectively. In the ATtRACT database only motif AGGUAG was found to interact with 4 (AGO1, DAZAP1, MSI1 and ZRANB2) RBPs. However, if we modify the motifs to AGGUA and GGAGU, 11 (AGO1, DAZAP1, HNRNPA1, HRB27C, KHSRP, MSI1, PTBP1, RBM5, SFPQ, SRSF1, and ZRANB2) and 5 (FMR1, HNRNPDL, LIN28A, MEX-3, and SRSF2) RBPs could be found for AGGUA and GGAGU, respectively. Unfortunately, none of these RBPs were found in the proteomics data. StarBaseV3 identified 18 RBPs binding with 3079 circRNA host genes by similarity. Among them, 19 circRNAs were differentially expressed between the muscle tissues of SCWrabs and ZIKArasb while 2 RBPs were identified in our dataset.

After simple interactions were filtered, a total of 2 miRNAs, 9 circRNAs, 2 RBPs, and 24 genes were subjected to Cytoscape for visualization (Figure 4). Two miRNAs ocu-N-miR-10-5p and ocu-miR-146b-5p shared 6 targets (3 protein coding genes and 3 circRNAs) while the two RBPs (HNRNPA2B1 and HNRNPK) shared no binding circRNAs. The interactions in noncoding RNAs, protein coding genes, and RBPs require further experiments to be explored in rabbits.

Interaction network of differentially expressed circRNAs, miRNAs, their target genes, and binding proteins.

 3.5 Gene Variants and miRNA Isoforms

Using the small RNA sequencing data, we found that the dominant sequences of 10 miRNAs (6 conserved and 4 novel miRNAs) differed in SCWrabs and ZIKArabs (Table 4). Among them, ocu-N-miR-26-5p was only detected in ZIKA rabbits. Most of these alternations in miRNA sequences could be produced by the cut site preference of DICER enzyme in the cytoplasm because the sequence difference occurred at the end of 3′ arm of miRNA. However, we found a base “A” added to the 5’ arm of ocu-miR-411-5p in SCWrabs, compared to ZIKArabs.

miRNA isoforms identified in the muscle tissues of SCWrabs and ZIKArabs.

miRNA SCWrabs Count a ZIKArabs Count a
ocu-miR-29a-3p CTA​GCA​CCA​TTT​GAA​ATC​AGT​G 96,129,85 CTA​GCA​CCA​TTT​GAA​ATC​AGT 73,85,90
ocu-miR-29d-3p CTA​GCA​CCA​TTT​GAA​ATC​AGT​G 96,129,85 CTA​GCA​CCA​TTT​GAA​ATC​AGT 73,85,90
ocu-miR-340b-3p TCC​GTC​TCA​GTT​ACT​TTA​TAG​CC 1,6,9 TCC​GTC​TCA​GTT​ACT​TTA​TAG​C 6,8,5
ocu-miR-411-5p ATA​GTA​GAC​CGT​ATA​GCG​TAC​G 117,244,252 TAG​TAG​ACC​GTA​TAG​CGT​ACG 119,318,365
ocu-miR-487-3p AAT​CGT​ACA​GGG​TCA​TCC​ACT​T 13,48,48 AAT​CGT​ACA​GGG​TCA​TCC​ACT 26,35,49
ocu-miR-532-3p CCT​CCC​ACA​CCC​AAG​GCT​TGC​A 54,81,78 CCT​CCC​ACA​CCC​AAG​GCT​TGC 42,80,98
ocu-N-miR-26-5p - 0,0,0 TTG​ACC​TAT​GAA​ATG​ACA​GAT​G 6,7,8
ocu-N-miR-48-5p TTG​CTC​TCT​CTC​TAA​CTC​TGT 41,66,43 TTG​CTC​TCT​CTC​TAA​CTC​TGT​C 32,92,59
ocu-N-miR-79-3p TCA​GCC​CCT​GGG​GTA​GTT​CTG​A 14,22,30 TCA​GCC​CCT​GGG​GTA​GTT​CTG 17,24,60
ocu-N-miR-80-3p TCA​GCC​CCT​GGG​GTA​GTT​CTG​A 14,22,30 TCA​GCC​CCT​GGG​GTA​GTT​CTG 17,24,60

Raw read counts of the miRNA sequence in the three biological replicates.

Next, we analyzed the variants in protein coding genes using the transcriptome sequencing data. Initially, 121,588 and 116,844 variants were identified in SCWrabs and ZIKA, respectively. After these variants were annotated by ANNOVAR, we found that 6605 exonic variants were shared while 6241 (from 4375 genes) and 6650 (from 4672 genes) exonic variants were specific to SCWrabs and ZIKArabs, respectively (Figure 5A; Supplementary Table 5). It was shown that the top 20 genes with most variants in SCWrabs and ZIKArabs were different (Figure 5B). Interestingly, HMCN1, COL15A1, NEB, VWA2, and AKAP6 were more frequently mutated in SCWrabs, while ENSOCUG00000027987, SYNE1, VPS13C, and BIVM-ERCC5 were more frequently mutated in ZIKArabs. Then, we annotated the specifically mutated genes for SCWrabs and ZIKArabs using DAVID Bioinformatics Resource. After shared terms were filtered, specific pathways and gene ontology terms were shown in Figure 5C. It is notable that muscle development related pathways and biological processes were involved by these specific genes, such as “zinc ion binding,” “regulation of actin cytoskeleton,” “metal ion binding,” “protein serine/threonine kinase activity,” and “calcium ion binding.”

Gene variants identified in the muscle tissues of SCWrabs and ZIKArabs. (A) Number of total variants and exonic variants identified in the muscle tissues of SCWrabs and ZIKArabs. (B) Number of mutation events in the top 20 genes. (C) Functional annotation for the mutated genes specifically identified in the muscle tissues of SCWrabs and ZIKArabs. Values in the blocks represent the numbers of genes involved in the pathways and processes.

 4 Discussion

In this study, we identified 3159 circRNAs, 356 miRNAs, and 755 proteins in the longissimus dorsi tissues of two rabbit breeds (SCWrabs and ZIKArabs). They will provide a valuable resource for future research in rabbits. The circRNAs were mainly derived from chromosomes chr13, chr1, and chr2, which agreed with a previous study (Zhang et al., 2018). The features of miRNA sequences identified in this study, such as the first nucleotide preference and length distribution, also agreed with the miRNA database miRBase (Grey et al., 2005). This evidence first supported that the noncoding RNAs identified in this study might be credible. Compared to SCWrabs we identified 267 circRNAs (120 up-regulated and 147 down-regulated), 3 miRNAs (down-regulated), and 29 proteins (19 up-regulated and 10 down-regulated) differentially expressed in the muscle tissues of ZIKArabs and they will be further discussed.

Myosin proteins are motor proteins that convert chemical energy in the form of ATP to mechanical energy, thus generating force and movement (Cooper, 2000). In our study, the transcriptome sequencing and iTRAQ proteomics identified 4 myosin genes (e.g., MYH4, MYH13, MYLK4, MYL6B) (Table 1) and 3 myosin proteins (MYH7B, MYH8, and MYH13) (Table 3), respectively, up-regulated in the muscle tissues of SCWrabs compared to ZIKArabs. In addition, myosin genes (MYH1, MYH7, and MYH13) derived circRNAs were dysregulated in the muscle tissues of SCWrabs compared to ZIKArabs. Also, we found that differentially expressed miRNAs ocu-miR-182-5p and ocu-N-miR-10-5p share the target gene MYH13 and that MYH4 is another target gene of ocu-N-miR-10-5p (Figure 2F). These indicate that myosin and myosin related noncoding RNAs might be one of the reasons for the differences in muscle fiber types and meat qualities of SCWrabs and ZIKArabs. Further, MYH13 is found to be preferentially expressed in most orbital and global fibers in the central innervation zone of rabbit extraocular muscles (Briggs and Schachat, 2002) and it has been proved to colocalize with the MYH4 in global fibers of rabbits (Briggs and Schachat, 2002). In synergy with MYH4 and an increased level of sarcoplasmic reticulum, the localized expression of MYH13 can result in a differentially responsive and very fast contracting region of extraocular muscle (Briggs and Schachat, 2002). The presence of MYH13 and MYH4 in the longissimus dorsi might be related to the good meat quality of SCWrabs. Myostatin (MSTN), also known as growth differentiation factor 8 (GDF-8), is a candidate gene related to the muscle growth. It is a secreted growth factor predominately expressed in skeletal muscle that negatively regulates skeletal muscle mass via the inhibition of myogenesis: muscle cell growth and differentiation (Lee, 2004). It has been shown that the inhibition of MSTN signaling in skeletal muscle can increase the lean mass, reduce the fat mass, and improve the glucose metabolism (Guo et al., 2009). Compared to ZIKArabs, we found a significant upregulation of the MSTN gene in the longissimus dorsi of SCWrabs (Supplementary Table S3). This might explain that SCwrabs have less weight, height, fat mass, and better taste than ZIKArabs.

circRNAs and miRNAs are regulatory RNAs related to the expression of multiple target genes and proteins (Qu et al., 2015; Chen et al., 2019b). They have been shown as key regulators in the muscle development and fiber types. For example, in sheep skeletal muscle 886 circRNAs were identified and many of them (e.g., circ776) interacted with muscle-specific miRNAs which were involved in growth and development of muscle (Cao et al., 2018). In Guizhou miniature pig (Sus scrofa), 149 circRNAs were identified to be associated with the muscle growth and the circRNA host genes were found to regulate the muscle development, contraction, chromatin modification, cation homeostasis, and ATP hydrolysis-coupled proton transport (Liang et al., 2017). Further, some circRNA host-genes were found to be related to the myogenesis in the longissimus dorsi of Landrace and Lantang pigs (Sun et al., 2017). In our study, we identified 267 differentially expressed circRNAs in the longissimus dorsi of SCWrabs and ZIKArabs. They might be responsible of the muscle growth and meat quality. In addition, the host genes of some circRNAs were predicted to be involved in the muscle development and fiber type (Figure 1E). However, more experiments are required to validate the association between circRNAs and muscle growth/fiber types in rabbit.

Some miRNAs have been shown to be highly enriched in cardiac and/or skeletal muscle, such as miR-1, miR-133a, miR-133b, miR-206, miR-208, miR-208b, miR-486, and miR-499 (Guller and Russell, 2010). These miRNAs were identified in the current study except miR-486 (Supplementary Table S4). miR-1, miR-133a, and miR-206 were the top three highly expressed miRNAs while miR-208 was lowly expressed (average TPM < 2) in the longissimus dorsi of both SCWrabs and ZIKArabs (Supplementary Table S4). miR-1 can stimulate the myoblast differentiation through the inhibition of HDAC4 (histone deacetylase 4) (Chen et al., 2006); miR-133a increases the myoblast proliferation via the repression of SRF (serum response factor) (Chen et al., 2006); and miR-206 has the influence of the differentiation program via the indirect down-regulation of the helix–loop–helix protein Id, a repressor of MyoD (myogenic differentiation 1) (Kim et al., 2006). We identified three differentially expressed miRNAs in the muscle tissues of SCWrabs compared to ZIKArabs, such as ocu-miR-182-5p, ocu-miR-146b, and ocu-N-miR-10-5p. By targeting NF-B expression, miR-146b promotes myogenic differentiation and modulates multiple gene targets in muscle cells and its homolog miR-146a regulates the maturation and differentiation of vascular smooth muscle cells by targeting NF-κB expression (Dong et al., 2013; Khanna et al., 2014). Interestingly, miR-182 has been studied to inhibit the muscle cell differentiation through the regulation of myogenin (Zhang et al., 2016). It has also been shown to attenuate the atrophy-related gene expression by targeting FoxO3 in skeletal muscle (Hudson et al., 2014). Considering that many miRNAs are unknown for rabbit and the high expression of miRNAs in rabbit longissimus dorsi, it is important to study their associations with the muscle development and fiber type.

In this study, we also identified exonic variants in the muscle tissues of SCWrabs and ZIKArabs (Supplementary Table S5). It is notable that SCWrabs and ZIKArabs had distinct mutation frequency in some genes (Figure 5B). For example, HMCN1, COL15A1, NEB, VWA2, and AKAP6 were more frequently mutated in SCWrabs than ZIKArabs, while ENSOCUG00000027987, SYNE1, VPS13C, and BIVM-ERCC5 were more frequently mutated in ZIKArabs than SCWrabs. Among them, HMCN1 (hemicentin-1), encoding a fibrillar extracellular matrix glycoprotein, is up-regulated in the mice with chronic skeletal muscle fibrosis compared to wild type and up-regulated in the human skeletal muscle in response to oral shilajit supplementation (Das et al., 2016; Chapman et al., 2017), suggesting that HMCN1 might be related to the muscle fibrosis. Also, we found the sequences of mature miRNAs in SCWrabs and ZIKArabs were different (Table 4). Although most of them had the difference at the 3′ arms of miRNA sequences, ocu-miR-411-5p had a base “A” added to the 5’ arm in SCWrabs, compared to ZIKArabs. These miRNA isoforms may result in different target genes in the formation of rabbit longissimus dorsi.

Although this study is the first time to study the profiles of circRNAs, miRNAs, and proteins in the muscle tissues of SCWrabs and ZIKArabs, we have to admit the existence of some limitations and gaps that could be filled by our future research. For example, more samples and timepoints can be applied to compare the difference of muscle development of SCWrabs and ZIKArabs; more data can be generated to support the noncoding RNAs, proteins, and variants that might be related to the muscle development and meat quality of rabbit; and more stringent statistical cut-offs can be used to narrow down the gene candidates.

 5 Conclusion

In summary, we identified 3159 circRNAs, 356 miRNAs, and 755 proteins in the muscle tissues of two rabbit breeds (SCWrabs and ZIKArabs), of which 267 circRNAs (120 up-regulated and 147 down-regulated), 3 miRNAs (down-regulated), and 29 proteins (19 up-regulated and 10 down-regulated) were differentially expressed. We next performed the interaction network analysis to identify the regulation networks between differentially expressed noncoding RNAs, mRNAs, and proteins in the rabbit muscle tissues. Interestingly, we found 10 miRNAs which had different dominant sequences in the muscle tissues of SCWrabs and ZIKArabs. Further, we analyzed the gene variants which might help explain the genetic background of the two rabbit lines. This is the first time to study noncoding RNA and protein profiles in the muscle tissues of SCWrabs and ZIKArabs. The noncoding RNAs identified in this work can be a valuable resource for future rabbit studies. More importantly, our findings will provide the basis of understanding the molecular mechanisms in the muscle development and differentiation and will benefit the rabbit breeding program.

 Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

 Ethics Statement

The animal study was reviewed and approved by the Medical Ethics Committee on Animal Research of the Sichuan Academy of Animal Science. Written informed consent was obtained from the owners for the participation of their animals in this study.

 Author Contributions

Conceptualization, XZ and XX; methodology, XZ, ZG, and CZ; software, XZ, YJ, and CY; validation, LK, JZ, and LT; formal analysis, XZ, XD, and DH; investigation, XZ; resources, ML, YR, and GW; data curation, CL; writing—original draft preparation, XZ; writing—review and editing, GW and XX; visualization, XZ; supervision, XX; project administration, XX; funding acquisition, XX. All authors have read and agreed to the published version of the manuscript.

 Funding

This research was supported by Livestock Breeding Project of Sichuan Province (2021YFYZ0033), Special Financial Operation of Sichuan Province (SASA202105), National Modern Agricultural Industrial Technology System of China (CARS-43-D-1), and Transformation Scientific Institution Achievements Project of Sichuan Province (2019JDZH0023). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

 Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

We acknowledge Zhenhua Zhuang from Chengdu Baseline Ltd. for his support on the bioinformatics analysis and suggestions.

 Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.777232/full#supplementary-material

References Bailey T. L. (2011). DREME: Motif Discovery in Transcription Factor ChIP-Seq Data. Bioinformatics 27 (12), 1653–1659. 10.1093/bioinformatics/btr261 Briggs M. M. Schachat F. (2002). The Superfast Extraocular Myosin (MYH13) Is Localized to the Innervation Zone in Both the Global and Orbital Layers of Rabbit Extraocular Muscle. J. Exp. Biol. 205 (Pt 20), 3133–3142. 10.1242/jeb.205.20.3133 Cao Y. You S. Yao Y. Liu Z.-J. Hazi W. Li C.-Y. (2018). Expression Profiles of Circular RNAs in Sheep Skeletal Muscle. Asian-australas J. Anim. Sci. 31 (10), 1550–1557. 10.5713/ajas.17.0563 Chapman M. A. Mukund K. Subramaniam S. Brenner D. Lieber R. L. (2017). Three Distinct Cell Populations Express Extracellular Matrix Proteins and Increase in Number during Skeletal Muscle Fibrosis. Am. J. Physiology-Cell Physiol. 312 (2), C131–C143. 10.1152/ajpcell.00226.2016 Chen J.-F. Mandel E. M. Thomson J. M. Wu Q. Callis T. E. Hammond S. M. (2006). The Role of microRNA-1 and microRNA-133 in Skeletal Muscle Proliferation and Differentiation. Nat. Genet. 38 (2), 228–233. 10.1038/ng1725 Chen M. Mithraprabhu S. Ramachandran M. Choi K. Khong T. Spencer A. (2019a). Utility of Circulating Cell-free RNA Analysis for the Characterization of Global Transcriptome Profiles of Multiple Myeloma Patients. Cancers 11 (6), 887. 10.3390/cancers11060887 Chen M. Xu R. Ji H. Greening D. W. Rai A. Izumikawa K. (2016). Transcriptome and Long Noncoding RNA Sequencing of Three Extracellular Vesicle Subtypes Released from the Human colon Cancer LIM1863 Cell Line. Sci. Rep. 6, 38397. 10.1038/srep38397 Chen M. Xu R. Rai A. Suwakulsiri W. Izumikawa K. Ishikawa H. (2019b). Distinct Shed Microvesicle and Exosome microRNA Signatures Reveal Diagnostic Markers for Colorectal Cancer. PLoS One 14 (1), e0210003. 10.1371/journal.pone.0210003 Chen Y. Chen Y. Shi C. Huang Z. Zhang Y. Li S. (2018). SOAPnuke: a MapReduce Acceleration-Supported Software for Integrated Quality Control and Preprocessing of High-Throughput Sequencing Data. Gigascience 7 (1), 1–6. 10.1093/gigascience/gix120 Cooper G. M. (2000). “Actin, Myosin, and Cell Movement,” in The Cell: A Molecular Approach. 2nd ed (Sunderland (MA): Sinauer Associates). Dalle Zotte A. Szendrő Z. (2011). The Role of Rabbit Meat as Functional Food. Meat Sci. 88 (3), 319–331. 10.1016/j.meatsci.2011.02.017 Das A. Datta S. Rhea B. Sinha M. Veeraragavan M. Gordillo G. (2016). The Human Skeletal Muscle Transcriptome in Response to Oral Shilajit Supplementation. J. Med. Food 19 (7), 701–709. 10.1089/jmf.2016.0010 Dong S. Xiong W. Yuan J. Li J. Liu J. Xu X. (2013). MiRNA-146a Regulates the Maturation and Differentiation of Vascular Smooth Muscle Cells by Targeting NF-Κb Expression. Mol. Med. Rep. 8 (2), 407–412. 10.3892/mmr.2013.1538 Giudice G. Sánchez-Cabo F. Torroja C. Lara-Pezzi E. (2016). ATtRACT-A Database of RNA-Binding Proteins and Associated Motifs. Database (Oxford) 2016, baw035. 10.1093/database/baw035 Greco S. Cardinali B. Falcone G. Martelli F. (2018). Circular RNAs in Muscle Function and Disease. Ijms 19 (11), 3454. 10.3390/ijms19113454 Grey F. Antoniewicz A. Allen E. Saugstad J. McShea A. Carrington J. C. (2005). Identification and Characterization of Human Cytomegalovirus-Encoded microRNAs. J. Virol. 79 (18), 12095–12099. 10.1128/JVI.79.18.12095-12099.2005 Güller I. Russell A. P. (2010). MicroRNAs in Skeletal Muscle: Their Role and Regulation in Development, Disease and Function. J. Physiol. 588 (Pt 21), 4075–4087. 10.1113/jphysiol.2010.194175 Guo T. Jou W. Chanturiya T. Portas J. Gavrilova O. McPherron A. C. (2009). Myostatin Inhibition in Muscle, but Not Adipose Tissue, Decreases Fat Mass and Improves Insulin Sensitivity. PLoS One 4 (3), e4937. 10.1371/journal.pone.0004937 Holcomb M. Ding Y.-H. Dai D. McDonald R. J. McDonald J. S. Kallmes D. F. (2015). RNA-sequencing Analysis of Messenger RNA/MicroRNA in a Rabbit Aneurysm Model Identifies Pathways and Genes of Interest. AJNR Am. J. Neuroradiol 36 (9), 1710–1715. 10.3174/ajnr.A4390 Hudson M. B. Rahnert J. A. Zheng B. Woodworth-Hobbs M. E. Franch H. A. Russ Price S. (2014). miR-182 Attenuates Atrophy-Related Gene Expression by Targeting FoxO3 in Skeletal Muscle. Am. J. Physiology-Cell Physiol. 307 (4), C314–C319. 10.1152/ajpcell.00395.2013 John B. Enright A. J. Aravin A. Tuschl T. Sander C. Marks D. S. (2004). Human MicroRNA Targets. Plos Biol. 2 (11), e363. 10.1371/journal.pbio.0020363 Kee H. J. Park E. W. Lee C. K. (2008). Characterization of Beef Transcripts Correlated with Tenderness and Moisture. Mol. Cell 25 (3), 428–437. Khanna N. Ge Y. Chen J. (2014). MicroRNA-146b Promotes Myogenic Differentiation and Modulates Multiple Gene Targets in Muscle Cells. PLoS One 9 (6), e100657. 10.1371/journal.pone.0100657 Kim H. K. Lee Y. S. Sivaprasad U. Malhotra A. Dutta A. (2006). Muscle-specific microRNA miR-206 Promotes Muscle Differentiation. J. Cell Biol 174 (5), 677–687. 10.1083/jcb.200603008 Kim S. Scheffler K. Halpern A. L. Bekritsky M. A. Noh E. Källberg M. (2018). Strelka2: Fast and Accurate Calling of Germline and Somatic Variants. Nat. Methods 15 (8), 591–594. 10.1038/s41592-018-0051-x Korfage J. A. M. Helmers R. de Goüyon Matignon M. van Wessel T. Langenbach G. E. J. van Eijden T. M. G. J. (2009). Postnatal Development of Fiber Type Composition in Rabbit Jaw and Leg Muscles. Cells Tissues Organs 190 (1), 42–52. 10.1159/000155226 Kuang L. Lei M. Li C. Zhang X. Ren Y. Zheng J. (2018). Identification of Long Non-coding RNAs Related to Skeletal Muscle Development in Two Rabbit Breeds with Different Growth Rate. Ijms 19 (7), 2046. 10.3390/ijms19072046 Kuang L. Xie X. Zhang X. Lei M. Li C. Ren Y. (2014). Expression Profiles of Myostatin, Myogenin, and Myosin Heavy Chain in Skeletal Muscles of Two Rabbit Breeds Differing in Growth Rate. Anim. Biotechnol. 25 (4), 223–233. 10.1080/10495398.2013.865639 Lee S.-J. (2004). Regulation of Muscle Mass by Myostatin. Annu. Rev. Cell Dev. Biol. 20, 61–86. 10.1146/annurev.cellbio.20.012103.135836 Lehnert S. A. Reverter A. Byrne K. A. Wang Y. Nattrass G. S. Hudson N. J. (2007). Gene Expression Studies of Developing Bovine Longissimus Muscle from Two Different Beef Cattle Breeds. BMC Dev. Biol. 7, 95. 10.1186/1471-213X-7-95 Li J.-H. Liu S. Zhou H. Qu L.-H. Yang J.-H. (2014). starBase v2.0: Decoding miRNA-ceRNA, miRNA-ncRNA and Protein-RNA Interaction Networks from Large-Scale CLIP-Seq Data. Nucl. Acids Res. 42 (Database issue), D92–D97. 10.1093/nar/gkt1248 Li R. Yu C. Li Y. Lam T.-W. Yiu S.-M. Kristiansen K. (2009). SOAP2: an Improved Ultrafast Tool for Short Read Alignment. Bioinformatics 25 (15), 1966–1967. 10.1093/bioinformatics/btp336 Li S. Zeng W. Li R. Hoffman L. C. He Z. Sun Q. (2018a). Rabbit Meat Production and Processing in China. Meat Sci. 145, 320–328. 10.1016/j.meatsci.2018.06.037 Li X. Yang L. Chen L.-L. (2018b). The Biogenesis, Functions, and Challenges of Circular RNAs. Mol. Cell 71 (3), 428–442. 10.1016/j.molcel.2018.06.034 Liang G. Yang Y. Niu G. Tang Z. Li K. (2017). Genome-wide Profiling of Sus scrofa Circular RNAs across Nine Organs and Three Developmental Stages. DNA Res. 24 (5), 523–535. 10.1093/dnares/dsx022 Liu Z. Huang J. Huo Y. Gong J. Zhang Y. Wei C. (2015). Identification of Proteins Implicated in the Increased Heart Rate in ShenSongYangXin-Treated Bradycardia Rabbits by iTRAQ-Based Quantitative Proteomics. Evidence-Based Complement. Altern. Med. 2015, 1–7. 10.1155/2015/385953 McKoy G. Léger M.-E. Bacou F. Goldspink G. (1998). Differential Expression of Myosin Heavy Chain mRNA and Protein Isoforms in Four Functionally Diverse Rabbit Skeletal Muscles during Pre- and Postnatal Development. Dev. Dyn. 211 (3), 193–203. 10.1002/(SICI)1097-0177(199803)211:3<193:AID-AJA1>3.0.CO;2-C Memczak S. Jens M. Elefsinioti A. Torti F. Krueger J. Rybak A. (2013). Circular RNAs Are a Large Class of Animal RNAs with Regulatory Potency. Nature 495 (7441), 333–338. 10.1038/nature11928 Noce A. Cardoso T. F. Manunza A. Martínez A. Cánovas A. Pons A. (2018). Expression Patterns and Genetic Variation of the Ovine Skeletal Muscle Transcriptome of Sheep from Five Spanish Meat Breeds. Sci. Rep. 8 (1), 10486. 10.1038/s41598-018-28760-9 Paci G. Preziuso G. D’Agata M. Russo C. Dalle Zotte A. (2013). Effect of Stocking Density and Group Size on Growth Performance, Carcass Traits and Meat Quality of Outdoor-Reared Rabbits. Meat Sci. 93 (2), 162–166. 10.1016/j.meatsci.2012.08.012 Park S. K. Gunawan A. M. Scheffler T. L. Grant A. L. Gerrard D. E. (2009). Myosin Heavy Chain Isoform Content and Energy Metabolism Can Be Uncoupled in Pig Skeletal Muscle. J. Anim. Sci. 87 (2), 522–531. 10.2527/jas.2008-1269 Pertea M. Kim D. Pertea G. M. Leek J. T. Salzberg S. L. (2016). Transcript-level Expression Analysis of RNA-Seq Experiments with HISAT, StringTie and Ballgown. Nat. Protoc. 11 (9), 1650–1667. 10.1038/nprot.2016.095 Qiu X. Gong L. Sun X. Guo J. Chodara A. M. (2011). Efficacy of Acupuncture and Identification of Tear Protein Expression Changes Using iTRAQ Quantitative Proteomics in Rabbits. Curr. Eye Res. 36 (10), 886–894. 10.3109/02713683.2011.601843 Qu S. Yang X. Li X. Wang J. Gao Y. Shang R. (2015). Circular RNA: A New star of Noncoding RNAs. Cancer Lett. 365 (2), 141–148. 10.1016/j.canlet.2015.06.003 Robinson M. D. McCarthy D. J. Smyth G. K. (2010). edgeR: a Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics 26 (1), 139–140. 10.1093/bioinformatics/btp616 Ruusunen M. Puolanne E. (2004). Histochemical Properties of Fibre Types in Muscles of Wild and Domestic Pigs and the Effect of Growth Rate on Muscle Fibre Properties. Meat Sci. 67 (3), 533–539. 10.1016/j.meatsci.2003.12.008 Sun J. Xie M. Huang Z. Li H. Chen T. Sun R. (2017). Integrated Analysis of Non-coding RNA and mRNA Expression Profiles of 2 Pig Breeds Differing in Muscle Traits. J. Anim. Sci. 95 (3), 1092–1103. 10.2527/jas.2016.086710.2527/jas2016.0867 Sun L. Bai M. Xiang L. Zhang G. Ma W. Jiang H. (2016). Comparative Transcriptome Profiling of Longissimus Muscle Tissues from Qianhua Mutton Merino and Small Tail Han Sheep. Sci. Rep. 6, 33586. 10.1038/srep33586 Sweeney H. L. Hammers D. W. (2018). Muscle Contraction. Cold Spring Harb Perspect. Biol. 10 (2), a023200. 10.1101/cshperspect.a023200 van Harten S. Cardoso L. A. (2010). Feed Restriction and Genetic Selection on the Expression and Activity of Metabolism Regulatory Enzymes in Rabbits. Animal 4 (11), 1873–1883. 10.1017/S1751731110001047 Wang J. Elzo M. A. Yan L. Chen S. Jia X. Zhang M. (2016). A Single Nucleotide Polymorphism in CAST Gene Is Associated with Meat Quality Traits in Rabbits. Anim. Sci. Pap. Rep. 24 (3), 269–278. Wang J. Hu Y. Elzo M. A. Shi Y. Jia X. Chen S. (2017a). Genetic Effect of Myf5 Gene in Rabbit Meat Quality Traits. J. Genet. 96 (4), 673–679. 10.1007/s12041-017-0822-7 Wang J. Li G. Elzo M. A. Yan L. Chen S. Jia X. (2015). A Novel Single Nucleotide Polymorphism of the POU1F1 Gene Associated with Meat Quality Traits in Rabbits. Ann. Anim. Sci. 15 (3), 611–620. 10.1515/aoas-2015-0015 Wang J. Shi Y. Elzo M. A. Su Y. Jia X. Chen S. (2017b). Myopalladin Gene Polymorphism Is Associated with Rabbit Meat Quality Traits. Ital. J. Anim. Sci. 16 (3), 400–404. 10.1080/1828051X.2017.1296333 Wang K. Li M. Hakonarson H. (2010). ANNOVAR: Functional Annotation of Genetic Variants from High-Throughput Sequencing Data. Nucleic Acids Res. 38 (16), e164. 10.1093/nar/gkq603 Xing J. Jing W. Jiang Y. (2013). Molecular Characterization and Expression Analysis of Fat Mass and Obesity-Associated Gene in Rabbit. J. Genet. 92 (3), 481–488. 10.1007/s12041-013-0298-z Yang J. Yang L. Li B. Zhou W. Zhong S. Zhuang Z. (2016). iTRAQ-Based Proteomics Identification of Serum Biomarkers of Two Chronic Hepatitis B Subtypes Diagnosed by Traditional Chinese Medicine. Biomed. Res. Int. 2016 (3290260), 1–29. 10.1155/2016/3290260 Zhang D. Li Y. Yao X. Wang H. Zhao L. Jiang H. (2016). miR-182 Regulates Metabolic Homeostasis by Modulating Glucose Utilization in Muscle. Cell Rep. 16 (3), 757–768. 10.1016/j.celrep.2016.06.040 Zhang F. Zhang R. Zhang X. Wu Y. Li X. Zhang S. (2018). Comprehensive Analysis of circRNA Expression Pattern and circRNA-miRNA-mRNA Network in the Pathogenesis of Atherosclerosis in Rabbits. Aging 10 (9), 2266–2283. 10.18632/aging.101541 Zhou L. Wei R. Zhao P. Koh S. K. Beuerman R. W. Ding C. (2013). Proteomic Analysis Revealed the Altered Tear Protein Profile in a Rabbit Model of Sjögren's Syndrome-Associated Dry Eye. Proteomics 13 (16), 2469–2481. 10.1002/pmic.201200230