Next Article in Journal
The Psychoemotional Stress-Induced Changes in the Abundance of SatIII (1q12) and Telomere Repeats, but Not Ribosomal DNA, in Human Leukocytes
Next Article in Special Issue
Novel Single Nucleotide Polymorphisms and Haplotype of MYF5 Gene Are Associated with Body Measurements and Ultrasound Traits in Grassland Short-Tailed Sheep
Previous Article in Journal
The First Russian Patient with Native American Myopathy
Previous Article in Special Issue
BMP7 Functions to Regulate Proliferation of Dermal Papilla Cells in Hu Sheep
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Integrated Hair Follicle Profiles of microRNAs and mRNAs to Reveal the Pattern Formation of Hu Sheep Lambskin

1
Joint International Research Laboratory of Agriculture and Agri-Product Safety of Ministry of Education of China, Yangzhou University, Yangzhou 225009, China
2
College of Animal Science and Technology, Yangzhou University, Yangzhou 225009, China
3
International Centre for Agricultural Research in the Dry Areas, Addis Ababa 999047, Ethiopia
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Genes 2022, 13(2), 342; https://doi.org/10.3390/genes13020342
Submission received: 14 January 2022 / Revised: 5 February 2022 / Accepted: 10 February 2022 / Published: 14 February 2022

Abstract

:
Hair follicle development is closely associated with wool curvature. Current studies reveal the crucial role of microRNAs (miRNAs) in hair follicle growth and development. However, few studies are known regarding their role in wool curvature. To reveal the potential roles of miRNAs in Hu sheep lambskin with different patterns, a total of 37 differentially expressed (DE) miRNAs were identified in hair follicles between small waves (SM) and straight wool (ST) groups using RNA-seq. Through functional enrichment and miRNA-mRNA co-expression analysis, some key miRNAs (oar-miR-143, oar-miR-200b, oar-miR-10a, oar-miR-181a, oar-miR-10b, oar-miR-125b, etc.) and miRNA-mRNA pairs (miR-125b target CD34, miR-181a target FGF12, LMO3, miR-200b target ZNF536, etc.) were identified. Though direct or indirect ways affecting hair follicle development, these miRNAs and mRNAs may have possible effects on wool curvature, and this study thus provides valuable insight on potential pattern formation.

1. Introduction

Hu sheep is a unique species of lambskin in China, characterized by precocity and high productivity. It is well known for its gentle lambskin and beautiful water-wave pattern. Hu lambskin thus holds economic value because of its popularity among certain consumers and its quality is based on the visible wave pattern. The wave pattern determines the quality of lambskin, and the quality of small waves is excellent, while straight wool is poor. The type of wave pattern is determined by many factors, such as curvature, density, and fineness of the wool. Wool curvature is the key factor of the wave pattern of lambskin, which is decided by hair follicle growth and development [1]. Hair follicle growth and development is a complex process, which participated in delicate molecular regulatory mechanisms. Although its genetic basis in French breeds has recently been unraveled as an insertion/deletion polymorphism (InDel) in the 3′UTR part of the IRF2BP2 gene [2], it remains unclear whether the same mechanism exists for Hu sheep. If so, a Hu sheep population of high-quality lambskin might be bred via molecular breeding, thus allowing for the conservation of the Hu sheep lambskin trait and increasing the production of high-quality lambskin.
Recently, many reports have indicated that multiple genes may be involved in hair follicle growth and development, such as Wnt10a [3], lef1 [4], Sox9 [5], BMP4 [6]. Additionally, some signaling pathways, such as BMP [7], Eda [8], Shh [9], and TGF-β [10] signaling, were revealed to promote or suppress the process of hair follicle growth and development. With the advancement of genomics, environment–gene interactions have become a research hotspot, and epistatic regulation may represent a bridge between the environment and genetics. Studies have shown that rearing cashmere goats under controlled lighting conditions may enhance cashmere production during non-cashmere seasons. It was also suggested that these goats may have acquired epigenetic memory, as they continued with this production pattern in the following year (after experiments have ceased). Wu et al. noted a significant increase in melatonin-related Hoxc13 expression which also positively triggered hair follicle development-relevant genes. It is, however, known that melatonin changes seasonally and environmental factors may thus influence genetic factors that are involved in hair follicles’ development [11].
MicroRNAs (miRNAs) are small endogenous RNAs (consisting of approximately 22 nucleotides) that are considered to play important regulatory roles in biological processes by targeting mRNAs or repressing translation [12]. In addition to being associated with many diseases, studies have linked miRNAs to the regulation of cell development, proliferation, and apoptosis, especially associated with many diseases [13]. With the further research of miRNA, the hair follicle growth and development process may also be controlled by epigenetic regulator factor miRNAs, except for multiple genes and signaling pathways. MiRNA and its target gene form diverse regulatory networks, playing an important role in hair follicle growth and development. Current research has found that miRNAs are expressed in a variety of different hair follicle cells, such as hair follicle stem cells, matrix cells, outer root sheath cells, inner root sheath cells, and lots of miRNAs are verified to be specifically expressed in hair follicle cells [14]. For example, both miR-125b [15] and miR-205 [16] are expressed in hair follicle stem cells, but they have the opposite effect. MiR-125b may act as a repressor to suppress stem cell differentiation, while miR-205 has a positive effect on stem cells and the proliferation of their progenies. Similarly, miR-214 [17] and miR-31 [18] are expressed in the outer root sheath and hair follicle matrix cells where they induce different hair follicle functions. MiR-24 is richly expressed in the inner root sheath and competitively binds to TCF-3 to promote hair follicle keratinocyte differentiation. However, it is unclear in which cells miR-137 is richly expressed, but miR-137 is associated with coat color.
MiRNA is normally about 20–25 nucleotides with a seed sequence of 6–8 nucleotides. Various studies have shown that miRNAs could competitively bind to 3′UTR of mRNAs to affect the expression of the target genes [19]. Of course, the targeting relationship between miRNA and mRNA is not one by one. Years of research have indicated that the targeting relationship between miRNA and mRNA is not one-on-one, but that miRNA may regulate the expression of multiple or even hundreds of genes [20]. For example, a previous study on the changes in hair follicle morphogenesis indicated that the latter had been related to the expression of miR-125b and around 800 transcripts with which miR-125b may be targeted [15]. Similarly, miR-205 has been identified to competitively bind to 3′UTR of INPP1b, INPP4b, PHLDA2, and FRK, whereas miR-31 also has four validated target genes, KRT16, KRT17, FGF10, and DLX3 [14]. It has also been suggested that, under certain circumstances, the biological function of miRNAs may be predicted through their target genes.
The keratin (KRT) family is an important components of hair follicle. More specifically, KRT27, KRT35, and KRT85 were found to be involved in the construction of the initially differentiated cortical structure, which ultimately leads to fibrillar bending [21], while, KRT17 was involved in the regulation of the hair follicle cycle and acted as a positive regulator of hair follicle growth [22]. Since miRNAs can target genes to regulate traits, we hypothesize that when either miRNAs or target genes are associated with hair follicle development, it is highly likely that the corresponding miRNAs or target genes might also be involved in the same process. For example, since KRT17 is involved in the hair follicle cycle, then miR-31 may be a candidate miRNA to investigate its role in hair follicle development.
Much work has been performed on the hair follicles of Hu sheep and we previously reported a total of 114 DE circular RNA (cirRNAs) [23] and 72 DE long non-coding RNAs (lncRNAs) [24] between small waves and straight wool after completing whole transcriptome sequencing. In the same study, 137 DE mRNAs between large waves and small waves, and BMP7, MMP2, SNAI1, CDKNIC, MT3, and POU1F1 were selected to act as candidate genes [25]. Furthermore, to distinguish between small, medium, and large waves, miRNA sequencing was performed and DE-miRNAs (miR-143, miR-10a, let-7a, etc.) were screened association analysis with hair follicle metrics (diameters and numbers) [26]. Since the consideration of small differences between the three pattern types, the DE miRNAs were not accurate enough, so this experiment further mined DE miRNAs from straight wool and small waves.
In summary, miRNA plays an extremely important role in hair follicle development. By thoroughly screening the DE miRNA and its target genes between small waves and straight wool of hair follicle in Hu sheep lambskin. It may provide the basis for further explorations. Specifically, understanding miRNA regulation in the formation of different patterns in Hu sheep lambskin might allow for the breeding of new Hu sheep strains of with high-quality lambskin.

2. Materials and Methods

2.1. Sample Collection

Animal research proposals were approved by the Animal Care and Use Committee at Yangzhou University. All animal experimental procedures were implemented in strict obedience of the guidelines for the Animal Care and Use Committee at Yangzhou University to minimize the suffering of animals.
A total of 6 two-day-old healthy Hu lambs were selected from the Suzhou Stud Farm in Jiangsu Province, China. Three pairs of full-sib individuals were divided into small waves groups (SM1, SM2, SM3) and straight wool groups (ST1, ST2, ST3) with three replicates (Figure 1). The phenotype of Lambskin was classified as small waves and straight wool according to the width of the pattern that is the length between the bulges at the ends of a single pattern. A pattern with a width of less than 1.25 cm was deemed a “small” pattern, whereas wool that barely curved was deemed as “straight wool”. Approximately 1 cm of hair root was cut off from the dorsal side of Hu lambs. A total of 1/3 volume of the tube was collected and stored in Drikold.

2.2. Library Preparation for Sequencing

Total RNA was extracted from the hair roots of small waves and straight wool groups using RNAiso Plus Kit (Takara, Kusatsu, Shiga, Japan) according to the manufacturer’s instructions. After total RNA isolation, the quality and concentration of RNA were measured to meet the sequencing requirements. The mRNA libraries were generated with 3 μg high-quality total RNA per sample through NEBNext® UltraTM RNA Library Prep Kit Illumina® (NEB, Ipswich, MA, USA) following the manufacturer’s recommendations. Firstly, mRNA was purified with ploy-T oligo-attached magnetic beads. First-strand cDNA synthesized by M-MuLV Reverse Transcriptase and random hexamer primer. Then, RNase H and polymerase I were used to synthesize the second cDNA. The cDNA fragments of 150–200 bp were selected by AMPure XP system (Beckman Coulter, Beverly, MA, USA) and ligated to NEBNext Adaptor. The quality of PCR product was assessed on the Agilent Bioanalyzer 2100 system. The library was sequenced (150 bp paired-end) using Illumina Hiseq 2500 platform.
Simultaneously, a total of 3 μg total RNA per sample was applied to construct the small RNA library through NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB, Ipswich, MA, USA) following the manufacturer’s recommendations. The 3′ ends and 5′ ends adaptor were directly and specifically ligated to the 3′ end and 5′ ends of miRNA, siRNA, and piRNA. Then, the first-strand of cDNA was synthesized using M-MuLV Reverse Transcriptase, and PCR amplification was implemented. DNA fragments corresponding to 140–160 bp were recovered and dissolved in the elution buffer. After the cDNA library was assessed on the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips, the libraries were sequenced on an Illumina Hiseq 2500 platform (50 bp single-end).

2.3. Quality Control, Reference Genome Alignment, and Assembly

Raw reads were obtained from sequencing and the quality control was carried out through custom perl and python scripts. Clean reads were purified by removing low-quality reads and reads containing ploy-N, polyA/T/G/C, with 5′ adapter contaminants, without 3′ adapter or the insert tag. Subsequently, error rates, Q20, Q30, and GC content were calculated [1,27]. Then, clean reads were used for all the downstream analyses. The reference genome (Ovis aries v4.0) and gene annotation files were downloaded from the genome website (http://www.ensembl.org/index.html, accessed on 9 May 2021). Hisat2 (http://ccb.jhu.edu/software/hisat2, accessed on 9 May 2021) and Stringtie (http://ccb.jhu.edu/software/stringtie/, accessed on 9 May 2021) were used to map mRNA to reference sequence and to splice the transcript.
For miRNAs, the reads mapped to reference sequence were compared with the sequences in the specified range in miRBase [28]. To remove tags formed by rRNA, tRNA, snRNA, snoRNA, repeat sequences, and mRNA, small RNA tags were mapped to Rfam database and RepeatMasker. At the same time, miREvo [29] and mirdeep2 software [30] were integrated to predict novel miRNA through exploring the minimum free energy, Dicer cleavage site, and secondary structure of the unannotated small RNA tags, the miRNA first nucleotide bias with different lengths and miRNA nucleotide bias at each position, respectively. According to the priority rule: known miRNA > rRNA > tRNA > snRNA > snoRNA > repeat > novel miRNA; each small RNA was mapped to only one annotation.

2.4. Quantification and Differential Expression of miRNA and mRNA

In this experiment, miRNA and mRNA expression levels were estimated through transcript per million (TPM) and fragments per kilobase per million mapped reads (FPKM) [31], respectively. To ensure the accuracy of the subsequent screening of DE miRNAs, we firstly carried out Pearson correlation analysis of samples.
Due to small waves and straight wool groups with three biological replicates, DESeq2 R package (1.8.1) was used for differential expression analysis of these two groups. Flog change and corrected p-values that were adjusted using the Benjamini and Hochberg approach were used to assess the significant difference of miRNAs and mRNAs expression levels in different groups. Corrected p-values < 0.05 were considered to be significantly differentially expressed [32]. Additionally, volcano plots and heat maps were adopted to screen differentially expressed miRNAs and mRNAs.

2.5. Target Gene Prediction and Functional Enrichment Analysis of miRNAs

Predicting the target gene of known and novel miRNA was performed by miRanda and RNAhybrid software for straight wool and small waves groups [33]. The correspondence relationship between miRNA and target gene was obtained through the intersection of predicted results.
According to the results of predicting the target gene of DE miRNA, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were used on the target gene candidates of differentially expressed miRNAs. GO enrichment analysis was implemented using the GOseq method based on Wallenius non-central hyper-geometric distribution [34]. KEGG database resource is a tool for understanding high-level functions and utilities of the biological system. The statistical enrichment of the target gene candidates was estimated by KOBAS software [35]. The clusterProfiler R package was applied to analyze GO and KEGG pathway enrichment of the target gene candidates. It (p < 0.05) was considered significantly enriched.

2.6. Integrated Analysis of miRNA and mRNA Expression Profiles

To further analyze the relationship between DE miRNAs and DE mRNAs related to hair follicle development, a network of miRNA-mRNA needed to be constructed based on their expression levels and functions through Cytoscape_2.8.3 software.

2.7. Verification of RNA-seq Results

To verify the reliability of the RNA-seq results, 6 DE miRNAs and 6 DE mRNAs were randomly selected with U6 and GAPDH as reference genes, respectively. The detailed information on relevant primers is shown in Supplementary Table S1. Each sample with three replicates was tested to establish a standard curve. The Real-time Quantitative Polymerase Chain Reaction (RT-qPCR) amplification was performed under the following procedure: initial denaturation at 95 °C for 5 min, circular reaction at 95 °C for 10 s, and 60 °C for 30 s with 40 cycles. The relative expression levels of DE miRNAs and mRNAs were calculated by 2−ΔΔCt method. A dependent sample t-test was performed to calculate relative expression and SPSS13.0 was used to analyze the significant difference. All the results were shown as mean ± SD.

3. Results

3.1. The Expression Profile of miRNA and mRNA

The expression profiles of miRNAs and mRNAs were investigated in hair follicles between straight wool and small waves using RNA-seq. Regarding miRNA, a total of approximately 32,724,357 and 38,996,977 raw reads were generated in SM and ST, respectively, and 31,173,567 and 37,529,529 clean reads were filtered for subsequent analysis by RNA-seq (Supplementary Table S2). The length interval of sRNA was 18 to 35 nt by length distribution statistics, and miRNAs were concentrated in 21 to 22 nt (Figure 2A). Additionally, miRNAs were distributed on chromosome 1 to X (Supplementary Figure S1), and most miRNAs were located on the chromosome 1, 2, 3. Moreover, known miRNA, novel miRNA, tRNA, rRNA, etc. were identified by alignment and annotation, but the known miRNA and novel miRNA account for only a small fraction of non-coding RNA (Figure 2B,C). Base preference analysis of miRNAs with different lengths showed that first base distribution of the miRNAs with 20 to 24 nt were mainly U (Figure 2D–G), which were consistent with the characteristics of miRNA.
According to the results of Pearson correlation analysis between samples, the correlation coefficient between biological replicates was higher in the small waves group, and the expression pattern was similar, as well as in the straight wool group (Supplementary Figure S2). According to the TPM of miRNAs estimating expression levels, 595 miRNAs were identified in hair follicles between SM vs. ST, containing 37 DE miRNAs (18 known miRNAs and 19 novel miRNAs) consisting of 18 up-regulated and 19 down-regulated miRNAs (Figure 3A and Supplementary Table S3). Furthermore, the expression profile of DE miRNAs was obviously presented to different expression patterns in Figure 3B.
Regarding mRNA, a total of approximately 330,913,612 and 296,509,262 raw reads were generated in SM and ST, respectively, and 325,445,988 and 290,032,082 clean reads were filtered for subsequent analysis (Supplementary Table S4). According to the FPKM of mRNAs estimating gene expression levels, 17,786 mRNAs were identified in hair follicle, suggesting that these data were reasonable. The mRNA difference significance analysis showed that a total of 489 DE mRNAs were identified between SM vs. ST (Figure 4A and Supplementary Table S5), of which 222 were up-regulated and 267 down-regulated DE mRNAs between SM vs. ST. Additionally, the expression profile of DE mRNAs was obviously presented to different expression patterns in Figure 4B.

3.2. Functional Enrichment Analysis of Target Genes

To better understand the biological functions of DE miRNAs, the Gene Ontology (GO) and KEGG pathway analysis were performed on their target genes using the target gene prediction method. The results showed that 37 DE miRNAs corresponded to 533 target genes, and GO enrichment analysis of target genes, a total of 249 GO terms were significantly enriched, containing biological process (BP), cellular component (CC), and molecular function (MF) (Supplementary Table S6). From Figure 4, the top enriched GO term is intracellular organelle (GO:0043229), cellular component organization or biogenesis (GO:0071840), and GTP binding (GO:0005525) in CC, BP, and MF, respectively. Moreover, some GO terms, such as epithelial cell differentiation (GO:0030855), regulation of epithelial cell differentiation (GO:0030856), and epidermal growth factor catabolic process (GO:0007174) in BP, suggest the crucial role in epithelial cell proliferation and differentiation, to a certain extent, an effect on hair growth. Some GO terms related to the preparation of transcription include DNA-directed RNA polymerase II, holoenzyme (GO:0016591) and nuclear DNA-directed RNA polymerase complex (GO:0055029) in CC. Some GO terms involved in energy supply and enzyme synthesis include nucleoside triphosphate adenylate kinase activity (GO:0046899) and tRNA-specific ribonuclease activity (GO:0004549) in MF (Figure 5A–C and Supplementary Table S6).
KEGG pathway analysis revealed that these target genes were enriched in 52 KEGG terms, with the most enriched pathway being that of metabolic pathways (oas01100) (Figure 5D and Supplementary Table S7). Specifically, Ras signaling pathway (oas04014), MAPK signaling pathway (oas04010), and PI3K-Akt signaling pathway (oas04151) were found, which are known to be associated with hair follicle development [36,37,38].

3.3. Co-Expression Analysis of Integrated DE miRNA-mRNA

To further understand the potential role of miRNAs in hair follicle development, an interaction network was constructed using DE miRNAs and mRNAs. In fact, target gene prediction of 37 DE miRNAs was performed to target 533 genes (Supplementary Table S8). Furthermore, we performed intersection analysis of 489 DE mRNAs and 533 predicted target genes; 33 overlapped genes were obtained (Figure 6A). Based on the principle of opposite miRNA and mRNA expression levels, a miRNA-mRNA co-expression network was constructed, where 8 DE miRNAs were targeted to 16 genes. The co-expression network revealed that 7 up-regulated miRNAs targeted 15 down-regulated genes, while 1 down-regulated miRNA targeted 1 up-regulated gene (Figure 6B).

3.4. Data Validation

To verify the accuracy of RNA-seq, RT-qPCR was performed. The latter showed similar expression patterns for both miRNAs and mRNAs to their respective RNA-seq results, which suggests reliability in the sequencing (Figure 6). In addition, the expression patterns of miR-125b and CD34, miR-181a and FGF12, LMO3, miR-200b, and ZNF536 displayed an opposite trend in SM and ST, coinciding with the miRNA-mRNA co-expression network result (Figure 7 and Figure 8).

4. Discussion

4.1. miRNAs and Hair Follicle Development

Hair follicle development is controlled by the changes of multiple cellular activities in the hair follicle which, in turn, are regulated by various transcription factors, signaling pathways, and epigenetic regulation. Among them, miRNA-mRNA regulatory networks participate in hair follicle cellular activity and thereby influence hair follicle development [14]. Some studies have suggested a tight association between the latter and wool curvature or, more specifically, that dermal papilla and other hair follicle cell activities are responsible for wool curvature via direct or indirect ways [39]. By far, the studies show that a lack of the miRNA processors Dier may lead to abnormal development, suggesting the indispensable role of miRNAs in hair follicle growth and development [40]. Previously, miR-125b and miR-205 were confirmed to respectively act as an activator and repressor of hair follicle stem cell differentiation [15,16], while miR-214 and miR-31 have been found to be highly expressed in proliferation hair matrix cells [17,18]. As such, it is not only necessary but also feasible to study hair follicle development from a miRNA perspective.

4.2. DE miRNAs Analysis in Different Degree of Curliness

In a previous study, we identified 13 DE-miRNAs between small and large wave patterns in Hu sheep [41]. In the current study, as more DE-miRNAs were screened (seeing that the difference between SM and ST was greater than that between small and large waves), a total of 37 DE-miRNAs were identified between SM vs. ST via RNA-seq. Fortunately for the sake of comparison, three DE-miRNAs (oar-miR-10a, oar-miR-143, oar-let-7c) were found in both studies [41]. To date, more studies on the miRNA expression profiles in sheep or goats are being published (with some also revealing profiles during different stages of wool production). For example, Liu et al. [42] investigated the miRNA expression profile of Chinese tan sheep hair follicles with different degrees of curliness, and reported the identification of 49 DE miRNAs. Interestingly, some of these DE-miRNAs (oar-miR-99a, oar-let-7c, and oar-miR-125b) have been associated with differing curliness in the hair follicles of both Hu sheep and Chinese tan sheep. Li et al. [43], on the other hand, explored 1354 DE miRNAs in the telogen skin between fine-wool sheep and Liaoning cashmere goats and reported the highest abundance of oar-miR-143 and oar-let-7c for these two groups. Similarly, Liu et al. [44] indicated that oar-let-7 and oar-miR-200 families might be the critical miRNAs in sheep hair follicle development.
Thus, in terms of available literature, many previously reported DE-miRNAs such as oar-miR-99a, oar-let-7c, oar-miR-125b, oar-miR-143, and oar-miR-200 were also present in our study. We, therefore, hypothesize that these miRNAs are likely to be involved in the growth and development of hair follicles in Hu sheep and influence the pattern formation of lambskin.

4.3. Functional Enrichment Analysis of DE miRNAs between SM vs. ST

To better understand the function of miRNAs, we performed functional annotation of miRNAs target genes. Research on hair follicle growth and development revealed that hair follicle development depended on a series of signaling molecules that participated in dermal and epithelial cell interactions—leading to the ordered proliferation and differentiation of the two cell populations, and ultimately resulting in the formation of intact hair follicles [45]. Fortunately, some target genes were enriched in an epithelial cell process, including the epithelial cell differentiation (GO:0030855), regulation of epithelial cell differentiation (GO:0030856), and epidermal growth factor catabolic process (GO:0007174), which were closely associated with hair follicle development.
The KEGG pathway analysis of target genes between SM vs. ST, some signaling pathways including MAPK signaling pathway, Ras signaling pathway, and PI3k-Akt signaling pathway were enriched. Recent studies suggested that the MAPK signaling pathway plays a positional role in hair follicle development and that activated MAPK signaling positively affected dermal papilla cell proliferation [38,46]. Additionally, the periodic changes of the important regulators in MAPK signaling pathway were shown to influence other key regulatory factors or signalings of Notch and FGF signaling pathways to regulate the development and differentiation of the hair follicle [47,48]. Surprisingly, activated PI3k-Akt signaling promoted the expression level of the crucial gene TGF-β2 that forms part of the TGF-β signaling pathway thus indirectly affected hair regeneration [49]. Especially, activating PI3k-Akt signaling through TPA treatment could also increase the differentiation rate of hair follicle stem cells [37], whereas the Ras signaling pathway could also regulate hair follicle morphogenesis via acting as a role of regulating the expression of Shh signaling pathway [36].

4.4. Co-Expression Analysis of miRNA-mRNA

To explore the regulatory mechanism of miRNA, co-expression analyses of integrated miRNA-mRNA were performed to construct the interaction network by DE miRNAs and target genes. Regarding SM vs. ST, a total of 16 overlapped genes were targeted by 8 DE miRNAs.
Studies have revealed the consistent role of miR-125b in controlling the activity of various skin stem cells and its inhibition of stem cell proliferation by acting on the TGF-β/BMP signaling pathway [50,51,52,53]. It’s target gene, vitamin D receptor (VDR), could have a positive role in hair follicle differentiation [15], while in our study, miR-125b targeted CD34 and ZNF687. CD34 was the marker of hair follicle-associated-pluripotent stem cells, regenerating part of hair follicle [54]. ZNF687, on the other hand, is a member of the zinc finger family and zinc finger proteins (structural proteins) have been described to affect the hair follicle cycle. Among them, the zinc finger transcription factor Trps1 is specifically expressed in the mesenchymal cell nucleus during hair follicle morphogenesis [55,56] and ZNF687 may be similarly involved [57]. MiR-200b was predicted to target 14 genes, including MAP4K5, CTNND2, PAK3, SLC4A4, ZNF76, ZNF831, etc. Previously, the miR-200 family has been reported as being highly expressed in epithelial cells where they contribute to the regulation of proliferation in hair morphogenesis [58]. As such, there is a possibility that its target genes may also be involved in hair follicle morphogenesis. It has also been found that miR-29a/b1 could target Bmpr1a and Lrb6 to inhibit hair follicle stem cell differentiation via suppressing BMP and WNT signaling [59]. Considering that NAV3 was the target gene of miR-29a, the role of NAV3 in hair follicle morphogenesis was also worth exploring. In addition, miR-431 targeted CD34, FGF12, and ZNF831, which might be participated in hair follicle morphogenesis [60,61]. Overall, these key miRNAs may therefore be either directly or indirectly associated with hair follicle activity, and further research is needed to explore their function(s) regarding wool curvature.
In summary, the current study provided a Hu sheep hair follicle miRNA expression profile and unraveled its miRNA-mRNA interaction network. Several DE miRNAs, including miR-125b, miR-10a, miR-10b, miR-200b, miR-143, miR-29b, etc., were identified while miRNAs-mRNAs interaction relationships, which may be involved in hair follicle morphogenesis (miR-200b targeting MAP4K5, miR-125b targeting CD34, and miR-431 targeting FGF12), were highlighted. These results will provide a great basis for further exploration of wool curvature—such as the use of molecular breeding to produce new strains of Hu sheep with high-quality lambskin.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/genes13020342/s1, Figure S1: Statistics analysis of miRNA chromosome distribution; Figure S2: Pearson correlation between samples; Table S1: Primers information of miRNAs and mRNAs; Table S2: Overview of the quality control of miRNA reads generated from hair follicle tissues; Table S3: DE miRNAs from SM vs. ST in hair follicle; Table S4: Overview of the quality control of mRNA reads generated from hair follicle tissues; Table S5: DE mRNAs from SM vs. ST in hair follicle; Table S6: GO enrichment analysis of target genes; Table S7: KEGG enrichment analysis of target genes; Table S8: Target genes prediction analysis of DE miRNAs.

Author Contributions

X.L. and W.S. designed this study; X.L. and W.C. performed the experiments; S.W., X.C. and Z.Y. contributed to materials and data collection in this study; X.L. and W.C. analyzed the data; X.L. and W.C. wrote the manuscript; T.G., J.M.M. and A.H. revised the final manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (31872333, 32172689), National Natural Science Foundation of China—CGIAR (32061143036), Major New Varieties of Agricultural Projects in Jiangsu Province (PZCZ201739), Key Research and Development Plan (modern agriculture) in Jiangsu Province (BE2018354) and Jiangsu Agricultural Science and Technology Innovation Fund (CX(18)2003).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The sequence data from this article have been submitted to SRA with the bioproject numbers PRJNA647222 and PRJNA632702.

Conflicts of Interest

All the authors declare no competing interests.

References

  1. Jiang, L.; Schlesinger, F.; Davis, C.A.; Zhang, Y.; Li, R.; Salit, M.; Gingeras, T.R.; Oliver, B. Synthetic spike-in standards for RNA-seq experiments. Genome Res. 2011, 21, 1543–1551. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Demars, J.; Cano, M.; Drouilhet, L.; Plisson-Petit, F.; Bardou, P.; Fabre, S.; Servin, B.; Sarry, J.; Woloszyn, F.; Mulsant, P.; et al. Genome-Wide Identification of the Mutation Underlying Fleece Variation and Discriminating Ancestral Hairy Species from Modern Woolly Sheep. Mol. Biol. Evol. 2017, 34, 1722–1729. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Reddy, S.; Andl, T.; Bagasra, A.; Lu, M.M.; Epstein, D.J.; Morrisey, E.E.; Millar, S.E. Characterization of Wnt gene expression in developing and postnatal hair follicles and identification of Wnt5a as a target of Sonic hedgehog in hair follicle morphogenesis. Mech. Dev. 2001, 107, 69–82. [Google Scholar] [CrossRef]
  4. Zhou, P.; Byrne, C.; Jacobs, J.; Fuchs, E. Lymphoid enhancer factor 1 directs hair follicle patterning and epithelial cell fate. Genes Dev. 1995, 9, 700–713. [Google Scholar] [CrossRef] [Green Version]
  5. Ouspenskaia, T.; Matos, I.; Mertz, A.F.; Fiore, V.F.; Fuchs, E. WNT-SHH Antagonism Specifies and Expands Stem Cells prior to Niche Formation. Cell 2016, 164, 156–169. [Google Scholar] [CrossRef] [Green Version]
  6. Chen, D.M.; Jarrell, A.; Guo, C.T.; Lang, R.; Atit, R. Dermal β-catenin activity in response to epidermal Wnt ligands is required for fibroblast proliferation and hair follicle initiation. Development 2012, 139, 1522–1533. [Google Scholar] [CrossRef] [Green Version]
  7. Li, A.; Lai, Y.C.; Figueroa, S.; Yang, T.; Widelitz, R.B.; Kobielak, K.; Nie, Q.; Chuong, C.M. Deciphering principles of morphogenesis from temporal and spatial patterns on the integument. Dev. Dyn. 2015, 244, 905–920. [Google Scholar] [CrossRef] [Green Version]
  8. Zhang, Y.; Tomann, P.; Andl, T.; Gallant, N.M.; Huelsken, J.; Jerchow, B.; Birchmeier, W.; Paus, R.; Piccolo, S.; Mikkola, M.L.; et al. Reciprocal requirements for EDA/EDAR/NF-kappaB and Wnt/β-catenin signaling pathways in hair follicle induction. Dev. Cell 2009, 17, 49–61. [Google Scholar] [CrossRef] [Green Version]
  9. Cetera, M.; Leybova, L.; Joyce, B.; Devenport, D. Counter-rotational cell flows drive morphological and cell fate asymmetries in mammalian hair follicles. Nat. Cell Biol. 2018, 20, 541–552. [Google Scholar] [CrossRef]
  10. Qiu, W.; Li, X.; Tang, H.; Huang, A.S.; Panteleyev, A.A.; Owens, D.M.; Su, G.H. Conditional activin receptor type 1B (Acvr1b) knockout mice reveal hair loss abnormality. J. Invest. Derm. 2011, 131, 1067–1076. [Google Scholar] [CrossRef] [Green Version]
  11. Wu, J.H.; Zhang, Y.J.; Zhang, J.X.; Chang, Z.L.; Li, J.Q.; Yan, Z.W.; Zhang, W.G. Hoxc13/β-catenin Correlation with Hair Follicle Activity in Cashmere Goat. J. Integr. Agric. 2012, 11, 1159–1166. [Google Scholar] [CrossRef]
  12. Cai, Y.; Yu, X.; Hu, S.; Yu, J. A brief review on the mechanisms of miRNA regulation. Genom. Proteom. Bioinform. 2009, 7, 147–154. [Google Scholar] [CrossRef] [Green Version]
  13. Vishnoi, A.; Rani, S. MiRNA Biogenesis and Regulation of Diseases: An Overview. Methods Mol. Biol 2017, 1509, 1–10. [Google Scholar] [CrossRef] [PubMed]
  14. Andl, T.; Botchkareva, N.V. MicroRNAs (miRNAs) in the control of HF development and cycling: The next frontiers in hair research. Exp. Derm. 2015, 24, 821–826. [Google Scholar] [CrossRef] [Green Version]
  15. Zhang, L.; Stokes, N.; Polak, L.; Fuchs, E. Specific microRNAs are preferentially expressed by skin stem cells to balance self-renewal and early lineage commitment. Cell Stem Cell 2011, 8, 294–308. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  16. Wang, D.; Zhang, Z.; O’Loughlin, E.; Wang, L.; Fan, X.; Lai, E.C.; Yi, R. MicroRNA-205 controls neonatal expansion of skin stem cells by modulating the PI(3)K pathway. Nat. Cell Biol. 2013, 15, 1153–1163. [Google Scholar] [CrossRef] [PubMed]
  17. Ahmed, M.I.; Alam, M.; Emelianov, V.U.; Poterlowicz, K.; Patel, A.; Sharov, A.A.; Mardaryev, A.N.; Botchkareva, N.V. MicroRNA-214 controls skin and hair follicle development by modulating the activity of the Wnt pathway. J. Cell Biol. 2014, 207, 549–567. [Google Scholar] [CrossRef] [Green Version]
  18. Mardaryev, A.N.; Ahmed, M.I.; Vlahov, N.V.; Fessing, M.Y.; Gill, J.H.; Sharov, A.A.; Botchkareva, N.V. Micro-RNA-31 controls hair cycle-associated changes in gene expression programs of the skin and hair follicle. FASEB J. 2010, 24, 3869–3881. [Google Scholar] [CrossRef] [Green Version]
  19. Baek, D.; Villen, J.; Shin, C.; Camargo, F.D.; Gygi, S.P.; Bartel, D.P. The impact of microRNAs on protein output. Nature 2008, 455, 64–71. [Google Scholar] [CrossRef] [Green Version]
  20. Stark, A.; Brennecke, J.; Russell, R.B.; Cohen, S.M. Identification of Drosophila MicroRNA targets. PLoS Biol. 2003, 1, E60. [Google Scholar] [CrossRef]
  21. Yu, Z.; Gordon, S.W.; Nixon, A.J.; Bawden, C.S.; Rogers, M.A.; Wildermoth, J.E.; Maqbool, N.J.; Pearson, A.J. Expression patterns of keratin intermediate filament and keratin associated protein genes in wool follicles. Differentiation 2009, 77, 307–316. [Google Scholar] [CrossRef] [PubMed]
  22. Liu, G.; Liu, R.; Tang, X.; Cao, J.; Zhao, S.; Yu, M. Expression profiling reveals genes involved in the regulation of wool follicle bulb regression and regeneration in sheep. Int. J. Mol. Sci. 2015, 16, 9152–9166. [Google Scholar] [CrossRef] [Green Version]
  23. Lv, X.; Chen, W.; Sun, W.; Hussain, Z.; Chen, L.; Wang, S.; Wang, J. Expression profile analysis to identify circular RNA expression signatures in hair follicle of Hu sheep lambskin. Genomics 2020, 112, 4454–4462. [Google Scholar] [CrossRef] [PubMed]
  24. Lv, X.; Chen, W.; Sun, W.; Hussain, Z.; Chen, L.; Wang, S.; Wang, J. Analysis of lncRNAs Expression Profiles in Hair Follicle of Hu Sheep Lambskin. Animals 2020, 10, 1035. [Google Scholar] [CrossRef] [PubMed]
  25. Sun, W.; Ni, R.; Yin, J.F.; Musa, H.H.; Ding, J.T.; Chen, L. Genome array of hair follicle genes in lambskin with different patterns. PLoS ONE 2013, 8, e68840. [Google Scholar] [CrossRef] [Green Version]
  26. Gao, W.; Sun, W.; Yin, J.; Lv, X.; Bao, J.; Yu, J.; Wang, L.; Jin, C.; Hu, L. Screening candidate microRNAs (miRNAs) in different lambskin hair follicles in Hu sheep. PLoS ONE 2017, 12, e0176532. [Google Scholar] [CrossRef]
  27. Erlich, Y.; Mitra, P.P.; delaBastide, M.; McCombie, W.R.; Hannon, G.J. Alta-Cyclic: A self-optimizing base caller for next-generation sequencing. Nat. Methods 2008, 5, 679–682. [Google Scholar] [CrossRef]
  28. Langmead, B.; Trapnell, C.; Pop, M.; Salzberg, S.L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009, 10, R25. [Google Scholar] [CrossRef] [Green Version]
  29. Wen, M.; Shen, Y.; Shi, S.; Tang, T. miREvo: An integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinform. 2012, 13, 140. [Google Scholar] [CrossRef] [Green Version]
  30. Friedlander, M.R.; Mackowiak, S.D.; Li, N.; Chen, W.; Rajewsky, N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012, 40, 37–52. [Google Scholar] [CrossRef]
  31. Zhou, L.; Chen, J.; Li, Z.; Li, X.; Hu, X.; Huang, Y.; Zhao, X.; Liang, C.; Wang, Y.; Sun, L.; et al. Integrated profiling of microRNAs and mRNAs: microRNAs located on Xq27.3 associate with clear cell renal cell carcinoma. PLoS ONE 2010, 5, e15224. [Google Scholar] [CrossRef]
  32. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Enright, A.J.; John, B.; Gaul, U.; Tuschl, T.; Sander, C.; Marks, D.S. MicroRNA targets in Drosophila. Genome Biol. 2003, 4, 1–27. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Chen, C.M.; Lu, Y.L.; Sio, C.P.; Wu, G.C.; Tzou, W.S.; Pai, T.W. Gene Ontology based housekeeping gene selection for RNA-seq normalization. Methods 2014, 67, 354–363. [Google Scholar] [CrossRef] [PubMed]
  35. Mao, X.; Cai, T.; Olyarchuk, J.G.; Wei, L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics 2005, 21, 3787–3793. [Google Scholar] [CrossRef]
  36. Mukhopadhyay, A.; Krishnaswami, S.R.; Cowing-Zitron, C.; Hung, N.J.; Reilly-Rhoten, H.; Burns, J.; Yu, B.D. Negative regulation of Shh levels by Kras and Fgfr2 during hair follicle development. Dev. Biol. 2013, 373, 373–382. [Google Scholar] [CrossRef] [Green Version]
  37. Qiu, W.; Lei, M.; Zhou, L.; Bai, X.; Lai, X.; Yu, Y.; Yang, T.; Lian, X. Hair follicle stem cell proliferation, Akt and Wnt signaling activation in TPA-induced hair regeneration. Histochem. Cell Biol. 2017, 147, 749–758. [Google Scholar] [CrossRef]
  38. Xiao, S.; Wang, J.; Chen, Q.; Miao, Y.; Hu, Z. The mechanism of activated platelet-rich plasma supernatant promotion of hair growth by cultured dermal papilla cells. J. Cosmet. Derm. 2019, 18, 1711–1716. [Google Scholar] [CrossRef] [PubMed]
  39. Nissimov, J.N.; Das Chaudhuri, A.B. Hair curvature: A natural dialectic and review. Biol. Rev. Camb. Philos. Soc. 2014, 89, 723–766. [Google Scholar] [CrossRef]
  40. Yi, R.; Pasolli, H.A.; Landthaler, M.; Hafner, M.; Ojo, T.; Sheridan, R.; Sander, C.; O’Carroll, D.; Stoffel, M.; Tuschl, T.; et al. DGCR8-dependent microRNA biogenesis is essential for skin development. Proc. Natl. Acad. Sci. USA 2009, 106, 498–502. [Google Scholar] [CrossRef] [Green Version]
  41. Lv, X.; Sun, W.; Yin, J.; Ni, R.; Su, R.; Wang, Q.; Gao, W.; Bao, J.; Yu, J.; Wang, L.; et al. An Integrated Analysis of MicroRNA and mRNA Expression Profiles to Identify RNA Expression Signatures in Lambskin Hair Follicles in Hu Sheep. PLoS ONE 2016, 11, e0157463. [Google Scholar] [CrossRef] [PubMed]
  42. Liu, Y.; Zhang, J.; Xu, Q.; Kang, X.; Wang, K.; Wu, K.; Fang, M. Integrated miRNA-mRNA analysis reveals regulatory pathways underlying the curly fleece trait in Chinese tan sheep. BMC Genom. 2018, 19, 360. [Google Scholar] [CrossRef] [PubMed]
  43. Li, J.; Qu, H.; Jiang, H.; Zhao, Z.; Zhang, Q. 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. DNA Cell Biol. 2016, 35, 696–705. [Google Scholar] [CrossRef]
  44. Liu, Y.; Wang, L.; Li, X.; Han, W.; Yang, K.; Wang, H.; Zhang, Y.; Su, R.; Liu, Z.; Wang, R.; et al. High-throughput sequencing of hair follicle development-related micrornas in cashmere goat at various fetal periods. Saudi J. Biol. Sci. 2018, 25, 1494–1508. [Google Scholar] [CrossRef] [PubMed]
  45. Saxena, N.; Mok, K.W.; Rendl, M. An updated classification of hair follicle morphogenesis. Exp. Dermatol. 2019, 28, 332–344. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  46. Sohn, K.M.; Jeong, K.H.; Kim, J.E.; Park, Y.M.; Kang, H. Hair growth-promotion effects of different alternating current parameter settings are mediated by the activation of Wnt/β-catenin and MAPK pathway. Exp. Dermatol. 2015, 24, 958–963. [Google Scholar] [CrossRef]
  47. Sawada, A.; Shinya, M.; Jiang, Y.J.; Kawakami, A.; Kuroiwa, A.; Takeda, H. Fgf/MAPK signalling is a crucial positional cue in somite boundary formation. Development 2001, 128, 4873–4880. [Google Scholar] [CrossRef]
  48. Jiang, L.; Xu, J.; Jin, R.; Bai, H.; Zhang, M.; Yang, S.; Zhang, X.; Zhang, X.; Han, Z.; Zeng, S. Transcriptomic analysis of chicken cochleae after gentamicin damage and the involvement of four signaling pathways (Notch, FGF, Wnt and BMP) in hair cell regeneration. Hear. Res. 2018, 361, 66–79. [Google Scholar] [CrossRef]
  49. Zhu, H.L.; Gao, Y.H.; Yang, J.Q.; Li, J.B.; Gao, J. Serenoa repens extracts promote hair regeneration and repair of hair loss mouse models by activating TGF-β and mitochondrial signaling pathway. Eur. Rev. Med. Pharm. Sci. 2018, 22, 4000–4008. [Google Scholar] [CrossRef]
  50. Battista, M.; Musto, A.; Navarra, A.; Minopoli, G.; Russo, T.; Parisi, S. miR-125b Regulates the Early Steps of ESC Differentiation through Dies1 in a TGF-Independent Manner. Int. J. Mol. Sci. 2013, 14, 13482–13496. [Google Scholar] [CrossRef] [Green Version]
  51. Emmrich, S.; Rasche, M.; Schoning, J.; Reimer, C.; Keihani, S.; Maroz, A.; Xie, Y.; Li, Z.; Schambach, A.; Reinhardt, D.; et al. miR-99a/100 similar to 125b tricistrons regulate hematopoietic stem and progenitor cell homeostasis by shifting the balance between TGF β and Wnt signaling. Gene Dev. 2014, 28, 858–874. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  52. Le, M.T.; Shyh-Chang, N.; Khaw, S.L.; Chin, L.; Teh, C.; Tay, J.; O’Day, E.; Korzh, V.; Yang, H.; Lal, A.; et al. Conserved regulation of p53 network dosage by microRNA-125b occurs through evolving miRNA-target gene pairs. PLoS Genet. 2011, 7, e1002242. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  53. Ooi, A.G.; Sahoo, D.; Adorno, M.; Wang, Y.; Weissman, I.L.; Park, C.Y. MicroRNA-125b expands hematopoietic stem cells and enriches for the lymphoid-balanced and lymphoid-biased subsets. Proc. Natl. Acad. Sci. USA 2010, 107, 21505–21510. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  54. Amoh, Y.; Hoffman, R.M. Hair follicle-associated-pluripotent (HAP) stem cells. Cell Cycle 2017, 16, 2169–2175. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  55. Fantauzzo, K.A.; Bazzi, H.; Jahoda, C.A.; Christiano, A.M. Dynamic expression of the zinc-finger transcription factor Trps1 during hair follicle morphogenesis and cycling. Gene Expr. Patterns 2008, 8, 51–57. [Google Scholar] [CrossRef]
  56. Oliver, C.H.; Nichols, J.; Watson, C.J. The KRAB domain zinc finger protein, Zfp157, is expressed in multiple tissues during mouse embryogenesis and in specific cells in adult mammary gland and skin. Genesis 2013, 51, 179–186. [Google Scholar] [CrossRef]
  57. Su, R.; Fan, Y.; Qiao, X.; Li, X.; Zhang, L.; Li, C.; Li, J. Transcriptomic analysis reveals critical genes for the hair follicle of Inner Mongolia cashmere goat from catagen to telogen. PLoS ONE 2018, 13, e0204404. [Google Scholar] [CrossRef] [Green Version]
  58. Hoefert, J.E.; Bjerke, G.A.; Wang, D.; Yi, R. The microRNA-200 family coordinately regulates cell adhesion and proliferation in hair morphogenesis. J. Cell Biol. 2018, 217, 2185–2204. [Google Scholar] [CrossRef] [Green Version]
  59. Ge, M.; Liu, C.; Li, L.; Lan, M.; Yu, Y.; Gu, L.; Su, Y.; Zhang, K.; Zhang, Y.; Wang, T.; et al. miR-29a/b1 Inhibits Hair Follicle Stem Cell Lineage Progression by Spatiotemporally Suppressing WNT and BMP Signaling. Cell Rep. 2019, 29, 2489–2504. [Google Scholar] [CrossRef] [Green Version]
  60. Hillege, M.M.G.; Galli Caro, R.A.; Offringa, C.; de Wit, G.M.J.; Jaspers, R.T.; Hoogaars, W.M.H. TGF-β Regulates Collagen Type I Expression in Myoblasts and Myotubes via Transient Ctgf and Fgf-2 Expression. Cells 2020, 9, 375. [Google Scholar] [CrossRef] [Green Version]
  61. Mossahebi-Mohammadi, M.; Quan, M.; Zhang, J.S.; Li, X. FGF Signaling Pathway: A Key Regulator of Stem Cell Pluripotency. Front. Cell Dev. Biol. 2020, 8, 79. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Figure 1. Identification of two pattern phenotype in Hu sheep. (A) Small waves. (B) Straight wool.
Figure 1. Identification of two pattern phenotype in Hu sheep. (A) Small waves. (B) Straight wool.
Genes 13 00342 g001
Figure 2. Identification and characterization of miRNA. (A) Length distribution of total sRNA fragment. (B) Classified identification of all sRNA in small waves. (C) Classified identification of all sRNA in straight wool. (D) First base preference analysis of known miRNA in small waves. (E) First base preference analysis of known miRNA in straight wool. (F) First base preference analysis of novel miRNA in small waves. (G) First base preference analysis of novel miRNA in straight wool.
Figure 2. Identification and characterization of miRNA. (A) Length distribution of total sRNA fragment. (B) Classified identification of all sRNA in small waves. (C) Classified identification of all sRNA in straight wool. (D) First base preference analysis of known miRNA in small waves. (E) First base preference analysis of known miRNA in straight wool. (F) First base preference analysis of novel miRNA in small waves. (G) First base preference analysis of novel miRNA in straight wool.
Genes 13 00342 g002
Figure 3. Significant difference analysis of miRNAs. (A) Volcano plot of identified miRNAs between SM vs. ST, where blue and red represent down-regulated and up-regulated, respectively. (B) Heat maps revealing the differentially expressed (DE) miRNA in hair follicle between SM vs. ST.
Figure 3. Significant difference analysis of miRNAs. (A) Volcano plot of identified miRNAs between SM vs. ST, where blue and red represent down-regulated and up-regulated, respectively. (B) Heat maps revealing the differentially expressed (DE) miRNA in hair follicle between SM vs. ST.
Genes 13 00342 g003
Figure 4. Significant difference analysis of mRNAs. (A) Volcano plot of identified mRNAs between SM vs. ST, where blue and red represent down-regulated and up-regulated, respectively. (B) Heat maps revealing the differentially expressed (DE) mRNA in hair follicle between SM vs. ST.
Figure 4. Significant difference analysis of mRNAs. (A) Volcano plot of identified mRNAs between SM vs. ST, where blue and red represent down-regulated and up-regulated, respectively. (B) Heat maps revealing the differentially expressed (DE) mRNA in hair follicle between SM vs. ST.
Genes 13 00342 g004
Figure 5. GO and KEGG functional enrichment analysis of target genes. (A) Top 30 GO term with biological process in SM vs. ST. (B) Top 30 GO term with cellular component in SM vs. ST. (C) Top 30 GO term with molecular function in SM vs. ST. (D) Top 30 KEGG pathways in SM vs. ST.
Figure 5. GO and KEGG functional enrichment analysis of target genes. (A) Top 30 GO term with biological process in SM vs. ST. (B) Top 30 GO term with cellular component in SM vs. ST. (C) Top 30 GO term with molecular function in SM vs. ST. (D) Top 30 KEGG pathways in SM vs. ST.
Genes 13 00342 g005
Figure 6. Co-expression network of miRNA-mRNA. (A) Overlapped genes between DE genes and target genes. (B) The interaction network of miRNAs and overlapped genes in SM vs. ST. “v” sharp represents DE miRNAs and circle represents overlapped genes. Red represents up-regulated miRNAs or genes and blue represents down-regulated miRNAs or genes.
Figure 6. Co-expression network of miRNA-mRNA. (A) Overlapped genes between DE genes and target genes. (B) The interaction network of miRNAs and overlapped genes in SM vs. ST. “v” sharp represents DE miRNAs and circle represents overlapped genes. Red represents up-regulated miRNAs or genes and blue represents down-regulated miRNAs or genes.
Genes 13 00342 g006
Figure 7. Verification of miRNAs sequencing by RT-qPCR between SM and ST. (A) Relative expression of let-7c. (B) Relative expression of miR-200b. (C) Relative expression of miR-181a. (D) Relative expression of miR-143. (E) Relative expression of miR-125b. (F) Relative expression of miR-10a. * represents statistical significance (p < 0.05) and ** represents extremely significant difference (p < 0.01).
Figure 7. Verification of miRNAs sequencing by RT-qPCR between SM and ST. (A) Relative expression of let-7c. (B) Relative expression of miR-200b. (C) Relative expression of miR-181a. (D) Relative expression of miR-143. (E) Relative expression of miR-125b. (F) Relative expression of miR-10a. * represents statistical significance (p < 0.05) and ** represents extremely significant difference (p < 0.01).
Genes 13 00342 g007
Figure 8. Verification of mRNAs sequencing by RT-qPCR between SM and ST. (A) Relative expression of SOBP mRNA. (B) Relative expression of LMO3 mRNA. (C) Relative expression of GABRA1 mRNA. (D) Relative expression of FGF12 mRNA. (E) Relative expression of CD34 mRNA. (F) Relative expression of ZNF536 mRNA. * represents statistical significance (p < 0.05) and ** represents extremely significant difference (p < 0.01).
Figure 8. Verification of mRNAs sequencing by RT-qPCR between SM and ST. (A) Relative expression of SOBP mRNA. (B) Relative expression of LMO3 mRNA. (C) Relative expression of GABRA1 mRNA. (D) Relative expression of FGF12 mRNA. (E) Relative expression of CD34 mRNA. (F) Relative expression of ZNF536 mRNA. * represents statistical significance (p < 0.05) and ** represents extremely significant difference (p < 0.01).
Genes 13 00342 g008
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Lv, X.; Chen, W.; Wang, S.; Cao, X.; Yuan, Z.; Getachew, T.; Mwacharo, J.M.; Haile, A.; Sun, W. Integrated Hair Follicle Profiles of microRNAs and mRNAs to Reveal the Pattern Formation of Hu Sheep Lambskin. Genes 2022, 13, 342. https://doi.org/10.3390/genes13020342

AMA Style

Lv X, Chen W, Wang S, Cao X, Yuan Z, Getachew T, Mwacharo JM, Haile A, Sun W. Integrated Hair Follicle Profiles of microRNAs and mRNAs to Reveal the Pattern Formation of Hu Sheep Lambskin. Genes. 2022; 13(2):342. https://doi.org/10.3390/genes13020342

Chicago/Turabian Style

Lv, Xiaoyang, Weihao Chen, Shanhe Wang, Xiukai Cao, Zehu Yuan, Tesfaye Getachew, Joram M. Mwacharo, Aynalem Haile, and Wei Sun. 2022. "Integrated Hair Follicle Profiles of microRNAs and mRNAs to Reveal the Pattern Formation of Hu Sheep Lambskin" Genes 13, no. 2: 342. https://doi.org/10.3390/genes13020342

APA Style

Lv, X., Chen, W., Wang, S., Cao, X., Yuan, Z., Getachew, T., Mwacharo, J. M., Haile, A., & Sun, W. (2022). Integrated Hair Follicle Profiles of microRNAs and mRNAs to Reveal the Pattern Formation of Hu Sheep Lambskin. Genes, 13(2), 342. https://doi.org/10.3390/genes13020342

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop