Transcriptome-Wide Comparative Analysis of microRNA Profiles in the Telogen Skins of Liaoning Cashmere Goats ( Capra hircus ) and Fine-Wool Sheep ( Ovis aries ) by Solexa Deep Sequencing

Abstract

Compare the microRNA (miRNA) trancriptomes of goat and sheep skin using Solexa sequencing to understand the development of skin and hair follicles (HFs). miRNA expression patterns vary in the two small RNA libraries from goat (G library) and sheep (S library) telogen skin samples. Analysis of the size distribution of 25.32 million clean reads revealed that most are 21–23 nucleotides. A total of 1910 known miRNAs and 2261 novel mature miRNAs were identified in this study. Among them, 107 novel miRNAs and 1246 known miRNAs were differentially expressed in the two libraries; 10 of the known miRNAs were identified using stem-loop quantitative real-time PCR. Furthermore, GO and KEGG pathway analysis of predicted miRNA targets illustrated the roles of these differentially expressed miRNAs in telogen HF development and growth. This study provides important information about the role of miRNAs in the regulation of HF development and their function in the telogen phase. This observation may help future investigations of the regulation of miRNAs during wool quality improvement.

Introduction

MiRNAs, are an extensive class of endogenous, small single-stranded and noncoding RNA molecules. miRNAs play an important role in diverse biological processes by regulating post-transcriptional gene expression via pairing to the 3′-untranslated regions (UTRs) of target mRNAs (Wenguang et al., 2007; Ghelani et al., 2012). miRNAs play key roles in a broad range of biological functions, such as development, differentiation, tissue morphogenesis, and even some diseases. However, the regulatory mechanism between miRNAs and hair follicle (HF) development has not been thoroughly examined. Solexa sequencing, a high-throughput sequencing technology, provide effective tools to study miRNAs. Thus, our research may provide a clue for further unraveling mechanism between miRNAs and HF development.

HF development result from extensive epithelial–mesenchymal interactions (Schmidt-Ullrich and Paus, 2005; Mikkola and Millar, 2006; Ohyama et al., 2010). Morphologically, the process can be divided into three distinct stages (induction, organogenesis, and cytodifferentiation) (Schmidt-Ullrich and Paus, 2005). The histological structures of mature HF consist of the hair bulb (Hardy, 1992), inner root sheath (IRS), hair shaft (Sperling, 1991), outer root sheath, and connective tissue sheath (Stenn and Paus, 2001). Thickening and elongation of epidermal keratinocytes (hair placodes) are known to be the HFs (Paus et al., 1999). The fine adjustment of the process relies on a series of signal factors and receptors. Emerging evidences indicates that proper small RNA processing is essential for normal HF development (Andl et al., 2006; Yi et al., 2008; Yi and Fuchs, 2010). For example, the expression of miR-31 markedly increases during anagen, whereas it decreases during catagen and telogen. Additionally, it can negatively regulate the expression of Fgf10 and Dlx3 transcription factors (Mardaryev et al., 2010). Bone morphogenetic proteins (BMPs) play essential roles in the control of skin development, postnatal tissue remodeling, and tumorigenesis. Recently, research has demonstrated that miR-21 prevents the inhibitory effects of BMP4 on cell proliferation and migration (Ahmed et al., 2011). In addition, miR-7-5p inhibits melanoma cell migration and invasion by inhibiting IRS-2 (a target of miR-7-5p) expression and oncogenic Akt signaling (Giles et al., 2013). Furthermore, the self-renewal property of the growth cycle of HFs, including anagen, telogen, and catagen, is known. Although the relationship between anagen HFs and small RNA (sRNA) has been reported widely (Liu et al., 2012), few studies have investigated telogen and catagen HFs. Therefore, in this study, areas of skin where the fur grows in telogen were sampled.

Cashmere goats and fine-wool sheep provide high quality wool, primary source of materials in the animal fiber industry (Wenguang et al., 2007). The hairy coat of Liaoning Cashmere goats is primarily composed of a double coat structure while fine-wool sheep has a single coat structure. Due to differences in physical properties of their hair, the market prices between them differ significantly. Along with the discovery of sRNAs, increasing attention has been paid to the regulatory roles of miRNAs in HF development. Therefore, Cashmere goats and fine-wool sheep were chosen as appropriate animal models for our research. To further clarify the function of miRNAs in skin telogen follicles, we identified and comparatively analyzed the miRNA profiles in Cashmere goat and fine-wool sheep using Solexa deep sequencing. Our findings provide a solid basis for improving the quality of wool and increasing the economic value, possibly even optimizing breeding.

Materials and Methods

Ethics statement

Liaoning cashmere goats (Capra hircus) and fine-wool sheep (Ovis aries) were obtained from a livestock farm located in the College of Animal Science and Technology, JiLin University, Changchun, China. Experiments were performed according to the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China; revised in June 2004) and approved by the Ethics Committee of Jilin University in accordance with the guidelines on animal care and use.

Skin collection and total RNA isolation

Three cashmere goats (one female and two male) and three fine-wool sheep (one female and two male) aged 1.5 years were used in this study. All the animals were fed under the same conditions and were sacrificed humanely to minimize suffering. All samples collection protocols were performed in accordance with the guide for the Care and Use of Laboratory Animals published by the U.S. National Institutes of Health. A skin sample was obtained from the scapular region of each animal in June 2011 during the telogen phase of HF growth and immediately stored in liquid nitrogen.

Total RNA was extracted using TRIzol reagent (TransGen Biotech) according to the manufacturer's protocol. Equal volumes of total RNA (2.5 μg from each skin sample) were detected by agarose (1%) gel-electrophoresis, and only two bands were clearly visible. Then, RNA concentration, purity and integrity were measured using the Agilent Technologies 2100 Bioanalyzer (RNA Nano Chip, Agilent). The ratios of OD260/OD280 were both above 1.8. We pooled the six samples into two groups (goat group and sheep group) with the same concentration and stored these at −80°C until further analysis.

Small RNA library construction and high-throughput sequencing

Equal quantities of total RNA isolated from the skin of three cashmere goats (hereafter referred to as the goat library or marked G) and three fine-wool sheep (hereafter referred to as the sheep library or marked S) were pooled, respectively, and prepared according to the Solexa EAS Small RNA Sample Prep Protocol. Solexa deep sequencing (Illumina Genome Analyzer) was performed at the Beijing Genomics Institute (BGI). Approximately, 10 μg of total RNA representing each sample set was provided.

In short, solexa sequencing was performed as follows: small RNA fragments (14–30 bases) of the total RNA pool were selected and enriched by 15% polyacrylamide/urea gel-electrophoresis and then ligated with proprietary adaptors to the 5′-monophosphate and 3′-hydroxyl termini of the RNA with T4 RNA ligase (Ambion). The 5′ ligation products (40–60 bases in length) were enriched and purified by 15% TBE-Urea gel and RNA fragments with adaptors at both ends (70–90 bases in length) were also enriched and purified by 10% TBE-Urea gel. The purified adaptor-ligated samples were then used as templates for cDNA synthesis. cDNA was amplified with 15 cycles using Illumina's small RNA primer sets. The purified cDNA was quantified on an Agilent Technologies 2100 Bioanalyzer, and sequencing libraries were constructed on the Illumina/Solexa sequencing platform with SBS-sequencing by synthesis method (Malone et al., 2012), which can decrease the loss of nucleotides caused by the secondary structure.

Analysis of sequencing data

Basic analysis

Raw sequence reads from HiSeq sequencing were first subjected to data cleaning according to the BGI small RNA reads pipeline as previously reported (Wei et al., 2009). This process includes removing the adaptor reads, several types of contaminants and all low-quality reads. The remaining clean reads in 18-30nt were obtained, and the unique sequences with their associated read counts were mapped onto the genome assembly using the SOAP program. The subsequent processes, including evaluation of sequencing quality and calculating length distribution of clean reads, were then performed. So far, the basic analysis process had been finished.

Advanced analysis

Standard bioinformatics analysis annotated the clean reads into different categories including miRNA, siRNA, piRNA, rRNA, tRNA, snRNA, snoRNA, repeat associate sRNA, exon, or intron RNA, by comparison with the known classes in small RNA databases (e.g., Rfam [ftp.sanger.ac.uk/pub/databases/Rfam], NCBI GenBank database) and by using software (e.g., Tag2repeat and Overlap developed by BGI). To avoid repeat annotations, we traversed these unique sequences in order (rRNAetc [in which GenBank>Rfam] >known miRNAs> repeat-associated small RNAs>exon>intron). To identify known miRNAs, unique sequences were aligned with all of animal precursor/mature miRNAs sequences from miRBase 18.0 (http://mirbase.org/). Then, the conservations of known miRNAs were compared with other species.

For those clean tags that were not annotated to any category, the novel miRNAs and their potential to form characteristic hairpin structures were predicted using MIREAP (https://sourceforge.net/projects/mireap/). The criteria of predicted miRNAs were as described previously (Meyers et al., 2008). Additionally, for comparing the differential expression of known miRNAs in two libraries, a traditional method was applied (Zhang et al., 2013b).

Furthermore, miRNA analysis was performed using the MIREAP software for target gene prediction and DAVID 6.7 tool (http://david.abcc.ncifcrf.gov/) for gene functional annotation. The whole process is shown in Supplementary Figure S1 (Supplementary Data are available online at www.liebertpub.com/dna).

miRNA validation through stem-loop quantitative real-time PCR

To verify the accuracy of the result for solexa sequencing, we randomly chose ten differential expressions of known miRNAs for stem-loop RT-qPCR as described previously (Kramer, 2011). According to these mature miRNA sequences in miRBase18.0, we designed the specific primers of the reverse transcription and quantitative PCR (showed in Supplementary Table S1). Primers were synthesized by SBS Biotechnology Co. The process of stem-loop RT-qPCR was described as follows: First, equal amounts of total RNA, 2 μg from each sample, were reverse-transcribed to cDNA using the EasyScript First-Strand cDNA Synthesis SuperMix Kit for RT-PCR (TransGen Biotech Co). The mixtures were incubated at 42°C for 30 min and 85°C for 5 min. Then, for detection and quantification, the MyiQ Real-time PCR detection system (Bio-Rad) and TransStart Top Green qPCR SuperMix (TransGen Biotech Co) were used. Each reaction (25 μL reaction mixtures) contained primers (1 μL), cDNA (1 μL), and the supplied enzyme mixture (12.5 μL). Next, 0.5 μL of the DNA double-strand-specific SYBR Green I dye was added for the detection of qPCR products. The PCR was performed using the following conditions: The mixtures were incubated at 95°C for 10 min, followed by 40 cycles of 95°C for 15 s and 60°C for 30 s in a 96-well optical plate. All samples were assayed in triplicate. The PCR products were verified by melting curve analysis and agarose gel electrophoresis. U6 snRNA was used as a reference gene for normalization. The difference in co-efficiencies was less than 0.1, indicating similar amplification efficiencies for all of the cDNA. The data were analyzed using the ΔΔCt method (Livak and Schmittgen, 2001), and the fold change in miRNAs expression were expressed as 2^−ΔΔCt.

Results

Overview of two libraries' small RNA composition

A total of 12,317,481 raw reads from three Cashmere goat skins (G) and 16,592,081 raw reads from three fine-wool sheep skins (S) were obtained by Solexa deep sequencing. All sequencing data in each small RNA library were classified in Table 1. The clean reads without adaptor, contaminants, and low-quality reads were accounting for 11,827,856 in G and 13,488,550 in S (with length 18-30nt). And 539,776 unique sequences were obtained in G and 1,214,068 unique sequences in S. Furthermore, the number of unique reads in the sheep library were over 2.2-fold (1214068/539776) higher than that in the goat library, while the average value of each unique read, representing the counts of clean total reads in the goat library, is much lower than that in the sheep library (11.11%-G/21.91%-S).

Table 1.

Small RNA Sequence Statistics for Liaoning Cashmere Goats (G) and Fine-Wool Sheep (S) Libraries

	Total reads	Percent of total reads (%)	Unique reads	Percent of unique reads (%)
G library
Raw reads	12,317,481	100.00
Low-quality reads	45,591	0.37
High-quality reads	12,271,890	99.63
High-quality reads	12,271,890	100.00
3′ adaptor sequence	31,868	0.26
5′ adaptor contaminants	39,821	0.32
Sequences smaller than 18nt^a	200,847	1.64
Insert sequence^b	171,174	1.39
Ploy A	324	2.64E-5
Clean reads ≥18nt	11,827,856	96.38	539,776
Clean reads ≥18nt	11,827,856	100.00	539,776	100.00
Mapping to sheep genome	7,396,049	62.53	129,582	24.01
Clean reads ≥18nt	11,827,856	100.00	539,776	100.00
Repeat-associated small RNAs	23,644	0.20	13,103	2.43
Noncoding RNA(rRNA etc.)^c	1,575,646	13.32	155,309	28.77
MiRNAs	8,876,920	75.05	51,197	9.48
Unannotated small RNAs^d	1,351,646	11.43	320,167	59.31
S library
Raw reads	16,592,081	100.00
Low-quality reads	64,413	0.39
High-quality reads	16,527,668	99.61
High-quality reads	16,527,668	100.00
3′ adaptor sequence	28,533	0.17
5′ adaptor contaminants	2,143,450	12.97
Sequences smaller than 18nt^a	372,423	2.25
Insert sequence^b	493,962	2.99
Ploy A	750	4.54E-5
Clean reads ≥18nt	13,488,550	81.61	1,214,068
Clean reads ≥18nt	13,488,550	100.00	1,214,068	100.00
Mapping to sheep genome	8,155,384	60.46	706,786	58.22
Clean reads ≥18nt	13,488,550	100.00	1,214,068	100.00
Repeat-associated small RNAs	811,883	6.02	249,032	20.51
Noncoding RNA(rRNA etc.)^c	1,610,379	11.94	192,130	15.83
MiRNAs	8,871,675	65.77	79,245	6.53
Unannotated small RNAs^d	2,194,613	16.27	693,650	57.13

Number of reads less than 18nt. Generally, small RNA tags are between 18-30nt long. So too short tags should be removed from data for further analysis.

Number of reads with no insertion.

rRNA.etc represents the total reads of four noncoding RNAs (rRNA, tRNA, snRNA, scRNA, snoRNA, and srpRNA).

Unannotated small RNAs represent the remaining small RNAs that failed to be mapped to any of them (rRNA.etc, known skin miRNAs and repeat-associated small RNAs).

miRNA, microRNA.

To further identify the disparity of library compositions, we compared the total reads and unique reads of two libraries. We found the common reads of two libraries occupied the most of total reads (82.31%) and the S library accounted for the major of unique reads (67.93%) (Supplementary Fig. S2). We next mapped these clean reads onto the genome assembly by the soap program (Li et al., 2008) to further acknowledge the possible mechanisms underlying the diversity of the small RNA. The 129,582 (G) and 706,786 (S) unique reads represented the total reads (Table 1). The total reads mapping to parts of the genome account for more than 60% of clean reads in two libraries. Analysis of the unique reads showed that the unique read content of the G library (24.01%) was lower than that of the S library (58.22%). In G library, the number of isoRNA species is more than that in S library, whereas the number of sRNA species is lower than that in S library. In addition, the recent distribution of sRNAs in sheep chromosomes, which includes 26 pairs of autosomes and the X chromosome, were displayed in the sheep genome database (Supplementary Fig. S3). By contrast, in the number of miRNA on each chromosome, the G library was more uneven than that of the S library. Furthermore, chromosomes 1, 15, 2, 3, and 5 had the extrude distribution state.

The small RNA length distribution and abundance were different between the two libraries (Fig. 1). In S library, the majority of the clean reads (60.55%) were 21–23 nucleotides long, whereas the G sRNAs was characterized by a bimodal distribution, including a weak peak at 19–21 nucleotides and a strong peak at 21–23 nucleotides.

FIG. 1.

The length and count distribution of sRNAs in two libraries.

Homology analysis

Because there are no goat miRNAs available in the miRBase 18.0 database, we analyzed the distribution of these known miRNAs in 85 species including Caenorhabditis elegans, Drosophila, and others. Then, we filtered 278 miRNA families that were commonly expressed in two libraries. After collating these data, we found over 100 miRNA families that were mapped to each of nine species, including cattle, horses, and pigs, which are close homology with goats and sheep (Supplementary Table S2). In addition, we compared 278 miRNA families to the known Ovis aries miRBase, and we found that 21 matched. To explore the conservation of these 21 miRNA families in other species, we further compared them to the nine species mentioned above (Supplementary Table S2). These results revealed that 7 miRNA families were conserved in all these species and among them, including miR-127, -136, -154, -299, -368, undef, and miR-368. Among them, the undef miRNAs referred the paired miRNAs, which were located on a different chromosome and miR-368 family was newly discovered.

Because of the highest homology between cattle and goats or sheep, we mapped the 21 miRNAs to Bos taurus genome and found that 19 families were divided into two clusters, except for the miR-368 family, which was not expressed in cattle (Fig. 2).

FIG. 2.

Twenty-one miRNA families cluster on the chr21 of bovine.

Known miRNAs and miRNA*s

A total of 1238 known mammalian miRNAs in G library and 1574 miRNAs in S library were identified according to the miRBase 18.0 (www.mirbase.org). Finally, a total of 1910 known miRNAs were obtained in two libraries. Furthermore, miRNA*s, a strand of miRNA duplex that generated from premicroRNA (premiRNA), was identified (Mah et al., 2010). One hundred seventy-eight miRNA*s, including 110 duplex-like miRNA:miRNA* pairs, were obtained in G library while 266 miRNA*s, including 164 pairs, were found in S library. Next, we chose 10 known miRNAs of the highest abundance from each library (74.92% of total reads-G and 52.28%-S) to contrast their counts (Fig. 3) and their most abundant sequences (Supplementary Table. S3). Among them, 7 miRNAs (miR-143, -203, and let-7, -7b, -7, -5p, -7f, and -7c) were found in both libraries. And these miRNAs, except for miR-203, expression levels were higher in G library than these in S library. In addition, miR-378, -378e and -378d were only detected in G library whereas miR-21, -101, and -101c were only expressed in S library.

FIG. 3.

The expressed pattern of 10 highest abundance known miRNAs from each library. ** means a statistically significant difference (**P < 0.01). miRNA, microRNA.

Small RNA classification and annotation

Excluding the known small RNAs, the remaining clean reads were further classified as mappable sequences (Material and Methods section). The number and proportions of each small RNA type in the two libraries were different (Table 2). The results showed that 75.8% of total reads in G library and 67.66% in S library were matched with protein-coding RNA; however, the percentage of unique reads was lower (18.48% in G and 13.79% in S). By contrast, a lower percentage of total reads (12.77% in G and 16.07% in S) mapped to nonprotein-coding RNA with a higher ratio of unique reads (22.21% in G and 29.10% in S). In addition, the repeat-associated small RNAs made up a larger proportion in S library (6.02% of total reads, representing 20.51% of unique reads) than in G library (0.2% of total reads, representing 2.43% of unique reads). The highest proportion of unique sequences (59.31% and 57.13%, respectively) were attributed to unannotated small RNAs in both libraries, while a smaller proportion of total reads were observed and might cover other classes of small RNAs such as piRNAs. In both libraries of miRNA, the smaller proportion of unique reads (9.48% G and 6.53% S) contained a larger proportion of total reads (>60%).

Table 2.

Reads Abundance of Various Classes of Small RNAs in G and S Library

Annotation	Total reads	Percent of total reads	Unique reads	Percent of unique reads
G library
Clean reads	11,827,856	100.00	539,776	100.00
Protein-coding RNA	8,965,519	75.80	99,751	18.48
Exon sRNA	71,108	0.60	40,319	7.47
Intron sRNA	17,491	0.15	8235	1.52
miRNA	8,876,920	75.05	51,197	9.48
Nonprotein-coding RNA	1,510,655	12.77	119,858	22.21
rRNA	1,257,982	10.64	83,125	15.40
repeat sRNA	23,644	0.20	13,103	2.43
scRNA	1331	0.01	127	0.02
snRNA	39,033	0.33	2578	0.48
snoRNA	44,432	0.38	2416	0.45
srpRNA	4588	0.04	680	0.13
tRNA	139,645	1.18	17,829	3.30
Unannotation	1,351,646	11.43	320,167	59.31
S library
Clean reads	13,488,550	100.00	1,214,068	100.00
Protein-coding RNA	9,126,975	67.66	167,110	13.76
Exon sRNA	87,312	0.64	28,440	2.34
Intron sRNA	167,988	1.25	59,425	4.90
miRNA	8,871,675	65.77	79,245	6.53
Nonprotein-coding RNA	2,166,962	16.07	353,297	29.10
rRNA	1,086,028	8.05	67,987	5.60
repeat sRNA	811,883	6.02	249,032	20.51
scRNA	3370	0.02	220	0.02
snRNA	41,393	0.31	6203	0.51
snoRNA	54,062	0.40	3608	0.30
srpRNA	5175	0.04	759	0.06
tRNA	165,051	1.22	25,488	2.10
Unannotation	2,194,613	16.27	693,650	57.13

Novel miRNA prediction

Novel candidate miRNAs in the unannotated database of our libraries were identified by MIREAP software based on their characteristic hairpin structure (Schwab et al., 2005; Allen et al., 2007). Sixty novel candidate miRNAs in G library and 2425 ones in S library were identified (Supplementary Table S4). These predicted premiRNA sequences of the novel candidate miRNAs were compared with the mature sequences in miRBase 18.0. Among 60 and 2425 potential novel miRNAs identified, 50 and 2217 novel mature miRNAs were predicted in the G and S library, respectively. Overall, the count of novel miRNA in G library was less than that in S library (Fig. 4B). In addition, the length of the novel miRNA sequences ranged from 18 to 26 nt.

FIG. 4.

Comparison of expression pattern of novel miRNAs and known miRNAs identified between goats and sheep. Differentially expressed novel miRNAs and known miRNAs in G library and S library. (A, C) Venn diagram displaying the distribution of unique novel (A) miRNAs and known miRNAs (C) in the G library (left, blue circle) and S library (right, violet circle). The overlapping region indicates co-expressed unique miRNAs; dashed circles indicate differentially expressed unique miRNAs in two libraries. (B, D) The scatter plot of differentially expressed novel miRNAs (B) and known miRNAs (D) (Control: X-axis; Treatment: Y-axis). The X axis and Y axis show expression level of miRNAs in the two species respectively. Red points represent miRNAs with ratio > 2; Blue points represent miRNA with 1/2 < ratio ≤ 2; Green points represent miRNAs with ratio ≤ 1/2. Ratio, normalized expression of the treatment/ normalized expression of the control. Color images available online at www.liebertpub.com/dna

We next compared the 50 novel mature miRNAs in G library with that 2217 in S library (Fig. 4 and Supplementary Table S5). As results, 11 and 96 differential expression novel miRNAs were specific distinguished in G library and S library respectively. Simultaneously, there were six differentially expressed novel miRNAs (miR-21, -26, -27, -41, -14, and -55) commonly enriched in both libraries (Fig. 4A and Supplementary Table S5). Among them, only the expression levels of miR-21, -14, and -55 in G library were lower than that in S library, especially miR-55.

Analysis and validation of differentially expressed miRNAs

The biogenesis of miRNAs must be highly regulated by various space-times. Thus, comparing differentially expressed miRNAs between libraries provides insight into molecular events related to the skin HFs of different species. We first normalized the reads and then used a p-value < 0.001 and an absolute value of log₂-ratio fold-change >1.0 as a threshold to judge the statistical significance of miRNA expression. Then, we filtered 1246 known miRNAs (out of 1910, 65.24%) that were significantly differentially expressed in the two libraries. Of those, 842 miRNAs (266 Sheep-specific, 576 co-expressed) were significantly downregulated, and 404 miRNAs (152 Goat-specific, 252 co-expressed) were significantly upregulated in G library when compared with the S library (Supplementary Table S6 and Fig. 4C). The number of the specific differentially expressed miRNAs in S library (21.35%) was much higher than that in G library (12.20%) (Fig. 4D). The result was similar to that of the differentially expressed novel miRNAs (Fig. 4A, B).

To validate the solexa sequencing results, we selected 10 differentially expressed miRNAs randomly from the two libraries for stem-loop quantitative real-time PCR (let-7a, miR-184-3p, miR-101a, miR-200a, miR-145, miR-17-5p, miR-451, miR-30b, miR-29b, and miR-125b-5p). The expression patterns of these 10 miRNAs are shown in Figure 5. Most of them showed similar expression patterns to those identified by Solexa sequencing analysis.

FIG. 5.

Validation of the sequencing results by qRT-PCR.

Known miRNAs and novel miRNAs target gene prediction and functional annotation analysis

To gain insight into the differences in regulatory mechanisms for known and novel miRNAs between the two libraries, MIREAP software was used to predict target genes of the miRNA. In G library, 2,195,185 miRNA-mRNA interaction sites for 1238 known miRNAs and 101,984 miRNA-mRNA interaction sites for 60 novel miRNAs were predicted. In S library, 2,880,071 miRNA-mRNA interaction sites for 1574 known miRNAs and 4,381,137 miRNA-mRNA interaction sites for 2425 novel miRNAs were predicted. These results illustrate that one target gene can be regulated by more than one novel miRNA and one miRNA can regulate several target genes simultaneously (such as shown in Supplementary Table S7). To further explore the differences of target genes between the two libraries, we carried out a GO enrichment and KEGG pathway analysis using the DAVID gene annotation tool. GO enrichment analysis for target genes based on the cellular component showed that 8706 target genes were involved in cellular components. Molecular function analysis showed that 8444 target genes were assigned different functions, specifically 85.2% of genes were related to copper ion binding functions. In addition, 8073 target genes were involved in biological processes. Not every Go-term had a direct effect on HF development, but they would likely take part in the regulation process (we filtered several major classes of Go-terms, shown in Supplementary Table S8).

KEGG pathway annotation showed that 8163 target genes were involved in 276 biological functions. (Supplementary Table S8 as examples). The most enriched pathway was the MAPK (2.52%) signaling, followed by the Wnt (1.83%), cytokine-cytokine receptor (1.70%), neuroactive ligand–receptor interaction (1.65%), and so on.

Discussion

miRNAs are involved in various cell processes including cell differentiation, proliferation, migration and apoptosis. There is growing evidence that miRNAs take part in the regulation of HF development. A large number of studies have proven that miRNAs play important roles in the regulation of HF development (Ning and Andl, 2013). However, few studies have focused on the difference of the miRNA profiles in the telogen skins of goat and sheep. Our study tried to explore the gene expression profiles during HF development in cashmere goats and fine-wool sheep.

In S library, the majority of the clean reads (60.55%) were 21–23 nucleotides long, which is within the typical size range for animal miRNAs (Blaszczyk et al., 2001). The abundance of 22nt small RNAs (32.18%), which is the typical size of Dicer-derived products (Xie et al., 2011), was more than the sum of 21nt (13.75%) and 23nt (14.62%) small RNAs. This result agrees with the recent deep sequencing in Ovis aries (Zhang et al., 2013a). However, the size of the goat sRNAs was characterized by a notable bimodal distribution, including a weak peak (∼19–21 nucleotides, accounting for 34.19%) and a strong peak (∼21–23 nucleotides, accounting for 55.08%). The 22-nt small RNAs were the most abundant (33.88%) and the next larger fractions were the 20-nt (20.29%), 21-nt (10.91%), and 23-nt (10.29%) small RNAs. This result is consistent with that previously reported for other Cashmere goats (Liu et al., 2012). Overall, the abundant 21–23-nt size classes may primarily represent miRNAs and were similar in both libraries.

In this study, a total of 1238 known mammalian miRNAs in G library and 1574 ones in S library were identified. These miRNAs in G library were more uneven than those in S library, especially located in chromosomes 1, 15, 2, 3, and 5, which means that the number of miRNAs on positive strands was not equal to that on negative strands in both libraries. We found seven miRNA families were conserved in 85 species. Among them, miR-127, miR-136, miR-154, miR-299, and miR-368, had been reported (Liu et al., 2012) and undef and miR-368 families were newly discovered. This outcome indicates that these seven miRNA families are necessary for the development of HFs in livestock. In the highest abundant 10 known miRNAs, 7 miRNAs (miR-143, -203, let-7, -7-5p, -7b, -7c, and -7f) expressed in both libraries and their total reads occupied more than a half. As expected, these miRNAs may play an important role in HF development. For example, miR-143 controls smooth muscle and endothelial cell function (Deng et al., 2015), miR-203 regulates keratinocytes differentiation (Deppe et al., 2016), and let-7 family regulates cell differentiation (Deppe et al., 2016). On the other hand, miR-378, −378e and −378d were only detected in G library, which promote angiogenesis (Joyce et al., 2011). And miR-21, −101, and 101c were only found in S library, which promote HF-derived neural crest stem cell differentiation (Ni et al., 2014). This may be the main cause of the diversity of compositions in these two small RNA libraries.

In addition, 60 novel candidate miRNAs in G library and 2425 ones in S library were identified. We predicted novel miRNAs via characteristic hairpin structure (Schwab et al., 2005; Allen et al., 2007). Fifty and 2217 novel mature miRNAs were predicted in G library and the S library. And the expression levels of novel mature miRNAs were much lower than conserved miRNAs. This result provides further evidence that nonconserved miRNAs are often expressed at a lower level than conserved miRNAs. The low abundance of novel miRNAs might provide a special role for HF development and growth during various stages. These differences in 50 novel miRNAs we analyzed in the two species might reflect the characteristics of miRNA expression in a spatial- or tissue-specific manner. Thus, whether these low abundance miRNAs are expressed at higher levels in other tissues and organs, or whether their regulation mechanism is related to environmental and others factors needs further investigation.

The differences of miRNAs target genes between the two libraries were explored via a GO term and KEGG pathway annotation. As expected, the Notch signaling pathway, which is a primary signaling pathway in telogen and is universally required for regulating HF differentiation and maintaining hair cycle homeostasis, was also included in our results (Ezratty et al., 2011; Lin et al., 2011; Aubin-Houzelstein, 2012) (Supplementary Fig. S4). In addition, several signaling pathways that play a role in controlling HF development were found by KEGG analysis, such as the TGF-beta signaling pathway (Jamora et al., 2004; Lin and Yang, 2013), Jak-STAT signaling pathway (Nishio et al., 2001; Schmitt et al., 2012), Wnt signaling pathway (Andl et al., 2002; Millar, 2002), VEGF signaling pathway (Li et al., 2012a, 2012b), MAPK signaling pathway (Dickinson et al., 2011), and ErbB signaling pathway (Wieduwilt and Moasser, 2008). In contrast, a number of enriched pathways did not play a major role in HF development and its related functions, such as neuroactive ligand–receptor interaction, cytokine–cytokine receptor interaction, and Axon guidance. In summary, the GO term and KEGG pathway annotations for the predicted miRNA targets further supported the significant roles of these novel miRNAs during HF development and growth.

Conclusion

Overall, this study provides important information about the role of miRNAs in the regulation of HF development and their function in the telogen phase. A total of 107 novel miRNAs and 1246 known miRNAs were differentially expressed in goat library and sheep library. Further identifying their targets will illustrate the roles of these differentially expressed miRNA in telogen HF development and growth. This observation may help future investigations of the regulation of miRNAs during wool quality improvement.

Footnotes

Acknowledgments

This research was supported by the Cashmere goat breeding center of Liaoning province. This work was supported by grants from the National Natural Science Foundation of China (NSFC) (Nos. 30800807 and 31072097), the Special Foundation for Postdoctor of China Ministry of Education (No. 20100471261) and Special Funds for Scientific Research on Public Causes (201303119).

Disclosure Statement

No competing financial interests exist.

References

Ahmed

M.I.

, Mardaryev

A.N.

, Lewis

C.J.

, Sharov

A.A.

, and Botchkareva

N.V.

(2011). MicroRNA-21 is an important downstream component of BMP signalling in epidermal keratinocytes. J Cell Sci, 124, 3399–3404.

Allen

, Xie

, Gustafson

A.I.

, and Carrington

J.C.

(2007). microRNA-directed phasing during trans-acting siRNA biogenesis in plants. Cell, 131, 74–88.

Andl

, Murchison

E.P.

, Liu

, Zhang

, Yunta-Gonzalez

, Tobias

J.W.

, et al. (2006). The miRNA-processing enzyme dicer is essential for the morphogenesis and maintenance of hair follicles. Curr Biol, 16, 1041–1049.

Andl

, Reddy

S.T.

, Gaddapara

, and Millar

S.E.

(2002). WNT signals are required for the initiation of hair follicle development. Dev Cell, 2, 643–653.

Aubin-Houzelstein

(2012). Notch signaling and the developing hair follicle. Adv Exp Med Biol, 727, 142–160.

Blaszczyk

, Tropea

J.E.

, Bubunenko

, Routzahn

K.M.

, Waugh

D.S.

, Court

D.L.

, et al. (2001). Crystallographic and modeling studies of RNase III suggest a mechanism for double-stranded RNA cleavage. Structure, 9, 1225–1236.

Deng

, Blanco

F.J.

, Stevens

, Lu

, Caudrillier

, McBride

, et al. (2015). MicroRNA-143 activation regulates smooth muscle and endothelial cell crosstalk in pulmonary arterial hypertension. Circ Res, 117, 870–883.

Deppe

, Steinritz

, Santovito

, Egea

, Schmidt

, Weber

, et al. (2016). Upregulation of miR-203 and miR-210 affect growth and differentiation of keratinocytes after exposure to sulfur mustard in normoxia and hypoxia. Toxicol Lett, 244, 81–87.

Dickinson

S.E.

, Olson

E.R.

, Zhang

, Cooper

S.J.

, Melton

, Criswell

P.J.

, et al. (2011). p38 MAP kinase plays a functional role in UVB-induced mouse skin carcinogenesis. Mol Carcinog, 50, 469–478.

10.

Ezratty

E.J.

, Stokes

, Chai

, Shah

A.S.

, Williams

S.E.

, and Fuchs

(2011). A role for the primary cilium in notch signaling and epidermal differentiation during skin development. Cell, 145, 1129–1141.

11.

Ghelani

H.S.

, Rachchh

M.A.

, and Gokani

R.H.

(2012). MicroRNAs as newer therapeutic targets: a big hope from a tiny player. J Pharmacol Pharmacother, 3, 217–227.

12.

Giles

K.M.

, Brown

RAM

, Epis

M.R.

, Kalinowski

F.C.

, and Leedman

P.J.

(2013). miRNA-7-5p inhibits melanoma cell migration and invasion. Biochem Biophys Res Commun, 430, 706–710.

13.

Hardy

M.H.

(1992). The secret life of the hair follicle. Trends Genet, 8, 55–61.

14.

Jamora

, Lee

, Kocieniewski

, Azhar

, Hosokawa

, Chai

, et al. (2004). A signaling pathway involving TGF-β2 and snail in hair follicle morphogenesis. Plos Biol, 3, e11.

15.

Joyce

C.E.

, Zhou

, Xia

, Ryan

, Thrash

, Menter

, et al. (2011). Deep sequencing of small RNAs from human skin reveals major alterations in the psoriasis miRNAome. Hum Mol Genet, 20, 4025–4040.

16.

Kramer

M.F.

(2011). Stem-loop RT-qPCR for miRNAs. Curr Protoc Mol Biol. Chapter 15, Unit 15.10. DOI: 10.1002/0471142727.mb1510s95.

17.

, Li

, Kristiansen

, and Wang

(2008). SOAP: short oligonucleotide alignment program. Bioinformatics, 24, 713–714.

18.

, Lu

Z.-F.

, Man

X.-Y.

, Li

C.-M.

, Zhou

, Chen

J.-Q.

, et al. (2012a). VEGF upregulates VEGF receptor-2 on human outer root sheath cells and stimulates proliferation through ERK pathway. Mol Biol Rep, 39, 8687–8694.

19.

, Man

X.-Y.

, Li

C.-M.

, Chen

J.-Q.

, Zhou

, Cai

S.-Q.

, et al. (2012b). VEGF induces proliferation of human hair follicle dermal papilla cells through VEGFR-2-mediated activation of ERK. Exp Cell Res, 318, 1633–1640.

20.

Lin

H.-Y.

, Kao

C.-H.

, Lin

K.M.-C.

, Kaartinen

, and Yang

L.-T.

(2011). Notch signaling regulates late-stage epidermal differentiation and maintains postnatal hair cycle homeostasis. PLoS One 6, e15842.

21.

Lin

H.-Y.

, and Yang

L.-T.

(2013). Differential response of epithelial stem cell populations in hair follicles to TGF-beta signaling. Dev Biol, 373, 394–406.

22.

Liu

, Xiao

, Li

, Zhao

, Lai

, Yu

, et al. (2012). Identification of conserved and novel microRNAs in cashmere goat skin by deep sequencing. PLoS One, 7, e50001.

23.

Livak

K.J.

, and Schmittgen

T.D.

(2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2(T)(-Delta Delta C) method. Methods, 25, 402–408.

24.

Mah

S.M.

, Buske

, Humphries

R.K.

, and Kuchenbauer

(2010). miRNA*: a passenger stranded in RNA-induced silencing complex?. Crit Rev Eukaryot Gene Expr, 20, 141–148.

25.

Malone

, Brennecke

, Czech

, Aravin

, and Hannon

G.J.

(2012). Preparation of small RNA libraries for high-throughput sequencing. Cold Spring Harbor Protoc, 2012, 1067–1077.

26.

Mardaryev

A.N.

, Ahmed

M.I.

, Vlahov

N.V.

, Fessing

M.Y.

, Gill

J.H.

, Sharov

A.A.

, et al. (2010). Micro-RNA-31 controls hair cycle-associated changes in gene expression programs of the skin and hair follicle. FASEB J, 24, 3869–3881.

27.

Meyers

B.C.

, Axtell

M.J.

, Bartel

D.P.

, Baulcombe

, Bowman

J.L.

, et al. (2008). Criteria for annotation of plant microRNAs. Plant Cell, 20, 3186–3190.

28.

Mikkola

M.L.

, and Millar

S.E.

(2006). The mammary bud as a skin appendage: unique and shared aspects of development. J Mammary Gland Biol Neoplasia, 11, 187–203.

29.

Millar

S.E.

(2002). Molecular mechanisms regulating hair follicle development. J Invest Dermatol, 118, 216–225.

30.

, Zhang

, Liu

, Yang

, Wang

, Fu

, et al. (2014). miR-21 promotes the differentiation of hair follicle-derived neural crest stem cells into Schwann cells. Neural Regen Res, 9, 828–836.

31.

Ning

M.S.

, and Andl

(2013). Control by a hair's breadth: the role of microRNAs in the skin. Cell Mol Life Sci, 70, 1149–1169.

32.

Nishio

, Matsui

, Tsuji

, Tamura

, and Suzuki

(2001). Immunolocalisation of the janus kinases (JAK) - signal transducers and activators of transcription (STAT) pathway in human epidermis. J Anat, 198, 581–589.

33.

Ohyama

, Zheng

, Paus

, and Stenn

K.S.

(2010). The mesenchymal component of hair follicle neogenesis: background, methods and molecular characterization. Exp Dermatol, 19, 89–99.

34.

Paus

, Muller-Rover

, van der Veen

, Maurer

, Eichmuller

, Ling

, et al. (1999). A comprehensive guide for the recognition and classification of distinct stages of hair follicle morphogenesis. J Invest Dermatol, 113, 523–532.

35.

Schmidt-Ullrich

, and Paus

(2005). Molecular principles of hair follicle induction and morphogenesis. Bioessays, 27, 247–261.

36.

Schmitt

M.J.

, Philippidou

, Reinsbach

S.E.

, Margue

, Wienecke-Baldacchino

, Nashan

, et al. (2012). Interferon-γ-induced activation of signal transducer and activator of transcription 1 (STAT1) up-regulates the tumor suppressing microRNA-29 family in melanoma cells. Cell Commun Signal 10, 41.

37.

Schwab

, Palatnik

J.F.

, Riester

, Schommer

, Schmid

, and Weigel

(2005). Specific effects of MicroRNAs on the plant transcriptome. Dev Cell, 8, 517–527.

38.

Sperling

L.C.

(1991). Hair anatomy for the clinician. J Am Acad Dermatol, 25, 1–17.

39.

Stenn

K.S.

, and Paus

(2001). Controls of hair follicle cycling. Physiol Rev, 81, 449–494.

40.

Wei

, Chen

, Yang

, Ma

, and Kang

(2009). Characterization and comparative profiling of the small RNA transcriptomes in two phases of locust. Genome Biol 10, R6.

41.

Wenguang

, Jianghong

, Jinquan

, and Yashizawa

(2007). A subset of skin-expressed microRNAs with possible roles in goat and sheep hair growth based on expression profiling of mammalian microRNAs. OMICS, 11, 385–396.

42.

Wieduwilt

M.J.

, and Moasser

M.M.

(2008). The epidermal growth factor receptor family: biology driving targeted therapeutics. Cell Mol Life Sci, 65, 1566–1584.

43.

Xie

S.-S.

, Li

X.-Y.

, Liu

, Cao

J.-H.

, Zhong

, and Zhao

S.-H.

(2011). Discovery of porcine microRNAs in multiple tissues by a solexa deep sequencing approach. PLoS One, 6, e16235.

44.

, and Fuchs

(2010). MicroRNA-mediated control in the skin. Cell Death Differ, 17, 229–235.

45.

, Poy

M.N.

, Stoffel

, and Fuchs

(2008). A skin microRNA promotes differentiation by repressing 'stemness'. Nature, 452, 225–229.

46.

Zhang

, Zhao

, Wei

, Sheng

, Ren

, Xu

, et al. (2013a). Identification and characterization of the miRNA transcriptome of Ovis aries. PLoS One, 8, e58905.

47.

Zhang

X.-D.

, Zhang

Y.-H.

, Ling

Y.-H.

, Liu

, Cao

H.-G.

, Yin

Z.-J.

, et al. (2013b). Characterization and differential expression of microRNAs in the ovaries of pregnant and non-pregnant goats (Capra hircus). BMC Genomics, 14, 157.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.19 MB

0.03 MB

0.02 MB

0.06 MB

0.02 MB

0.07 MB

0.02 MB

0.03 MB

0.08 MB

0.00 MB

0.23 MB

0.03 MB

0.02 MB