Expression Pro ling and Functional Analysis of Circular RNAs in Inner Mongolian Cashmere Goats Hair Follicles

Fangzheng Shang (  1540850474@qq.com ) Inner Mongolia Agriculture University: Inner Mongolia Agricultural University Yu Wang Inner Mongolia Agriculture University: Inner Mongolia Agricultural University Rong Ma Inner Mongolia Agriculture University: Inner Mongolia Agricultural University Zhengyang Di Inner Mongolia Agriculture University: Inner Mongolia Agricultural University Zhihong Wu Inner Mongolia Agriculture University: Inner Mongolia Agricultural University Erhan Hai Inner Mongolia Agriculture University: Inner Mongolia Agricultural University Youjun Rong Inner Mongolia Agriculture University: Inner Mongolia Agricultural University Jianfeng Pan Inner Mongolia Agriculture University: Inner Mongolia Agricultural University Lili Liang Inner Mongolia Agriculture University: Inner Mongolia Agricultural University Zhiying Wang Inner Mongolia Agriculture University: Inner Mongolia Agricultural University Ruijun Wang Inner Mongolia Agriculture University: Inner Mongolia Agricultural University Rui Su Inner Mongolia Agriculture University: Inner Mongolia Agricultural University Zhihong Liu Inner Yanhong Zhao Inner Mongolia Agriculture University: Inner Mongolia Agricultural University Zhixin Wang Inner Mongolia Agriculture University: Inner Mongolia Agricultural University Jinquan Li Inner Mongolia Agriculture University: Inner Mongolia Agricultural University Yanjun Zhang Inner Mongolia Agriculture University: Inner Mongolia Agricultural University


Background
There are two kinds of hair follicles in the skin of cashmere goats, namely primary hair follicles and secondary hair follicles. Primary hair follicles produce coarse hairs, and secondary hair follicles produce cashmere. The structural characteristics of the skin and hair follicles of cashmere goats are not only important correlates of their biological characteristics, but also have a direct and important effect on the yield and quality of cashmere.The hair follicle is a skin accessory organ with complex morphology and structure;it controls the growth of hair, and its most prominent feature is regeneration. The initiation of morphogenesis of primary and secondary hair follicles in cashmere goats occurs at different stages of embryonic development, and the initiation of primary hair follicles is earlier than that of secondary hair follicles. At the embryonic stage of 45-55 days, the skin forms a complete epidermal structure, and the hair follicles have not yet appeared; at 55-65 days, the primary hair follicles begin to develop in various parts of the fetus, and the keratinocytes in the basal layer of the epithelium are arranged in a fence to form hair buds, but the formation of primary hair follicles in the lateral part of the body is later than that in other parts (such as the top of the head, shoulder and neck). At 65 days, obvious primary follicle hair buds are observed on the sides of the body. At 65-75 days, the primordial bodies of secondary hair follicles are observed in various parts of the fetus, and secondary hair follicles begin to occur and grow from the epidermis near the primary hair follicles. Similar to the primary hair follicles, the formation of secondary hair follicles in the lateral part of the body is later than that in other parts. At 75 days, obvious secondary hair follicle hair buds are observed on the side of the body(Additional les 1: Figure S1) [1,2].
The skin follicle traits of cashmere goats have a direct and important effect on the yield and quality of cashmere. The ultimate goal of the research in the eld of hair follicle growth and development in cashmere goats is to reveal the mechanism of cashmere growth and nd the important genes related to cashmere growth. The genesis and development of hair follicles may be related to some protein coding genes.
At present, it is considered that most of the signal molecules regulating hair follicle morphogenesis belong to the Wnt pathway [3], tumor necrosis factor(TNF) family, broblast growth factor (FGF)family [4], bone morphogenetic protein (BMP)family [5], Sonic hedgehog (SHH)conduction pathway [6], transforming growth factor(TGF)family [7], and NOTCH conduction pathway [8]. Some of these coding genes are stimulants of hair follicle development and some are inhibitors, which are repeatedly used to regulate each other.
CircRNA is a special kind of RNA, which has no free 5'cap structure and 3'poly (A) structure, and is insensitive to nuclease [9]. CircRNAs can be divided into four types according to its origin:intron circRNA, exon circRNA, exon-intron circRNA,and intergenic circRNA [10]. The main mechanisms of action of circRNAs include: regulating the expression of parental genes [11,12]; interacting with RNA binding proteins [13]; translating proteins [14]; acting as competitive endogenous RNA to regulate the expression of genes [15][16][17]. A current research focus is the action of circRNA,through competitive binding to miRNA ,in regulating gene expression to complete the regulation of life activities has become a research hotspot. Studies have shown that circLMO7 can enhance the HDAC4 expression of the miR-378a-3p target gene through competitive binding of miR-378a-3p, promote muscle cell proliferation and inhibit the differentiation of bovine myoblasts [18]; CircARF3 adsorbs miR-103, to alleviate the targeted inhibition of miR-103 on TRAF3 and alleviate adipose in ammation by promoting mitochondrial autophagy [19]. CircRNA3669 as ceRNA adsorbs miR-26a, and removes the down-regulation of RCN2 by miR-26a in dairy goat endometrial epithelial cells [20]. However, circRNA studies on the regulation of hair follicle development in the embryonic stage of cashmere goats are still scarce.
On the basis of the previous joint analysis of the research group on the regulation of miRNA-mRNA in the development of fetal hair follicles in cashmere goats [21]. In this study, in order to explore the pattern of expression pattern and functional role of circRNAs in the development of fetal skin hair follicles of cashmere goats, we rst used a high-throughput sequencing technique to construct the circRNA expression pro les of cashmere goats during the fetal period (45,55, 65 and 75 days), and identi ed the differentially expressed circRNAs in different control groups. At the same time, the host genes of differentially expressed circRNAs were analyzed by GO and the KEGG pathaway. Following this,the co-expression network of circRNA-miRNA-mRNA was constructed by combining miRNA and mRNA databases, and the binding of circRNA-miRNA was veri ed by dual-luciferase reporter assays. This study has laid a foundation for further exploration of the regulatory role of circRNAs as ceRNAs in the hair follicles of cashmere goat skin, and has also provided a new direction for the study of villus growth and development.

Animals And Samples
In this experiment,three-year-old ewes were selected for estrus synchronization treatment at Jinlai Animal Husbandry, Inner Mongolia, all animal experiments were performed in accordance with the 'Guidelines for Experimental Animals' of the Ministry of Science and Technology (Beijing, China). All surgery was performed according to recommendations proposed by the European Commission (1997), and was approved by the experimental animal ethics committee of Inner Mongolia Agricultural University(hohhot,China).And their mating time was recorded.Skin samples were taken from three sheep at gestational periods of 45, 55, 65and 75 days, immediately treated with DPEC water and placed in liquid nitrogen.The samples were then stored in a refrigerator at-80℃ for RNA-seq and RT-qPCR tests.

RNA library construction and sequencing
The study sequenced 12 samples from the lateral skin of the fetuses of cashmere goats at 45, 55, 65,75 days of. Total RNA was isolated and puri ed using Trizol reagent (Invitrogen, Carlsbad, CA, USA) ,following the manufacturer's procedure. The amount of RNA and purity of each sample were quanti ed using NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA). The RNA integrity was assessed using Agilent 2100. Approximately 5 ug of total RNA was used to deplete ribosomal RNA according to the instructions of the Ribo-Zero™ rRNA Removal Kit (Illumina, San Diego, USA),and the remaining RNA fragments were reverse transcribed using an RNA-seq Library Preparation Kit (Illumina) to form the nal cDNA.Finally, we performed the paired-end sequencing on an Illumina Hiseq 4000 (LC Bio, Hangzhou, Zhejiang,China),following the vendor's recommended protocol.

Identi cation of Transcripts
According to the characteristics of circRNA structure and splicing sequence, and combined with literature reports, we use CIRCExploter2 [22,23] and CIRI [24,25] to predict circRNAs, and integrate the results of the two software programs according to the starting and ending positions of circRNA. The following is the circRNA identi cation standard: 1.mismatch ≤ 2;2.back-spliced junctions reads ≥ 1;3.two splice sites comprise less than 100 kb of the genome. According to the above identi cation and screening criteria, circRNA was identi ed more accurately.

Differential Expression Analysis
When analyzing the differential expression of the identi ed circRNAs, differential multiples (fold change) and p-values were used by default to screen differential circRNAs;that is, genes that simultaneously satisfy the absolute value of log2(fold change) greater than or equal to 1 and p-value less than or equal to 0.05 are marked as yes,otherwise, they are marked as no.

GO and KEGG Pathway Enrichment Analysis
A GO enrichment and KEGG path analysis of the host genes of the DE circRNAs was performed.The GO terms or KEGG pathways with pvalues < 0.05 or the top 20 terms are shown in the gures and tables.

Construction of circRNAs Regulatory Networks
CeRNA was rst proposed by the Pandol team at Harvard Medical School in Cell in 2011. It is pointed out that there are competitive endogenous RNA (ceRNA) molecules in cells. These ceRNA molecules (including lncRNA, circRNA, mRNA, pseudogenes, etc.) can compete through miRNA response elements (MRE) for combination with the same miRNA,in order to regulate each other's expression level [26]. In this study, two software programs, Targetscan (http://www.targetscan.org/mamm_31/) and miRanda (http://www.microrna.org/microrna/home.do), were used to predict the miRNA targeted by circRNA.The target genes targeting miRNA were predicted, and the circRNA-miRNA-mRNA regulatory network was constructed preliminarily. Finally, Cytoscape was used to visualize it.

qRT-PCR
In accordance with the manufacturer's instructions, we used Trizol reagent (Takara, Dalian, Liaoning, China) to extract the total RNA from 12 skin samples representing the fetal period of cashmere goats. Subsequently, we used a PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Dalian, Liaoning, China) to reverse transcribe RNA to cDNA.Following this,qRT-PCR was performed on a LightCycler®96 Real-Time PCR system(Roche,Basel, Switzerland) using TB GreenPremix Ex Taq (Takara, Dalian, Liaoning, China).The qRT-PCR conditions were: 95 °C for 30 s and then 40 cycles at 95 °C for 10 s, 60 °C for 30 s, followed by 72 °C for 10 s. β-actin was used as the reference gene.The validated primers used for RT-qPCR are listed in(Additional les 2: Table S1).The expression levels were calculated using the 2 −△△CT method [27].

Dual-Luciferase Reporter Assays
Chi-miR-27b-3p mimics and chi-miR-16a-3p mimics were synthesized by Hanbio Biotechnology Company(Shanghai, China). The miRNA mimics were transfected into HEK 293T cells useing the LipoFiter transfection reagent according to the manufacturer's instructions. The psiCHECK2-circRNA3236-WT construct was generated by inserting the circRNA3236 fragments containing the miRNA binding sequence into the psiCHECK-2 vector (Promega) at the 3 end of the Renilla luciferase gene. The mutant psiCHECK2-circRNA3236-MUT construct was generated by mutating the miRNA-binding sequence to the complementary sequence using overlapping extension PCR. For circRNA3236 luciferase assays, the HEK 293T cells were transfected with miRNA mimics and either the psiCHECK2-circRNA3236-WT or mutated psiCHECK2-circRNA3236-Mut reporter plasmid. At 48 h post-transfection, luciferase activity was measured using a dual-luciferase reporter assay system (Promega) according to the manufacturer's instructions. The relative luciferase activities were calculated by comparing the Fire y/Renilla luciferase activity ratio.

Identi cation 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 identi ed 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.

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 identi ed 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 les 3: Figure S2).

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 les 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 les 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 identi ed 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 les 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 speci c 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 signi cantly 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.

Discussion
MiRNA is an non-coding RNA studied earlier in skin hair follicles.The miRNA in animal'bodies binds to the untranslated region (UTR) of the 3' end of mRNA to achieve the function of inhibiting mRNA translation into protein,without reducing the number of mRNA molecules [31].
Each miRNA can regulate multiple target genes, and speci c target genes can also be regulated by multiple miRNAs at the same time. In terms of the regulation of hair follicles in the skin of cashmere goats, the researchers identi ed 22 new miRNA and 316 conservative miRNA in adult Inner Mongolia cashmere goats. It is speculated that they may play an important role in the growth of skin and hair follicles [32].
Among the 399 conserved miRNAs screened from northern Shaanxi white cashmere goats, 326 miRNAs were expressed in the growth, retrogression and resting periods, while 3, 12 and 11 miRNAs were speci cally expressed in the growth, retrogression and resting periods, respectively [33]. In Shanbei white cashmere goats, it was found that miRNA-206 regulated the periodic changes of hair follicles by affecting the expression of genes related to the growth and development of hair follicles [34]. LncRNA is a kind of non-coding RNA that is more than 200 bp long and is widespresd in the cytoplasm and nuclei of the cells of animals, plants and microorganisms. In recent years, the role of lncRNA in skin hair follicles has been gradually investigated by researchers. The main mechanism of lncRNA is that ceRNA directly or indirectly interacts with miRNA to inhibit the binding of miRNA to the target gene 3' UTR, or lncRNA may bind directly to the coding gene to achieve the effect of regulating gene regulation. Previous studies have identi ed 403 new lncRNA transcripts and 350TUCP transcripts;3500 differentially encoded mRNA transcripts (a total of 3357 genes) and 172 lncRNA transcripts;at the same time, they identi ed 411 known miRNAs, predicted 307 new miRNAs, and identi ed 72 differentially expressed miRNAs [35]. The expression of lncRNA, microRNA and mRNA in the skin of cashmere goats was systematically identi ed by a multi-group method in two key stages of hair follicle growth (retrogression and growth). The results showed a synergistic effect of lncRNA and miRNA in the process of hair follicle growth transformation;it was found that the retrogression inducing factors TGFβ1 and BDNF were regulated by miRNA-873 and lnc108635596 in lncRNA-miRNA-mRNA network [36]. In addition,regulatory networks such as miR-221-5p-lnc_000679-WNT3,miR-34a-lnc_000181-GATA3,and miR-214-3p-lnc_000344-SMAD3 were constructed, and the interaction between miRNA-mRNA and lncRNA-miRNA was identi ed to illustrate their role in the hair follicle biology of cashmere goats [37]. CircRNA is mostly located in the cytoplasm, and insights that circRNA to mediates the occurrence and development of liver cancer through a ceRNA mechanism [38]. It was found that circ_PUM1 could compete with miRNA-136, resulting in up-regulation of NOTCH3 expression, thus promoting the occurrence and development of endometrial carcinoma [39]. It was identi ed that hsa_circ_0000467 plays a regulatory role in gastric cancer by regulating the level of miRNA-326-3p, and that circRNA may be a potential diagnostic marker and therapeutic target in gastric cancer [40]. It is veri ed that circFUT8 plays an inhibitory role in bladder cancer cells by targeting miRNA-570-3p/KLF10 [41]. In non-human animal research, Zhang et al found that circRNA-006258 regulates the growth and lactation of goat mammary epithelial cells through a circRNA-006258-miR-574-5p-EVI5L regulatory network [42]. Wang et al constructed a circRNA expression pro le related to milk fat metabolism under heat stress and established the related ceRNA regulatory network [43]. Li et al analyzed the differentially expressed circRNA, between fast contractile muscle and slow contractile muscle of porcine skeletal muscle and established a ceRNA network by combining the differentially expressed circRNA with miRNA and mRNA databases, which was preliminarily veri ed by double luciferase [44]. However, there are no reports on circRNA related to fetal hair follicle development in cashmere goats. In this study, a total of 21784 circRNA were identi ed in four stages of fetal skin hair follicles of cashmere goats. The differentially expressed circRNA of each comparison group was screened and the differential circRNAs were combined with miRNA and mRNA using the mode of upregulation-downregulation-upregulation or downregulation-upregulation-downregulation to construct the ceRNA network. Some circRNAs and differentially expressed circRNAs at different stages have been identi ed previously in the hair follicle development of previous Angora rabbits [45],but the total number of circRNAs, and differentially expressed circRNAs were greater and a ceRNA network was established in this study. It is suggested that circRNAs may play an irreplaceable role in the development of hair follicles in cashmere goats.
Studies have shown that circRNAs may affect biological function by regulating the level of linear mRNA expression [46]. In this study, the parental genes of differential circRNAs were analyzed by GO and the KEGG pathaway. The biological process of GO enrichment includes hair follicle morphogenesis, hair follicle maturation and cell growth; the pathway of KEGG enrichment includes the Notch signaling pathway and NF-kappaB signaling pathway. Previous studies have shown that there is a direct relationship between the Notch signal pathway and hair follicle morphogenesis [47], high expression of Notch1 and Notch2 can accelerate the formation of mouse hair substrate and, at the same time, it can inhibit the cells around hair substrate to form substrate [48]. Ocu-miR-205 can promote the transition of Rex rabbit hair follicles from the growing stage to a degenerative quiescent stage by regulating the expression of related genes and proteins in Notch, BMP and other signal pathways,changing the hair density [49]. Krieger et al.used mice model to study the regulation of the NF-kb signal pathway in the hair follicle cycle, and found that the NF-kb signal pathway is essential for the growth and activation of hair follicle stem cells [50].
The study of the NF-kappaB signal pathway in mice by Schmid et al showed that NF-kappaB signal pathway was activated downstream of EdaA1 and EDAR, thus playing an important role in the development of hair follicles [51]. Zhang et al.found that Wnt/ β-catenin signal transduction is a necessary signal for NF-kB activation, and EDAR is the direct target gene of Wnt. The initial activation of the Wnt/ βcatenin signaling pathway depends on the activity of EDA/EDAR/NF-kB in the prohair substrate of primary hair follicles. The complex interaction and interdependence of Wnt/ β-catenin and EDA/EDAR/NF-kB signaling pathways in the initiation and maintenance of primary hair follicle substrate were revealed [52]. Therefore, we speculate that circRNA may play a regulatory role in the primary stage of hair follicle development via the Notch signaling pathway and the NF-kappa B signaling pathway .
In order to explore further the functional mechanism of circRNA in the hair follicle development of cashmere goats, we predicted the miRNA targeted by circRNA and the target genes targeted by miRNA, which were differentially expressed by d75vsd45, and thus constructed the ceRNA network.Among them, there were 16 circRNAs targeted 6 miRNAs,and 6 miRNAs targeting 23 target genes in the up-regulation-down-

Conclusion
We constructed the expression pro les of circRNAs in the hair follicles of Inner Mongolia cashmere goats at different embryonic stages (45,55, 65, and 75 days), and a total of 21784 circRNA were identi ed. The differentially expressed circRNA in six control groups were screened.
Six differentially expressed circRNAs were randomly selected for a qRT-PCR experiment, which showed that the sequencing results were reliable. The results of GO and KEGG analysis of parental genes showed that there may be some relationship between circRNA and its parental genes, and circRNA may play a role in the biological process of hair follicle growth and development in cashmere goats through the Notch signaling pathway and NF-kappaB signaling pathway.At the same time, the co-expression network of circRNA-miRNA-mRNA was constructed, and the interaction between 102 pairs of circRNA-miRNA and 126 pairs of miRNA-mRNA was studied. The results of the dual luciferase assay showed that circRNA3236 and chi-miR-27b-3p,and circRNA3236 and chi-miR-16b-3p have a targeted binding relationship,.The speci c binding sites were veri ed,which provides an important basis for exploring the molecular mechanism of circRNA as ceRNA in the genesis and development of hair follicles, and also provides important information for studying the mechanism of action of circRNA in human hair follicles.

Abbreviations
CircRNA:circular RNA;ceRNA:competing endogenous RNA;qRT-PCR:quantitative RT-PCR Declarations Figure 1 Identi cation results of circRNAs.The purple columns represent the number of circRNAs and the green columns represents the number of circRNA-hosting genes.   Differential expression of circRNAs in different groups.The red columns represent up-regulated circRNAs and the green columns represents down-regulated circRNAs.

Figure 5
The expression quantity and expression trend of circRNA in different periods. 0.8< |Rs|<1 indicates a strong correlation; 0.6< |Rs|<0.8 indicates a strong correlation; 0.4< |Rs|<0.6 indicates a moderate correlation; 0.2< |Rs|<0.4 indicates a weak correlation; 0< |Rs|<0.2 indicates no correlation; and the degree of proximity between |Rs| and 1 represents the degree of closeness and correlation between two variables.

Figure 6
CircRNA-miRNA-mRNA regulatory network analysis in Cashmere goat hair follicle. a Upregulated circRNA networks and b downregulated circRNAs networks for d75vsd45 .