High-throughput sequencing of hair follicle development-related micrornas in cashmere goat at various fetal periods.

Inner Mongolia cashmere goat marks a precious gerplasm genetic resource due to its excellent cashmere traits. Therefore, it is of crucial importance to investigate the cashmere development mechanism of cashmere goat and to search for the important cashmere growth-related candidate genes. Fetal skin samples at 10 different periods of cashmere goat were collected in this research. Moreover, high-throughput sequencing was conducted on RNA samples from side skin of cashmere goat fetuses collected at three critical periods of skin hair follicle initiation, growth and development (namely, 45, 55 and 65 days) after balanced mix in line with the previous research results. Meanwhile, 3 samples at corresponding periods were used as the biological duplications. Data regarding microRNA and mRNA expression in skin and hair follicles of cashmere goats at various fetal periods were obtained using the high-throughput sequencing method. The results indicated that microRNAs in the oar-let-7 and oar-miR-200 families in 55 days and 66 days of pregnancy samples had been notably up-regulated relative to those in 45 days of pregnancy samples. This revealed that they might be the critical microRNAs in hair follicle development.

Introduction

Research background

Cashmere goat industry is an important component of animal husbandry in China. Cashmere, which is famous at home and abroad for its pure white color, softness, slender texture and good gloss, is an important animal by-product for foreign exchange earning in China (Zhang et al., 2005a, Zhang et al., 2005b). Cashmere is the derivative of goat skin, with primary hair follicle growing coarse wool while secondary hair follicle growing cashmere. The secondary hair follicle in cashmere goat is one of the reproducible hairs, which shows cyclic variation and controls the growth and loss of cashmere. Cyclic variation regulation of hair follicle is a complicated physiological process, which is subject to the interaction of numerous signal molecules and signaling pathways (Zhang, 2007). A series of signal regulatory molecules are extensively involved in all links of hair follicle morphogenesis through signaling pathways such as Wnt, Notch, Shh, TGF-β, EGF, FGF and BMP. They have played vital roles in the development and morphogenesis of all types of hair follicles (Yuan, 2014). The cyclic development of hair follicle is controlled by the complicated network structure of signaling pathway. Different signaling pathways can balance hair follicle development through inhibiting or promoting the hair follicle development signal molecules. Eventually, they can transfer the information to the hair follicles and peri-hair follicle tissues, thus regulating the initiation and cease of hair follicle development and activities as well as the rate of hair follicle mitosis (Zhang et al., 2014, Gu et al., 2013).

Growth trait and density of hair follicles can partly affect the yield and quality of cashmere. As a result, research on the molecular regulatory mechanism of the hair follicle growth variation in cashmere goat has become a research focus, which has also become the major task that we encounter (Jiang, 2015, Zhang et al., 2014a, Zhang et al., 2014b, Memili, 1998). Currently, some researchers at home and abroad have concentrated on carrying out molecular regulatory mechanism research on the hair follicle cyclic growth in cashmere goat. However, molecular regulatory mechanism research on the initiation, growth and development process of skin hair follicles in cashmere goat at the fetal period is lacking. Moreover, reports regarding regulating the cyclic development of skin hair follicles at fetal period by microRNA are not available (Baley and Li, 2012a, Shuhong et al., 2006, Cutting et al., 2011, Sontakke et al., 2014, Weng et al., 2014). As a kind of skin derivative, hair follicle has complicated structure and is constituted by multi-layer of special cells (Zhihong et al., 2012, Gebremedhn et al., 2014, Zielak-Steciwko et al., 2014, Shyh-Chang et al., 2013, Baley and Li, 2012a). Hair follicle is the only organ that shows cyclic growth throughout the life of mammal.

Literature review

Hair is the product of hair follicle cell proliferation and differentiation (Rong et al., 2011). Morphogenesis of hair follicle depends on the balance between keratinocyte proliferation and apoptosis. In addition, the induction and development of embryonic hair follicle is regulated by interstitial cell and fibroblast locating in hair basilar part, which is the interaction consequence of epithelial-interstitial cells. 3 signals are involved in the hair follicle morphogenesis, which are epithelial thickening, hair papilla formation and deepening of hair follicle epithelium to the dermis. The interaction between epithelial cells on skin hair follicles and interstitial cells on dermis has promoted cell population proliferation and differentiation, thus hair follicle development can be completed eventually. The continuous and complicated interaction between epithelial cell and hypodermal cell regulates the hair follicle morphogenesis and the subsequent hair cycle (Hossain et al., 2013, Chao et al., 2013).

Firstly, the hypodermal cell sends the initial signal, which induces the formation of hair bud on the epidermis. The hair bud can thus release certain factors to induce the formation of hair papilla by dermal fibroblast. Subsequently, the hair papilla releases the second signal to stimulate epithelial cell proliferation and differentiation, thus forming the complete hair follicle structure (Wenguang et al., 2007, Su, 2011, Jiang, 2012). Placode differentiation is the first step of hair follicle formation, which is formed by epidermal invagination. Later, local basal layer cell protrudes to interior dermis to form the hair bud (Wang, 2014). The hypodermal cells also aggregate into a cluster under the hair bud, which makes up the subsequent hair papilla (He et al., 2011). Then, the hair bud continues to stretch deep into the dermis, which forms a solid columnar structure. Its end has surrounded more dermal fibroblasts, which forms a cap-like structure. The further growth of columnar structure has rendered end expansion, which forms a hair bulb (Ouji et al., 2013). The hair bulb has recessed, thus part of the surrounding fibroblasts and mesenchymal components have invaginated to form hair papilla, while the other part of fibroblasts has formed the connective tissue sheath (Zhang et al., 2012). The hair follicle grows downward; in the meantime, the hair matrix cells begin to rapidly divide and migrate upward, which have differentiated into IRS and hair shaft. The hair follicle can reach the deep dermis as it grows; meanwhile, its accessory structures such as sebaceous gland and arrector pili have completely appeared. Moreover, the differentiated hair shaft grows upward gradually, which has eventually penetrated the body surface (Hu and Wang, 2012).

Generally, the hair cycle includes three periods, namely, anagen, catagen and telogen (Cui and Ren, 2012, Ma et al., 2015, Botchkarev et al., 2002, Yang et al., 2013, Yu, 2015). New hair shaft will come into being while the aging hair shaft will eventually fall off in each hair cycle. Signal molecules regulating the interaction between epithelium and dermis include Wnt/wingless signaling pathway, hedgehog, fibroblast growth factor (FGF), tumor necrosis factor (TNF), transforming growth factor-b/bone morphogenetic protein (TGF-b/BMP) and Shh signal pathways. The different combinations of these signals may determine the development of hair follicle (Zhang et al., 2014).

Wnt/β-catenin signaling pathway: Wnt ligand protein is a lipid-modified extracellular glycoprotein, which can be secreted in the activation of wntless protein (Guo, 2014). The extracellular Wnt protein binds with Frizzled receptor and low-density lipoprotein (LRP) on target cell membrane, and accumulates β-catenin in cytoplasm through the Disheveled protein. Later, β-catenin enters cell nucleus, acts with the TCF/LEF transcription factor family, and promotes expression of specific genes (Bai et al., 2013, Cai, 2014, Luo, 2014). β-catenin is an important adhesion molecule between epithelial cells, which is necessary to induce the production of hair follicle placode and adult animal hair follicle keratinocytes (He et al., 2015). β-catenin mutation at embryonic development period will lead to blocked hair follicle placode formation. Stem cell cannot differentiate into hair follicle keratinocyte in the absence of β-catenin (Teh and Blaydon, 2007).

TGF- b/BMP signal: BMP, which belongs to the TGF superfamily, can realize signal transduction through binding with the transmembrane receptor (bone morphogenetic protein receptor) on cell membrane. The BMP family members mainly play inhibitory roles in hair follicle morphogenesis (Qu, 2014). TGF-b2 can activate Smad2/3, which prompts the tissue regeneration of hair follicle stem cells; in contrast, hair follicle regeneration is remarkably delayed when hair follicle stem cells cannot sense the TGF-b 2 signal. TGF-b2 signaling pathway acts against the BMP pathway in hair follicle stem cell, which can activate hair follicle stem cell through the direct target gene Tmeff1 of TGF-b 2/Smad 2/3. Abnormal BMP2 and BMP4 expression has suppressed the formation of hair follicle embryo; while noggin and chordin are inhibitors of BMP4. In noggin-deletion mouse embryo, induction of secondary hair follicle formation has been completely ceased, while the morphogenesis of primary hair follicle remains unaffected (Su et al., 2013).

FGF: fibroblast growth factor (FGF) has participated in proliferation and differentiation processes of numerous cells and tissues as a key participator. FGF5 is a key inducible factor at catagen, the deletion of which will result in prolonged hair follicle anagen in mouse (Gao et al., 2014). FGF18-knockout mouse has notably shortened telogen, while FGF18 in wild-type mouse can negatively regulate hair follicle growth at anagen (Li, 2011).

The miR-199 family is highly expressed in hair follicle, while it cannot be detected in epidermis, and miR-200 family member deriving from two genetic loci tends to be expressed in both hair follicle and epidermis. Many microRNAs have been extensively expressed in other tissues, such as let-7 family, miR-21 family and miR-17 family. Moreover, some microRNAs enrich in epidermal tissues, including miR-203 (Li et al., 2015). Millar and his colleagues determined the hair follicle specific microRNAs using the previously constructed Dkk1 overexpression transgenic mouse model. It could be seen in the model that; hair development had almost been eliminated since the Wnt signaling pathway was greatly suppressed. Some microRNAs, such as miR-200b and miR-196a, have been considered as the critical candidate microRNAs for hair follicle development, since their expression levels have been down-regulated in skin of Dkk1 transgenic mouse (Philpott and Sanders, 1994).

This research is carried out based on the new generation high-throughput sequencing platform, and the differentially expressed microRNAs in skin hair follicle of fetal cashmere goat are screened out. In the meantime, combined analysis of microRNA and mRNA sequencing data is carried out by means of the bioinformatic method, to search for the microRNAs related to cyclic development of skin hair follicle in fetal cashmere goat as well as the target genes (Zhang et al., 2005a, Zhang et al., 2005b). Moreover, the functions of microRNAs in the genesis and development of hair follicle in fetal cashmere goat are analyzed, with an aim to provide theoretical foundation for realizing artificial regulation of cashmere growth and accelerating the genetic breeding of cashmere goat.

Achievements concerning the changes in cyclic development of hair follicle have been attained in human and mouse. In addition, some basic data with important reference value have been accumulated from research regarding the molecular regulatory mechanism of cyclic growth and development of cashmere goat secondary hair follicle. Firstly, genes related to the growth and development of cashmere goat secondary hair follicle are selected based on methods such as gene library construction and subtractive hybridization inhibition. Secondly, genes related to growth and development of cashmere goat secondary hair follicle are cloned and their expression patterns in skin tissues are analyzed, to reveal the roles of genes in cyclic variation of secondary hair follicle. However, the molecular regulatory mechanism of the growth and development of skin hair follicle in fetal cashmere goat remains unclear (Sam et al., 2015).

microRNA is a class of endogenous non-coding single-strand RNA, which is about 22 nucleotides in length. It binds with the specific site of target gene mRNA through base complementation; meanwhile, it inhibits protein synthesis or induces mRNA degradation, thus regulating target gene expression at post-transcription level. Moreover, it plays an important role in numerous processes such as organism growth and development (Janssen et al., 2013).

Andl et al. and Yi et al. had reported respectively that many microRNAs had been found in skin of fetal mouse, which played vital roles in hair follicle development and maintenance (Pritchard et al., 2012). Mardaryev et al. further discovered that microRNA-31 played an important regulatory role in expression of hair follicle cyclic variation-related genes in mouse skin. In addition, it was involved in establishing the balanced expression of hair follicle genes required for normal hair growth and hair fibrillation. Yin Jun et al. constructed the small molecular library of 70-day-old fetal cashmere goat skin and discovered that 4 microRNAs were expressed in cashmere goat and mouse skin. Furthermore, high expression could be seen in mouse hair follicle, which was speculated to be directly related to hair follicle development. Lateral skin tissues from 70-day-old fetus and 2-week-old lamb were screened using microRNA chip. It was found that 4 microRNAs were up-regulated in lateral skin of 70-day-old fetus compared with that of 2-week-old lamb, while 64 were up-regulated in lateral skin of 2-week-old lamb relative to that of 70-day-old fetus. Liu et al. had determined the microRNA expression in skin from cashmere goat secondary hair follicle at anagen. 22 out of the identified 316 microRNAs were new microRNAs, which shared the same positions on chromosome with sheep microRNA. They were speculated to play important regulatory roles in skin hair follicle, but the precise regulatory network and target gene remained unclear. Therefore, targeted experiment should be designed, to verify the functions of related microRNAs in the cyclic development of skin hair follicle at fetal period at a deeper level.

Let-7 is one of the microRNAs that have been first studied. The sequence and function of let-7 are highly conserved among various species. A total of 12 genes have been discovered to be coded by let-7 in human genome, which have constituted the let-7 family, including let7a, let7b and let7c. let-7 family is one of the most thoroughly studied microRNAs, which plays an important regulatory role in cell proliferation, differentiation and tumor formation by post-transcriptional regulatory method (Xue et al., 2013, Pan et al., 2014). Deletion or low expression of miR-let7 may result in the high expression of oncogenes Fos and Sis, thus promoting tumor formation. However, high miR-let7 expression will inhibit expression of multiple oncogenes, such as Ras and C-myc, thus suppressing tumor formation.

Besides, it has been discovered that there are genes expressing mouse skin and hair follicle in microRNA. It is pointed out in literature to isolate and clone numerous mouse E17.5 epidermis and hair follicle RNA, thus the microRNA abundance can be calculated through the clone efficiency. Similar results can be obtained by a distinct method. Subsequently, Wnt inhibitor Dkk1 hair follicle developer gene microRNA-expressed newly born mouse skin matrixes are compared and analyzed. The control group and Dkk1 treatment group microRNA expression profiles are constructed, which has further verified that microRNA is an important component of skin and hair follicle. Expression of miR-203 can hardly be detected in E13.5 skin (only a simple epithelial cell is present), but E15.5 miR-203 is the microRNA with the highest abundance. This suggests that it is expressed in skin layering and differentiation processes. Immunohistochemically detection on adult mouse skin has indicated that miR-203 is highly expressed in specific cells, including skin base layer or hair follicular root sheath or hair follicular base layer. It demonstrates that miR-203 is suppressed in zebrafish and outer skin, thus suggesting its conserved function among species. miR-203 transgene is expressed in the basal layer of epidermis of mice soon after birth, while death and histological analysis finds that the mice have thin epidermis, deletion of basal cell protein and excessive expression of microRNA. All these have reduced the cell proliferation in vitro and cell cycle arrest abilities of miR-203 (Zheng et al., 2014). Research results have indicated that epidermal keratinocyte is the most important target gene of miR-203 (Heo et al., 2012). miR-203 is an important regulatory factor in transcription factor stratified epithelium of different species. It can result in defective loss of all stratified epithelial cells in the meantime of inhibiting miR-203 performance as well as expression pattern. Bioinformatic analysis and experimental data have proved that the direct inhibition of miR-203 expression renders the cycle of basal layer cells into the bottom layer (Boon et al., 2013). Meanwhile, it can regulate the cycle of keratinocyte through controlling miR-34 family members (Boudreau et al., 2014). microRNA plays an important regulatory role in mature embryonic stem cells, which can predict and regulate skin stem cells. A recent study has indicated a large amount of miR-125b expression after differentiation and gradually decreased skin stem cell differentiation (Lin et al., 2012). Consequently, it is proposed that gene miR-125b can inhibit stem cell differentiation.

Mir-31 has been proved to be an extremely important factor regulating hair follicle cycle, which shows high expression in hair follicle at anagen while low expression at telogen. Mice receiving injection of miR-31 inhibitor at the early and moderate anagen display accelerated hair growth as well as altered matrix differential of hair keratinocytes and hair shaft formation. Analysis of the number of growth regulatory molecules and cytoskeletal proteins reveals that miR-31 in hair follicle is related to the balance required for hair growth and hair fibration (Åkerblom et al. 2013). In addition, the miRNA-200b and miRNA-196a WNT signaling pathways can control hair follicle growth through regulating the potential target genes. Expression quantities of miRNA-200b and miRNA-196a will be notably reduced in epidermis of transgenic mouse with Dkk1 overexpression.

Cashmere goat epidermis contains the primary and secondary hair follicles, which correspond to the products of wool and cashmere, respectively. The cycle of primary hair follicle is different from that of secondary hair follicle in cashmere goat. The primary hair follicle cycle lasts for a month, which can be classified into cell apoptosis and renewal. It displays little seasonal variation and obvious seasonal cycle compared with the secondary hair follicle. Moreover, it enters the anagen, catagen and telogen in one cycle at a time (Yu et al., 2015, Guang-Jun et al., 2013, Aritro et al., 2014, Song et al., 2013, Chandra et al., 2015). MicroNA plays an important role in cashmere goat skin. It is found in research that the lack of cap-like structure in RNA will render RNA death at the early embryonic development period in mice. However, removal of the cap-like structure at epidermal anagen will not lead to death of newly born mice at anagen. Histological examination has found that newly born mice with cap-like structure develop dysplasia defect in hair follicle, which will manifest as hair follicle cell proliferation, apoptosis and renewal, thus entering in the dermal tissue.

The technical route is shown as Fig. 1.

Fig. 1

Technical route.

Growth trait and density of hair follicles can partly affect the yield and quality of cashmere. As a result, research on the molecular regulatory mechanism of the hair follicle growth variation in cashmere goat has become a research focus, which has also become the major task that we encounter. Currently, some researchers at home and abroad have concentrated on carrying out molecular regulatory mechanism research on the hair follicle cyclic growth in cashmere goat. However, molecular regulatory mechanism research on the initiation, growth and development process of skin hair follicles in fetal cashmere goat is lacking. Moreover, reports regarding regulating the cyclic development of skin hair follicles at fetal period by microRNA are not available. Research on skin hair follicle morphogenesis, growth and development rules of fetal cashmere goat had been carried out previously. Meanwhile, skin sample library at fetal periods (55, 65, 75, 85, 95, 105, 115, 125 and 135 days, with 3 fetuses at each period) was constructed. Differentially expressed microRNAs in lateral skin hair follicles of cashmere goats at various fetal periods were analyzed using the high-throughput solexa sequencing method based on the previous research results. Moreover, a microRNA expression database was established. (The following periods had been determined based on previous research on hair follicle morphogenesis process and development rule of cashmere goat. 45 days of embryonic period: dermal signal initiation; 55 days of embryonic period: primary hair follicle genesis; 65 days of embryonic period: secondary hair follicle genesis; 105 days of embryonic period: rapidest growth of primary hair follicle; 115 days of embryonic period: rapidest growth of secondary hair follicle.) microRNAs related to skin hair follicle genesis and development of cashmere goat at fetal periods were screened out. The functions of microRNAs in hair follicle genesis and development of cashmere goat at fetal periods were analyzed, which could provide theoretical foundation for realizing artificial regulation of cashmere growth and accelerating the genetic breeding of cashmere goat.

Selection of micrornas related to hair follicle genesis and development in cashmere goat at fetal periods through high-throughput sequencing

Cashmere goat skin contains the primary and secondary hair follicles. The cycle of primary hair follicle is different from that of secondary hair follicle. The primary hair follicle cycle exists in each month, which can be classified into cell apoptosis and renewal. It displays little seasonal variation and shows obvious seasonal cycle compared with the secondary hair follicle (Song et al., 2013, Chandra et al., 2015, Aberdam and Candi, 2006). The production periodicity and complicated structure of secondary hair follicle in cashmere goat has controlled the growth and development of cashmere, which is also closely related to wool quality. However, research on microRNA expression changes during the cyclic variation of skin hair follicle in cashmere goat at fetal periods is rarely reported (Alonso and Fuchs, 2006, Botchkareva, 2012, Butcher and Denkers, 2002). Therefore, the small RNA library at anagen of skin hair follicle in cashmere goat at fetal periods was constructed in this experiment. The conserved microRNAs on cashmere goat skin at different periods were evaluated using the second-generation sequencing technique. Meanwhile, the new RNA was predicted, which was used for the subsequent hair follicle cyclic variation of cashmere goat at fetal periods as well as microRNA test and screening.

Experimental materials

The mother ewes aged 3–4 years were selected from Inner Mongolia cashmere goat farm as the experiment population. Ewes at various gestation periods were conducted caesarean section from March 2014 based on the breeding record. 1 cm2 lateral skin samples were collected at 10 periods from cashmere goat fetuses, namely, 45d, 55d, 65d, 75d, 85d, 95d, 105d, 115d, 125d and 135d. 3 fetuses were selected at each period as the biological duplications. Skin samples were washed with DEPC-treated water, numbered and placed rapidly into liquid nitrogen for quick-freezing. They were then placed in a refrigerator at −80 °C for use.

Results and analysis

Total RNA was extracted using the one-step acid guanidine thiocyanate-phenol-chloroform method. The plant tissues were ground in liquid nitrogen and mixed with the ice-cold guanidine thiocyanate; followed by 5 min of incubation on ice and centrifugation. Later, the mixture was incubated at 4 °C, and the RNA-containing supernatant was transferred to a new tube, followed by ethanol precipitation, washing and dissolving in the DEPC-treated water. RNA was further treated with phenol-chloroform and washed over twice (Chen et al., 2012, Huelsken et al., 2001, Kimura-Ueki et al., 2012, Oshimori and Fuchs, 2012). Purity and integrity of RNA were detected using agarose gel electrophoresis; meanwhile, the absorbance was detected at A260.

A total of 3 mg total RNA was used for the preparation of cDNA library of small molecule RNA. The RNA was connected from the 3′-end to 5′ in succession and was reversely transcribed into cDNA for PCR amplification. The total library was conducted 10% polyacrylamide gel electrophoresis, and the corresponding microRNA strips were extracted and eluted. The purified small RNA library was performed fluorescent quantitation after ethanol precipitation and washing, and was analyzed using Illumina Nextseq 500 single-end sequencing.

The electrophoresis detection is shown as Fig. 2. The small RNA cDNA library construction is shown as Fig. 3.

Fig. 2

Electrophoresis detection.

Fig. 3

cDNA library construction.

Result analysis: strips adjacent to 150 bp were collected through gel extraction

The Sample description (Samples are cashmere goats selected on Day 45, 55 and 65 of gestation) is shown as Table 1.

Table 1

Sample description (Samples are cashmere goats selected on Day 45, 55 and 65 of gestation).

Sample ID	Name	Description
Goat_45	Goat_45	Day 45 of gestation
Goat_55	Goat_55	Day 55 of gestation
Goat_65	Goat_65	Day 65 of gestation

Library construction and high-throughput sequencing had been successfully carried out on samples. The library construction results were sequenced using the Illumina Nextseq 500 sequencing platform, which had obtained the high-quality transcriptome data.

Method description: double-end sequencing was performed by means of the NextSeq500 sequencing platform in this project, with the single-end read length of 151nt. Clean reads referred to the clean reads after removing the adapter interloping.

Extraction strategy of clean reads

end clean strategy

Cut the adapter of sRNA and maintained reads with the length of over 17;

Got rid of basic groups with the mass of lower than 20 and maintained those with the length of 16;

Removed reads with the 30% basic group mass of lower than 20;

Got rid of the NN.

end clean strategy:

Cut the adapter of mRNA and maintained reads with the length of over 20;

Cut 1 basic group at the end;

Got rid of basic groups with the mass of lower than 20 and maintained those with the length of 16;

Removed reads with the 30% basic group mass of lower than 20;

Got rid of the NN.

The result presentation is shown as Table 2

Table 2

Obtain the high-quality clean reads.

Sample	Raw data	Clean data	Unique reads
Goat_45	5,946,209	5,253,789(88.36%)	901,652(15.16%)
Goat_55	9,252,829	7,998,727(86.45%)	1,178,643(12.74%)
Goat_65	6,049,887	5,172,144(85.49%)	1,012,503(16.74%)

Data analysis

The RAW reads of every sample have reached over 5.94 million.

The proportion of CLEAN was over 85–88%, which shows that the database is good enough;

The total sequencing was higher than 7.08 million reads/sample.

The sRNA-seq workflow is shown as Fig. 4.

Fig. 4

sRNA-seq workflow.

The clean reads sequencing quality analysis is shown as Fig. 5, Fig. 6, Fig. 7.

Fig. 5

Clean reads sequencing quality analysis (Cashmere_goat_45_fastqc).

Fig. 6

Clean reads sequencing quality analysis (Cashmere_goat_55_fastqc).

Fig. 7

Clean reads sequencing quality analysis (Cashmere_goat_65_fastqc).

Sequencing quality distribution

γ-axis represented the sequencing quality score, which could be calculated through the following formula:Note: p stood for the probability of incorrect corresponding basic group, while x-axis represented the location of each read. p = .001582385 when Q = 28.

Sequencing length: After removing the index sequence, using the random basic group in a balanced database, and cutting out basic groups of poor quality in the following part, the longest effective length. obtained was 147nt. Since this project is building the database for MicroRNA, sequences over 30nt-length belong to the adapter sequences

Sequencing quality analysis: The 100%Q scores of all sequences are over 30 from the 1st to 30th place, suggesting the excellence of database establishment and sequencing figure.

Conclusion: Excellent sequencing quality had been achieved in this project.

The effective length statistics is shown as Fig. 8, Fig. 9, Fig. 10.

Fig. 8

Effective length statistics (Cashmere_goat_45_fastqc).

Fig. 9

Effective length statistics (Cashmere_goat_55_fastqc).

Fig. 10

Effective length statistics (Cashmere_goat_65_fastqc).

Method description

The index sequence was removed, balance of the constructed library was tested using random basic group, and the low quality basic groups in the latter part were cut out. Finally, the obtained clean reads were used for effective length analysis.

Result presentation

Data analysis: Clean reads, refers to the basic groups after removing the mass value lower than 20. The effective sequencing length of a majority in the samples of this project ranges from 20 to 24nt, which suggests that only a small number of basic groups is in low quality, and that the effective insert in database establishment is eligible.

Representative analysis of cDNA library

Method description:

cDNA construction was processed to conduct random readation of the retrieved polyA-mRNA/ncRNA, which was added adapter and carried out RT-PCR. Excessive PCR amplification would lead to a large amount of reads with completely same sequences in the cDNA library. Reads with completely same sequence derived from excessive PCR amplification; in addition, abundant RNA would produce a large amount of 5′end same reads at the time of readation. PCR amplification cycles were strictly controlled in library construction, to guarantee that cDNA library could better represent the original mRNA/ncRNA expression quantity.

Calculation of PCR duplication level: 0.2 million reads were randomly selected from the sequencing data as the Total Reads, and calculated based on the following formula:

The result presentation is shown as Fig. 11, Fig. 12, Fig. 13.

Fig. 11

Result presentation (Cashmere_goat_45_fastqc), Sample of Day 45 cashmere goat PCR duplication level〉94.53%.

Fig. 12

Result presentation (Cashmere_goat_55_fastqc), Sample of Day 55 cashmere goat PCR duplication level >95.87%.

Fig. 13

Result presentation (Cashmere_goat_65_fastqc), Sample of Day 65 cashmere goat PCR duplication level >93.82%.

Data analysis: In this project, the PCR duplication level of all samples was around 93–95%.

The Comparison of Reads with reference genome is shown as Table 3.

Table 3

Mapping of clean reads on the reference genome.

Sample	Input	Total	Unique	Multiple
Cashmere_goat_45	5,253,789	4,960,403(94.42%)	3,533,470(71.23%)	1,426,933(28.77%)
Cashmere_goat_55	7,998,727	7,494,511(93.70%)	5,521,486(73.67%)	1,973,025(26.33%)
Cashmere_goat_65	5,172,144	4,750,883(91.86%)	3,548,127(74.68%)	1,202,756(25.32%)

Method description

Bowtie was used in this part with tolerance to mismatch of 1 basic group, to compare the effective sequencing data on the reference genome, which could be downloaded at http://goat.kiz.ac.cn/GGD/download.htm.

Result presentation

Data analysis: There are around 91–94% of reads can be mapped on the reference genome, among which 71–74% of reads are unique sequence alignment.

The reference genome information is shown as Table 4.

Table 4

Length distribution across Genomic Regions.

Species	CDS	Intergenic	Intron
Goat	29986545.00(1.10%)	2004425745.00(75.3%)	626904551.00(23.6%)

Data analysis: 1.10% of reference genome is distributed on CDS, 75.3% on intergenic, while 23.6% is on intron (3) Distribution of Reads across genomic regions, shown as Table 5, Fig. 14, Fig. 15, Fig. 16.

Table 5

Reads distribution across reference Genomic Regions.

Sample	CDS	Intergenic	Intron
Cashmere_goat_45	69358.94(1.40%)	3667378.71(73.93%)	1223665.35(24.67%)
Cashmere_goat_55	107152.99(1.43%)	5648481.36(75.37%)	1738876.66(23.20%)
Cashmere_goat_65	98700.90(2.08%)	3587757.63(75.52%)	1064424.48(22.40%)

Fig. 14

Cashmere_goat_45 Reads distribution across reference Genomic Regions.

Fig. 15

Cashmere_goat_55 Reads distribution across reference Genomic Regions.

Fig. 16

Cashmere_goat_65 Reads distribution across reference Genomic Regions.

Method description: The Clean Reads were compared on the reference genome, and distribution in all regions was counted.

Result presentation. The reads distribution across reference Genomic Regions is shown as Table 5.

The Reads distribution across reference Genomic Regions (Cashmere_goat_45) is shown as Fig. 14. The Reads distribution across reference Genomic Regions (Cashmere_goat_55) is shown as Fig. 15. The Reads distribution across reference Genomic Regions (Cashmere_goat_65) is shown as Fig. 16.

Data analysis: From the above figures, it is shown that Reads are mainly distributed in intergenic and intron regions, and secondly in CDS region.

Method description

Clean Reads of all samples were compared to Rfam database using the bowtie procedure; moreover, the reads that could be compared were classified based on RNA type.

Result presentation

The Proportion of clean reads mapped to Rfam is shown as Table 6. The Rfam classification is shown as Table 7. The Goat_45 Rfam classification is shown as Fig. 17. The Goat_55 Rfam classification is shown as Fig. 18. The Goat_65 Rfam classification is shown as Fig. 19.

Table 6

Proportion of clean reads mapped to Rfam.

Sample	Total input reads	Total mapped reads	Unipue mapped reads	Multiple mapped reads
Cashmere_goat_45	5,253,789	4,613,121(87.81%)	1,056,830(22.91%)	355,6291(77.09%)
Cashmere_goat_55	7,998,727	7,038,175(87.99%)	1,820,594(25.87%)	5,217,581(74.13%)
Cashmere_goat_65	5,172,144	4,417,209(85.40%)	968,614(21.93%)	3,448,595(78.07%)

Table 7

Rfam classification.

Sample	Cis-reg	Others	lncRNA	microRNA	rRNA	sRNA	snRNA	tRNA
Cashmere_goat_45	15,142	69,910	4725	4,067,689	146,483	3459	95,435	210,278
	(0.33%)	(1.52%)	(0.10%)	(88.18%)	(3.18%)	(0.07%)	(2.07%)	(4.56%)
Cashmere_goat_55	18,276	88,877	4339	6,248,739	222,632	4518	133,314	317,480
	(0.26%)	(1.26%)	(0.06%)	(88.78%)	(3.16%)	(0.06%)	(1.89%)	(4.51%)
Cashmere_goat_65	12,292	57,603	2535	387,9923	199,966	2512	107,063	155,315
	(0.28%)	(1.30%)	(0.06%)	(87.84%)	(4.53%)	(0.06%)	(2.42%)	(3.52%)

Fig. 17

Goat_45 Rfam classification.

Fig. 18

Goat_55 Rfam classification.

Fig. 19

Goat_65 Rfam classification.

Data analysis: From the above figures, it is shown that microRNA accounts for 88%, and tRNA and rRNA are in the second place.

Method description: The obtained clean reads were compared to mature microRNA downloaded from miRBase as well as our predicted microRNA library using bowtie, mismatching of 1 basic group was tolerable, and two optimum matchings were exported at most. The Clean reads mapping against mature microRNAs is shown as Table 8.

Table 8

Clean reads mapping against mature microRNAs.

Sample	Input	Total	Uniq	Multiple
Goat_45	5253789	2,455,491(46.74%)	2,254,982(91.83%)	200,509(8.17%)
Goat_55	7998727	4,125,742(51.58%)	3,804,662(92.22%)	321,080(7.78%)
Goat_65	5172144	2,344,798(45.34%)	2,239,021(95.49%)	105,777(4.51%)

Data analysis: From the above figures, it is shown that around 45–51% of reads can be compared to identified microRNAs, among which 91–95% were unique position matching.

Method description: clean reads that could not be compared to the mature microRNA were compared to the pre-microRNA library in cashmere goat miRBase using the bowtie procedure, with tolerance to 1 mismatching. The Reads distribution on the annotated pre-microRNA Regions is shown as Table 8 (see Table 9).

Table 9

Reads distribution on the annotated pre-microRNA Regions.

Sample	Input	Total	Uniq	Multiple
Goat_45	2,798,298	1,104,844(39.48%)	1,070,503(96.89%)	34,341(3.11%)
Goat_55	3,872,985	1,306,499(33.73%)	1,290,825(98.80%)	15,674(1.20%)
Goat_65	2,827,346	1,016,026(35.94%)	1,006,224(99.04%)	9802(0.96%)

Data analysis: 33–39% of reads can be compared to the pre-microRNA database.

Transcripts Per Million (TPM) was used to represent the standardized value of microRNA expression quantity, which revealed the number of reads contained in each microRNA in every 1 million reads.

Statistical analysis of expressed mature microRNAs is shown as Table 10. The Statistical analysis of expressed pre-microRNAs is shown as Table 11.

Table 10

Statistical analysis of expressed mature microRNAs.

Sample	TPM > 0	TPM > 10
Goat_45	628	385
Goat_55	634	353
Goat_65	611	340

Table 11

Statistical analysis of expressed pre-microRNAs.

Sample	TPM > 0	TPM > 10
mouse_assay	571	310
Goat_55	555	304
Goat_65	523	262

The Standard DEmicroRNA workflow is shown as Fig. 20

Fig. 20

DEG workflow.

Method description: First, correlation analysis between two samples was conducted, which meant to examine the correlation of the same microRNA expression between different samples. Similar expression quantities of most microRNAs suggested high uniformization quality of data. However, most microRNA expression quantities of a sample being remarkably higher than those of another demonstrated problem in data quality or statistically significant difference in microRNA expression pattern. The Sample correlation analysis is shown as Fig. 21.

Fig. 21

Sample correlation analysis.

Data analysis: Seen from the figures above, result suggests that the expression patterns of microRNA are similar in developmental phase of both primary and secondary hair follicle, when the samples of Day 55 and Day 65 of gestation are put together.

Method description: DEmicroRNAs were analyzed using the fisher test.

Two important parameters

FC: absolute ratio of expression change.

P-value: statistical confidence level of expression change.

Result presentation (fold change ≥2 or ≤0.5, p-value ≤.01)

The Total number of differentially expressed microRNA is shown as Table 12.

Table 12

Total number of differentially expressed microRNA.

DEmicroRNA sample	Total microRNA	Up regulated miR NO	Down regulated miR NO	Total DEmiR
Goat_55_vs_Goat_45	649	35	71	106
Goat_65_vs_Goat_45	649	37	82	119
Goat_65_vs_Goat_55	649	37	59	96

Data analysis: Seen from the above, it is obviously shown that the quantity of up-regulated microRNAs is larger than that of down-regulated microRNA when the later one compared to the former one.

Discussion

Research on hair follicle microRNA mainly concentrates on mouse and human based on hair removal. However, little research about microRNA in skin hair follicle of cashmere goat is available (Halin et al., 2017, Abd Samad et al., 2017). Wen-guang Zhang (2007) discovered expression of 159 microRNAs in goat skin and sheep skin through gene microarray; in addition, 9 microRNAs were specifically expressed in lateral skin of goats (cashmere growth region). Of them, oar-let-7b/d/f and oar-miR-200b/c has been detected in our experiment. A total of 649 miRNAs have been detected in cashmere goat samples in this project, including 151 known ora-miRNA recorded in miRBase database as well as 498 newly predicted miRNAs. The newly predicted miRNAs can be served as candidate miRNAs for experimental verification and differential expression analysis. Liu (2012) had detected microRNA expression in cashmere goat skin at catagen and discovered 316 conserved microRNAs and 22 new microRNA. In this project, 498 new conserved mature microRNAs have been discovered through sheep genomic sequence comparison (Razali and Said, 2017, Abdul Halim and Phang, 2017, Shamsudin et al., 2017). This may be related to differences in reference genome and goat variety (all experimental animals are produced by Inner Mongolia cashmere goats).

Summary

Among the reported miRNAs, oar-let-7b/d/f and oar-miR-200b/c in oar-miR-377-3p, oar-let-7 and oar-miR-200 families have been markedly up-regulated in 55 and 65 days of pregnancy samples relative to those in 45 days of pregnancy samples, while oar-miR-106a and oar-miR-218a have been evidently down-regulated. It suggests that they may be the critical miRNAs in hair follicle development.

Conclusion and innovations

Conclusion

microRNA and mRNA expression data in skin hair follicles at different fetal periods of cashmere goats are obtained through high-throughput sequencing. microRNAs in oar-let-7 and oar-miR-200 families have been distinctly up-regulated in 55 and 65 days of pregnancy samples compared with those in 45 days of pregnancy sample. It suggests that they may be the critical microRNAs in hair follicle development.

Innovations

At present, research achievements in terms of cashmere goat hair follicle in China has reached the international advanced level. High-throughput sequencing analysis at catagen of cashmere goats during fetal periods is a good way for obtaining secondary hair follicle genesis and development-related microRNAs. To the best of our knowledge, no related research has been reported at home and abroad, which marks the important innovation of this project.