Comparison of hair follicle density in high and low wool production rabbits
To characterize the hair follicle density, the follicle densities of the backs, abdomens, sides, and hips of Wan strain Angora rabbits with HWP and LWP were compared (Fig. 1). A morphological analysis showed that the hair follicle densities of backs, abdomens, sides, and hips were higher in the HWP group (Fig. 1a, b, c, d) than in the LWP group (Fig. 1e, f, g, h). The results demonstrate that a high hair follicle density leads to high wool production in Wan strain Angora rabbits.
Sequencing and assembly
Eight libraries of the HWP groups (H1, H2, H3, and H4) and LWP groups (L1, L2, L3, and L4) were constructed. For the HWP and LWP libraries, above 84,456,770 and 94,769,312 clean reads per sample were obtained, respectively (Table 1). Above 89.19% and 89.02% of the reads were aligned with the rabbit reference genome uniquely located by above 77.55% and 75.67% of the clean reads for the HWP and LWP libraries, respectively. Above 17,380,601 (52.78%) and 21,898,377 (46.66%) reads, respectively, were identified as protein-coding mRNAs of the HWP and LWP groups (Additional file 1: Table S1). The other types of reads amounted to 12,349,910 (36.07%) and 16,461,100 (39.19%) for HWP and LWP groups, respectively, and these reads may include lncRNAs (Additional file 1: Table S1).
Table 1. The analyses of reads mapped to the rabbit reference genome
Sample name
|
H1
|
H2
|
H3
|
H4
|
L1
|
L2
|
L3
|
L4
|
Total reads
|
95,426,664
|
84,456,770
|
149,301,000
|
89,038,716
|
123,865,064
|
94,769,312
|
109,570,388
|
130,245,954
|
Total mapped
|
85,943,449 (90.06%)
|
75,328,473 (89.19%)
|
133,972,487 (89.73%)
|
79,823,704 (89.65%)
|
110,263,589 (89.02%)
|
84,819,962 (89.5%)
|
99,272,110 (90.6%)
|
117,162,667 (89.95%)
|
Multiple mapped
|
11,938,587 (12.51%)
|
9,315,870 (11.03%)
|
13,827,931 (9.26%)
|
10,659,938 (11.97%)
|
16,534,984 (13.35%)
|
9,777,660 (10.32%)
|
14,500,967 (13.23%)
|
11,066,700 (8.5%)
|
Uniquely mapped
|
74,004,862 (77.55%)
|
66,012,603 (78.16%)
|
120,144,556 (80.47%)
|
69,163,766 (77.68%)
|
93,728,605 (75.67%)
|
75,042,302 (79.18%)
|
84,771,143 (77.37%)
|
106,095,967 (81.46%)
|
Reads map to '+'
|
36,872,804 (38.64%)
|
32,897,302 (38.95%)
|
60,000,375 (40.19%)
|
34,636,589 (38.9%)
|
46,454,860 (37.5%)
|
37,454,686 (39.52%)
|
42,357,439 (38.66%)
|
52,928,464 (40.64%)
|
Reads map to '-'
|
37,132,058 (38.91%)
|
33,115,301 (39.21%)
|
60,144,181 (40.28%)
|
34,527,177 (38.78%)
|
47,273,745 (38.17%)
|
37,587,616 (39.66%)
|
42,413,704 (38.71%)
|
53,167,503 (40.82%)
|
Non-splice reads
|
56,778,516 (59.5%)
|
51,943,073 (61.5%)
|
89,625,539 (60.03%)
|
51,395,189 (57.72%)
|
74,099,166 (59.82%)
|
57,787,882 (60.98%)
|
65,280,836 (59.58%)
|
79,655,006 (61.16%)
|
Splice reads
|
17,226,346 (18.05%)
|
14,069,530 (16.66%)
|
30,519,017 (20.44%)
|
17,768,577 (19.96%)
|
19,629,439 (15.85%)
|
17,254,420 (18.21%)
|
19,490,307 (17.79%)
|
26,440,961 (20.3%)
|
Reads mapped in proper pairs
|
69,994,050 (73.35%)
|
62,138,530 (73.57%)
|
113,304,356 (75.89%)
|
65,315,796 (73.36%)
|
88,424,106 (71.39%)
|
70,592,926 (74.49%)
|
80,480,580 (73.45%)
|
100,326,142 (77.03%)
|
Characterization of lncRNAs in rabbit skin tissue
The RNA-seq analysis produced 22,136 lncRNAs (Additional file 2: Table S2). The lncRNA transcripts included 10,692 lincRNAs (48.3%), 2,612 antisense lncRNAs (11.8%), and 8,832 intronic lncRNAs (39.9%) (Fig. 2a). The average length of the novel lncRNAs was considerably shorter than the mRNAs, but longer than the known lncRNAs (Fig. 2b). The exon numbers of the novel lncRNAs were less than the mRNAs while greater than the known lncRNAs (Fig. 2c). In addition, ORF size in novel lncRNAs was longer than that in annotated lncRNAs, but shorter than that in protein-coding genes (Fig. 2d).
Long noncoding RNAs and mRNAs expression profiles in rabbit skin tissue
The results showed that the expression levels of mRNAs were higher than those of lncRNAs (Additional file 3: Figure S1). 50 and 38 differentially expressed (DE) lncRNAs and genes, respectively, were screened in the LWP and HWP groups (Additional file 4: Table S3, Table S4). Of these lncRNAs and genes, 15 lncRNAs and 21 genes were upregulated, and 35 lncRNAs and 17 genes were downregulated in the LWP group. Hierarchical cluster analysis of lncRNA and mRNA expression levels between LWP and HWP groups revealed distinct expression patterns (Fig. 3).
Long noncoding RNA target prediction and functional analysis
The potential target genes of lncRNAs were predicted accordingly their position (co-location) and expression correlation (co-expression) with the protein-coding genes. Gene ontology (GO) analysis was applied to investigate the potential functions of the lncRNAs’ co-location and co-expression mRNAs on the regulation of hair follicle development and wool production (Fig. 4). The significance of enrichment of each GO term was assessed by P-value < 0.05, and then the GO terms were filtered by the enrichment scores (-Lg P-value). The GO enrichment analysis showed that the lncRNAs’ co-location mRNAs were significantly enriched in phospholipid, lipid metabolic, and epithelial cell apoptotic processes in the biological process category (Fig. 4a), while co-expression mRNAs were significantly enriched in the cellular metabolic, lipoprotein, lipid biosynthetic, lipid, and fatty acid transport processes (Fig. 4b). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis offered a reliable way of elucidating the candidate biological pathways that the integrated target genes were enriching. The cytokine-cytokine receptor interaction, chemokine signalling pathway and JAK-STAT signalling pathway were significantly involved in lncRNAs’ co-location mRNAs (Fig. 5a). In addition, pathways related to the biosynthesis of amino acids, arginine and proline metabolism, ether lipid metabolism, and the hedgehog signalling pathway were highly enriched by lncRNAs’ co-expression mRNAs (Fig. 5b). Therefore, the target genes of the DE lncRNAs between the LWP and HWP groups were related to lipid metabolism, amino acid synthesis, JAK-STAT, and the hedgehog signalling pathway. According to the functional enrichment analyses, five DE lncRNAs (LNC_002171, LNC_000797, LNC_005567, LNC_013595, and LNC_020367) were selected to construct regulatory networks (Fig. 6). LNC_002171 and LNC_000797 were involved in JAK-STAT and the hedgehog signalling pathway.
Validation of DE lncRNAs and mRNAs with quantitative real-time polymerase chain reaction
To validate the RNA-Seq results, LNC_000797, LNC_013595, LNC_020367, KRTAP15-1, TCHHL1, and ALOX15B were selected and their expression patterns in the LWP and HWP groups were examined by q-PCR. The results showed that the three DE lncRNAs and mRNAs were differentially expressed in the LWP and HWP groups. In addition, they exhibited a similar trend in the results of the RNA-seq and the q-PCR (Fig. 7). Therefore, the fragments per kilobase of transcript per million mapped reads (FPKM) obtained from RNA-seq could be reliably used to determine lncRNA and mRNA expression in the LWP and HWP groups.