Analysis of histology and long noncoding RNAs involved in the rabbit hair follicle density using RNA sequencing

Background Hair follicle density influences wool fibre production, which is one of the most important traits of the Wan Strain Angora rabbit. However, molecular mechanisms regulating hair follicle density have remained elusive. Results In this study, hair follicle density at different body sites of Wan Strain Angora rabbits with high and low wool production (HWP and LWP) was investigated by histological analysis. Haematoxylin-eosin staining showed a higher hair follicle density in the skin of the HWP rabbits. The long noncoding RNA (lncRNA) profile was investigated by RNA sequencing, and 50 and 38 differentially expressed (DE) lncRNAs and genes, respectively, were screened between the HWP and LWP groups. A gene ontology analysis revealed that phospholipid, lipid metabolic, apoptotic, lipid biosynthetic, and lipid and fatty acid transport processes were significantly enriched. Potential functional lncRNAs that regulate lipid metabolism, amino acid synthesis, as well as the Janus kinase (JAK)-signal transducer and activator of transcription (STAT) and hedgehog signalling pathways, were identified. Consequently, five lncRNAs (LNC_002171, LNC_000797, LNC_005567, LNC_013595, and LNC_020367) were considered to be potential regulators of hair follicle density and development. Three DE lncRNAs and genes were validated by quantitative real-time polymerase chain reaction (q-PCR). Conclusions LncRNA profiles provide information on lncRNA expression to improve the understanding of molecular mechanisms involved in the regulation of hair follicle density. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-07398-4.


Background
The Angora rabbit is an economically important livestock breed in several countries, especially in China and France. Wool production is one of the most important traits in Angora rabbits. The fur quality of rabbits is largely dependent on hair density, and hair follicle density determines hair density [1,2]. The wool production of Wan strain Angora rabbits is in uenced by genetics and environment.
For Angora rabbits under the same environmental conditions, gender, body site, and the month of age are closely related to wool bre production [3]. Genetic factors that in uence wool bre production are bre diameter, length, neness, and the bre density [3][4][5][6][7]. The mean hair follicle density depends on the skin area. The development of wool follicles occurs during prenatal life and no new hair follicles are formed after birth, implying that hair density in an adult rabbit will depend on how much that particular body part grows after the formation of the hair follicles [7,8]. Correspondingly, hair follicle density and other wool characteristics are highly variable over the human and rabbit body [7,9,10]. The molecular mechanism underlying hair follicle density in rabbit skin and hair follicle development remains unclear.
Hair follicle development is a complex morphogenetic process and undergoes periodic stages of growth (anagen), regression (catagen), and relative quiescence (telogen) [11][12][13]. The process of hair follicle formation and differentiation relies on many regulating molecules including messenger RNAs (mRNAs) and micro RNAs (miRNAs) [14][15][16], as well as a variety of signalling systems, such as the Wnt, Notch, bone morphogenetic protein (BMP), and broblast growth factor (FGF) pathways [17][18][19][20][21][22]. LncRNAs are RNA transcripts longer than 200 nucleotides that lack open reading frames (ORF) and protein-coding capabilities [23]. They regulate protein-coding gene expression at posttranscriptional and transcriptional levels [24,25]. It is generally known that lncRNAs are also involved in the regulation of the hair follicle development and skin homeostasis [26][27][28]. RP11-766N7.3, H19 and HOTAIR are speci c lncRNAs that are involved in Wnt signalling to regulate hair follicle development [29]. Strand-speci c RNA sequencing (ssRNA-seq) also showed that lncRNAs may be considered as potential candidate markers for further study on the molecular mechanisms of hair follicle initiation [30]. However, hair follicle density-related lncRNAs in rabbits have not been pro led so far.
In this study, the RNA-seq based approach was used to determine lncRNA expression levels in Angora rabbits with high wool production (HWP) and low wool production (LWP) after hair follicle density analysis. The results should provide fundamental resources to reveal the regulatory function of lncRNAs in hair follicle density in rabbits, as well as supply information for understanding human hair disorders such as hypotrichosis.

Results
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. 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 identi ed as protein-coding mRNAs of the HWP and LWP groups (Additional le 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 le 1: Table S1).  (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 pro les in rabbit skin tissue
The results showed that the expression levels of mRNAs were higher than those of lncRNAs (Additional le 2: Figure S1). 50 and 38 differentially expressed (DE) lncRNAs and genes, respectively, were screened in the LWP and HWP groups (Supplementary Data 2). 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 (coexpression) 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 signi cance of enrichment of each GO term was assessed by P-value < 0.05, and then the GO terms were ltered by the enrichment scores (-LgP). The GO enrichment analysis showed that the lncRNAs' co-location mRNAs were signi cantly enriched in phospholipid, lipid metabolic, and epithelial cell apoptotic processes (Fig. 4a), while co-expression mRNAs were signi cantly 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 signi cantly 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' coexpression 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, ve 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.

Discussion
Wool density is one of the most important indices to evaluate the quality of the fur of the Wan strain Angora rabbit [6]. Hair follicle density determines wool density [1]. The quality of fur is associated mainly with the traits of the hair follicles [31]. To characterize the hair follicle density, the follicle density of the backs, abdomens, sides, and hips of Wan Strain Angora rabbits with HWP and LWP was compared (Fig. 1). A morphological analysis showed that the hair follicle density of backs, abdomens, sides, and hips of the HWP group was higher compared to the LWP group (Fig. 1). The results demonstrated that high hair follicle density contributed to high wool production in the Wan strain Angora rabbit. In French Angora rabbits, divergent selection of total eece weight led to a positive difference of 0.55 genetic standard deviation for secondary to primary follicle ratio (S/P), although a low genetic correlation existed between them [32].
The formation of hair follicles is divided into prenatal hair morphogenesis and the postnatal hair cycle [33]. Once established during embryogenesis, hair follicle density is permanently xed in postnatal life, and the hair follicle location eventually becomes xed as a result of anchoring in the subcutis [34]. LncRNAs are widely involved in various biological processes, including the hair follicle cycle [35,36]. The lncRNA and mRNA expression pro les were compared in the dorsal skin of LWP and HWP rabbits, and 50 and 38 DE lncRNAs and genes were obtained, respectively. These lncRNAs and genes might play crucial roles in regulating hair follicle density, and their differential expression might be the reason for differences in hair follicle density and wool production between HWP and LWP rabbits. Liu et al. (2020) analysed the miRNA effect on hair follicle density in the Rex rabbit [37], but lncRNA related to hair density in rabbits has only been done in the present study.
The GO analysis showed that the DE lncRNAs are potential regulators of phospholipid, lipid metabolic, epithelial cell apoptotic, lipid biosynthetic, and lipid and fatty acid transport processes. Keratin-associated proteins (KRTAPs) play a critical role in cross-linking the keratin intermediate laments to build a hair shaft [38]. KRTAP7-1, KRTAP8-1, and KRTAP15-1 were predicted as the targets of LNC_005567 in this study. KRTAP7-1 is involved in supporting the mechanical strength and shape of hair [38]. KRTAP15-1 is expressed in secondary follicles in the skin and associated with bre diameter [39]. COL3A1 and LOXL4 were the target genes of LNC_013595 and LNC_020367, respectively. COL3A1 is one of collagens forming different extracellular matrix (ECM) components [40]. The lysyl oxidase like 4 (LOLX4) enzyme is responsible for initiating covalent cross-linking in collagen brils and is involved in providing additional mechanical strength to the ECM [41,42]. The amount of ECM per cell contributes to the volume of the dermal papilla [43]. Hedgehog and JAK-STAT signalling pathway were signi cantly enriched by target genes ENSOCUG00000021211 and ENSOCUG00000023782 of LNC_002171 and LNC_000797, respectively. The hedgehog signalling pathway is correlated to the initiation of hair follicle formation and is a pivotal growth signal for dermal papilla maturation and growth [34,44,45]. The JAK-STAT signalling pathway is involved in maintaining the quiescence of hair follicles during telogen [46], and JAK-STAT inhibition contributes to the promotion of hair growth and the activation of hair follicle stem cells [47]. These ndings demonstrate that LNC_002171, LNC_000797, LNC_005567, LNC_013595, and LNC_020367 are potentially important regulators of hair follicle density and development.
TCHHL1 is a hair-speci c protein given its high expression in scalp and chin skin [48]. TCHHL1 was identi ed in a genome-wide association study (GWAS) to have a signi cant association with hair shape within the top-associated single nucleotide polymorphisms (SNPs) (rs17646946), and showed nominally signi cant association with hair curliness [49]. ALOX15B is restricted to terminally differentiating keratinocytes (in particular the stratum granulosum) and 8(S)-lipoxygenase activity seems to be involved in terminal differentiation of mouse epidermis [50]. Clements et al. (2012) identi ed reduced expression of ALOX15B gene in ankyloblepharon-ectodermal defects-clefting (AEC) syndrome skin, with downregulated genes (KRT25 and KRT27) encoding keratins involved in the morphogenesis of hair follicles [51]. Thus, in combination with the current research, three genes may participate in the regulation of hair follicle density in Angora rabbits. The results of q-PCR of LNC_000797, LNC_013595, LNC_020367, KRTAP15-1, TCHHL1, and ALOX15B showed similar expression patterns between RNA-Seq and q-PCR, demonstrating the reliability of these data.

Conclusions
In conclusion, differences in the histology and lncRNA pro les of skin were identi ed in HWP and LWP rabbits. The histological analysis showed a higher hair follicle density in HWP rabbits. The analyses of lncRNA pro les identi ed candidate lncRNAs involved in lipid metabolism, apoptosis, and hair follicle development. Further studies are required to investigate the roles of candidate lncRNAs in hair follicle density to improve rabbit breeding programmes.

Animals
All rabbits were procured from the rabbit farm and acquired an approval from the farm owner in the Animal Husbandry and Veterinary Medicine Institute of Anhui Academy of Agriculture Sciences, Hefei, Anhui, China. 60 Wan Strain Angora rabbits (about one year old) were reared in the same conditions with regular pellets and water ad libtum. The wool weight of ve successive collections in one year from adult rabbits were determined. The sixty rabbits were divided into two populations designated as high wool production (HWP) and low wool production (LWP) according to wool production. The average wool weights showed remarkable difference (HWP: 401.3 ± 36.5 g vs LWP: 314.4 ± 29.2 g, P < 0.001). Finally, four rabbits with high and low wool production (430.1 ± 16.5 g vs 291.6 ± 13.3 g, P < 0.0001) were selected for the present study, respectively.

Sample collection, preparation for histological examination
The eight rabbits selected (four rabbits with high wool production, four rabbits with low wool production) were given anesthesia through an ear vein injection of 0.7% pentobarbital sodium (6 ml/kg) before sampling. Skin tissue samples (1 cm 2 ) from the backs, abdomens, sides and hips were collected at the fourth week after plucking for histological analysis. The skin samples were xed with 4% formaldehyde (40% formaldehyde solution and distilled water mixed at a 1:9 ratio) solution. Cross sections of the xed samples were washed with running water, dehydrated using an ethyl alcohol series, cleaned in xylene, and embedded in para n wax. The specimens were sectioned to a thickness of 4 μm using a Leica RM2235 microtome (Leica, Wetzlar, Germany). Cross-sections of the xed and para n-embedded samples were stained with HE, examined and photographed using an Olympus BX51 biomicroscope (Olympus Optical Company, Tokyo, Japan). The iodine solution was smeared on the resultant lesion to prevent bacterial infection.
After the experiment, the rabbits were retained in the rabbit farm and reared and protected from external stimuli.

cDNA library construction and sequencing
Under anesthesia, skin samples from the back of the eight rabbits selected (four rabbits with high wool production and four rabbits with low wool production; 430.1 ± 16.5 g vs 291.6 ± 13.3 g, P < 0.0001) were collected at the fourth week after plucking for RNA-seq. The skin samples were rstly frozen in liquid nitrogen immediatelly after cutting and then stored at -80℃ before RNA extration. The total RNA was extracted from the skin samples of HWP rabbits (designated as H1, H2, H3, and H4) and LWP rabbits (designated as L1, L2, L3, and L4) using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). The samples were sent to the Beijing Novogene Co., LTD. The quality of RNA (purity and integrity) was evaluated using the NanoPhotometer spectrophotometer (Implen, CA, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). The RNA was puri ed by removing rRNA, fragmented randomly, and converted to double Mapping, assembling and screening Impurity data were removed from the raw reads and more than 12 Gb clean reads per sample were generated. The clean reads with high quality were then aligned using HISAT2 to the rabbit reference genome (https://www.ncbi.nlm.nih.gov/genome/?term=Rabbit) sequence. The mapped reads of each sample were assembled by StringTie (v2.0.4) [52]. The candidate lncRNAs were distinguished according to its sequence charecteristics (length>200 nt and noncoding potential) and meantime transcripts predicted with coding potential were ltered out by multiple tools. Conservative and comparative analyses were done between lncRNAs and mRNAs, and classi cation of lncRNAs was also analyzed.

Quanti cation, target prediction and function analysis
The expression levels of the lncRNAs and mRNAs in each sample were calculated by fragments per kilobase of transcript per million fragments mapped (FPKM) using StringTie-eB software [14]. Differential expression between LWP and HWP groups was analyzed by using cuffdiff (https://www.genepattern.org/modules/docs/Cuffdiff/7), and the threshold was set as |log2 (Fold Change)| ≥ 1 and P value < 0.05. Target prediction was conducted by searching coding genes 100 kb up-and down-stream of lncRNAs (co-location) and analyzing co-expression relationship (pearson correlation) of mRNAs to lncRNAs. Then, GO and KEGG enrichment analyses were performed on targets and the function of key lncRNAs were predicted. GOseq R package [53] and KOBAS The relative expression level of each gene was estimated by the 2 -ΔΔCT method.
Student's t-test with two-sided was used in statistical comparisons in wool weight between HWP and LWP groups and RNA expression. Error bars represent the mean ± standard deviation (SD) as determined using GraphPad Prism 5 (GraphPad Sofware, Inc., Availability of data and materials The data was presented in the manuscript and the supporting materials.
Authors´ contributions HLZ, DWH, and HSD conceived the study. HLZ, DWH, and XFW performed sample collection and total RNA preparation. XWZ and YXQ performed the q-PCR validation. HSD and DWH conducted the data analysis and prepared gures and tables. HSD and DWH wrote the manuscript. All authors read and approved the nal manuscript.
Ethics approval and consent to participate The present study was carried out in strict accordance with relevant guidelines and regulations by the Ministry of Agriculture of the People's Republic of China. All experimental protocols were approved by the Ethics Committee of Anhui Academy of Agricultural Sciences.

Consent for publication
Not applicable.

Figure 1
The histological observation of skin tissue from different body sites of Wan strain Angora rabbits with HWP and LWP. a, b, c, d Transverse section of the backs, abdomens, sides, and hips of Wan strain Angora rabbits with HWP. e, f, g, h Transverse section of the backs, abdomens, sides, and hips of Wan strain Angora rabbits with LWP. HWP, High wool production; LWP, Low wool production.  Heatmaps of differentially expressed lncRNAs and mRNAs between HWP and LWP rabbits. a lncRNAs. b mRNAs. "L" and "H" represent low wool production and high wool production groups, respectively. "Red" and "blue" indicate up-regulated and downregulated transcripts, respectively. HWP, High wool production; LWP, Low wool production.

Figure 4
GO enrichment analysis of cis-regulated target genes. a GO analysis of lncRNA co-location mRNAs according to biological process b GO analysis of lncRNA co-expression mRNAs according to biological process. The signi cance of enrichment of each GO term was assessed by P-value < 0.05, and GO terms were subsequently ltered by the enrichment scores (-LgP).