The intragenic mRNA-microRNA regulatory network during telogen-anagen hair follicle transition in the cashmere goat
It is widely accepted that the periodic cycle of hair follicles is controlled by the biological clock, but the molecular regulatory mechanisms of the hair follicle cycle have not been thoroughly studied. The secondary hair follicle of the cashmere goat is characterized by seasonal periodic changes throughout life. In the hair follicle cycle, the initiation of hair follicles is of great significance for hair follicle regeneration. To provide a reference for hair follicle research, our study compared differences in mRNA expression and microRNA expression during the growth and repose stages of cashmere goat skin samples. Through microRNA and mRNA association analysis, we found microRNAs and target genes that play major regulatory roles in hair follicle initiation. We further constructed an mRNA-microRNA interaction network and found that hair follicle initiation and development were related to MiR-195 and the genes CHP1, SMAD2, FZD6 and SIAH1.
Cashmere goats generate wool and cashmere, which are produced in the primary hair follicles and secondary follicles in the skin1. The growth of hair follicles is cyclical throughout the cashmere goat’s life2, undergoing cyclic transformations from stages consisting of anagen (rapid growth) to catagen (apoptosis-driven regression) and back to anagen3. Paus R proposed that a biological “clock” drives hair follicle cycling4. If so, what are the key players, and how many key players are molecular controls? The transformation and growth of cashmere are regulated by various complex factors in the skin5,6,7, and each follicle growth period in the skin has a specific activated/silenced gene expression pattern8,9.
MicroRNAs are a class of endogenous non-coding RNAs with lengths of approximately 20–24 nt10. These sequences combine with the sequences of mRNA-specific regions to complement target gene expression at the level of translation and play a role in the regulation of gene expression11,12,13. MicroRNAs play a role in regulating the cashmere growth cycle by targeting various signalling pathways and transcription factors14. MiR-203 regulates the differentiation of epithelial keratinocytes and directly inhibits the expression of p6312,15. MiR-200b and MiR-196a are associated with hair follicle development16,17. MiR-31 is associated with hair matrix differentiation and hair stem formation18, and its expression increases significantly during villus growth and decreases during degenerative and regenerative periods19. Hence, changes in microRNA expression patterns in skin cells are closely related to the cashmere growth cycle20.
The regulation of gene expression in life is relatively complex. Through analysis of microRNA sequencing and gene transcriptome sequencing data21, we can study the molecular mechanisms of microRNA-regulated gene expression systematically and obtain a more accurate understanding of the regulation of gene expression. In 2012, Liu performed microRNA sequencing of cashmere goat skin at various stages to determine whether microRNA expression varies in the skin or hair. Differences in the expression of microRNAs may occur through key gene regulation signalling pathways that modulate hair follicle growth22, and microRNA-mRNA interactions may regulate the development of hair follicles and their growth cycle23. In the present study, to clarify the regulatory mechanism of the cashmere growth cycle and gain insight into the related gene regulatory network, we characterized the interactions between microRNAs and mRNAs in goat skin, aiming to identify the key microRNAs regulating the growth and development of hair follicles and targeted mRNAs. We also characterized the gene regulatory network of differentially expressed microRNAs and their target genes, performed an annotated analysis of target genes related to differentially expressed genes, and explored key genes associated with hair follicle initiation and the regulation of microRNAs.
Quality control of RNA and sequencing
Total RNA samples from six Inner Mongolian Arbas cashmere goats (Fig. 1a) were analysed. Total RNA was not degraded, and there were no stray DNA bands. The OD260/OD280 ranged from 1.8 to 2.1, and the RNA integrity value exceeded 8.6. The RNA integrity met the sequencing requirements (Fig. 1b). The base distribution and mass fluctuation analysis of each circle sequence (Fig. 1c) showed that A, T, C, G, and N began to fluctuate at each location and then tended to stabilize. (The percentages of each base vary by species.) The base distribution of this study was uniform, the N% proportion was stable and low, and the base quality of the library was satisfactory (Fig. 1d).
RNA quality inspection. (a) Agarose electrophoresis of total RNA from cashmere goat skin. (b) RNA electrophoretic and RIN values of the skin. (c) Base distribution of the original data; the abscissa is the reads base coordinate, and the ordinate is the percentage of A, T, C, G, and N bases among all reads. (d) Original data quality distribution; the abscissa coordinates are reads base coordinates, and the ordinate is the base mass of reads (Solexa scale: 40 = highest, −15 = lowest).
Full size image
Transcriptome and MicroRNA data annotation
The amounts of transcriptome data in the telogen (T-1, T-2, and T-3) and anagen (A-1, A-2, and A-3) stages were 2.6 Gb, 2.3 Gb, 2.4 Gb, 2.4 Gb, 3.2 Gb and 2.4 Gb, respectively. After filtering out the low-quality sequences, we used Trinity software for de novo splicing to obtain the following amounts of transcriptome data: 198,542 bp, 219,454 bp, 219,411 bp, 247,695 bp, 250,081 bp, 247,695 bp and 203,451 bp. The total number of genes was 55,541, the number of alternative splices was 105,854, the average sequence length was 1965.15 bp, the length of the longest sequence was 23,311 bp, and the length of the shortest sequence was 351 bp.
The genes that control skin development differed significantly in the course of one year (Table 1). The first significant difference in gene number occurred from February to March, and the second occurred from March to April; thus, in March, there was a period of significant change for the genes in the skin, and these changes were relatively independent. The third significant change was from June to July, and the fourth significant change was October to November. After November and until February of the following year, changes in skin genes were almost absent.
Table 1 Statistical data for the differentially expressed genes between each pair of months.
Full size table
From sRNA sequencing of skin tissue samples from the repose and growth stages, we obtained telogen (T-1, T-2, T-3) and anagen data (A-1, A-2, A-3) consisting of 12.2 Mb, 11.0 Mb, 10.4 Mb, 11.3 Mb, 10.2 Mb and 10.9 Mb reads, respectively. After the removal of low-quality reads, we obtained 12.1 Mb, 10.6 Mb, 9.8 Mb, 11.2 Mb, 9.7 Mb and 10.8 Mb reads (Table 2) that could be used for subsequent analysis.
Table 2 Statistical results for sRNA sequencing data.
Full size table
We also performed statistical analysis of the length distribution of the high-quality sRNA reads (Fig. 2b), and the results showed that total reads were primarily 22 nt long, accounting for 31.46% of reads, and that the proportions of reads that were 20, 21, 22, and 23 nt in length among distinct reads were 14.4%, 14.01%, 13.54% and 12.43%, respectively.
Sequencing data processing and analysis of differences. (a) Gene differential expression statistics; the ordinate is the Log_fold_change, the abscissa is the absolute expression level, that is, logConc. (b) SRNA high-quality reads length distribution statistics. (c) Comparison of sequence data with MiRbase database.
Full size image
Comparing our high-quality sequencing data to known data in the MiRbase database (Fig. 2c), we obtained 437,517 known microRNA sequences, accounting for 47.41% of sequences, and 4,852,979 unknown microRNA sequences, accounting for 52.29% of sequences, which included rRNAs, tRNAs and snoRNAs, among others. There were 333 precursor microRNAs and 470 mature body microRNAs, which were classified into 150 families. We carried out a forecast of these new miRNAs and found that there were 92 precursor microRNAs and 99 mature body microRNAs.
Differential analysis of the transcriptome and miRNAs
From a differential analysis of the differentially expressed genes that were related to villus growth in March and July (Table 1), we obtained 12,865 differentially expressed genes in the telogen and anagen stages, of which there were 7,664 upregulated and 5,201 downregulated genes (Fig. 2a). A differential analysis of mature microRNAs expressed during the growth and rest periods revealed 35 microRNA genes that were upregulated in the telogen stage compared to their expression in the anagen stage, and there were 16 genes with more than 100-fold higher expression in the anagen stage. We also found 9 downregulated microRNA genes in the anagen stage relative to the telogen stage, and there were 6 genes with more than 100-fold higher expression in the telogen stage: MiR-148a, MiR-34b, MiR-195, MiR-335, MiR-7g*-1, and MiR-101*-1. Thus, we successfully identified the major microRNAs and their target genes.
Analysis of microRNA-mRNA interactions related to the initiation of cashmere production
Using transcriptome data as the target gene database, we obtained 12,927 differentially expressed microRNAs, among which there were 6,965 positive regulatory relationships and 5,114 negative regulatory relationships between the miRNA and the transcriptome. The minimum number of mRNAs with a target relationship with a single microRNA was 51, and the maximum number of mRNAs was 1682. The average number of target genes for a single microRNA was 783 (Table 3). Therefore, microRNA-mRNA interactions in cashmere goat skin were complex, and regulatory relationships can be expressed as one-to-one, one-to-many and many-to-many.
Table 3 MicroRNA and mRNA correlation analysis results.
Full size table
Mutual analysis of microRNAs and mRNAs
The relationship between microRNA and mRNA is complex: one microRNA can target multiple mRNAs, and vice versa. A target gene can be regulated by multiple microRNAs simultaneously, and the relationship between microRNAs and mRNAs is not always limited to negative regulation. The mechanism by which microRNAs regulate mRNAs is not yet clear. This study explored the expression correlations between microRNAs and mRNAs through linear programming and determined the values of the differential multipliers for microRNAs and mRNAs (Fig. 3). For the first quadrant (upper right corner), there were 1,315 mRNAs and 43 microRNAs with higher expression in the anagen stage. For the second quadrant (upper left corner), there were 736 mRNAs with higher expression and 43 microRNAs with lower expression in the anagen stage. For the third quadrant (lower left corner), there were 449 mRNAs and 13 microRNAs with lower expression in the anagen stage. For the fourth quadrant (lower right corner), there were 890 mRNAs with lower expression and 13 microRNAs with higher expression in the anagen stage. Above all, the points in the second quadrant and the fourth quadrant were negatively associated with miRNA and mRNA. These miRNAs that were negatively related to mRNA will be our next area of study.
Correlation analysis of the positive and negative regulation of differential microRNA-mRNA expression. The abscissa is the log2 value of the multiple difference of target genes, and the longitudinal coordinate is the log2 value of the multiple difference of microRNAs. The red dashed line is 1.5 times the difference of the ascending line. The blue dashed line is 1.5 times the difference of the demarcation line.
Full size image
Functional enrichment analysis of differentially expressed microRNA target genes
We also performed GO annotation and KEGG pathway analysis. As the results show, the GO database (Fig. 4a) included 8,409 differentially expressed mRNAs in the telogen relative to the anagen stage, of which there were 3,205 differentially expressed mRNAs upregulating microRNAs, and the number of differentially expressed mRNAs downregulating microRNAs was 4,723. The number of microRNAs negatively regulating their target mRNAs was 2,025. The KEGG database (Fig. 4b) included 5,432 differentially expressed mRNAs between the telogen and anagen stage; among these, there were 1,322 differentially expressed mRNAs upregulated by microRNAs.
Enrichment analyses of GO orthology functions and KEGG pathways for differentially expressed microRNA target genes. (a) Function enrichment and analysis of target gene GO orthology for differentially expressed MicroRNAs. (b) Enrichment analysis of target gene KEGG pathways for differentially expressed MicroRNAs. The P value is the critical threshold for differential microRNA enrichment; the Rich factor is calculated for S/B, S is the mRNA with microRNA target relations annotated to a certain function or path, and B is used to annotate the mRNA of a function or path, based on the mRNA ratio. In other words, the Rich factor represents the ratio and is a percentage; the dot size represents the mRNA number, that is, S.
Full size image
A collection of the most significant top 8 results for GO and KEGG enrichment is presented in Fig. 5.
The first 8 results obtained for functional enrichment analysis of the differential expression of microRNA target genes. The length of the bar is the significantly enriched P value of the value (10 base); the smaller the P value, the more significant the enrichment, and the longer the length of the bar.
Full size image
The GO annotations of the most concentrated microRNA target genes were as follows, from high to low: integral component of the membrane, DNA binding, nucleus, catalytic activity, nucleotide binding, transcription factor activity, transmembrane transport and the oxidation-reduction process. These processes were related to cell development, indicating that a large number of cell proliferation processes occurred during the initiation of hair follicles.
The results of target gene KEGG annotation from high to low were as follows: insulin signalling pathway, focal adhesion, cancer pathways, Wnt signalling pathway, phagosome, bacterial invasion of epithelial cells, arrhythmogenic right ventricular cardiomyopathy and regulation of the actin cytoskeleton. Similar to the results obtained with KEGG pathway enrichment, we found that the concentration of insulin signalling pathways was the highest; there were 220 genes annotated to this pathway. Although there have been no reports of the role of insulin signalling in hair follicle development to date, the insulin signalling pathway is closely related to cell growth, and its potential involvement in hair follicle development will be further studied. We also analysed the signalling pathways related to hair follicles, including the Jak-STAT signalling pathway, Notch signalling pathway and Wnt signalling pathway, which were previously reported. The corresponding numbers of genes identified by annotation were 123, 348 and 262, respectively. The function enrichment P value of the Wnt signalling pathway was 7.29E-64, which indicated a closer relationship with villus development, as shown by the large number of genes in the signalling pathway in skin tissue exhibiting targeted regulatory relationships with the differentially expressed microRNAs. Therefore, Wnt signalling plays an important role in the development of hair follicles.
mRNA-microRNA network analysis
The telogen and anagen stages in cashmere goat skin showed differential expression of microRNAs and their negatively regulated target genes according to the construction of a data relation network based on GO annotation and network construction based on KEGG relationships (Fig. 6). Genes differentially expressed in the telogen and anagen stages have many functions in cashmere goats (Fig. 6a), comprising a complex regulatory network with microRNAs. With GO annotation, there were 350 upregulated genes and 13 microRNAs clustered in the lower part of the network. Additionally, 233 downregulated genes and 33 microRNAs were clustered above the network. Among GO network relationships, genes related to cashmere growth in the Wnt signalling pathway were FZD6, LEF1, FZD3, WNT5A and TCF7, and related microRNAs were MiR-195, MiR-148a, MiR-4206 and gamma.
Network analysis of mRNA-microRNA regulation. (a) Analysis results for the microRNA-mRNA enriched GO function network. (b) MicroRNA-mRNA KEGG pathway enrichment network. The dark dots are signalling pathways, the small green dots are differentially expressed genes, and the red dots are microRNAs with target regulation relationships.
Full size image
The KEGG pathway network diagram was similar to the GO annotation network diagram (Fig. 6b). The genes were clearly divided into two groups, with 184 upregulated genes and 13 microRNAs clustered in the lower part of the network. Fourteen genes were annotated in the Wnt signalling pathway, including calcineurin B homologous protein 1 (CHP1), E3 ubiquitin-protein ligase (SIAH1), FZD6, WNT5A and mothers against decapentaplegic homolog 2 (SMAD2), as well as MiR-195, MiR-335star and other microRNA genes. One hundred twenty-eight downregulated genes and 31 microRNAs were clustered above the network. The genes and microRNAs annotated in the Wnt signalling pathway were 2A5G, CUL1, LRP5, FOSL1, MiR-17astar, MiR-2285m, MiR-136 and MiR-34astar. Therefore, there is a complex signalling pathway between the genes expressed during rest and the genes downregulated during the growth stage.
Up- and downregulated target genes as well as the negative regulatory effects of microRNAs during the period of transition from the telogen to the anagen stage were clearly clustered into two groups. The brown dots in the network map indicate the functional annotation of target genes. Some functions in the map were annotated only for upregulated genes, some functions were annotated only for downregulated genes, and others had functions in both up- and downregulated transcription. This finding indicates that the thresholds for gene balance or gene expression and for differentially expressed genes were relatively important for a specific function. The interactions between microRNAs and their target genes played roles in the initiation and growth of cashmere. This finding further indicates that the initiation of hair follicle growth was controlled by multiple genes and that microRNAs play a role in the post-transcriptional negative regulation of target genes in hair follicle initiation.
Construction of a control network for the Wnt signalling pathway
Multiple genes and microRNAs were involved in the Wnt signalling pathway during goat hair follicle development (Fig. 7a). In addition, the genes and microRNAs in the Wnt signalling pathway also played roles in other related pathways and interacted with other genes. In this study, we used a “biogrid” database to download the genes that interacted with target genes and then merged the genes from the Wnt signalling pathway to construct the Wnt signalling network (Fig. 7b).
Wnt signal path annotation results and construction of the regulatory network. (a) Wnt signal pathway51,52,53 (http://www.kegg.jp/kegg/kegg1.html.). Rectangular nodes represent gene products (such as enzymes or some RNA regulator), gene products all belong to the KEGG ONTOLOGY classification system indicated with a blue background (KO) (some have highly similar sequences, and the same pathway proteins with similar functions were classified as a group KO), and gene products on the white background are not part of the KO classification system. Sequencing of the genes annotated in red revealed these non-KO gene products (nodes indicate gene products with the same or similar functions); circular nodes represent compounds (i.e., a substrate or product); a white background with a rounded rectangle indicates another pathway. The arrow shows the direction of the enzyme reaction or the direction of information transfer, the solid line indicates a direct solution, and the dotted line indicates an indirect solution. (b) Wnt signalling pathway regulation network.
Full size image
Using the GO regulatory network, KEGG regulatory network and the Wnt signalling pathway regulatory gene interaction network in the goat hair follicle resting phase, we screened for genes with upregulated expression in the Wnt signalling pathway from the telogen to the anagen stage, including CHP1, FZD6, SIAH1 and other genes. Keratin and keratin-associated proteins were constitutive proteins during cashmere production. MiR-195 was predicted using target gene prediction and found to negatively regulate KAP24-1, KAP3-1, KAT8 and the important regulatory factor EGFR, which constitute the main components of cashmere. The genes corresponding to KAP24-1, KAP3-1, KAT8 and EGFR were CHP1, FZD6, SIAH1 and SMAD2. Therefore, we chose CHP1, FZD6, SIAH1 and FZD6 as the regulatory genes during the initiation of villus growth for further verification.
Expression of CHP1, FZD6, SIAH1, and SMAD2
The expression of CHP1, FZD6, SIAH, and SMAD2 was verified via fluorescence quantitative PCR, using β-actin as a reference gene. As shown in Fig. 8, compared to the expression of microRNA-195 (Fig. 8a), the expression level of the CHP1 gene (Fig. 8d) was higher in the anagen stage than in the telogen stage, and the difference was relatively significant (P < 0.01), with a descending trend from the early growth stage to the late growth stage that was not significant (P > 0.05). The expression level of the Frizzled-6 (FZD6) gene (Fig. 8e) increased from the early stage to the telogen stage, but the difference was not significant (P > 0.05). However, the expression level increased significantly from the early growth stage to the late growth stage (P < 0.01). The expression of the gene SIAH1 (Fig. 8c) increased significantly from the telogen to the early growth stage (P < 0.01), but there was no significant difference from the early growth stage to the late growth stage (P > 0.05). The SMAD2 gene (Fig. 8b) showed no significant difference from the telogen to the anagen stage (P > 0.05), but its expression increased from the telogen to the late growth stage. In the Wnt signalling network, two of the three key genes were significantly expressed from the telogen to the anagen stage.
Fluorescence quantitative expression results. **Indicates the difference between representatives is very significant, and * indicates the difference is significant. (a) Relative expression of microRNA-195. (b) Relative expression of SMAD2. (c) Relative expression of SIAH1. (d) Relative expression of CHP1. (e) Relative expression of FZD6.
Full size image
Farazi T A et al. showed that a computer-aided algorithm could be used to integrate miRNA and mRNA paired expression spectra to characterize miRNA-mRNA interactions and that this approach identified miRNA targets with high accuracy. This approach effectively avoids certain software predictions that do not consider specific cases of miRNA and mRNA expression24,25,26. Chen et al. published an article describing a model for predicting disease associations by mixing microRNAs into computations (HAMDA), which significantly reduced experimental time and cost27,28,29. With HAMDA, the functions of microRNAs can be correlated with disease characteristics, and the role of new microRNAs in disease can be revealed more clearly. Additionally, the authors proposed a method of using Laplacian Regularized Sparse Subspace Learning for MicroRNA-Disease Association (LRSSLMDA), which projects a microRNA and a disease into a common subspace to interpret the function of the microRNA subspace27,30,31,32. Based on these methods and on a joint analysis of miRNA-mRNA, we conducted a targeted prediction of miRNA-mRNA interactions and used reverse locations of mRNAs to predict miRNA target genes. After obtaining positive target genes, we used RT-qPCR to further verify the target genes and determine the interactions between target genes regulated by miRNAs.
In this study, target genes of differentially expressed microRNAs were analysed using GO orthology enrichment analysis and KEGG pathway enrichment33, as described in Chen L’s study. The most commonly annotated GO function was integral component of the membrane within cellular components, indicating that microRNAs can regulate villus initiation and the villus growth cycle by regulating cell membrane-associated gene expression. The KEGG pathway results showed significant enrichment of the Jak-STAT signalling pathway, Notch signalling pathway, p53 signalling pathway and Wnt signalling pathway. These cashmere growth-related pathways were found via analysis of differentially expressed microRNA target genes in the KEGG pathway map, and the Wnt signalling pathway was found to play a particularly important role. Therefore, this study constructed the Wnt signalling pathway regulatory network during hair follicle initiation.
The American sociologist Granovetter first posited the concept of relationship intensity and divided relationships into strong and weak relationships34, with strong ties maintaining the relationships between groups and organizations and weak ties between groups and organizations. Information obtained from strong relationships is often highly repeatable, while weak relationships provide more information and other resources than strong relationships. Based on the ideas provided by social networks, we hypothesize that strong and weak relationships exist between genes. A gene associated with more genes in a network is probably an important gene or key node gene that controls or influences a trait35. Thus far, efforts such as pathway analysis have raised awareness of the functional contributions of gene mutations and DNA copy number variations to cancer development, progression and metastasis36. Wang E commented on the challenges associated with studying cancer omic data using an integrative network approach and suggested possible research directions. As a result, most current cancer genome sequencing work has mapped mutations onto biological pathways, including signalling pathways37,38. Mcgee S R indicated that by studying the distributions of such network motifs, insight into cancer signalling regulatory molecular mechanisms of tumorigenesis can be gained, and the identification of these loops has practical implications for predicting prognosis and the clinical outcomes of cancer patients. Collectively, network motifs and modules are critical for cancer signalling and are associated with clinical outcomes39. Additionally, as stated by Han P, analysis of GRNs can identify regional subnetworks for certain biological processes; in-cloud regulatory structures between genes, key regulators, and cancer hallmark subnetworks; network dynamics for network rewiring; and network motifs. Furthermore, such results may reveal the molecular mechanisms of signalling pathways associated with cancer hallmarks and cancer patient outcomes38,40,41,42. We found a complex network between genes in the Wnt signalling pathway and other genes and signalling pathways during hair follicle development from rest to growth, indicating that hair follicle initiation is determined not by single genes but by differences in the number of interactions between genes. One or more key genes are present in this gene regulatory network. In the Wnt signalling pathway, gene interactions are centred on the CHP1, SIAH1, and SMAD2 genes; however, many inter-gene interactions and multiple signalling pathways are involved.
As indicated by previous studies, microRNAs have negative regulatory roles43,44,45,46. In this study, quantitative analysis of MiR-195 expression during each period of chorionic villus initiation was performed. MiR-195 showed a trend of declining expression during the periods of repose, prophase and growth. In our experimental results, we found that two genes, SMAD2 and FZD6, were negatively correlated during the periods of repose, prophase and growth. Therefore, we speculate that MiR-195 continues to inhibit the expression of SMAD2 and FZD6 during the entire initiation process. At present, no information is available regarding the regulation of SMAD2 and FZD6 targeting in by-195. However, a study by Mo J showed that MiR-195 regulates cell proliferation (Mo J. et al.47). Zheng R., Du J. et al. showed that downregulated MiR-195 was targeted to promote the expression of SMAD7 in two-leaf aortic valve calcification48,49. MiR-195 may play an important role in the regulation of SMAD family genes. However, two genes, SIAH1 and CHP1, showed a trend that first increased and then decreased. Presumably, MiR-195-mediated inhibition of these two genes after the initiation of cashmere production is affected by other microRNAs and is relieved or reduced. A study by Zhang X showed that MiR-195 affects the proliferation of colon cancer cells and regulation of Wnt/β-catenin pathway protein-specific MiR-195 by targeting FGF2 and regulation of the Wnt/β-catenin signalling pathway, consistent with studies showing the effects of MicroRNA-195 on Wnt signalling47,50.
Materials and Methods
The study samples were collected from an Inner Mongolian Arbas white cashmere goat breeding farm. All animal experiments were performed in accordance with the ‘Guidelines for Experimental Animals’ of the Ministry of Science and Technology (Beijing, China). All surgeries were performed according to recommendations proposed by the European Commission (1997) and were approved by the experimental animal ethics committee of Inner Mongolia Agricultural University. Six cashmere goats were selected from the same growth environment and had equal body weights, unrelated relationships, were of the same sex, and had good growth conditions. The sampling times occurred during the telogen and anagen stages. The sampling position was the side body or the middle of the body at a distance of 10–15 cm from the scapula, and the sample size was 2 cm2. The samples were stored at −80 °C in a freezer.
Extraction of total RNA with liquid nitrogen lapping and RNAiso Plus (Takara, China), with separation on a Bioanalyzer 2100 (Agilent, USA), was used to detect RNA and microRNA quality. A cDNA library was constructed for the transcriptional group using a TruSeqTM RNA sample preparation kit (Illumina, USA) according to the manufacturer’s instructions. A microRNA sequencing library was constructed according to the instructions of the Illumina (USA) reference kit (USA). The transcriptome was constructed, and Solexa sequencing of the microRNA library was performed by Shanghai M-image Biological Medicine Science and Technology Co., Ltd.
Original data processing and differential expression analysis
Seqprep, sickle, and condetri_v2.0.pl software were used to evaluate the original data for transcriptional sequencing, and Trinity software was used for de novo splicing and open reading frame (ORF) prediction. We predicted ORF sequences using BLAST (BLAST Version 2.2.25, with an E-value parameter value setting of less than 10−5) in the NCBI NR database with BLASTP database identification of the string database and KEGG database without predicted sequences or NR, string, or gene library ORF BLASTX comparison notes. The expression of each gene in each sample was counted based on a comparison between the sequencing sample and the reference genome, and gene expression levels were calculated using the maximum likelihood method in RSEM software. EdgeR software was used to find significant differences in the expression of all transcripts in the sample at a threshold of P = 0.05.
Short oligonucleotide alignment program (SOAP) software was used to locate sRNAs in the genome. “Tag2repeat” software was used for repeated sRNA comparisons. Using overlap software, sRNAs were compared to the exons and introns of mRNAs, and sRNAs derived from mRNA degradation were found. The GenBank and Rfam (9.1) databases (database link) were used to annotate the data and remove rRNAs, scRNAs, snoRNAs, snRNAs, and tRNAs from the sRNAs as much as possible. The sequences were compared to microRNA precursors and the mature bodies of cattle and sheep in MiRBase21 (http://www.MiRbase.org/) to obtain known microRNAs. We analysed whether there were significant differences in the expression of known microRNAs (P = 0.05) and compared microRNA expression patterns using a log 2-ratio and scatter plot.
Association analysisTarget gene prediction for the initiation of related microRNAs
Targetscan (parameter: context score percentile of ≥90) and MiRanda (parameter: max energy of ≤−20) software were used to predict the differentially expressed microRNA target genes in the telogen and anagen stages, and the target regulatory relationships between microRNAs and mRNAs were determined based on sequence complementarity, sequence conservatism, thermo-kinetic factors, site-binding ability and UTR base distribution. The relationships between microRNAs and mRNAs are complex. A microRNA can target multiple mRNAs simultaneously. In contrast, a target gene mRNA can be regulated by multiple microRNAs simultaneously. We used ACGT101-CORR (1.1 version) software to extract and perform statistical analysis for data describing target regulation relationships, specifically positive and negative relationships between the expression of microRNAs and their target genes expressed as the difference multiplier of microRNA and the difference multiplier of mRNA (R correlation analysis).
Analysis of genes related to chorionic cycle initiation
This study used GO enrichment and KEGG functions for gene analysis. The target functions of all genes or annotations among the differentially expressed miRNAs were counted, and hypergeometric tests were used to determine the functions or pathways that were significantly enriched among differentially expressed miRNA-mRNA relationships. The formula for calculating P value significance was
If the microRNA target gene of the functional annotation satisfied this condition, we defined the result as distinct, significant expression. In the formula, TB is the total number of mRNAs with functional annotations or path annotations, TS is the number of mRNAs corresponding to the differential expression of miRNA in the TB, B is the number of mRNAs annotated for a particular single function or a particular path, and S is the number of mRNAs annotated for a specific single function or a specific miRNA in a particular path. According to the first 8 GO functions and KEGG pathways with the smallest P values, we sorted the results and performed an enrichment analysis.
Among KEGG pathways, the pathway that is most closely related to hair follicle development is the Wnt signalling pathway. We separated microRNAs from their target genes in the Wnt pathway using the biogrid database and downloaded the network interaction, microRNAs, target genes, Wnt genes and signalling pathway interactions, which were distinguished by various colours, characteristics and attributes such as size, using Cytoscape software to obtain the functions in the Wnt pathway network analysis chart.
Analysis of main genes and the expression of microRNAs acting on the chorionic cycle
Real-time fluorescent quantitative PCR was performed using the FAST SYBR Green Master Mix Kit (Applied Biosystems, USA). Fluorescent quantitative primers were designed with primer5 (Table 4). PCR reactions were performed using an Agilent 3000XP Real Time PCR amplification instrument; the PCR reaction program was 95 °C for 20 sec, 1 cycle; 95 °C for 5 sec; 60 °C for 30 sec, 40 cycles; 95 °C for 1 min, 55 °C for 30 sec; and 95 °C for 30 sec, 1 cycle. Each sample was tested 3 times, and 3 blank controls were used for each primer. We used the 2−ΔΔCt method for relative quantitation between samples, and the reference gene was β-actin. Using the statistical analysis software SPSS 17.0, we tested differences in the relative expression of keratin-associated protein genes in Arbas cashmere goat skin in different months using two methods: LSD and Duncan’s test. The results were expressed as the average value ± standard deviation.
Table 4 Primers used for PCR validation.
Full size table
Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Li, H., Ruan, J. & Durbin, R. Mapping short DNA sequencing reads and calling variants using mapping quality scores. Genome research 18, 1851–1858, https://doi.org/10.1101/gr.078212.108 (2008).
CASArticlePubMedPubMed Central Google Scholar2.
Ingolia, N. T., Ghaemmaghami, S., Newman, J. R. & Weissman, J. S. Genome-wide analysis in vivo of translation with nucleotide resolution using ribosome profiling. Science 324, 218–223, https://doi.org/10.1126/science.1168978 (2009).
ADSCASArticlePubMedPubMed Central Google Scholar3.
Cloonan, N. et al. Stem cell transcriptome profiling via massive-scale mRNA sequencing. Nature methods 5, 613–619, https://doi.org/10.1038/nmeth.1223 (2008).
CASArticlePubMed Google Scholar4.
Paus, R. & Foitzik, K. In search of the “hair cycle clock”: a guided tour. Differentiation; research in biological diversity 72, 489–511, https://doi.org/10.1111/j.1432-0436.2004.07209004.x (2004).
CASArticlePubMed Google Scholar5.
Costa, V., Angelini, C., De Feis, I. & Ciccodicola, A. Uncovering the complexity of transcriptomes with RNA-Seq. Journal of biomedicine & biotechnology 2010, 853916, https://doi.org/10.1155/2010/853916 (2010).
CASArticle Google Scholar6.
Margulies, M. et al. Genome sequencing in microfabricated high-density picolitre reactors. Nature 437, 376–380, https://doi.org/10.1038/nature03959 (2005).
ADSArticlePubMedPubMed Central Google Scholar7.
Zhang, Y. Y., Zan, L. S. & Wang, H. B. Genome array on differentially expressed genes of muscle tissue in intact male and castrated Qinchuan cattle. Yi chuan = Hereditas 32, 1166–1174 (2010).
CASPubMed Google Scholar8.
Mardis, E. R. Next-generation DNA sequencing methods. Annual review of genomics and human genetics 9, 387–402, https://doi.org/10.1146/annurev.genom.9.081307.164359 (2008).
CASArticlePubMed Google Scholar9.
Bentley, D. R. et al. Accurate whole human genome sequencing using reversible terminator chemistry. Nature 456, 53–59, https://doi.org/10.1038/nature07517 (2008).
ADSCASArticlePubMedPubMed Central Google Scholar10.
Hudson, M. B. Assays for micro-rna-182 as a biomarker for muscle atrophy and therapeutic applications. (2014).
Mardis, E. R. A decade’s perspective on DNA sequencing technology. Nature 470, 198–203, https://doi.org/10.1038/nature09796 (2011).
ADSCASArticlePubMed Google Scholar12.
Ozsolak, F. & Milos, P. M. RNA sequencing: advances, challenges and opportunities. Nature reviews. Genetics 12, 87–98, https://doi.org/10.1038/nrg2934 (2011).
CASArticlePubMed Google Scholar13.
Maher, C. A. et al. Chimeric transcript discovery by paired-end transcriptome sequencing. Proceedings of the National Academy of Sciences of the United States of America 106, 12353–12358, https://doi.org/10.1073/pnas.0904720106 (2009).
ADSArticlePubMedPubMed Central Google Scholar14.
Rana, T. M. Illuminating the silence: understanding the structure and function of small RNAs. Nature reviews. Molecular cell biology 8, 23–36, https://doi.org/10.1038/nrm2085 (2007).
CASArticlePubMed Google Scholar15.
Chen, H. L. et al. Galectin-7 Regulates Keratinocyte Proliferation and Differentiation through JNK-miR-203-p63 Signaling. The Journal of investigative dermatology 136, 182–191, https://doi.org/10.1038/JID.2015.366 (2016).
CASArticlePubMedPubMed Central Google Scholar16.
Shen, F. et al. Genetic variants in miR-196a2 and miR-499 are associated with susceptibility to esophageal squamous cell carcinoma in Chinese Han population. Tumour biology: the journal of the International Society for Oncodevelopmental Biology and Medicine 37, 4777–4784, https://doi.org/10.1007/s13277-015-4268-3 (2016).
CASArticle Google Scholar17.
Au, K. F., Jiang, H., Lin, L., Xing, Y. & Wong, W. H. Detection of splice junctions from paired-end RNA-seq data by SpliceMap. Nucleic acids research 38, 4570–4578, https://doi.org/10.1093/nar/gkq211 (2010).
CASArticlePubMedPubMed Central Google Scholar18.
Weilner, S. et al. Secreted microvesicular miR-31 inhibits osteogenic differentiation of mesenchymal stem cells. Aging cell 15, 744–754, https://doi.org/10.1111/acel.12484 (2016).
CASArticlePubMedPubMed Central Google Scholar19.
Edgren, H. et al. Identification of fusion genes in breast cancer by paired-end RNA-sequencing. Genome biology 12, R6, https://doi.org/10.1186/gb-2011-12-1-r6 (2011).
CASArticlePubMedPubMed Central Google Scholar20.
Fu, S., Zhao, H., Zheng, Z., Li, J. & Zhang, W. Melatonin regulating the expression of microRNAs involved in hair follicle cycle of cashmere goats skin. Yi chuan = Hereditas 36, 1235–1242, https://doi.org/10.3724/SP.J.1005.2014.1235 (2014).
CASArticlePubMed Google Scholar21.
Jang, I. et al. miRseqViewer: multi-panel visualization of sequence, structure and expression for analysis of microRNA sequencing data. Bioinformatics 31, 596–598, https://doi.org/10.1093/bioinformatics/btu676 (2015).
CASArticlePubMed Google Scholar22.
Liu, Z. et al. Identification of conserved and novel microRNAs in cashmere goat skin by deep sequencing. PloS one 7, e50001, https://doi.org/10.1371/journal.pone.0050001 (2012).
ADSCASArticlePubMedPubMed Central Google Scholar23.
Song, M. K. et al. Analysis of microRNA and mRNA expression profiles highlights alterations in modulation of the MAPK pathway under octanal exposure. Environmental toxicology and pharmacology 37, 84–94, https://doi.org/10.1016/j.etap.2013.11.005 (2014).
CASArticlePubMed Google Scholar24.
Farazi, T. A., Spitzer, J. I., Morozov, P. & Tuschl, T. microRNAs in human cancer. The Journal of pathology 223, 102–115, https://doi.org/10.1002/path.2806 (2011).
CASArticlePubMed Google Scholar25.
Zhao, L. et al. Integration analysis of microRNA and mRNA paired expression profiling identifies deregulated microRNA-transcription factor-gene regulatory networks in ovarian endometriosis. Reproductive biology and endocrinology: RB&E 16, 4, https://doi.org/10.1186/s12958-017-0319-5 (2018).
Article Google Scholar26.
Mavrakis, K. J. & Wendel, H. G. TargetScreen: an unbiased approach to identify functionally important microRNA targets. Cell cycle 9, 2080–2084, https://doi.org/10.4161/cc.9.11.11807 (2010).
CASArticlePubMed Google Scholar27.
Chen, X. et al. WBSMDA: Within and Between Score for MicroRNA-Disease Association prediction. Scientific reports 6, 21106, https://doi.org/10.1038/srep21106 (2016).
ADSCASArticlePubMedPubMed Central Google Scholar28.
You, Z. H. et al. PBMDA: A novel and effective path-based computational model for microRNA-disease association prediction. PLoS computational biology 13, e1005455, https://doi.org/10.1371/journal.pcbi.1005455 (2017).
CASArticlePubMedPubMed Central Google Scholar29.
Chen, X., Yan, C. C., Zhang, X. & You, Z. H. Long non-coding RNAs and complex diseases: from experimental results to computational models. Briefings in bioinformatics 18, 558–576, https://doi.org/10.1093/bib/bbw060 (2017).
ArticlePubMed Google Scholar30.
Chen, X. & Huang, L. LRSSLMDA: Laplacian Regularized Sparse Subspace Learning for MicroRNA-Disease Association prediction. PLoS computational biology 13, e1005912, https://doi.org/10.1371/journal.pcbi.1005912 (2017).
ADSCASArticlePubMedPubMed Central Google Scholar31.
Chen, X., Huang, L., Xie, D. & Zhao, Q. EGBMMDA: Extreme Gradient Boosting Machine for MicroRNA-Disease Association prediction. Cell death & disease 9, 3, https://doi.org/10.1038/s41419-017-0003-x (2018).
CASArticle Google Scholar32.
Chen, X., Niu, Y. W., Wang, G. H. & Yan, G. Y. HAMDA: Hybrid Approach for MicroRNA-Disease Association prediction. Journal of biomedical informatics 76, 50–58, https://doi.org/10.1016/j.jbi.2017.10.014 (2017).
CASArticlePubMed Google Scholar33.
Chen, L. et al. Gene Ontology and KEGG Pathway Enrichment Analysis of a Drug Target-Based Classification System. PloS one 10, e0126492, https://doi.org/10.1371/journal.pone.0126492 (2015).
CASArticlePubMedPubMed Central Google Scholar34.
Carrasco-Garcia, E. et al. Paradoxical role of SOX2 in gastric cancer. American journal of cancer research 6, 701–713 (2016).
CASPubMedPubMed Central Google Scholar35.
Ma, L., Ballantyne, C., Brautbar, A. & Keinan, A. Analysis of multiple association studies provides evidence of an expression QTL hub in gene-gene interaction network affecting HDL cholesterol levels. PloS one 9, e92469, https://doi.org/10.1371/journal.pone.0092469 (2014).
ADSCASArticlePubMedPubMed Central Google Scholar36.
Li, L. et al. The human phosphotyrosine signaling network: evolution and hotspots of hijacking in cancer. Genome research 22, 1222–1230, https://doi.org/10.1101/gr.128819.111 (2012).
CASArticlePubMedPubMed Central Google Scholar37.
Wang, E. Understanding genomic alterations in cancer genomes using an integrative network approach. Cancer letters 340, 261–269, https://doi.org/10.1016/j.canlet.2012.11.050 (2013).
CASArticlePubMed Google Scholar38.
Wang, E. et al. Cancer systems biology in the genome sequencing era: part 2, evolutionary dynamics of tumor clonal networks and drug resistance. Seminars in cancer biology 23, 286–292, https://doi.org/10.1016/j.semcancer.2013.06.001 (2013).
CASArticlePubMed Google Scholar39.
McGee, S. R., Tibiche, C., Trifiro, M. & Wang, E. Network Analysis Reveals A Signaling Regulatory Loop in the PIK3CA-mutated Breast Cancer Predicting Survival Outcome. Genomics, proteomics & bioinformatics 15, 121–129, https://doi.org/10.1016/j.gpb.2017.02.002 (2017).
Article Google Scholar40.
Han, P., Gopalakrishnan, C., Yu, H. & Wang, E. Gene Regulatory Network Rewiring in the Immune Cells Associated with Cancer. Genes 8, https://doi.org/10.3390/genes8110308 (2017).
Article Google Scholar41.
Cloutier, M. & Wang, E. Dynamic modeling and analysis of cancer cellular network motifs. Integrative biology: quantitative biosciences from nano to macro 3, 724–732, https://doi.org/10.1039/c0ib00145g (2011).
CASArticle Google Scholar42.
Wang, E. et al. Cancer systems biology in the genome sequencing era: part 1, dissecting and modeling of tumor clones and their networks. Seminars in cancer biology 23, 279–285, https://doi.org/10.1016/j.semcancer.2013.06.002 (2013).
CASArticlePubMed Google Scholar43.
Yang, Y. et al. MicroRNA-195 acts as a tumor suppressor by directly targeting Wnt3a in HepG2 hepatocellular carcinoma cells. Molecular medicine reports 10, 2643–2648, https://doi.org/10.3892/mmr.2014.2526 (2014).
CASArticlePubMed Google Scholar44.
Xie, Z. R., Yang, H. T., Liu, W. C. & Hwang, M. J. The role of microRNA in the delayed negative feedback regulation of gene expression. Biochemical and biophysical research communications 358, 722–726, https://doi.org/10.1016/j.bbrc.2007.04.207 (2007).
CASArticlePubMed Google Scholar45.
Tang, W. F. et al. Host MicroRNA miR-197 Plays a Negative Regulatory Role in the Enterovirus 71 Infectious Cycle by Targeting the RAN Protein. Journal of virology 90, 1424–1438, https://doi.org/10.1128/JVI.02143-15 (2016).
CASArticlePubMedPubMed Central Google Scholar46.
Xiao, B. et al. Induction of microRNA-155 during Helicobacter pylori infection and its negative regulatory role in the inflammatory response. The Journal of infectious diseases 200, 916–925, https://doi.org/10.1086/605443 (2009).
CASArticlePubMed Google Scholar47.
Mo, J., Zhang, D. & Yang, R. MicroRNA-195 regulates proliferation, migration, angiogenesis and autophagy of endothelial progenitor cells by targeting GABARAPL1. Bioscience reports 36, https://doi.org/10.1042/BSR20160139 (2016).
Article Google Scholar48.
Du, J. et al. Downregulated MicroRNA-195 in the Bicuspid Aortic Valve Promotes Calcification of Valve Interstitial Cells via Targeting SMAD7. Cellular physiology and biochemistry: international journal of experimental cellular physiology, biochemistry, and pharmacology 44, 884–896, https://doi.org/10.1159/000485356 (2017).
CASArticle Google Scholar49.
Song, L. Y. et al. MicroRNA-195 Activates Hepatic Stellate Cells In Vitro by Targeting Smad7. BioMed research international 2017, 1945631, https://doi.org/10.1155/2017/1945631 (2017).
ArticlePubMedPubMed Central Google Scholar50.
Zhang, X. et al. MicroRNA-195 suppresses colorectal cancer cells proliferation via targeting FGF2 and regulating Wnt/beta-catenin pathway. American journal of cancer research 6, 2631–2640 (2016).
CASPubMedPubMed Central Google Scholar51.
Kanehisa, M., Furumichi, M., Tanabe, M., Sato, Y. & Morishima, K. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic acids research 45, D353–D361, https://doi.org/10.1093/nar/gkw1092 (2017).
CASArticlePubMed Google Scholar52.
Kanehisa, M. & Goto, S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic acids research 28, 27–30 (2000).
CASArticle Google Scholar53.
Kanehisa, M., Sato, Y., Kawashima, M., Furumichi, M. & Tanabe, M. KEGG as a reference resource for gene and protein annotation. Nucleic acids research 44, D457–462, https://doi.org/10.1093/nar/gkv1070 (2016).
CASArticlePubMed Google Scholar
The authors are thankful for the samples provided by the Aerbasi White Cashmere Goat Breeding Farm. Jinquan Li provided the test platform, whereas Zhihong Liu helped in designing and conducting the experiments and the analysis, evaluation, and interpretation of the results. This work has been funded under the National Natural Science Foundation of China (#31360537, #31660640), the National “863” Project (#2013AA102506), the Natural Science Foundation of Inner Mongolia (#2017MS0356), and the Major Projects of the Inner Mongolia Natural Science Foundation (#2016ZD02). The Major Projects of the Inner Mongolia Natural Science Foundation (#2016ZD02) and the Natural Science Foundation of Inner Mongolia (#2017MS0356) provided support for the design of the study and data collection.
Zhihong Liu, Feng Yang, Meng Zhao, Lina Ma and Haijun Li contributed equally.
AffiliationsCollege of Animal Science, Inner Mongolia Agricultural University, Hohhot, 010018, ChinaZhihong Liu, Feng Yang, Meng Zhao, Lina Ma, Yuchun Xie, Rile Nai, Tianyu Che, Rui Su, Yanjun Zhang, Ruijun Wang, Zhiying Wang & Jinquan LiKey Laboratory of Animal Genetics, Breeding and Reproduction, Inner Mongolia Autonomous Region, Hohhot, ChinaZhihong Liu, Feng Yang, Lina Ma, Yuchun Xie, Rui Su, Yanjun Zhang, Ruijun Wang, Zhiying Wang & Jinquan LiKey Laboratory of Mutton Sheep Genetics and Breeding, Ministry of Agriculture, Hohhot, ChinaZhihong Liu, Rile Nai & Jinquan LiEngineering Research Center for Goat Genetics and Breeding, Inner Mongolia Autonomous Region, Hohhot, ChinaMeng Zhao & Tianyu CheCollege of Veterinary Medicine, Inner Mongolia Agricultural University, Hohhot, 010018, ChinaHaijun Li
AuthorsSearch for Zhihong Liu in:Nature Research journals • PubMed • Google Scholar Search for Feng Yang in:Nature Research journals • PubMed • Google Scholar Search for Meng Zhao in:Nature Research journals • PubMed • Google Scholar Search for Lina Ma in:Nature Research journals • PubMed • Google Scholar Search for Haijun Li in:Nature Research journals • PubMed • Google Scholar Search for Yuchun Xie in:Nature Research journals • PubMed • Google Scholar Search for Rile Nai in:Nature Research journals • PubMed • Google Scholar Search for Tianyu Che in:Nature Research journals • PubMed • Google Scholar Search for Rui Su in:Nature Research journals • PubMed • Google Scholar Search for Yanjun Zhang in:Nature Research journals • PubMed • Google Scholar Search for Ruijun Wang in:Nature Research journals • PubMed • Google Scholar Search for Zhiying Wang in:Nature Research journals • PubMed • Google Scholar Search for Jinquan Li in:Nature Research journals • PubMed • Google Scholar
Z.L., N.R., F.Y. and J.L. made substantial contributions to the conception and design of the experiments. Conceived and designed the experiments: LM, MZ, YX, RW. Performed the experiments: Z.L., J.Q., M.Z., F.Y. Analysed the data: L.M., Y.Z., M.Z. Wrote the paper: R.S., F.Y., H.L. Critically revised the manuscript: F.Y., Z.L., H.L. All authors read and approved the final manuscript.
The authors declare no competing interests.
Correspondence to Jinquan Li.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
28 March 2018
23 August 2018
21 September 2018