Identification and Characterization of CircRNAs in Hair Follicles of Cashmere Goats
In order to explore the expression pattern of circRNA in fetal skin hair follicles of cashmere goats,in this study, high-throughput sequencing was performed at four stages of the fetal phase of Inner Mongolian cashmere goats (Alba type).First,we constructed 12 libraries that removed ribosomal RNA of cashmere goats during the fetal period, which were named D45-1, D45-2, D45-3, D55-1, D55-2, D55-3, D65-1, D65-2, D65-3, D75-1, D75-2 and D75-3,respectively. These libraries were applied to the IlluminaHiseq4000 platform for RNA sepuencing and 1063299566 raw reads were obtained from the 12 libraries(Table 1). In these original reads, 1023889360 effective reads were obtained by removing the reads with connector (Adaptor),the reads with N (N indicating that the base information cannot be determined) is greater than 5%, and the low-quality (the bases with a quality value Q ≤ 10 accounted for more than 20% of the total read) (Table 2). The Q20 (the proportion of bases with quality value ≥ 30, error rate < 0.001) of each library was 99.90% and the Q30 (the proportion of bases with quality value ≥ 20, error rate < 0.001) was above 98.1% (Table 1), indicating that the sequencing accuracy was high. Comparing the 1023889360 effective reads with the reference genome, the percentage of the number of reads on the reference genome as a percentage of valid reads was more than 94%, the percentage of the number of reads compared to the unique location of the reference genome as a percentage of valid reads was more than 77%, and the number of reads compared to the multiple locations of the reference genome as a percentage of valid reads was more than 17%. Therefore, the utilization rate of the data was normal, and the original data obtained met the requirements of the subsequent circRNA analysis (Table 2) in terms of quantity and quality. According to the characteristics of the circRNA structure and splicing sequence, we used CIRCExploter2 and CIRI to predict circRNAs. The results showed that the total number of circRNA identified in each library was more than 2800,and corresponding parental genes numbered more than 1700 (Fig. 1). Studies have found that most circRNAs contain two to four exons (Fig. 2a). At the same time, the exon length of circRNAs includes that the length of a single exon is longer than that of circRNAs composed of multiple exons (Fig. 2b). Chromosome distribution analysis showed that circRNAs were distributed on almost all chromosomes, and the number of circRNAs on chromosomes 1, 10 and 11 was higher than that on other chromosomes (Fig. 2c). Finally, exon circRNAs accounted for 92.01% (Fig. 2d) of all circRNAs in cashmere goat skin hair follicles.
Table 1
Data quality control statistics
Sample
|
Raw Data
|
|
Valid Data
|
|
ValidRatio(reads)
|
Q20%
|
Q30%
|
GC content%
|
|
Read
|
Base
|
Read
|
Base
|
|
|
|
|
d45_1
|
83377840
|
12.51G
|
80991260
|
12.15G
|
97.14
|
99.97
|
98.33
|
46
|
d45_2
|
80073650
|
12.01G
|
77655516
|
11.65G
|
96.98
|
99.97
|
98.33
|
45.50
|
d45_3
|
91499974
|
13.72G
|
87070648
|
13.06G
|
95.16
|
99.98
|
98.41
|
46
|
d55_1
|
100354248
|
15.05G
|
96199486
|
14.43G
|
95.86
|
99.98
|
98.45
|
46.50
|
d55_2
|
91779488
|
13.77G
|
88089256
|
13.21G
|
95.98
|
99.98
|
98.53
|
47
|
d55_3
|
94888486
|
14.23G
|
91438886
|
13.72G
|
96.36
|
99.98
|
98.48
|
47
|
d65_1
|
90684330
|
13.60G
|
87483330
|
13.12G
|
96.47
|
99.97
|
98.17
|
47
|
d65_2
|
102708880
|
15.41G
|
99047498
|
14.86G
|
96.44
|
99.98
|
98.47
|
46.50
|
d65_3
|
80732910
|
12.11G
|
77823688
|
11.67G
|
96.40
|
99.97
|
98.34
|
47
|
d75_1
|
81140600
|
12.17G
|
77983860
|
11.70G
|
96.11
|
99.98
|
98.45
|
45
|
d75_2
|
83141742
|
12.47G
|
80115068
|
12.02G
|
96.36
|
99.98
|
98.46
|
45
|
d75_3
|
82917418
|
12.44G
|
79990864
|
12.00G
|
96.47
|
99.97
|
98.33
|
46
|
Table 2
Reference genome alignment reads statistics
Sample
|
Valid reads
|
Mapped reads
|
Unique Mapped reads
|
Multi Mapped reads
|
d45_1
|
80991260
|
76343633(94.26%)
|
61948389(76.49%)
|
14395244(17.77%)
|
d45_2
|
77655516
|
74877365(96.42%)
|
61031456(78.59%)
|
13845909(17.83%)
|
d45_3
|
87070648
|
84828546(97.42%)
|
68283000(78.42%)
|
16545546(19.00%)
|
d55_1
|
96199486
|
93762003(97.47%)
|
76182223(79.19%)
|
17579780(18.27%)
|
d55_2
|
88089256
|
85888657(97.50%)
|
69725622(79.15%)
|
16163035(18.35%)
|
d55_3
|
91438886
|
89149678(97.50%)
|
71925770(78.66%)
|
17223908(18.84%)
|
d65_1
|
87483330
|
85122572(97.30%)
|
68170805(77.92%)
|
16951767(19.38%)
|
d65_2
|
99047498
|
96402904(97.33%)
|
79282917(80.05%)
|
17119987(17.28%)
|
d65_3
|
77823688
|
75757511(97.35%)
|
60872547(78.22%)
|
14884964(19.13%)
|
d75_1
|
77983860
|
75996848(97.45%)
|
62127737(79.67%)
|
13869111(17.78%)
|
d75_2
|
80115068
|
78137926(97.53%)
|
64459926(80.46%)
|
13678000(17.07%)
|
d75_3
|
79990864
|
77918087(97.41%)
|
63626735(79.54%)
|
14291352(17.87%)
|
Analysis of differences in circRNAs
In order to explore further the regulatory role of circRNAs in the early development of cashmere goat hair follicles,we divided the four stages into six comparison groups, and analyzed the differential expression by calculating the SRPBM value of circRNAs (Fig. 3a and Fig. 3b). When analyzing the differential expression of the identified circRNAs, the differential multiple (fold change) and p-value were used by default to screen differential circRNA: That is, simultaneously satisfying the absolute value of log2foldchange greater than or equal to 1 and p-value less than or equal to 0.05. The results showed: d75vsd45, circRNA up-regulated by 59 and down-regulated by 33; d75vsd55, circRNA up-regulated by 61 and down-regulated by 102; d75vsd65, circRNA up-regulated by 32 and down-regulated by 33; d65vsd55, circRNA up-regulated by 67 and down-regulated by 169; d65vsd45, circRNA up-regulated by 96 and down-regulated by 63; and d55vsd45, circRNA up-regulated by 76 and down-regulated by 42 (Fig. 4). By further exploring the law of differential expression of circRNA, it was found that the differential expression of circRNA of d65vsd55 was the greatest, while the differential expression of circRNA of d75vs65 was the lowest (Fig. 4 and Additional files 3:Figure S2).
Validation of circRNAs by qRT-PCR
To validate the accuracy of the circRNA sequencing results,the relative expression of six DE circRNAs(circRNA2049,circRNA3411,circRNA2225,circRNA5681,circRNA1604,and circRNA4351)(Additional files 4:Table S2),were measured by qRT-PCR (Fig. 5).The qRT-PCR results were consistent with the transcriptome sequencing data.
GO and KEGG Pathaway Analysis of Host Genes
Gene ontology analysis includes three domains describing the cellular and molecular roles of genes and gene products (biological process, cellular component, and molecular function) [28]. KEGG is a pathway database for the systematic analysis of gene function, linking genomic and functional information [29].In order to explore the regulatory role of host genes of DEcircRNAs in hair follicle genesis and development in cashmere goats, we analyzed the host genes of circRNA that were differentially expressed in different control groups by GO and KEGG pathway analyses. The results of GO enrichment showed that there was gene enrichment in the biological processes related to the growth and development of hair follicles, such as hair follicle morphogenesis, hair follicle maturation and cell development, and KEGG pathaway analysis showed that there was gene enrichment in the Notch signaling pathway, NF-kappaB signaling pathway, PI3K-AKt signaling pathway and other signal pathways(Additional files 5:Figure S3). Therefore, the parent genes corresponding to differentially expressed circRNAs may be involved in the process of hair follicle growth and development, and then play a regulatory role.
Functional Analysis of CircRNA as a
miRNA Sponge
The ceRNA hypothesis is a new model for post transcriptional gene regulation. According to the hypothesis, the expression of designated miRNAs is reduced by ceRNA [30]. To construct the ceRNA network of circRNA-miRNA-mRNA, we integrated our miRNA library data (unpublished data) and mRNA library data (unpublished data) to analyze the miRNA binding sites in the circRNAs and mRNA using miRanda and Targetscan. We selceted DE exon circRNA in d75vsd45 for construction of the ceRNA regulatory network .In the up-down-up regulation pattern, we predicted 46 circRNA-miRNA and 49 miRNA-mRNA interactions (Additional files 6:Table S3 and Fig. 6a). As shown in Fig. 6a, upregulated circRNA9106 may serve as a sponge for multiple miRNAs (chi-miR-1, chi-miR-18a-3p, and chi-miR-93-3p). Notably, three circRNAs, circRNA8058, circRNA6363 and circRNA8624,contained seed-targets of chi-miR-133a-5p, which were identified by searching for miRNA target sites. Moreover, these miRNAs can downregulate the expression of their target genes (chi-miR-133a-5p was predicted to bind with FLRT1, SOX5,and RBPJL genes). We further constructed a down-up-down co-expression network using miRanda and Targetscan with a strict model, for which 56 circRNA-miRNA and 77 miRNA-mRNA interactions were predicted (Additional files 6:Table S4 and Fig. 6b). In the down-up-down co-expression network,three downregulated circRNAs, including circRNA3382, circRNA1448and circRNA1896 contained seed-targets of chi-miR-26b-3p. chi-miR-26b-3p was predicted to bind with multiple target genes, including MCTP1,MEF2C,GPC6, and FZD5.At the same time, circRNA3236 was predicted to have two target microRNAs, ( chi-miR-27b-3p and chi-miR-16b-3p).
CircRNA3236 binds to chi-miR-27b-3p and chi-miR-16b-3p
CircRNA3236 was down-regulated in d75vsd45;Targetscan and miRanda software predicted that circRNA3236 was targeted to chi-miR-27b-3p and chi-miR-16b-3p, and there were two binding sites respectively. Therefore, two mutant vectors were constructed to verify the specific binding sites(Fig. 7a and Fig. 7b).The results showed that: compared with the NC group, chi-miR-27b-3p and chi-miR-16b-3p significantly decreased the expression of luciferase in circRNA3236 WT (p < 0.001). It shows that there is a binding effect between the two in this experiment. After mu1 mutation, chi-miR-27b-3p failed to down regulate the expression of luciferase in circRNA3236-mut1 (p > 0.05), indicating that the mutation was successful. Mu1 is the binding site of chi-miR-27b-3p and circRNA3236-wt. After mu4 mutation, chi-miR-16b-3p failed to down-regulated the expression of luciferase in circRNA236-mut4 (p > 0.05), indicating that the mutation was successful. Mu4 is the binding site of chi-miR-16b-3p and circRNA3236(Fig. 7c and Fig. 7d).