Identification of Long Non-Coding RNAs Related to Skeletal Muscle Development in Two Rabbit Breeds with Different Growth Rate

Skeletal muscle development plays an important role in muscle quality and yield, which decides the economic value of livestock. Long non-coding RNAs (lncRNAs) have been reported to be associated with skeletal muscle development. However, little is revealed about the function of lncRNAs in rabbits’ muscle development. LncRNAs and mRNAs in two rabbit breeds (ZIKA rabbits (ZKR) and Qixin rabbits (QXR)) with different growth rates at three developmental stages (0 day, 35 days, and 84 days after birth) were researched by transcriptome sequencing. Differentially expressed lncRNAs and mRNAs were identified for two rabbit breeds at the same stages by DESeq package. Co-expression correlation analysis of differentially expressed lncRNAs and mRNAs were performed to construct lncRNA–mRNA pairs. To explore the function of lncRNAs, Gene Ontology (GO) analysis of co-expression mRNAs in lncRNA–mRNA pairs were performed. In three comparisons, there were 128, 109, and 115 differentially expressed lncRNAs, respectively. LncRNAs TCONS_00013557 and XR_518424.2 differentially expressed in the two rabbit breeds might play important roles in skeletal muscle development, for their co-expressed mRNAs were significantly enriched in skeletal muscle development related GO terms. This study provides potentially functional lncRNAs in skeletal muscle development of two rabbit breeds and might be beneficial to the production of rabbits.


Introduction
The meat of rabbits as a functional food provides dietetic properties and remarkable nutritive value [1,2]. It is becoming more and more popular to people on account of its characteristics of rich protein, low cholesterol, and low fat. Thus, improving the yield and quality of rabbits muscle might be the central task for breeding rabbits.
Most long non-coding RNAs (lncRNAs) generate at certain stages of biological development in a specific manner of cell or tissue. Emerging research showed that lncRNAs participated in the development of skeletal muscle in livestock. For example, Ramayo-Caldas et al. identified 55 differentially expressed lncRNAs between high intramuscular fat tissues and low intramuscular fat tissues in pigs by RNA-Sequencing, suggesting that lncRNAs were related to the fat metabolism of muscle [3]. Billerey et al. found 418 intergenic lncRNAs in all nine muscle samples of Limousin calves by RNA-Sequencing and validated 13 intergenic lncRNAs by Real-Time Polymerase Chain Reaction (RT-PCR) [4]. Meanwhile, part intergenic lncRNAs were found located in meat quality traits related loci [4]. Novel lncRNAs identified from chicken skeletal muscle by transcriptome sequencing presented differential expression level in a variety of tissues, and overexpression of lncRNA Gallus gallus (gga)-lnc-0181 in skeletal muscle might play a significant role in the muscle development of chicken [5]. Using RNA-Sequencing, several lncRNAs and protein coding genes associated with muscle development were screened in sheep [6]. All researchers above indicated that lncRNAs play important roles in muscle development.
However, there is little research on rabbits' lncRNAs associated with muscle development. The expression patterns of lncRNAs in the rabbits' skeletal muscle development remain widely unknown. Thus, we detected the expression patterns of lncRNAs and mRNAs in two rabbit breeds differing in growth rate at three developmental stages (0 day, 35 days, and 84 days after birth). The potential lncRNAs related to muscle development in two different rabbit breeds were predicted according to the function of corresponding co-expressed mRNAs with the lncRNAs. The study will provide potential lncRNAs related to muscle development of rabbits. It will also provide important data for studying the molecular mechanism of different varieties feeding rabbits' growth difference and promoting the production of the meat rabbits.

Sample Information
The weight of two rabbit breeds at three developmental stages is shown in Figure 1. The weight of ZIKA rabbits (ZKR) was significantly higher than that of Qixin rabbits (QXR), suggesting that the two kinds of rabbits differed in growth rate and are suited for researching the molecular mechanism of muscle growth and development (all p < 0.05). The weight of ZIKA rabbits (ZKR) and Qixin rabbits (QXR) at three development stages. S1, S2, and S3 refer to the age of 0 day, 35 days, and 84 days after birth, respectively. * and ** refer to the statistically significant difference (p < 0.05) and extremely significant difference (p < 0.001), respectively.

Reads Filtering and Mapping
The filtering rate of each sample was greater than 90%. The Q30 was not less than 91.4%. After filtering, on average, 95,660,601, 94,297,177, 90,539,959, 97,386,913, 97,460,159, and 91,414,542 clean reads were obtained for three samples each of ZKR_S1, ZKR_S2, ZKR_S3, QXR_S1, QXR_S2, and QXR_S3, respectively, and more than 90% were mapped to the Oryctolagus cuniculus reference genome (OryCun2.0) ( Table 1). All these results suggested that the data of RNA-Sequencing were quietly credible. The weight of ZIKA rabbits (ZKR) and Qixin rabbits (QXR) at three development stages. S1, S2, and S3 refer to the age of 0 day, 35 days, and 84 days after birth, respectively. * and ** refer to the statistically significant difference (p < 0.05) and extremely significant difference (p < 0.001), respectively.

Reads Filtering and Mapping
The filtering rate of each sample was greater than 90%. The Q30 was not less than 91.4% .  After filtering, on average, 95,660,601, 94,297,177, 90,539,959, 97,386,913, 97,460,159, and 91,414,542  clean reads were obtained for three samples each of ZKR_S1, ZKR_S2, ZKR_S3, QXR_S1, QXR_S2, and QXR_S3, respectively, and more than 90% were mapped to the Oryctolagus cuniculus reference genome (OryCun2.0) ( Table 1). All these results suggested that the data of RNA-Sequencing were quietly credible.

Principal Component Analysis (PCA) and Differential Expression Analysis
Both the results of PCA for lncRNAs and mRNAs showed that the samples (ZKR and QXR) with the same stages (S1, S2, and S3) were more similar ( Figure 3A). The numbers of differentially expressed lncRNAs and mRNAs between ZKR and QXR at S1, S2, and S3 are shown in Figure 3B. A total of 128, 109, and 115 differentially expressed lncRNAs were identified between ZKR and QXR at S1, S2, and S3, respectively. Heatmaps of differentially expressed lncRNAs in the comparisons of ZKR_S1 vs. QXR_S1, ZKR_S2 vs. QXR_S2, and ZKR_S3 vs. QXR_S3 suggested that the samples of QXR and ZKR with the same stages were distinguished by clustering ( Figure 3C).

Principal Component Analysis (PCA) and Differential Expression Analysis
Both the results of PCA for lncRNAs and mRNAs showed that the samples (ZKR and QXR) with the same stages (S1, S2, and S3) were more similar ( Figure 3A). The numbers of differentially expressed lncRNAs and mRNAs between ZKR and QXR at S1, S2, and S3 are shown in Figure 3B. A total of 128, 109, and 115 differentially expressed lncRNAs were identified between ZKR and QXR at S1, S2, and S3, respectively. Heatmaps of differentially expressed lncRNAs in the comparisons of ZKR_S1 vs. QXR_S1, ZKR_S2 vs. QXR_S2, and ZKR_S3 vs. QXR_S3 suggested that the samples of QXR and ZKR with the same stages were distinguished by clustering ( Figure 3C). lncRNAs and mRNAs (B) and heatmaps of differentially expressed lncRNAs (C) in three comparisons. ZKR: ZIKA rabbits; QXR: Qixin rabbits; S1, S2, and S3 refer to the age of 0 day, 35 days, and 84 days after birth, respectively.

lncRNA-mRNA Co-Regulated Pairs
Co-expression correlations of differentially expressed lncRNA and differentially expressed mRNA from each comparison were performed according to the fragments per kilobase per million reads (FPKM). The co-expressed mRNAs in each lncRNA-mRNA co-regulated pair were selected to ZKR: ZIKA rabbits; QXR: Qixin rabbits; S1, S2, and S3 refer to the age of 0 day, 35 days, and 84 days after birth, respectively.

lncRNA-mRNA Co-Regulated Pairs
Co-expression correlations of differentially expressed lncRNA and differentially expressed mRNA from each comparison were performed according to the fragments per kilobase per million reads (FPKM). The co-expressed mRNAs in each lncRNA-mRNA co-regulated pair were selected to explore the main functional role of lncRNAs. Figure 4 shows the networks of lncRNAs TCONS_00013557 and XR_518424.2 with the corresponding co-expression mRNAs. explore the main functional role of lncRNAs. Figure 4 shows the networks of lncRNAs TCONS_00013557 and XR_518424.2 with the corresponding co-expression mRNAs.

Gene Ontology (GO) Analysis for Co-Expression mRNA of Each lncRNA
Based on the GO results, a total of 29 lncRNAs, whose co-expressed mRNAs had the most GO terms and the enriched mRNA ≥5, were selected for three comparisons with the same stages ( Table  2). In order to identify muscle development related lncRNAs in rabbits differing in growth rate, lncRNAs whose co-expressed mRNA enriched in the skeletal muscle development related GO terms (embryonic skeletal joint morphogenesis, embryonic skeletal system morphogenesis, and skeletal muscle tissue development) were selected from the 29 lncRNAs. For the comparison of ZKR_S1 vs. QXR_S1, lncRNA TCONS_00013557 was selected because of its co-expressed mRNAs enriching in the GO terms of collagen fibril organization, proteoglycan metabolic process, embryonic skeletal joint morphogenesis, extracellular matrix organization, collagen catabolic process, and embryonic skeletal system morphogenesis on biological process ( Figure 5A). As for cellular component, GO terms were dominantly composed of proteinaceous extracellular matrix and extracellular matrix. Within molecular function category, GO terms were significantly composed of extracellular matrix structural constituent conferring tensile strength, extracellular matrix structural constituent, and protein binding. Similarly, lncRNA XR_518424.2 was identified for the comparison of ZKR_S2 vs. QXR_S2 on account of its co-expressed mRNAs mainly enriching in the GO term of skeletal muscle tissue development ( Figure 5B). For the comparison of ZKR_S3 vs. QXR_S3, no lncRNAs were selected. The enriched mRNAs of the GO terms are shown in Table 3. Co-expression mRNAs [Odd-skipped related 2 (Osr2), Collagen type II α1 (Col2a1), and Collagen type XI α1 (Col11a1)] of TCONS_00013557 and co-expression mRNAs [Vestigial-like 2 (Vgll2), Caveolin 1 (Cav-1), and Hoxd10] of XR_518424.2 mainly enriched skeletal muscle development related GO terms (Table 3).

Gene Ontology (GO) Analysis for Co-Expression mRNA of Each lncRNA
Based on the GO results, a total of 29 lncRNAs, whose co-expressed mRNAs had the most GO terms and the enriched mRNA ≥5, were selected for three comparisons with the same stages ( Table 2). In order to identify muscle development related lncRNAs in rabbits differing in growth rate, lncRNAs whose co-expressed mRNA enriched in the skeletal muscle development related GO terms (embryonic skeletal joint morphogenesis, embryonic skeletal system morphogenesis, and skeletal muscle tissue development) were selected from the 29 lncRNAs. For the comparison of ZKR_S1 vs. QXR_S1, lncRNA TCONS_00013557 was selected because of its co-expressed mRNAs enriching in the GO terms of collagen fibril organization, proteoglycan metabolic process, embryonic skeletal joint morphogenesis, extracellular matrix organization, collagen catabolic process, and embryonic skeletal system morphogenesis on biological process ( Figure 5A). As for cellular component, GO terms were dominantly composed of proteinaceous extracellular matrix and extracellular matrix. Within molecular function category, GO terms were significantly composed of extracellular matrix structural constituent conferring tensile strength, extracellular matrix structural constituent, and protein binding. Similarly, lncRNA XR_518424.2 was identified for the comparison of ZKR_S2 vs. QXR_S2 on account of its co-expressed mRNAs mainly enriching in the GO term of skeletal muscle tissue development ( Figure 5B). For the comparison of ZKR_S3 vs. QXR_S3, no lncRNAs were selected. The enriched mRNAs of the GO terms are shown in Table 3. Co-expression mRNAs [Odd-skipped related 2 (Osr2), Collagen type II α1 (Col2a1), and Collagen type XI α1 (Col11a1)] of TCONS_00013557 and co-expression mRNAs [Vestigial-like 2 (Vgll2), Caveolin 1 (Cav-1), and Hoxd10] of XR_518424.2 mainly enriched skeletal muscle development related GO terms (Table 3).

Validation of the Selected lncRNAs and Co-Expression mRNAs
The RNA-Sequencing results of the selected lncRNAs and co-expression mRNAs are shown in Figure 6A,B. The expression levels of the selected genes were validated by RT-PCR. The RT-PCR results confirmed that lncRNA TCONS_00013557, Osr2, Col2a1, and Col11a1 were higher in ZKR than in QXR at S1 ( Figure 6C). The lncRNA XR_518424.2 and its co-expressed mRNAs (Vgll2 and Cav1) were all lower in ZKR than in QXR at S2, whereas co-expressed mRNA Hoxd10 was higher expressed in ZKR at S2 by RT-PCR ( Figure 6D). All the RT-PCR results were in agreement with the RNA-Sequencing results. Among these genes, TCONS_00013557, Col2a1, XR_518424.2, and Cav1 were significantly differentially expressed between ZKR and QXR. Correlation analysis showed significantly positive correlation in fold change data between RT-PCR and RNA-Sequencing (correlation coefficient R = 0.737, p < 0.05), confirming our transcriptome sequencing analysis ( Figure 6E).

Validation of the Selected lncRNAs and Co-Expression mRNAs
The RNA-Sequencing results of the selected lncRNAs and co-expression mRNAs are shown in Figure 6A,B. The expression levels of the selected genes were validated by RT-PCR. The RT-PCR results confirmed that lncRNA TCONS_00013557, Osr2, Col2a1, and Col11a1 were higher in ZKR than in QXR at S1 ( Figure 6C). The lncRNA XR_518424.2 and its co-expressed mRNAs (Vgll2 and Cav1) were all lower in ZKR than in QXR at S2, whereas co-expressed mRNA Hoxd10 was higher expressed in ZKR at S2 by RT-PCR ( Figure 6D). All the RT-PCR results were in agreement with the RNA-Sequencing results. Among these genes, TCONS_00013557, Col2a1, XR_518424.2, and Cav1 were significantly differentially expressed between ZKR and QXR. Correlation analysis showed significantly positive correlation in fold change data between RT-PCR and RNA-Sequencing (correlation coefficient R = 0.737, p < 0.05), confirming our transcriptome sequencing analysis ( Figure  6E). Figure 6. The expression levels of lncRNAs TCONS_00013557 and XR_518424.2 and the corresponding co-expression mRNAs by transcriptome sequencing and RT-PCR. The expression levels of lncRNA TCONS_00013557 and its co-expressed mRNAs (Osr2, Col2a1, and Col11a1) at the stage of S1 by RNA-Sequencing (A) and RT-PCR (C). The expression levels of XR_518424.2 and its co-expressed mRNAs (Vgll2, Cav1, and Hoxd10) at the stage of S2 by RNA-Sequencing (B) and RT-PCR (D). (E) Linear regression analysis of fold change (FC) data between RT-PCR and RNA-Sequencing. ZKR: ZIKA rabbits; QXR: Qixin rabbits; S1, S2, and S3 refer to the age of 0 day, 35 days, and 84 days after birth, respectively. * refers to the statistically significant difference (p < 0.05). Black dots represent log2 transformed FC values of a gene obtained from RT-PCR (X-axis) and RNA-Sequencing (Y-axis). R: correlation coefficient. Figure 6. The expression levels of lncRNAs TCONS_00013557 and XR_518424.2 and the corresponding co-expression mRNAs by transcriptome sequencing and RT-PCR. The expression levels of lncRNA TCONS_00013557 and its co-expressed mRNAs (Osr2, Col2a1, and Col11a1) at the stage of S1 by RNA-Sequencing (A) and RT-PCR (C). The expression levels of XR_518424.2 and its co-expressed mRNAs (Vgll2, Cav1, and Hoxd10) at the stage of S2 by RNA-Sequencing (B) and RT-PCR (D). (E) Linear regression analysis of fold change (FC) data between RT-PCR and RNA-Sequencing. ZKR: ZIKA rabbits; QXR: Qixin rabbits; S1, S2, and S3 refer to the age of 0 day, 35 days, and 84 days after birth, respectively. * refers to the statistically significant difference (p < 0.05). Black dots represent log2 transformed FC values of a gene obtained from RT-PCR (X-axis) and RNA-Sequencing (Y-axis). R: correlation coefficient.

Discussion
Emerging research suggests that lncRNAs play a significant role in the muscle development of pigs [3], bull calves [4], chickens [5], sheep [6], and humans [7]. Nevertheless, the lncRNAs analysis of rabbits' muscle has not been studied. In the present work, we identified lncRNAs and mRNAs in muscle tissues of two rabbit breeds by transcriptome sequencing. We explored the rabbits' muscle lncRNAs from structure and expression level. To reveal the functions, the co-expression correlations of differentially expressed lncRNA and differentially expressed mRNA from each comparison and GO analysis of co-expression mRNAs were performed.
In this study, a total of 1997 lncRNA transcripts were found. Most lncRNAs were longer than 2000 bp and contained 2 exons. The PCA results of mRNA and lncRNA showed higher similarity in the same development stage between the two rabbit species, implying that the comparison was reasonable. Comparing the QXR to ZKR, 128, 109, and 115 differentially expressed lncRNAs were identified between ZKR and QXR at S1, S2, and S3, respectively.
There is a vital interplay between muscle cells and the extracellular matrix in skeletal muscle development [8]. The extracellular matrix, which was mainly composed of collagens, proteoglycans, and glycoproteins, maintained the integrity of skeletal muscle [9]. While inhibiting collagen synthesis, myoblasts differentiation in vitro was blocked [10,11]. Proteoglycans can regulate collagen fibrillogenesis and suppress cell growth [12]. Melo et al. revealed that proteoglycans were essential for skeletal muscle differentiation [13]. Interestingly, co-expression mRNAs of TCONS_00013557 significantly enriched in the GO terms of collagen and proteoglycans related GO terms, indicating that co-expression mRNAs and the corresponding lncRNA TCONS_00013557 affect skeletal muscle development via altering the formation of collagens and proteoglycans in extracellular matrix in development of stage 1.
Co-expression correlation analysis combining GO analysis showed that co-expression mRNAs (Osr2, Col2a1, and Col11a1) of TCONS_00013557, which was differentially expressed between QXR_S1 and ZKQ_S1, significantly enriched in the GO terms such as collagen fibril organization, proteoglycan metabolic process, embryonic skeletal joint morphogenesis, and embryonic skeletal system morphogenesis. Osr2, a zinc-finger transcription factor, was expressed in numerous murine tissues including skeletal muscle tissues and had transcriptional activity involving in postnatal development [14,15]. Col2a1 was expressed in human rotator cuff-derived mesenchymal stem cells, which might be a cell source for muscle repair [16]. Its mRNA level was increased throughout the process of chondrogenic differentiation [17]. Col11a1 was regulated by a transcription activator FP9C, which was associated with cell differentiation in myoblasts and osteoblasts [18]. Consistent with the reports above, these co-expressed mRNAs with lncRNA TCONS_00013557 were all expressed higher in ZKR than in QXR at S1, suggesting that differentially expressed lncRNA TCONS_00013557 was likely involved in the skeletal muscle development of rabbits with different growth rates.
Co-expression mRNAs (Vgll2, Cav-1, and Hoxd10) of XR_518424.2 differentially expressed between QXR_S2 and ZKQ_S2 mainly enriched in skeletal muscle tissue development, in which Vgll2, Cav-1, and Hoxd10 were significantly enriched. Vgll2 played an important role in skeletal muscle differentiation and development myotubes. It was a cofactor of transcription enhancer factor 1 and myocyte enhancer factor 2 [19,20]. Vgll2 expression was skeletal muscle-specific in mammalian adult tissues and increased in differential myotubes [19]. Similarly, Vgll-2 was expressed in different sites of chick skeletal myogenesis and related to skeletal muscle differentiation as downstream gene of myogenic factor [21]. Vgll2 defecting resulted in an increase in fast-twitch fibers' numbers and Myh7 down-regulated in mice, suggesting that Vgll2 might be related to slow muscle fibers' programing [22]. Cav-1 was detected in various adult monkey tissues, including skeletal muscle, and was co-located with dystrophin on sarcolemma by immunohistochemistry [23]. Cav-1 was highly expressed in masticatory muscles of murine X-linked muscular dystrophy, which was with muscle injury and progressive muscle weakness caused by lack of dystrophin [24]. Hoxd10 was found differentially expressed in the back muscle of a mouse, absented Myf5, and a regulator for motor neuron development [25]. For all above, low expression of Vgll2 and Cav-1 might promote muscle development; differentially expressed Hoxd10 was related to muscle development. In this work, we found that the expression of Vgll2 and Cav-1 were lower in ZKR than in QXR at S2, whereas higher expression of Hoxd10 was found in ZKR than in QXR at S2. These mRNAs related to skeletal muscle development were co-expressed with the lncRNA XR_518424.2. Thus, XR_518424.2 probably participates in skeletal muscle development of rabbits with different growth rates.
In conclusion, we identified several lncRNAs and co-expressed genes related to skeletal muscle development in two rabbit breeds differing in growth rates. The co-expressed genes were mainly enriched in skeletal muscle development related GO terms. The lncRNAs (TCONS_00013557 and XR_518424.2) and co-expressed genes (Col2a1 and Cav-1) were validated to differentially expressed genes significantly by RT-PCR, confirming the important role of themselves and corresponding lncRNAs. This work provides candidate lncRNAs that may be used to explore the function of lncRNAs in the muscle development of rabbits. Further studies should be performed to validate the function and analyze the mechanism in detail.

Sample Collection
The meat rabbits of used in the experiments-German ZIKA rabbits (ZKR) and Sichuan native Qixin Rabbits (QXR)-were obtained from the rabbit farms of Sichuan Animal Sciences Academy in Chengdu, Sichuan, China. All rabbits used (all were male and belonged to the same family in each breed) were raised under the condition with the same diet and environmental temperature and given free access to water and food. The weight of each rabbit was recorded before longissimus muscle tissues were collected. The longissimus muscle tissues were collected from the ZKR and QXR at the age of 0 day (S1), 35 days (S2), and 84 days (S3) after birth (n = 3 for each stage and for each breed), respectively, and saved in liquid nitrogen immediately. The experiment was conducted according to the National Institutes of Health (NIH) Guidelines and National Research Council's publication "Guide for Care and Use of Laboratory Animals". 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 (the approval date: 22 March 2017).

RNA Isolation, Library Construction, and Sequencing
Total RNA of the longissimus muscle tissues were extracted with Trizol regent (Invitrogen, Carlsbad, CA, USA) and quality qualified RNA were treated with TruSeq Stranded Total RNA with Ribo-Zero Gold kit (Illumina, San Diego, CA, USA) to eliminate the ribosomal RNA. Strand-specific RNA-seq (ssRNA-seq) libraries were prepared according the manufacturer's instructions using the Illumina Standard RNA sample library preparation kit (Illumina, San Diego, CA, USA). After quantification using the Agilent 2100 bioanalyzer (Agilent, Santa Clara, CA, USA), the strand-specific libraries were sequenced on an Illumina HiSeq X ten instrument that generated paired-end reads of 150 nucleotides. Library construction and Illumina sequencing were performed by OE Biotech CO., LTD (Shanghai, China). The raw data have been deposited in the Sequence Read Archive database at the NCBI under the accession number SRP150254.

PCA and Differential Expression Analysis of lncRNAs and mRNAs
PCA was employed to explore the correlation among samples according to the expression level of lncRNAs and mRNAs, respectively. The differentially expressed lncRNA or mRNA for three comparisons (ZKR_S1 vs. QXR_S1, ZKR_S2 vs. QXR_S2, and ZKR_S3 vs. QXR_S3) were performed with the DESeq package (version 1.18.0, available online: http://www.bioconductor.org/ packages/release/bioc/html/DESeq.html), respectively. To control the false discovery rate, the p-value was adjusted by Benjamini and Hochberg's approach. The lncRNAs or mRNAs with the adjusted p-value < 0.05 and |log2(fold change)| > 1 were considered as differentially expressed genes.

Co-Expression Correlations of Differentially Expressed lncRNA and mRNA
To explore the functional role of lncRNA, the co-expression correlations of differentially expressed lncRNA and differentially expressed mRNA from each comparison were performed according to the FPKMs. Then the lncRNA-mRNA co-regulated pairs (Pearson's correlation coefficient >0.8 and p-value < 0.05) were screened for Gene Ontology (GO) analysis.

GO Enrichment Analysis
To explore the main functional role of lncRNAs in the muscle development of rabbits, the mRNAs in all lncRNA-mRNA co-regulated pairs were annotated by GO for differentially expressed lncRNAs. The GO terms with p-value < 0.05 were considered as significantly enriched. The top 10 lncRNAs whose co-expressed mRNAs had the most GO terms and the enriched mRNA ≥ 5 were screened. The GO enrichment graphs were drawn for the co-expressed mRNA of the selected lncRNAs in these selected lncRNAs.

RT-PCR
One microgram RNA was transcribed to cDNA. RT-PCR was determined using SYBR-Green PCR master mix kit (Applied Biosystems, Inc., Foster City, CA, USA) and performed on an ABI QuantStudio™ 6 Flex System (Applied Biosystems, Inc., Foster City, CA, USA) with the amplification conditions: one cycle of 95 • C for 10 min, followed by 45 cycles of 95 • C for 15 s, 60 • C for 60 s, and 95 • C for 15 s. The primers for amplication of genes and the internal control Gapdh are shown in Table 4. Three independent experiments were employed to detect the relative expression level. The relative expression level was calculated as below: relative quantification = 2 −∆∆Ct .

Statistical Analysis
The statistical significance was analyzed by the software of SPSS 21.0 (IBM Corp., Armonk, NY, USA). The experiment data was provided as mean value ± standard deviation. Difference between the groups was analyzed with one-way analysis of variance. p < 0.05 and p < 0.001 refer to the statistically significant difference (*) and extremely significant difference (**), respectively. The Pearson correlation analysis was performed to evaluate the fold change data between RT-PCR and RNA-Sequencing.
Author Contributions: L.K. conceived and designed the experiments, performed experiments, acquired and analyzed the data, wrote the article, and reviewed manuscript; M.L. performed experiments and sample processing, acquired and analyzed the data; C.L. performed experiments and sample processing, acquired and analyzed the data; X.Z. performed experiments and sample processing, acquired and analyzed the data; Y.R. performed experiments and sample processing, and analyzed the data; J.Z. performed experiments and sample processing, and analyzed the data; Z.G. performed experiments and sample processing, and reviewed manuscript; C.Z. performed experiments and sample processing, and reviewed manuscript; C.Y. performed experiments and sample processing, and reviewed manuscript; X.M. performed experiments and sample processing, and reviewed manuscript; M.F. performed experiments and sample processing, and reviewed manuscript; X.X. conceived and designed the experiments, analyzed the data, wrote the article, and reviewed manuscript.