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 five 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 cm2) from the backs, abdomens, sides and hips were collected at the fourth week after plucking for histological analysis. The skin samples were fixed with 4% formaldehyde (40% formaldehyde solution and distilled water mixed at a 1:9 ratio) solution. Cross sections of the fixed samples were washed with running water, dehydrated using an ethyl alcohol series, cleaned in xylene, and embedded in paraffin wax. The specimens were sectioned to a thickness of 4 μm using a Leica RM2235 microtome (Leica, Wetzlar, Germany). Cross-sections of the fixed and paraffin-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 firstly 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 purified by removing rRNA, fragmented randomly, and converted to double cDNA, then were ligated with NEBNext adaptors. Finally, the libraries were created by PCR using the NEBNext® Ultra™ Directional RNA Library Prep Kit ((NEB, USA) according to manufacturer’s recommendations. The libraries were quantified with Qubit2.0 and the insert size was determined using Agilent 2100. Eight libraries were sequenced on an Illumina HiSeq 2500 platform (Illumina, San Diego, CA, USA).
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 filtered out by multiple tools. Conservative and comparative analyses were done between lncRNAs and mRNAs, and classification of lncRNAs was also analyzed.
Quantification, 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 (http://www.genome.jp/kegg/) were used to conduct GO and KEGG enrichment analyses.
Quantitative real-time polymerase chain reaction
To verify the RNA-seq results, three candidate lncRNAs and mRNAs were selected from the list of DE lncRNAs and DEGs, and the relative expression level was determined by q-PCR on LightCycler 96 (Roche, Switzerland) using TransStart Green qPCR SuperMix (Transgen, Beijing, China). Before q-PCR, the first-strand cDNA was synthesised using EasyScript One-Step gDNA Removal and cDNA Synthesis SuperMix (Transgen, Beijing, China). The primers for q-PCR are listed in Additional file 3: Table S2. The amplification was carried out under the condictions as follow: 94 °C for 30 s, followed by 40 cycles of 94 °C for 5 s, 61 °C for 35 s, finally a melt curve analysis of 97 °C for 10 s, 65 °C for 1 min, and 97 °C for 1 s. The GAPDH gene was used as an internal control in these experiments. The q-PCR analysis was performed in triplicates for each sample (HWP group: H1, H2, H3, and H4; LWP group: L1, L2, L3, and L4). The relative expression level of each gene was estimated by the 2-ΔΔCT method.
Statistical analyses.
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., La Jolla, CA, USA). The results with a P value < 0.05 were considered to be statistically significant.