Secondary wool follicle growth process
Hematoxylin and eosin (H&E) staining at E90d showed primary and early secondary stage wool follicles (Figure 1a). From observing wool follicles at this stage, it is clear that primary wool follicles occur early, the bulbs are large, the wool follicles are long and have accessory structures such as sweat glands, sebaceous glands, and the arrector pili muscles Secondary wool follicles at this stage are smaller and grow nearer to the epidermis than the primary wool follicle (Figure 1b). At E120d, the secondary wool follicles are separated from and arranged in parallel with the primary wool follicles (Figure 1c, d). By birth, some of the secondary wool follicles have matured and their wool has passed through the body surface (Figure 1e, f).
Sequencing and mapping of the sheep skin tissue transcriptome
To examine the circRNAs expression profiles in sheep skin at different developmental stages, we performed RNA Integrity Number (RIN) tests on nine sheep skin tissue samples, three from each of the three developmental periods (embryonic day 90, embryonic day 120, and birth). The RIN values of the samples were 7.0, 7.3, and 8.7 for E90d, 6.9, 7.4, and 7.0 for E120d, and 7.8, 8.1, and 7.7 for samples at Birth. Results show that the RNA quality met the minimum requirements of sequencing. Library was thus constructed and the samples were sequenced. Raw reads were acquired via Illumina sequencing, which were then processed to remove rRNA, low-quality sequences, and junction contamination, among other processing. All subsequent analyses were based on these processed clean reads. We obtained the following clean reads from the three developmental stages of the sheep: 286,156,478 clean reads for E90d, 284,713,262 clean reads for E120d, and 261,413,764 clean reads for samples collected at birth. These reads were mapped to the sheep genome. A total of 8,753 candidate circRNAs and 3119 source genes were identified (Additional file 1: Table S1), 1,648 of which (18.8%) were expressed at all developmental stages (Figure 2a). The 30 highest-expressed circRNAs in each group are listed in Table 1. Based on their location in the genome, the 8,753 circRNAs were classified into six types: (1) Classic: when the formation site of the circRNA is exactly on the boundaries of exons (83.4%); (2) Alter-exon: when one end of the circRNA formation site is on the exon boundary, and the other end is inside the exon (8.6%); (3) Intron: when the formation site of the circRNA is completely in the intron region (1.2%); (4) Overlap-exon: when the formation site of the circRNA spans the exon region (5.5%); (5) Antisense: when the circRNA is formed by the antisense strand of the gene (0.3%); (6) Intergenic: when the formation site of circRNA is completely inside the intergenic region (1.0%) (Figure 2b). The circRNAs typically comprised of two to four exons (Figure 2c). In circRNAs with only one exon, the length of the exon was found to be significantly longer than that of a circRNAs comprised of multiple exons (Figure 2d). The peak gene density, based on the expression of circRNAs in all samples, was between 0.3 and 0.4 (Figure 2e).
Identification of differentially-expressed circRNAs
Based on the criterion of differentially-expressed circRNAs, clustering maps (Figure 3a) were used to illustrate their distribution. Significantly differentially-expressed circRNAs in the figure are in yellow (upregulated expression) or blue (downregulated expression). The number of differentially-expressed circRNAs in the three developmental stages are displayed in Figure 3b, c. We detected 377 differentially-expressed circRNAs and 314 source genes when comparing between Birth and E90d, 467 differentially-expressed circRNAs and 383 source genes when comparing between Birth and E120d, and 507 differentially-expressed circRNAs and 417 source genes when comparing between E120d and E90d (Additional file 2: Tables S2A, S2B, S2C).
Among the DEGs (Differentially expressed genes), circ_0004932 and circ_0004936 correspond to gene 13410 (TRPS1). It has been reported that Trps1 is involved in the growth and development of hair follicle cells [24]. There were also some circRNAs similar to circ_0004932 and circ_0004936, such as circ_0000997 and cir_0000999 that were mapped to source gene 851 (VAV3), and circ_0001520 and circ_0001524 that were mapped to source gene 3008 (TMEFF1), all associated with hair follicle growth [25, 26]. We also found that the expression level of circ_0006736 at E20d and Birth stages was significantly higher than at E90d. It might therefore play a role in the growth, development, and maturation of SF. Mapping showed that gene 20646 (SMAD1) is the source gene of circ_0006736. This gene can control the transformation of early hair follicle morphology by controlling the activity of stem cells [27]. The expression levels of circ_0005454 and circ_0005453 during E120d were significantly higher than at E90d, and SF grew significantly during the time between E90d and E120d, so we speculate that circ_0005454 and circ_0005453 participate in the growth of SF. Circ_0004116 is highly expressed at all three developmental stages, and might be active through the entire growth process of the wool follicle - growth of both PF and SF, the source gene of circ_0004116 is gene 11842 (RFX7). We hope to further study the function of RFX7 in AFWS wool follicles in the future.
Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses
The function of circRNA is reflected through their source gene. The function of circRNAs can be further studied by analyzing the Gene Ontology (GO) terms of their source genes. Based on differentially-expressed circRNAs and their source genes statistics (Additional file 2,Table S2), the top ten terms of candidate genes in each comparison group were selected for mapping (Figure 4 a-c), the detail information was listed in Additional file 3: Tables S3A, S3B, S3C. The most significantly enriched GO terms were: cellular component organization (GO: 0016043), regulation of primary metabolic process (GO: 0080090), intracellular part (GO: 0044424), intracellular organelle (GO: 0043229), membrane-bounded organelle (GO: :0043227), and protein binding (GO: 0005515).
To predict the pathways of the significantly enriched source genes, we performed an enrichment analysis using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (Figure 4d, Additional file 4: Table S4A, S4B, S4C). Results show that six significantly enriched pathways (endocytosis, lysine degradation, apoptosis, human papillomavirus infection, adherence junction, and tight junction) were identified. These involve 55 enriched source genes and their corresponding 255 circRNAs (Additional file 5: Table S5A). Among the 55 source genes, seven are associated with wool follicle growth. There were 35 circRNAs associated with these seven source genes (Additional file 5: Table S5B), of these 35 circRNAs, seven were found to be significantly differentially-expressed in our study: circ_0005720 from source gene 15869 (AKT3), circ_0001754 from source gene 3277 (TGFBR1), circ_0008036 from source gene 25354 (SMAD2), circ_0004032 from source gene 11746 (SOS2), circ_0005174 from source gene 13720 (RB1), circ_0005519 from source gene 15130 (EZH1), and circ_0007826 from source gene 24949 (FGFR2). A network describing the connections between the source genes and circRNAs was constructed (Figure 5).
Target miRNAs of differentially expressed circRNAs at the different developmental stages in sheep
To further understand the functions of circRNAs, the miRanda software was used to predict the interactions between the identified circRNAs and miRNAs. A total of 17 circRNAs and eight miRNAs were identified, and the relationship between them were constructed into a network (Figure 6, Table 2). For example, circ_0003042 is significantly differentially-expressed between Birth and E120d. It is predicted to interact with miR-432, and might act as a "miRNA sponge," binding all of the available miR-432 to prevent them from exerting their function.
Validation of circRNAs expression by qRT-PCR
To validate the expression levels of differentially-expressed circRNAs, we randomly selected seven highly expressed circRNAs and detected their expression levels by qRT-PCR (Additional file 6: Table S6). These results were consistent with the trends observed in the RNA-Seq data, the results for all circRNAs were r > 0.8, indicating that the RNA-Seq is reliable (Figure 7a-g). As can be seen in Figure 7 h, the circRNAs we selected could resist the digestion of RNase R, while the linear RNA in the sample (GAPDH) could not. After RNase R digestion, the expression of the seven circRNAs did not decrease significantly. On the contrary, most of them actually increased, we speculated that circRNAs were relatively enriched, and the efficiency during reverse transcription has relatively improved. The relative expression levels quantified by qRT-PCR have therefore also increased, RNase R digestion basically increased the purity of circRNAs. The results show that circRNAs can resist the digestion of RNase R, while linear RNAs cannot. Therefore, RNase R digestion is essential for the detection of circRNAs.