Differentially Expressed Circular RNA Profile Signatures Identified in Prolificacy Trait of Yunshang Black Goat Ovary at Estrus Cycle

CircRNAs acting as miRNA sponges play important roles in the growth process of animal individuals. The prolificacy trait of goats is involved in many pathways, however, the variation of circRNA expression profiles in the different phases of the estrus cycle at high and low fecundity groups is still unknown. Here, we analyzed the circRNA profiles of ovarian tissues among high and low fecundity groups in the follicular phase (HF vs LF), high and low fecundity groups in the luteal phase (HL vs LL), and high and low fecundity in the whole estrus cycle (HF vs HL and LF vs LL) using RNA-seq. A total of 283 (114 upregulated and 169 downregulated), 559 (299 upregulated and 260 downregulated), 449 (254 upregulated and 195 downregulated), and 314 (210 upregulated and 104 downregulated) differentially expressed (DE) circRNAs were screened in HF vs LF, HF vs HL, HL vs LL, and LF vs LL groups, respectively. Enrichment analysis suggested that the targeting of DE circRNAs was mainly enriched in oocyte meiosis, the GnRH signaling pathway, and estrogen signaling pathway. After integrating our previous study of miRNA-seq, there were 56 miRNAs that could target to 192 DE circRNAs, including the miR-133 family (including miR-133a-3p and miR-133b), miR-129-3p, and miR-21, which also had important influence on the prolificacy trait of goats. Then, 18 circRNAs with coding potential were obtained by four software predictions, and 9 circRNAs were validated by RT-qPCR. Together, circRNAs play a key role in the prolificacy trait and the transformation of the follicular phase to the luteal phase in the estrus cycle of goats.


INTRODUCTION
The kidding number of goats, an essential economic indicator of the goat industry, is mainly dependent on ovulation rate and ovarian growth (Notter, 2012). Many factors influence the kidding number and ovulation in goats. The reproductive performance of goats is different to that of sheep, and the genes conferring the major effects on prolificacy of goats are still unknown. Prolific and nonprolific breeds are characterized by differences in the estrus cycle. Determining the molecular mechanisms of follicular development has therefore been a focus of research in the genetic breeding of goats. In recent years, many studies of the mechanisms of prolificacy in goats have focused on candidate genes and their regulators, such as miRNAs, transcription factors, and varieties of regulated pathways. All these factors form a network of many interactions in vivo that regulate the goat reproduction trait Zhao et al., 2020). The genes GDF9, Smad2, BMP15, and BMP4, and the transcription factor Sox9 are involved in kidding number and increased ovulation rate, and are proven to be crucial for the reproductive performance of goats (Hanrahan et al., 2004;Lai et al., 2016;Chinmoy et al., 2017;He et al., 2020;Wang et al., 2021). In addition, many factors affected the reproductive ability of goats including follicle recruitment, development, selection, and ovulation (Luo et al., 2019). A few signaling pathways, such as the Wnt/β-catenin, nodal mammalian target of rapamycin (mTOR), and bone morphogenetic protein (BMP)/Smad pathways, are involved in follicular development .
CircRNAs are one type of long noncoding RNAs produced by pre-mRNA back-splicing that have been progressively identified and validated to form a covalently closed loop structure of singlestrand, non-polyadenylated circular molecules (Memczak et al., 2013;Ashwal-Fluss et al., 2014), which show a previously neglected regulatory method in the form of the circRNA-miRNA-mRNA axes (Hansen et al., 2013). Usually, circular RNAs act as a sponge for related miRNAs, and compared with their corresponding colinear mRNA isoforms, circRNAs are more stable (Chuang et al., 2018). Recently, many studies on circRNA function have shown that regulation of circRNA expression can affect a variety of molecular and physiological phenotypes, mainly related to the nervous system, innate immunity, and diseases (Xiao et al., 2020). CircRNAs also act as a template for translation (Kristensen et al., 2019). Many circRNAs were shown to play important roles in the process of reproduction by regulating the expression of genes at the transcriptional and post-transcriptional levels. In the study of reproduction trait, circRalGPS2s were lowly expressed in follicle growth and predominant in the ovarian stroma, suggesting that circularization of circRNAs may not have a relationship with the development of granulosa cells (Shen et al., 2019). A few studies showed that the expression of circRNA_103827 and circRNA_104816 was negatively correlated with embryonic development (Cheng et al., 2017). Although circRNAs have various functions and act through mechanisms, the regulatory effect of circRNAs on the goat reproduction trait is uncharacterized. Hence, it is crucial to reveal the molecular mechanism of circRNAs regulating the prolificacy trait and the transition from the follicular phase to the luteal phase in goats.
The Yunshang black goat is an excellent breed of meatproducing goat, which has many favorable characteristics including tender meat and unique flavor, and live in Yunnan Province, China. Yunshang black goat is suitable for large-scale breeding in all of regions of China. Many studies have examined the molecular mechanisms of resource conservation, strain selection, and trait formation of this native variety. One excellent feature of Yunshang black goat is prolificacy, making it suitable for the study of prolific goats. The molecular mechanisms underlying prolificacy trait regulation in the Yunshang black goat have not yet been confirmed. In order to reveal the key molecular mechanism networks of the prolificacy trait and follicular development in this variety, in this study, we screened for differentially expressed (DE) circRNAs among high and low fecundity goats in two phases of the estrus cycle (four groups: high fecundity in the follicular phase, high fecundity in the luteal phase, low fecundity in the follicular phase, and low fecundity in the luteal phase, respectively) by RNA-seq. The crucial correlated miRNA-circRNA and circRNA-mRNA regulator networks and pathways in the prolificacy trait and the transition from the follicular phase to the luteal phase were identified. Furthermore, the function and mechanism of DE circRNAs were predicted by functional enrichment and coding potential analysis. These analyses will provide a basis for the molecular mechanism of follicular development and the transition from the follicular phase to the luteal phase in goats. Our purpose was to clarify the molecular mechanisms of follicular phase regulation of goat ovaries and to provide a research basis for further study on the interactions between circRNAs, miRNAs, host genes, and signaling pathways for prolific traits and the transition from the follicular phase to the luteal phase in goats. This study also provides basic data for Yunshang black goat prolificacy traitrelated candidate gene screening and the conservation of native breed resources.

Sample Collection and Preparation
Twenty she-goats (2-3 years old), no significant differences in age and weight, were selected and grouped into four groups, including high fecundity goats in the follicular phase (HF) and in the luteal phase (HL) (n = 5, average kidding number ≥3) and low fecundity goats (n = 5, average kidding number ≤2, LH and LL groups) based on their kidding number records. All individuals were fed in the same environment and had free access to water in a Yunnan Province goat farm. The ovaries at the follicular phase of the estrous cycle were collected from animals with high (≥3) and low (≤ 2) numbers of kids, frozen immediately in liquid nitrogen, and stored at −80°C for conservation.

Preparation of RNA Library for Sequencing
According to the manufacturer's protocol, total RNA was isolated from ground ovarian tissue powder using TRIzol reagent (Invitrogen, Carlsbad, CA, United States ). The degree of degradation and contamination of RNA samples was monitored by 1% agarose gels, and the NanoPhotometer ® Finally, the Agilent Bioanalyzer 2,100 system was used to assesses library quality, and the AMPure XP system was used to purify the products.

Identification and Analysis of Differentially Expressed circRNAs
CircRNAs were identified by find_circ and CIRI2 (Memczak et al., 2013;Gao et al., 2017). The following are the basic principles of find_circ: according to Bowtie2 alignment results, find_circ extracted the 20 nt anchor sequences from the ends of the reads that were not aligned to the reference sequence, and aligned each pair of anchor sequences with the reference sequence again. If the 5′-end of the anchor sequence aligned with the reference sequence (start and end sites A3 and A4, respectively), the 3′-end of the anchor sequence was aligned upstream of the site (start and end sites A1 and A2, respectively), and there was a splice over the site (GT-AG) between the A2 and A3 of the reference sequence, then this read could be considered a candidate circRNA. Finally, the candidate circRNAs with read counts greater than or equal to 2 were identified as circRNAs. CIRI2 searched for a PCC (paired end mapping) signal and PEM (pair-end mapping) signal and junction reads of the shear site of paired chiastic clipping according to the comparison of the results of BWA (Li and Durbin, 2009). Then, based on the results of dynamic programming alignment, circRNA read support numbers and genome annotation information was used to filter the candidate circRNAs. Quantification was performed based on reads that mapped across the circularized junctions and were normalized using SRPBM (number of circular reads/number of mapped reads (units in billion)/read length) (Zheng et al., 2016). Differential circRNA expression was analyzed using DESeq2. CircRNAs with padj values of < 0.05 were considered to be significantly differentially expressed.

The GO and KEGG Enrichment of Host Genes With DE circRNAs
Functional annotation of host genes with DE circRNAs was based on the GO and KEGG annotations of the source genes. GO annotation was performed based on the corresponding genes in NCBI and their GO a nnotations. This information was held in the following database https://ftp.ncbi.nlm.nih.gov/gene/DATA/ gene2go.gz. For KEGG annotation, the KOBAS software was used to test the statistical enrichment of host genes with DE circRNAs (Xie et al., 2011).

Prediction Analysis of the Targeting miRNA of DE circRNAs
CircRNAs, a type of targeting molecule of miRNAs, are regulated by miRNA. miRanda was used to predict the miRNA interaction sites on circRNA for animals (Fahlgren and Carrington, 2010). The function and mechanism of action of circRNAs as miRNA sponges were investigated through the analysis of circRNA-miRNA interactions.

Coding Potential Prediction of circRNAs
CircRNAs have the potential to encode protein, and in this study IRESfinder software was used to calculate the IRES score for each circRNA. The score of >0.5 indicated the presence of IRES elements, with a higher score indicating higher accuracy. The coding potential of the circRNA sequences were further identified by the CPC2, CNCI, CPAT, and PLEK software programs.

Verification of Sequencing Results by RT-qPCR
Several DE circRNAs were randomly selected to verify the reliability of the sequencing data, and confirmed by RT-qPCR with RPL19 used as the reference gene. Accounting for circRNA characteristics, RT-qPCR primers were designed, and relevant primer information is shown in Table 1. Three duplicates per sample were analyzed, and a RocheLight Cycler ® 480 Ⅱ system (Roche Applied Science, Mannheim, Germany) was used to establish a standard curve. The following steps were applied for RT-qPCR: initial denaturation at 95°C for 5 min, followed by 40 cycles of denaturation at 95°C for 5 s, and annealing at 60°C for 30 s. The data were analyzed using the 2 -ΔΔCt method (Liu et al., 2021). The relative expression was calculated by a dependent sample t-test and SPSS 19.0 was used to analyze the significance of any differences.

PPI Analysis of Host Genes
As the goat database is not complete, we used the database of sheep, which has a high homology with goats, as a reference. The protein-protein interaction (PPI) analysis of host genes of DE circRNAs was conducted according to the STRING database (Organism: Ovis aries). We constructed the PPI network by extracting a list of host genes from the sheep database. In this study, the online tool STRING was used analyze the PPI of host genes, and PPIs with a score of >700 were selected as significant interactions.

Analysis of Transcriptome Sequencing of Ovaries in Yunshang Black Goats
Twenty cDNA libraries of ovaries from high and low fecundity goats in the follicular and luteal phases were constructed, and an Illumina Hiseq 2,500 sequencer was used to perform the sequencing. The samples were divided into four groups, and five replicates of each sample were considered: the high fecundity group in the follicular phase (HF-1, HF-2, HF-3, HF-4, HF-5), the high fecundity group in the luteal phase (HL-1, HL-2, HL-3, HL-4, HL-5), the low fecundity group in the follicular phase (LF-1, LF-2, LF-3, LF-4, LF-5), and the low fecundity group in the luteal phase (LL-1, LL-2, LL-3, LL-4, LL-5). In Table 2, the number of clean reads and GC contents for the four groups are shown. The average rate of clean reads that uniquely mapped to the Capra hircus genome was 90.80% (from 87.16 to 94.36%). After mapping the reference sequence, we identified 6,034 circRNAs. The spliced length distribution of circRNAs was primarily in the range from 200 bp to 10,000 bp ( Figure 1A). CircRNAs are mainly derived from exon splicing, with  smaller proportions arising from intron and intergenic splicing ( Figure 1B). The circRNA density statistics for each chromosome showed that circRNAs were distributed on chromosome 1 to 11, with the highest proportion of circRNAs located on chromosome 1, 2, and 3 (approximately 32%) (Figure 2 and Supplementary Figure  S1). In addition, circRNAs were named after the host genes (Supplementary Table S1). Frontiers in Physiology | www.frontiersin.org April 2022 | Volume 13 | Article 820459 5 circRNAs were randomly selected for RT-qPCR across the four groups. The RT-qPCR results accorded with the RNA-seq data ( Figure 4).

Enrichment Analysis of the Host Genes for DE circRNAs
GO and KEGG enrichment analyses of host genes for DE circRNAs were performed using the appropriate online databases. GO terms analysis showed that there were 558 host genes with DE circRNAs in the HF vs HL comparison, which were annotated into 4,607 functional subclasses with 814 significant GO items. In this comparison, the most significantly enriched GO terms were in cell and cellular process subclasses. In the LF vs LL comparison, a total of 3,684 functional subclasses and 492 terms were significant; most enrichment occurred in cellular component (CC) subclasses. Moreover, in the HF vs LF comparison, terms were related to cell part and cell terms. There were 3,643 enriched terms in this comparison, with 811 terms significantly enriched. In the HL vs LL comparison, we found 4,377 terms were enriched, 571 of which were significantly enriched. Of which, a lot of host genes for DE circRNAs were enriched in the cell (GO:0005623), cell part (GO:0044464), and cellular process (GO:0009987) subclasses, along with other subclasses ( Figure 5A). Our results suggested that some GO terms associated with cellular process showed the highest value in the HF vs HL comparison, followed by that in HL vs LL, LF vs LL, and HF vs LF. However, for almost all reproduction-related GO terms, such as reproductive process (GO:0022414), ovulation cycle (GO:0042698), oogenesis (GO:0048477) and sexual reproduction (GO:0019953), the enrichment in the HL vs HF comparison was higher than in the LL vs LF comparison. We speculated that the host gene enrichment in these GO terms was associated with the transformation of the follicular phase to the luteal phase. In the HF vs LF and HL vs LL comparisons, the GO terms were more focused on reproductive process (GO:0022414), hormone receptor binding (GO:0051427), steroid metabolic process (GO: 0008202), and ovulation cycle process (GO:0022602). These results showed that there were differences in the regulation of biological processes related to reproductive processes and hormone secretion in the different comparisons.
Comparison of the host genes and KEGG pathway analysis showed that the host genes were annotated into 154, 175, 157, and 183 pathways for the HF vs LF, HL vs LL, LF vs LL, and HF vs HL comparisons, respectively ( Figure 5B). Of these, progesterone-mediated oocyte maturation (ko04914), ovarian steroidogenesis (ko04913), oocyte meiosis (ko04114), the GnRH signaling pathway (ko04912), the estrogen signaling pathway (ko04915), and other pathways were involved in the development of ovarian and oocyte maturation.

Analysis of Targeting Relationship Between circRNAs and miRNAs
To further identify the key circRNAs in the prolificacy trait of goats, the binding sites of miRNA-bound circRNAs were predicted by miRanda software. We found that 56 miRNAs bind to 192 DE circRNAs. Of note, our result showed that the relationship of miRNA-circRNA showed the one-to-many and many-to-one characteristic ( Figure 6, Supplementary  Table S2). miR-133a-3p, miR-133b, miR-129-3p, and miR-21 were found, which mainly have an impact on ovary development, thus, their target circRNAs were likely to be associated with ovary development.

Identifying the Coding Potential of circRNAs
Among the four groups, in total, 4083 DE circRNAs were calculated with scores greater than 0.5 by IRESfinder software. Furthermore, four software programs were used to predict the coding potential of the above circRNAs, and 3,257 circRNAs were found. Finally, 18 circRNAs were identified with coding potential by the four software programs mentioned above ( Table 3).

PPI Network for Goat Prolificacy Trait
The PPI network of the host genes was constructed by Cytoscape using the analysis result of host genes from the STRING database ( Figure 7). For the four comparisons, the PPI network contained 69 protein-protein pairs. MAPK1, which has three nodes including MAPK1-PRS6KA1, MAPK1-PTPN11, and MAPK1-PLD1, is associated with ovarian development and diseases. Other pairs were known to influence the growth and development of the ovary and the prolificacy trait of mammals.

DISCUSSION
The process of reproduction is very complex. In particular, the study of the goat reproductive molecular mechanism is still in its infancy and there are still many unanswered questions. The reproduction traits in goats are governed by the precise regulation of genes, regulatory elements, and other factors, such as spatiotemporal expression, transcriptional regulation, and feedback mechanisms (Ye et al., 2018;Luo et al., 2019). Investigation of the ovary from goats in different estrus stages should provide information on potential follicle development and allow further improvement to kidding numbers. In a previous molecular biology study, hormone secretion and cell signaling transduction were found to play an important role in the mechanisms of ovarian and follicular development (Messina et al., 2016;Zou et al., 2020). Recently, circRNAs were   TRUE  TRUE  TRUE  TRUE  circSENP5  0.51  TRUE  TRUE  TRUE  TRUE  circUSP24  0.61  TRUE  TRUE  TRUE  TRUE  circASB1  0.53  TRUE  TRUE  TRUE  TRUE  circWASL  0.51  TRUE  TRUE  TRUE  TRUE  circCPSF6  0.57  TRUE  TRUE  TRUE  TRUE  circSLC12A2  0.72  TRUE  TRUE  TRUE  TRUE  circXPA  0.63  TRUE  TRUE  TRUE  TRUE  circSMC6  0.54  TRUE  TRUE  TRUE  TRUE  circPDS5B  0.58  TRUE  TRUE  TRUE  TRUE  circNBEA  0.57  TRUE  TRUE  TRUE  TRUE  circFAM188A  0.85  TRUE  TRUE  TRUE  TRUE  circISM1  0.51  TRUE  TRUE  TRUE  TRUE  circEIF3E  0.77  TRUE  TRUE  TRUE  TRUE  circZCCHC8  0.68  TRUE  TRUE  TRUE  TRUE  circANKRD27  0.52  TRUE  TRUE  TRUE  TRUE  circVRK1  0.52  TRUE  TRUE  TRUE  TRUE  circEIF4E3  0.59  TRUE  TRUE  TRUE  reported to have a range of extremely important effects on many kinds of biological processes (Memczak et al., 2013). A lot of studies have showed that circRNAs are widespread in plants and mammals, however, little is known about circRNAs in goats, especially those related to the prolificacy trait and ovarian development during the different phases of the estrus cycle (Rybak-Wolf et al., 2015;Ye et al., 2015).
Here, the high-throughput sequencing of circRNAs in the ovaries from goats with different kidding number records was performed in different estrus phases. The circRNAs were distributed widely among chromosomes 1-11 in the four comparisons, and which showed that circRNAs had complexity and functional diversity (Figure 2, Additional file Supplementary Figure S1). The majority of circRNAs in this study were <10,000 nt in length, which was longer than those in previous studies in chicken granulosa cells (200-300 nt) (Shen et al., 2019). In total, 53.8% of the circRNAs were shorter than 10,000 nt and longer than 1,000 nt, the proportion was also higher than that in the studies of goat endometrium circRNAs (23% > 1,000 nt) (Song et al., 2019). Species and tissue differences may be the main reason for the variation in circRNA length. In our study, the main source of circRNAs in the ovary was exon splicing, which was supported by previous research (Lv et al., 2020). Given all of the above features, it was concluded that there were species-and tissue-specific differences in circRNA expression. The changes in circRNA abundance in the four comparisons might be associated with their specific roles during the different estrus phases of the ovary and kidding numbers. We identified 283, 559, 449, and 314 DE circRNAs in HF vs LF, HF vs HL, HL vs LL, and LF vs LL comparisons, respectively ( Figure 3B, Supplementary Table S1). Then, these host genes of DE circRNAs were analyzed for GO terms and KEGG pathways. The majority of enriched GO terms were related to cell proliferation and reproductive process. There were 14 host genes of circRNAs related to the reproductive process; notably, mitogen-activated protein kinase 1 (MAPK1), which could mediate luteinizing hormone-induced communication disruption and oocyte maturation in ovarian follicles of rat (Gratao et al., 2008). MAPK1 was shown to play a basic role in regulating spindle assembly, microtubule organization, chromosome distribution, and meiosis resumption, and it could be hypothesized that activation of MAPK1 might improve meiosis resumption (Gordo et al., 2001). MAPK1 mRNA plays a crucial role in oocyte maturation in various species, including cattle and dogs, through its action on granulosa and cumulus cells (Quetglas et al., 2010;Lee et al., 2017). Moreover, some other reproduction-related pathways were also significantly enriched, including oocyte meiosis, the GnRH signaling pathway, and the estrogen signaling pathway.
In the HF vs LF comparison, circIGF2R, circMAPK8, circMYH10, and circPTEN were associated with the prolificacy trait. The IGF2R locus contains a highly abundant set of circIGF2R circRNAs, that showed multiple forms of different exon combinations in mice and plays significant roles in mammalian reproduction (Perera et al., 2018). MAPK8 might regulate the transcription of basal FSHb, and one study has shown that blocking its activation could inhibit the levels of major basal FSHb transcripts by half (Haisenleder et al., 2008). In the HL vs LL comparison, the host genes of circFSHR, circHMGCR, circSMAD3, and circFGFR2 were FSHR, HMGCR, SMAD3, and FGFR2, respectively. FSHR stimulates Gαs proteindependent activation of cAMP/PKA, leading to phosphorylation of cAMP response element binding protein and steroidogenic protein, which is required for the production of estrogens, which in turn is known to be a growth stimulant (Casarini and Crépieux, 2019). In a study of common genetic variants affecting dizygotic (DZ) twins, a GWAS conducted in 1980 mothers of spontaneous DZ twins and 12,953 controls found that SNPs similar to SMAD family member 3 (SMAD3) were significant associated with DZ twinning (Mbarek et al., 2016). Collectively, these results suggested that those circRNAs were involved in ovarian and oocyte development and thus influence the number of kidding.
Among the DE circRNAs of the comparison between the estrus phases (HF vs HL and LF vs LL), we found that the host genes of circLDLR, circMAPK1, circSMAD1, and circSMAD2 were LDLR, MAPK1, SMAD1, and SMAD2, respectively. SMAD1 and SMAD2 are members of the SMAD family, which is related to signal transduction in the TGF-β superfamily and the inhibition of BMP activity, and they have been thought to control follicular selection by balancing luteinizing hormone receptor (LHR) transcription (Yu et al., 2017;Xu et al., 2018). In addition, LDLR has been shown to take part in the biosynthesis of steroid hormones by affecting the absorption of cholesterol substrates via circulating lipoproteins (Ying et al., 2013). The MAPK1 (p42ERK2) isoform of MAPK is found in mammalian oocytes and the kinase is shown to be activated during meiotic maturation (Fan and Sun, 2004). Together, these circRNAs might be associated with the difference in the estrus phases. In previous studies, similar results were found for non-coding RNAs, including long noncoding RNA and microRNA enrichment analysis, suggesting that these circRNAs might act as sponges of endogenous RNAs (ceRNAs) to regulate gene expression (Kutzler et al., 2001). The results obtained from this study and previous reports suggested that the prolificacy trait and estrus induction are complex processes that includes numerous events, and that the circRNAs distributed in the two estrus phases in high and low fecundity goats are mainly involved in follicle growth, reproduction events, and other biological processes.
Prediction of circRNA-miRNA binding relationships was completed by the miRNA target gene prediction method. In total, binding was predicted for 56 miRNAs, including the miR-133 family (such as miR-133a-3p and miR-133b), miR-129-3p, and miR-21. Several of these miRNAs have been shown to be closely related to ovarian development and the prolificacy trait. Differential expression of miR-133 occurs during meiotic maturation of oocytes (Song et al., 2014). MiR-129-3p serves as a tumor suppressor by targeting BZW1 in ovarian cancer cells and the restoration of miR-129-3p might be a novel therapeutic strategy for ovarian cancer . MiR-21 promotes follicle development in the ovary during ovulation in mice (Carletti et al., 2010). In the chicken reproductive process, circRalGPS2 can compete with miRNA to modulate mRNA expression that regulates downstream target genes (Chiu et al., 2017). In our study, miR-129-3p was predicted to target chi-circC2H2orf76, chi-circDUS2, chi-circGPR137B, chi-circMYSM1, and chi-circSVEP1, miR-21-3p was predicted to targeting chi-circUTP20, and miR-133a-3p and miR-133b were predicted to targets to chi-circAKAP10. Overall, these differentially expressed circRNAs may regulate the function of the ovary, for estrus and different kidding numbers, by regulating miRNA, although specific mechanisms require further study.
As circRNAs have a loop-like structure consisting of exons, some circRNAs consist of an open reading frame and IRES. Therefore, researchers hypothesized that numerous circRNAs could be translated into protein in certain conditions (Nagarjuna et al., 2017). IRESfinder, CPC2, CNCI, CPAT, and PLEK software identified 18 circRNAs with coding potential. It should be noted that the biological function of these circRNAs, from the perspective of protein translation, should be explored further. In addition, we demonstrated the interactions between host genes in the ovary of Yunshang black goats. A total of 69 PPI pairs for the comparisons and several host genes were associated with the prolificacy trait. The relationships between these interactions require further study.

CONCLUSION
In summary, 283 (114 upregulated and 169 downregulated), 559 (299 upregulated and 260 downregulated), 449 (254 upregulated and 195 downregulated), and 314 (210 upregulated and 104 downregulated) DE circRNAs were identified in HF vs LF, HF vs HL, HL vs LL, and LF vs LL comparisons, respectively. Enrichment analysis showed that a number of host genes with DE circRNAs were enriched in oocyte maturation, ovarian steroidogenesis, oocyte meiosis, and the GnRH signaling pathway, which have important impacts on follicular growth and development. In addition, 192 DE circRNAs were found to target 56 miRNAs, 18 circRNAs had coding potential, and 69 PPI pairs were found. These results provide fundamental data on the molecular mechanisms underlying the prolific trait in goats.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The name of the repository and accession number can be found below: SRA, NCBI; PRJNA726268.

ETHICS STATEMENT
The animal study was reviewed and approved by All experimental procedures mentioned in the present study were approved by the Science Research Department (in charge of animal welfare issue)