An Integrated Analysis of MicroRNA and mRNA Expression Profiles to Identify RNA Expression Signatures in Lambskin Hair Follicles in Hu Sheep

Wave patterns in lambskin hair follicles are an important factor determining the quality of sheep’s wool. Hair follicles in lambskin from Hu sheep, a breed unique to China, have 3 types of waves, designated as large, medium, and small. The quality of wool from small wave follicles is excellent, while the quality of large waves is considered poor. Because no molecular and biological studies on hair follicles of these sheep have been conducted to date, the molecular mechanisms underlying the formation of different wave patterns is currently unknown. The aim of this article was to screen the candidate microRNAs (miRNA) and genes for the development of hair follicles in Hu sheep. Two-day-old Hu lambs were selected from full-sib individuals that showed large, medium, and small waves. Integrated analysis of microRNA and mRNA expression profiles employed high-throughout sequencing technology. Approximately 13, 24, and 18 differentially expressed miRNAs were found between small and large waves, small and medium waves, and medium and large waves, respectively. A total of 54, 190, and 81 differentially expressed genes were found between small and large waves, small and medium waves, and medium and large waves, respectively, by RNA sequencing (RNA-seq) analysis. Differentially expressed genes were classified using gene ontology and pathway analyses. They were found to be mainly involved in cell differentiation, proliferation, apoptosis, growth, immune response, and ion transport, and were associated with MAPK and the Notch signaling pathway. Reverse transcription-polymerase chain reaction (RT-PCR) analyses of differentially-expressed miRNA and genes were consistent with sequencing results. Integrated analysis of miRNA and mRNA expression indicated that, compared to small waves, large waves included 4 downregulated miRNAs that had regulatory effects on 8 upregulated genes and 3 upregulated miRNAs, which in turn influenced 13 downregulated genes. Compared to small waves, medium waves included 13 downregulated miRNAs that had regulatory effects on 64 upregulated genes and 4 upregulated miRNAs, which in turn had regulatory effects on 22 downregulated genes. Compared to medium waves, large waves consisted of 13 upregulated miRNAs that had regulatory effects on 48 downregulated genes. These differentially expressed miRNAs and genes may play a significant role in forming different patterns, and provide evidence for the molecular mechanisms underlying the formation of hair follicles of varying patterns.


Introduction
Persian lamb skin is one of the "three pillars" of the international fur market. Its trade volume is 11,000,000 to 13,000,000 tons, accounting for 15%-20% of the world's fur market in 2007. The Karakul breed of sheep is well known throughout the world, particularly for its lambskin that brand name is "Bukhara", which is mostly black and gray, and represents about 50% of the world's lambskin production. To increase the variety of colors in lambskin, black lambskin from Karakul sheep is usually decolorized and dyed with other colors, but the process of decolor can significantly affect its quality. The cultivation of sheep with high-quality white lambskin has been performed for centuries.
Hu sheep are a breed with white lambskin that is unique to China, and regarded as a protected breed by the Chinese government. Lambskin from Hu sheep is world famous for its wavy pattern and is known as a "soft gem" [1]. The production of Hu sheep lambskin has increased due to its increased market demand. However, in recent years, attention has been focused on meat characteristics rather than the quality of lambskin during the selection process, resulting in a gradual decrease in lambskin quality over time, which in turn has caused significant economic losses. Therefore, identifying, developing, and protecting unique germplasm resources to provide high-quality genetic material for breeding is of great economic value.
The quality of lambskin is affected by various factors such as types, visibility, and distribution area of pattern. These indices are generally used to evaluate the quality of lambskin. Fineness, density, and curvature of the hair follicles, in turn, determine the type of wave pattern [2]. Therefore, wool quality is based on pattern formation, and hair follicle characteristics form the basis for hair growth [3].
The hair follicle is a complex accessory organ of the skin that has a unique morphology and structure. It controls skin growth and plays a major function in skin regeneration. Hair follicles consist of epithelial and hypodermal layers, which include the connective tissue sheath, inner root sheath, outer root sheath, hair bulb, and hair shaft [4]. Hair follicles are divided into primary follicles and secondary follicles. The primary follicles form coarse wool, and secondary follicles produce the fine wool known as fluff [5]. The relative density of primary and secondary follicles determines the fineness and curvature of wool, thus affecting the quality of wool. Hair follicles in most mammals has a periodic regularity of development, including a growth stage, a regression stage, and a telogen phase [6,7]. Various types of cell participate in the biological processes for the regulation of this cyclical development [8,9], which is influenced by several regulatory pathways and molecules. Signaling pathways such as Wnt [10], TGF-β [11,12], MAPK [13], Shh [14], and Notch [15] are extensively involved in different components of hair follicle morphogenesis and its periodic variation. In addition, some cytokines also play an important role in hair follicle development and hair growth processes, including KGF [16], IGF [17], VEGF [18], EGF [19], and thymosin β4 [20], which promote the growth of hair follicles, and FGF-5 [21,22], which inhibits hair growth. MicroRNAs (miRNAs) are endogenously expressed, highly conserved small RNAs with 21-25 nucleotides that are derived from one arm of their precursor with a stem-loop structure that is expressed in specific cells and tissues [23]. Landgraf et al. [24] reported that miRNAs are differentially expressed in abnormal physiology and developmental processes. Viswanathan et al. [25] observed that miRNA can regulate RNA-induced silencing complex (RISC) by mRNA shearing and inhibition. It regulates the expression of target genes that participate in physiological and pathological processes by interacting with the target mRNA's 3'-untranslated region. In recent years, studies have shown that miRNAs are also expressed in hair follicles of mammals. Yi et al. [26] constructed a small RNA library from the skin and hair follicles of E17.5 fetal rats and reported that it predominantly consisted of miRNAs. In addition, by knocking out the key gene Dicer I that controls miRNA maturation, they discovered that mice fail to develop hair follicles. Mardaryev et al. [27] discovered that 200 miRNAs were expressed differentially in different growth stages of hair follicles. Yuan et al. [28] identified 399 known miRNAs and 172 new miRNAs in Cashmere goat skin. Around 326 miRNAs were expressed at different stages of hair follicle development and 26, 41, and 55 miRNAs were expressed during the growth stage, regression stage, and telogen phase, respectively. Liu et al. [29] identified differentially expressed miRNAs at 3 development stages of Tibetan sheep and found that miRNAs regulated the growth of hair follicles by the regulation of target genes in the MAPK and Wnt pathway.
There is currently very little information on miRNA expression during hair follicle development, and no relevant miRNA study on hair follicles has been conducted; therefore, the molecular mechanisms underlying the formation of hair follicles with different wave patterns are not currently known. Transcriptome sequencing was used to study the molecular mechanism underlying the formation of these types of hair follicles in Hu sheep and to screen differentially expressed miRNAs and genes that were related to the process of proliferation and differentiation, growth, and apoptosis of hair follicle cells. The present study may establish a theoretical basis for the molecular mechanism underlying the formation of different wave patterns in the hair follicles of Hu sheep.

Overview of miRNA and mRNA sequencing results
We sequenced the miRNA and mRNA from the total RNA samples including large-wave wool, medium-wave wool, and small-wave wool using an Illumina Hiseq 2500 sequencer. From the miRNA libraries, the number of clean reads in large-wave wool, medium-wave wool, and small-wave wool were 69,657,452, 6,609,614, and 6,422,619, respectively. From the mRNA libraries, the number of clean reads in large-wave wool, medium-wave wool, and small-wave wool were 67,172,566, 71,812,980, and 69,133,870, respectively. Detailed results of sequencing and assembly are shown in Tables 1 and 2. We obtained 522 differentially expressed miRNAs by analyzing the expression level of each miRNA in the gene library using DESeq, and then significant analysis of microarrays was used to select the significantly and differentially expressed miRNAs in the 3 types of pattern. The results showed that there were 24 significantly and differentially expressed miRNAs in small and medium waves, and compared to small waves, medium waves contained 18 downregulated miRNAs and 5 upregulated miRNAs. There were 18 significantly and differentially expressed miRNA in medium waves and large waves, including 3 downregulated miRNAs and 14 upregulated miRNAs compared to medium waves. There were 13 significantly and differentially expressed miRNAs in small waves and large waves, including 7 downregulated miRNAs and 4 upregulated miRNAs compared to small waves. These results showed that there were 36 novel miRNAs and 6 known miRNAs that were differentially expressed in lambskin ( Table 3).3333A total of 751, 745, and 764 differentially expressed genes were detected between small and large waves, medium and large waves, and small and medium waves, respectively. Analysis of differentially expressed genes between any two patterns indicated that there were 481 downregulated genes and 270 upregulated genes in large waves compared to that observed in small waves, 380 downregulated genes and 365 upregulated genes in medium waves compared to that observed in small waves, 493 downregulated genes and 271 upregulated genes in large waves compared to that in medium waves. (S1 Table). In addition, we determined that MAPK8 (XM_004021553.1) and ROCK2 (XM_004007209.1) were associated with the Wnt signaling pathway. These findings indicate that differentially expressed genes might be involved in hair follicle growth and development.

Cluster analysis of differentially expressed miRNAs
MiRNA with similar expression patterns usually have similar biological functions, and cluster analysis is generally used to analyze the expression levels of differentially expressed miRNAs in different patterns. The differentially expressed miRNAs from three types of patterns were analyzed by hierarchical cluster analysis (S1-S3 Figs). Expression pattern of miRNA were classified into the same cluster in 3 types of patterns, but the expression quantity of miRNAs between any 2 patterns had major differences; therefore, we speculated that there were various regulatory mechanisms that were involved in different patterns.

MiRNA target prediction and integration of miRNA and mRNA expression profiles
We ran the miRanda [30] and TargetScan [31] prediction software to study the biological functions of 42 differentially expressed miRNAs. The results showed that there were 43, 81, and 190 target genes in small and large waves, medium and large waves, and small and medium waves, respectively. We further integrated the sequencing data into the predicted miRNA-mRNA pairs to validate the target pairs. A total of 20, 52, 45 target pairs were found to be inversely expressed.

Functional annotation and classification
Gene ontology (GO) and biological pathway enrichment were done for the 130 genes [32,33]. Gene ontology analysis showed that 4, 25, and 7 GO items were significantly enriched with these genes in small and large waves, medium and large waves, and small and medium waves, respectively (Tables 4, 5 and 6). These GO items were enriched with biological processes, cellular component and molecular function. In addition, these genes were mainly involved in cell differentiation, proliferation, apoptosis, growth, immune response, and ion transport. Pathway analysis identified 5, 16, and 12 signaling pathways in small and large waves, medium and large waves and small and medium waves, respectively (Tables 7, 8 and 9). In addition, the mTOR signaling pathway and the MAPK signaling pathway were connected with hair follicle development. At the same time, many cancer pathways were overrepresented within these genes, such as the FoxO signaling pathway. Some genes were also enriched in pathways that were related to hair follicle growth and development such as the Wnt signaling pathway, the Notch signaling pathway, and the TGF-beta signaling pathway. Although the enrichment in these 3 pathways was not significant, these genes remain important topics for future investigations. Differentially expressed miRNAs and mRNA identified by RT-PCR analysis To confirm the small RNA and mRNA sequencing results, 5 miRNAs and mRNA were selected for RT-PCR verification. Although the results of sequencing and PR-PCR analyses indicated different fold changes, the expression trends were the same. This finding indicates that the sequencing results were reliable (Tables 10 and 11).

Regulatory network analysis of miRNA-mRNA
In organisms, miRNAs play a role in post-transcriptional regulation by suppressing or silencing specific target genes. Therefore, we looked for correlations between upregulated miRNAs and downregulated mRNAs, as well as between downregulated miRNAs and upregulated mRNAs. The results demonstrate the changes in the number of differentially expressed miR-NAs from small RNA sequencing and differentially expressed genes from transcriptome sequencing. Compared to small waves, there were 7 differentially expressed miRNAs, including 4 downregulated miRNAs and 3 upregulated miRNAs, and 20 differentially expressed genes, including 12 downregulated genes and 8 upregulated genes in large waves (S2 Table). Compared to small waves, there were 19 differentially expressed miRNAs, including 15 downregulated miRNAs and 4 upregulated miRNAs, and 52 differentially expressed genes, including 7 downregulated genes and 45 upregulated genes in medium waves (S3 Table). Compared to medium waves, there were 12 upregulated miRNAs and 41 downregulated genes. (S4 Table). Among these miRNAs and genes, NW_004080181.1_3961 and MT4 were selected for RT-PCR validation of miRNA and mRNA co-expression (Fig 1). The results showed that 1 upregulated

Discussion
Hu sheep lambskin is one of the most well-known in the world and is prized for several specific attributes, including its wave pattern, white color, soft hair, clarity, and beauty. It contains 3 kinds of waves, namely large, medium, and small, which are based on the peak distances in the wave patterns. The quality of small waves is considered far superior to that of large waves. In the present study, we used high-throughput sequencing to screen miRNAs and genes that were differentially expressed in large-wave, medium-wave, and small-wave follicles. To adjust for differences between individuals and environmental factors, we evaluated full-sibling sheep. Gene ontology was used to divide genes into 3 categories, namely cellular component, molecular function, and biological progress. The differentially expressed genes observed in the present study were mainly involved in development, cell differentiation, proliferation, apoptosis, immune response, and ion transport. Furthermore, pathway analysis indicated that specific genes were involved in the Wnt, TGF-β, and MAPK signaling pathways. These highly expressed genes may also be involved in hair growth and hair follicle development.

Predictive analysis of miRNAs and target genes
Studies have shown that miRNA accounts for 2%-5% of mammalian genes, and it regulates 60% of protein-coding genes [34]; therefore, miRNAs are the most widely distributed regulatory factor in animals [35]. Previous studies have found that miRNAs regulate protein expression through translation inhibition, which in turn influences various developmental processes, including muscle development, hair follicle development, and mammary gland development [36]. In the body, miRNAs exert their biological functions by targeting specific genes. In the present study, we determined that miRNAs regulate multiple target genes, whereas a target gene may be regulated by multiple miRNAs. The relationship between miRNAs and target genes plays an important role in the regulation of hair follicle growth and pattern formation, thus establishing it as a very good method for studying miRNAs [37]. In the present study, we identified 522 differentially expressed miRNAs and through significance analysis, 36 unknown miRNAs and 6 known miRNAs were identified. Among these miRNAs, miR-10 was determined to regulate the transcription of the Hox gene, as well as various activities involved in protein synthesis, whereas miR-10 possibly plays an important role in cell development and cell differentiation [38]. MiR-let-7, which includes miR-let-7i, miR-let-7c, and miR-let-7b, is expressed in the skin of goats, sheep, and mice, and its target genes were related to the growth of hair follicles and hair quality. Members of the MiR-200 family were expressed in the epidermis and hair follicles, whereas those of the miR-199 family were only expressed in hair follicles. MiR-199a-3p, which was a member of the miR-199 family, reduces cell proliferation in cancer and suppressed the occurrence of cancer [39], thus, we speculated that miR-200a and miR-199a-3p possibly play a role in hair follicle growth and development. In addition, Trakooljul et al. found that most target genes of miR-143 were related to cell proliferation, apoptosis, and tumorigenesis [40].

The association analysis of differentially expressed miRNAs and genes
MiRNAs have a negative regulatory relationship with target genes, and can inhibit and degrade the expression of target genes. We have analyzed the expression of miRNAs and genes between pairs of patterns in Hu sheep, including a combination of both differentially expressed miRNA and genes, and we have examined the regulatory network of differentially expressed miRNA and genes to identify additional target genes.
To study the effect of differentially expressed genes on different wave patterns in the hair of Hu sheep, association analysis was adopted to examine the expression and regulatory mechanism of differentially expressed miRNAs and genes during pattern formation, including upregulated miRNAs and downregulated genes as well as downregulated miRNAs and upregulated genes. Most miRNA-regulated target genes and differentially expressed genes are novel miRNAs that regulate various genes such as those expressed in small waves and large waves. For example, NW_004080189.1_7961 has a regulatory relationship with differentially expressed genes such as NFATC3, KPNB1, NRN1, and MEF2A. Gene ontology analysis showed that these genes were mainly involved in cell differentiation, proliferation, apoptosis, growth, immune response, and ion transport. LMO4, the target gene of miR-143, regulates EMT that is caused by TGF-β. EMT participates in various physiological and pathological processes such as the regeneration of damaged tissue, tissue fibrosis, tumorigenesis, and tumor metastasis [41]. FGF is involved in the process of cell migration and cell development, and it is a molecule that significantly promotes angiogenesis. FGF1, which is the target gene of NW_004080189.1_7961, is a member of the fibroblast growth factor (FGF) family, thus, it could also significantly promote angiogenesis and cell division, as well as participate in the growth and development of various tissues and organs. FGF1 plays a very important role in cell proliferation, migration, and differentiation [42,43,44]. IGF2BP2, the target gene of miRlet-7i, regulates the translation of IGF-2 [45]. IGF-2 is a member of the insulin-like growth factor system (IGFs), and IGF-2 stimulates cell proliferation and affects the growth and differentiation of tissues and organs to promote tumor formation. TCF4 is the target gene of NW_004080182.1_4070 and NW_004080177.1_2324. TCF4 is often combined with β-catenin of the Wnt signaling pathway and is recognized as the most important molecule of the T-cell factor [46]. TCF4 restrains the Wnt signaling pathway because some isomers of TCF4 lack the binding site for β-catenin [47], and yet some studies have shown that TCF4 activates the Wnt pathway, and may be related to the source of the cells or tissues [48]. Some studies have shown that MT3 is the encoding gene of metallothionein in cells; it has a unique molecular structure and high affinity for metal ions, and not only has the physiological function of other MTs, but also inhibits cell physiological activities, as well as acting as a suppressor gene [49]. Although it also belongs to the MT family, MT4 was determined not to be involved in the inhibition of cell proliferation. It could be also used as an important candidate gene because it is differentially expressed in hair follicles. KEGG pathway analysis showed that target genes were enriched in the cancer pathway, as well as other pathways related to hair follicle development such as the Wnt, TGF-β, MAPK, and mTOR signaling pathways. Su et al. [50] reported that miR-200a directly or indirectly affects the expression of β-catenin, as well as participates in the Wnt signaling pathway. In addition, miR-199a-3p influences the activity of the mTOR signaling pathway. Differentially expressed genes MAPK3 and MAPK8 are involved in the MAPK signaling pathway, and these have specific regulatory functions. In addition, TCF4 is regulated by internal signals in the Wnt pathway.

Ethics statement
The Institutional Animal Care and Use Committee (IACUC) of the government of Jiangsu Province (Permit Number: 45) and the Ministry of Agriculture of China (Permit Number: 39) approved the animal study proposal. All experimental procedures were conducted in strict compliance with the recommendations of the Guide for the Care and Use of Laboratory Animals of Jiangsu Province and of the Animal Care and Use Committee of the Chinese Ministry of Agriculture. All efforts were made to minimize animal suffering.

Experimental populations
Three pairs of full-sib individuals were selected at birth from among Hu sheep reared at a Suzhou stud farm in China and the farm owners gave permission for use of these sheep. Each pair contained 1 individual with predominantly large-wave, 1 individual medium-wave wool, and 1 with predominantly small-wave wool. About 1 cm of hair root was cut off and placed in a microtube surrounded by Drikold, and the quantity of hair root was collected 1/3 of the volume of the microtube, which was then collected into the freezing tube and stored in Drikold.

Total RNA extraction and cDNA library preparation
Total RNA was isolated using the Recover All Total Nucleic Acid Isolation Kit (Life Technologies, Carlsbad, CA, USA) for miRNA and mRNA sequencing, according to the manufacturer's instructions. Integrity of RNA was checked on an Agilent 2100 bioanalyzer. The sequencing library was prepared according to the standard protocol. Briefly, for miRNA sequencing, libraries were prepared by ligating different adaptors to the total RNA followed by reverse transcription and PCR amplification (RT-PCR). MiRNA libraries were sequenced on the Illumina Hiseq2500 system with 50-base-pair (bp) single-end reads, according to the manufacturer's standard protocol. All small sequencing raw data were deposited in the Sequence Read Archive database (SRR3099013, SRR3099014, and SRR3099015). For mRNA sequencing, total RNA was firstly poly-A-selected followed by fragmentation of RNA into small pieces. The cleaved RNA fragments were reverse transcribed to cDNA, end-repaired, and ligated with Illumina adapters using a Quick ligation TM kit (NEB) and DNA ligase. The libraries were then fractionated on agarose gel; short fragments (approximately 200 bp) were excised and amplified by PCR. The cDNA library was sequenced on the Illumina sequencing platform, and raw reads were generated using the Solexa pipeline according to the manufacturer's recommendations. All mRNA sequencing raw data were deposited in the Sequence Read Archive database (SRR3099016, SRR3099017, and SRR3099018).

Reads processing
For small sequencing reads, all low quality reads such as those with poly A/T and adapter contaminants were filtered. Sequences longer than 40 bp or shorter than 15 bp were removed. The high-quality clean reads were mapped to the Ovis aries genome using Bowtie software. Small RNA tags were aligned to the miRNA precursor/mature miRNA in miRBase. Furthermore, Rfam, RepBase and Genbank data were used to identify small RNA tags originating from miRNA, rRNA, tRNA, snRNA, and snoRNA.
For mRNA sequencing reads, adaptor sequences and low quality reads (threshold quality, 20; threshold length, 35 bp) were filtered. The matched reads were aligned to Ovis aries. The mRNA expression level was measured by FPKM (fragments per kilobase per million reads).

Real-time PCR verification
Differentially expressed miRNAs and genes were randomly selected, and U6 and GAPDH were used as reference genes. The SYBR Green I method was used for quantitative testing to verify the reliability of the sequencing results. Relevant information about genes that were assessed via RT-PCR analysis and the primers used in the assay are shown in S5 Table. The standard curve was established using a cDNA gradient dilution, and each sample was tested 3 times in a 7500 PCR instrument. The relative expression of the target gene was calculated according to the following formula [51]: Relative expression = 2 -ΔΔ Ct, ΔΔC t = (C, target gene-C, housekeeping gene) large waves −(C t, target gene-C, housekeeping gene) small waves . SPSS17.0 was used to compute for the Ct mean and standard deviation between replicate samples, and one-way ANOVA was adopted for analysis of significance.