Differential Expression and Functional Analysis of CircRNA in the Ovaries of Low and High Fecundity Hanper Sheep

Simple Summary Litter size is an important trait affecting reproductive capacity and breeding economics in meat sheep. Consequently, revealing its molecular mechanism helps us understand multiple lambs from the genetic perspective. In this study, we provide a genome-wide expression profile of circular RNAs (circRNAs) expression in Hanper sheep, which is a new breed of meat sheep raised by cross and self-group breeding for 15 years. In this study, ovarian circular RNAs and miRNAs associated with high and low fertility Hanper sheep are identified during the follicular and luteal phases of the estrous cycle, and their potential biological functions are predicted through Gene Ontology (GO), KEGG, GSEA, STEM, WGCNA analysis. Abstract Litter size is a considerable quality that determines the production efficiency of mutton sheep. Therefore, revealing the molecular regulation of high and low fertility may aid the breeding process to develop new varieties of mutton sheep. CircRNAs are the important factors regulating follicular development, but their mechanism role in the regulation of litter size in Hanper sheep is not clear. In the present study, ovarian tissues from the follicular (F) or luteal phase (L) of Hanper sheep that were either consecutive monotocous (M) or polytocous were collected. Then, we performed transcriptome sequencing to screen for differentially expressed circRNAs (DE-circRNAs) and elucidate their function. In total, 4256 circRNA derived from 2184 host genes were identified in which 183 (146 were upregulated, while 37 were downregulated) were differentially expressed in monotocous sheep in the follicular phase versus polytocous sheep in the follicular phase (MF vs. PF). Moreover, 34 circRNAs (14 were upregulated, while 20 were downregulated) were differentially expressed in monotocous sheep in the luteal phase versus polytocous sheep in the luteal sheep (ML vs. PL). This was achieved through DE-circRNAs function enrichment annotation analysis by GESA, GO, and KEGG, which function through the EGF-EGFR-RAS-JNK, TGF-β and thyroid hormone signaling pathway to affect the litter size of Hanper sheep in MF vs. PF and ML vs. PL. STEM results showed that MAPK signaling pathways play a key role in MF vs. PF and ML vs. PL. Through WGCNA analysis, AKT3 was a core gene in MF vs. PF and ML vs. PL. Moreover, competitive endogenous RNA (ceRNA) network analysis revealed the target binding sites for miRNA such as oar-miR-27a, oar-miR-16b, oar-miR-200a/b/c, oar-miR-181a, oar-miR-10a/b, and oar-miR-432 in the identified DE-cirRNAs.


Introduction
Litter size is an important trait affecting reproductive capacity and breeding economics in meat sheep. Consequently, revealing its molecular mechanism helps us understand multiple lambs from the genetic perspective. For example, FecB mutation, which is closely related to litter size, has a significant effect on small metabolic molecules in follicular fluid, making the follicular fluid stronger in its own antioxidant defense. From genomic, transcriptome, and proteomic levels, we found that FecB mutation promoted the expression time of the LHCGR gene on cumulus granulosa cells affecting the proliferation of cumulus granulosa cells [1]. Through miRNA and lncRNA analysis of FecB BB Small Tail Han sheep and FecB++ Dorset sheep, it was found the differentially expressed genes on cell apoptosis and proliferation were affected by gene mutation [2,3] and participated in progesterone mediated steroid biosynthesis and oocyte meiosis pathway [4]. Xist (loc101112291) and Gtl2 (loc101123329) were found to be highly expressed in Hanper sheep, suggesting the regulation of follicular development through methylation processes [5]. Mir-6404 and mir-29c are involved in the regulation of high and low fecundity in goats [6]. The expression of CTSB and CSSD in the ovarian tissue of Hu sheep with high and low fertility at the follicular stage was significantly up-regulated (p < 0.05) [7]. Comparative proteomics showed that taurine, sulfite metabolism, and apoptosis were involved in the regulation of ovarian function and closely related to follicular development and ovulation, revealing that FecB mutation was related to the high expression of mitochondrial oxidation related proteins [8]. In addition, the weighted gene co-expression network analysis showed that NR0B1, XLOC_041882, and MYH15 were mainly involved in the TGF-β signaling pathway, while NYAP1 and BCORL1 were significantly correlated with the oxytocin signaling pathway [9].
Ovary, the place where oocytes are produced and reproductive hormones are secreted, is directly related to follicular development and ovulation number. Previously, ncRNAs were regarded as a type of junk RNA that does not get translated. With the advancement of ultra-deep sequencing technology and bioinformatics, the functions of ncRNA became increasingly evident involving epigenetic and transcriptional regulations. Especially, the ncRNA mediated transcriptional regulation has gained tremendous attention. Different kinds of ncRNAs, namely the miRNA [10], LincRNA [11], and circRNA have been identified [12]. Among these, circRNAs are known to regulate various biological processes and organ function [13], such as pituitary [14], hypothalamus [15], muscle [16,17], spleen [18], and so on. However, the function of circRNAs in Hanper sheep ovaries is still unclear, in particular, their role in regulating the litter size is largely unknown. Cir-cRNAs are mainly formed by pre-mRNA reverse splicing [19] and are composed of an intron sequence [20], RNA binding protein [21], typical RNA splicing signal [22], and lasso precursor containing exons [23]. CircRNAs perform various biological functions, including acting as sponges for miRNAs [12], for example, testis-specific sex determination region Y (SRY) 9 sponges miR-138 [24]. CircDDXlO sponges miR-1301-3p or miR-4660 to control sirtuin 3 (SIRT3) expression, regulating the aging of human ovaries [25]. In bovine cumulus cells, circRNA_11396 functions as a molecular sponge for miR-187 to up-regulate the bone morphogenetic protein 2 receptor antibody (BMPR2) [26]. In porcine follicles, circINHA as the ceRNA of miR-10a-5p protects granulosa cells against apoptosis by upregulating the connective tissue growth factor CTGF [27]. Chi-circ_ 0008219 regulates ovine follicular development by combining miR-468-3p, miR-84c-5p, and miR-483a [28]. Likewise, there are numerous reports regarding the circRNAs mediated regulation of gene expression in several biological processes [29,30]. Moreover, circRNAs expression patterns could be tissue-specific [31] to transcriptionally regulate the parental genes [29]. Also, recent studies have demonstrated that circRNAs are translated in a cap-independent manner [32].
Hanper sheep is a new breed of meat sheep raised by Hebei Animal Husbandry and Veterinary Research Institute and Hebei Agricultural University experts, the population has been cross and self-group breeding for 15 years and purebred to 6-8 generations, and its genetic performance is stable. Therefore, in order to improve the economic benefits of breeding, it is particularly important to improve their fertility. However, the mean lamb Animals 2021, 11,1863 3 of 20 rate is 168%, which is higher than the 136% of Dorper sheep, but lower than the 260% of small-tailed Han sheep. The proportions of sheep that have 1, 2, or 3 or more lambs are 48.5%, 37.88%, and 13.62%, respectively. To explore effective breeding measures which could improve the regularity of fetal litter size in Hanper sheep, it is important to perform the differential expression analysis of ovarian circRNAs in the high-and low-fertility ewes. In this study, we collected follicles in the follicular phase and the corpus luteum tissues in the luteal phase from the ovaries of continuous monotocous (low reproduction) and polytocous (high reproduction) Hanper ewes respectively. Then, high-throughput sequencing was performed to screen the differentially expressed circRNAs (DE-circRNAs) and analysis functional enrichment of related genes. We aimed to reveal the candidate circRNAs and their target miRNAs that affect the fertility traits of Hanper sheep. Temporal and spatial expression models of circRNAs were constructed for better understanding and use of multi-lamb traits in sheep.

Ethics Approval
This study was authorized by the Animal Welfare and Ethical Committee of the Institute of Animal Sciences at the Hebei Agricultural University (IAS-CAAS, Baoding, China). All methods and procedures were performed on the experimental animals in accordance with the relevant guidelines (document approval number 2020070).

Selection of Experimental Animals and Sample Collection
Based on the lambing records of 1851 Hanper sheep (Hebei Liansheng Agricultural Development Co., Ltd. Zhuzhou, China), 12 unpregnant ewes aged 3-4 years were identified and divided into the monotocous group and the polytocous group, each of which consisted of six ewes. The ewes in the monotocous group lambed only one kid per litter within more than three consecutive litters. The ewes in the polytocous group lambed no less than three kids per litter within more than three consecutive litters. The natural estrus cycle of the 12 ewes was identified and the ewes were observed to have two normal natural estrus cycles in succession, then the bilateral ovaries of the ewes in the monotocous group (M) and the polytocous group (P) in the follicular phase (F) were collected at 36 h after the third natural estrus, denoted as MF and PF, respectively (three ewes in each group, six ewes in total). At 216 h (ninth day) after the third natural estrus, the bilateral ovaries of the ewes in the monotocous group (M) and the polytocous group (P) in the luteal phase (L) were collected, denoted as ML and PL, respectively (three ewes in each group, six ewes in total) ( Table 1). The mesangium around the ovary was removed by disinfected scissors. These follicles whose diameter was greater than 3 mm in the follicular phase ovaries from the same individual were separated and placed separately in 2 mL RNase-free cryopreserved tubes. The luteal samples in the luteal phase ovaries from the same individual were separated and put separately into 2 mL RNase-free cryopreserved tubes. The tubes with the samples were immediately frozen in liquid nitrogen before transferring to the laboratory. After that, the samples were moved and stored at −80 • C for later use.  3E202  3  3  5  6T01  3  4  2  6T06  3  3  4  3  CE104  3  3  3  3E161  3  3  4  CE103  3  4  3  4 Monotocous group

RNA Extraction and Library Construction
The follicles or luteum of 12 individual sheep tissues were frozen to grind into powder separately. Then, the total RNA in every sample was extracted with Trizol reagent (Thermo Fisher Scientific, Waltham, MA, USA). The concentration and purity of RNA were estimated by the RNA Nano 6000 Assay Kit with BioAnalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) and 1% agarose gel electrophoresis, respectively.
A total amount of 5 µg RNA per sample was used as input material for the RNA sample preparations. Firstly, ribosomal RNA was removed by Epicentre Ribozero™ rRNA Removal Kit (Epicentre, Madison, WI, USA), and rRNA free residue was cleaned up by ethanol precipitation. Subsequently, linear RNA was digested with 3 U of RNase R (Epicentre, Madison, WI, USA) per µg of RNA. Sequencing libraries were generated by NEBNext ® Ultra™ Directional RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) following the manufacturer's recommendations. Briefly, fragmentation was carried out using divalent cations under elevated temperature in the NEBNext First Strand Synthesis Reaction Buffer (5X). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNaseH-). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. In the reaction buffer, dNTPs with dTTP were replaced by dUTP. The remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3 ends of DNA fragments, the NEBNext Adaptor with hairpin loop structure was ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 150~200 bp length, the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, IL, USA). Then 3 µL USER Enzyme (NEB, Ipswich, MA, USA) was used with size-selected, adaptor-ligated cDNA at 37 • C for 15 min followed by 5 min at 95 • C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index Primer. At last, products were purified (AMPure XP system) and library quality was assessed using the Agilent Bioanalyzer 2100 system.

Sequence Mapping and CircRNA Prediction
The pure reads were mapped to the reference genome (Oar_v3.1) using the HiSAT2 alignment method which is a valid and quick tool to identify circRNA [33]. To ascertain the reliability of circRNAs, the BWA-MEM algorithm was applied for sequence segmentation and comparison, and the SAM file was used to find PCC (paired cut fragments), PEM (paired-end maps) sites, and GT-AG splicing signal [34]. Finally, the sequences with connection sites were reanalyzed using the dynamic programming algorithm, and circRNAs were annotated against the circBase database. CircRNAs without annotation were defined as novel circRNAs. Also, the type of exon chromosome and length distribution of the identified circRNAs were statistically analyzed.

Differential Expression Analysis of CircRNAs
Spliced Reads per Billion Mapping (SRPBM, defined as the number of circular reads/ number of mapped reads read × length) was used to estimate the expression of circR-NAs [35]. The DE-circRNAs were authenticated by the DEseq2 software package [36]. Since the chosen populations were identical sheep breeds and belonged to the same farm (Hebei Lian sheng Agricultural Development Co., Ltd. Baoding, Hebei, China), their growth characteristics were the same, only the number of litters was different. circRNAs with a fold change ≥ 1.5 and p-value < 0.05 between the two groups (MF vs. PF, ML vs. PL) were identified as the DE-circRNAs.

STEM Cluster and Weighted Gene Co-Expression Network Analysis (WGCNA) Analysis of CircRNAs
For clustering, comparison, and visualization, Short Time-series Expression Miner (STEM) software was used to display the trends in these two stages. All genes were divided into eight modules according to the expression pattern, and the p-value of each module was calculated using the Permutation Test method. The identified circRNAs with pvalue < 0.05 were considered enriched. Furthermore, WGCNA analysis was performed [37]. According to the formula Amn = [(1 + Smn)/2] β, the similarity matrix was converted into an adjacency matrix, and then the adjacency matrix was converted into a topological overlap matrix (TOM). Meanwhile, the function diss Tom = 1-TOM was employed to invert the topological overlap matrix to obtain the dissimilarity matrix. The purpose of this exercise was to eliminate errors from background noise and pseudo-correlation. Finally, the function "hclust" was used to perform hierarchical clustering of the dissimilar matrix, and the Dynamic Tree Cut method was used to generate a cluster tree cut. This process merged genes with similar expression patterns on the same branch, and each branch represented a co-expression module, while different colors represented different modules. Correlation analysis and clustering were performed according to the FPKM values of the distinct genes. Next, we determined the target gene module through feature vector analysis and performed cluster analysis on the representative genes in all modules-module eigengene (ME). A higher correlation between ME indicates the higher relevance of the module.

Functional Analysis of CircRNAs in Follicular and Luteal Period
The host genes enrichment annotation of DE-circRNA was determined by the David and Kyoto Encyclopedia of Genes and Genomes (KEGG) database. Using the GOseq R software package, ClueGO, and Gene Set Enrichment Analysis (GSEA), Gene Ontology (GO) enrichment analysis was performed on the DE-circRNAs of host genes [38]. p-value < 0.05 denoted the significance of GO and KEGG pathways and their functional enrichment.

Analysis of MiRNA Target Sites
To further inquire into the functional roles of miRNA, putative targets were identified by the combination of three software, Mireap, Miranda (v. 3.3a), and TargetScan (v. 7.0) [39]. Cytoscape (v. 3.5.1) used the circRNA-miRNA negative interactions to construct the circRNA-miRNA networks.

Validation of CircRNA Expression
Nine DE-circRNAs (FDR ≥ 1.5, p-value < 0.05) and three miRNAs were stochastically selected for the follow-up verification. The details of circRNA-specific reverse splicing and miRNA primer sequences are shown in Table 2. Quantitative real-time polymerase chain reaction (qRT-PCR) was performed in the Roche lightcycler 480ii system (Roche Applied Science, Mannheim, Germany) using SYBR green real-time PCR master mix (TaKaRa, Dalian, China) following the manufacturer's instructions. The reaction was executed at 95 • C for 3 min, followed by 40 cycles for 10 s at 95 • C and 30 s at 60 • C, respectively. GADPH and U6 were used as internal standards for circRNA and miRNA housekeeping genes, respectively. The relative change of gene expression among the control and experimental groups was calculated by the 2 −∆∆Ct method [40]. Also, the cyclization of circRNA was verified by Sanger sequencing [41].
We used Transcripts Per Million (TPM) map reads to normalize circRNA abundance. Among the normalized circRNAs, the transcript expression levels of each group (MF, PF, ML, PL) were kept similar. As shown in Supplementary Figure S1B, the TPM peak densities of MF, PF, ML, and PL ranged from 1-3 per million counts. Among the circRNAs in each group (MF, PF, ML, PL), the number of exons mainly ranged from 2 to 4, with a total of 3406 (Supplementary Figure S1C). Among them, the length of circRNA was mainly 300-400 nt (1415 in total) (Supplementary Figure S1D). The numbers in each column denote the expression level of circRNA in descending order. Red and green denote the higher and lower expression levels of circRNAs, respectively. Also, the lighter the color, the lower is the circRNA TPM value. Note: Normalized expression level in TPM = (read count × 1,000,000)/libsize (libsize is the sum of circRNA read count). MF: Monotocous sheep in the follicular phase, PF: Polytocous sheep in the follicular phase, ML: Monotocous sheep in the luteal phase, PL: Polytocous sheep in the luteal phase.

Diversity of CircRNA Isoforms and Circularization Forms
Since sheep data are not yet available in circBase databases, the identified circRNA can be considered new (Table S1). The 4256 identified circRNA were derived from 2184 loci, and most genes (1872, 80.11%) produced one to two circRNAs, 19.89% (465) produced 3-15 different circRNAs (Supplementary Figure S2A). For instance, the potential transforming growth factor binding protein 1 (LTBP1) can produce 15 different circRNAs (Supplementary Table S1). Notably, the circRNA spliced at the identical site contains a splicing signal, which usually consists of a GT dinucleotide (splice donor) at the 5 end of the main spliceosome and an AG dinucleotide (splice acceptor) at the 3 end. Given an AG acceptor splice site, there are usually multiple GT donor sites. In this study, transforming growth factor-beta receptor 3 (TGFBR3) produced three subtypes of circRNAs (novel_circ_0001803, novel_circ_0001805, and novel_circ_0001806; Figure S2B) with 3326, 8883, and 14,600 bp, respectively. Similarly, many circRNAs have different acceptors but share the same donor. In this study, the host gene AKT serine/threonine kinase 3 (AKT3) produced three subtypes of circRNAs that have three different acceptors but share a single donor NC_019469.2_31610710 (Supplementary Figure S2C). This showed that there were alternating and cross splicing events, where the parallel splicing site may join in multiple forward and reverse splicing reactions, or in coterminous exons distally placed exons.
In order to uncover the expression profiles of different circRNAs derived from the same host gene, we explored the expression abundance of circRNA subtypes at different developmental stages. Among the three subtypes of TGFBR3 circRNA, novel_circ_0001805 was highly expressed in each group, but novel_circ_0001803 and novel_circ_0001806 were missed in the PF group (Supplementary Figure S2D). Likewise, for AKT3 circRNA subtypes, novel_circ_0011336 was highly expressed in the MF and PL groups, while novel_circ_0011344 was highly expressed in the MF and PF groups ( Supplementary Figure S2E). In brief, these results implied that reverse splicing promoted multiformity and functional complexity in circRNAs by generating the subtypes.

Identification of DE-CircRNAs
The distribution and expression of circRNAs on each chromosome are shown by a chord diagram ( Figure 1A). In general, circRNAs are mainly distributed in chromosomes 1, 2, and 3. Ovarian DE-circRNAs in the polytocous and monotocous sheep were identified by  Table S2).

Identification of DE-CircRNAs
The distribution and expression of circRNAs on each chromosome are shown by chord diagram ( Figure 1A). In general, circRNAs are mainly distributed in chromosom 1, 2, and 3. Ovarian DE-circRNAs in the polytocous and monotocous sheep were iden fied by DEseq2, based on criteria FDR ≥ 1.5 and p-value < 0.05. We found that there we 146 up-and 37 down-regulated circRNAs in MF vs. PF groups. Likewise, there were up-and 20 down-regulated circRNAs in ML vs. PL groups, ( Figure 1B) (Supplementa Table S2).

Functional Analysis of CircRNAs in Sheep of Varying Fecundity in the Follicular and Luteal Phases
Based on the assumption that the functions of circRNA and host genes are know we analyzed the pathways enrichment of the host genes producing DE-circRNAs to p dict their potential functions during follicular development using GESA, GO, and KEG

Functional Analysis of CircRNAs in Sheep of Varying Fecundity in the Follicular and Luteal Phases
Based on the assumption that the functions of circRNA and host genes are known, we analyzed the pathways enrichment of the host genes producing DE-circRNAs to predict their potential functions during follicular development using GESA, GO, and KEGG. As shown in Figure 2A and Supplementary Figure S3, DE-circRNAs in MF vs. PF were mainly involved in the cellular response to endogen stimulus, phosphate-containing compound metabolic process, positive regulation of molecular function, protein complex subunit organization, tissue development, and metabolic process. The key cellular components were mitochondrion and vacuum, which performed molecular functions such as adenyl nucleotide binding, enzyme binding, transferase activity transferring phosphorus-containing groups, and carbohydrate derivative. The above enriched biological processes were mainly related to PTGDR, CYP11A1, INSR, SMAD1, ACVR2A, RPS6KA001, COL12A1, PGR, and LTBP1 involved in regulating the follicular development of Hanper sheep ( Figure 2B). As there were only 33 DE-circRNAs in ML vs. PL, GESA analysis of host genes could not be performed. Go Seq R enrichment analysis demonstrated that DE-circRNAs were mainly related to the biological processes include single-organism developmental process, acting on amino acids and derivatives, and cell surface composition (Supplementary Figure S4). LTBP1, ITGAV, SNX13, CAAP1, and COL14A1 are they key DE-circRNAs that are involved in the regulation of follicle development (Supplementary Table S3). As shown in Figure 2A and Supplementary Figure S3, DE-circRNAs in MF vs. PF were mainly involved in the cellular response to endogen stimulus, phosphate-containing compound metabolic process, positive regulation of molecular function, protein complex subunit organization, tissue development, and metabolic process. The key cellular components were mitochondrion and vacuum, which performed molecular functions such as adenyl nucleotide binding, enzyme binding, transferase activity transferring phosphoruscontaining groups, and carbohydrate derivative. The above enriched biological processes were mainly related to PTGDR, CYP11A1, INSR, SMAD1, ACVR2A, RPS6KA001, COL12A1, PGR, and LTBP1 involved in regulating the follicular development of Hanper sheep ( Figure 2B). As there were only 33 DE-circRNAs in ML vs. PL, GESA analysis of host genes could not be performed. Go Seq R enrichment analysis demonstrated that DE-circRNAs were mainly related to the biological processes include single-organism developmental process, acting on amino acids and derivatives, and cell surface composition (Supplementary Figure S4). LTBP1, ITGAV, SNX13, CAAP1, and COL14A1 are they key DE-circRNAs that are involved in the regulation of follicle development (Supplementary  Table S3). KEGG enrichment analysis of host genes of DE-circRNAs in MF vs. PF revealed the correlation between "Environmental Information Processing" and biological processes of "Organizational MAL Systems" (Figure 3A). Among these, the TGF-β signaling pathway, Rap1 signaling pathway, and cAMP signaling pathway involved in Environmental Information Processing were significantly enriched. Ovarian steroidogenesis, progesterone mediated oocyte maturation in Organizational MAL Systems, and the neurotrophin signaling pathway (Supplementary Table S4) were also significantly enriched. Finally, further analysis of the host genes involved in the enrichment pathway closely related to reproduction revealed that EGF-EGFR-RAS-JNK and TGF-β signaling pathways had important regulatory roles for sheep with different fertilities at the follicular stage (Figure KEGG enrichment analysis of host genes of DE-circRNAs in MF vs. PF revealed the correlation between "Environmental Information Processing" and biological processes of "Organizational MAL Systems" (Figure 3A). Among these, the TGF-β signaling pathway, Rap1 signaling pathway, and cAMP signaling pathway involved in Environmental Information Processing were significantly enriched. Ovarian steroidogenesis, progesterone mediated oocyte maturation in Organizational MAL Systems, and the neurotrophin signaling pathway (Supplementary Table S4) were also significantly enriched. Finally, further analysis of the host genes involved in the enrichment pathway closely related to reproduction revealed that EGF-EGFR-RAS-JNK and TGF-β signaling pathways had important regulatory roles for sheep with different fertilities at the follicular stage ( Figure 3B). KEGG enrichment analysis of DE-circRNAs host genes in ML vs. PL revealed that these were mainly enriched in the Organic Systems process ( Figure 3C). As shown in Table 4, ITGAV participated in thyroid hormone and PI3K-Akt signaling pathway regulating the fecundity of Hanper sheep in ML vs. PL. Also, COL14A1 mainly controlled the fecundity levels in Hanper sheep through protein digestion and absorption. In short, the enrichment processes in MF vs. PF and ML vs. PL are mainly related to the reproduction regulation genes involving COL14A1, ITGAV, CYP11A1, ACVR2A, and SMAD1. 3B). KEGG enrichment analysis of DE-circRNAs host genes in ML vs. PL revealed that these were mainly enriched in the Organic Systems process ( Figure 3C). As shown in Table 4, ITGAV participated in thyroid hormone and PI3K-Akt signaling pathway regulating the fecundity of Hanper sheep in ML vs. PL. Also, COL14A1 mainly controlled the fecundity levels in Hanper sheep through protein digestion and absorption. In short, the enrichment processes in MF vs. PF and ML vs. PL are mainly related to the reproduction regulation genes involving COL14A1, ITGAV, CYP11A1, ACVR2A, and SMAD1.

STEM Cluster Analysis of DE-CircRNAs
STEM software was used to perform trend analysis and 20 distribution chromatograms were constructed ( Figure 4A). Profile 1 and Profile 11 which contain 945 and 895 cirRNAs host genes are two prominent profiles ( Figure 4A). ClueGO results showed that host genes in profile 1 are mainly related to organelle organization, cellular macromolecule metabolic process, and anatomical structure morphogenesis biological processes (Supplementary Figure S5A, Supplementary Table S5). The genes of Profile 11 are mainly enriched in enzyme binding, and cellular macromolecule metabolic process ( Figure S5B, Table S6). In addition, KEGG enrichment analysis using GESA showed that profile 1 and profile 11 are mainly enriched in the MAPK signaling pathway ( Figure 4B). Moreover, Profile 1 mainly regulated follicular development in MF vs. PF and ML vs. PL through MAP3K5, SOS1, AKT3, and SOS2 in the MAPK signaling pathway. Profile 11 mainly functions through the MAPK signaling pathway involving IL1RL1, SOS2, AKT3, BRAF, MAPK1, and PAK2 to regulate follicle development in MF vs. PF and ML vs. PL ( Figure 4C). Figure S6 revealed that the MAPK signaling pathways had important regulatory effects on sheep with different fertility and different stages in different circRNAs host genes of profile 1 and profile 11.
STEM software was used to perform trend analysis and 20 distribution chromatograms were constructed ( Figure 4A). Profile 1 and Profile 11 which contain 945 and 895 cirRNAs host genes are two prominent profiles ( Figure 4A). ClueGO results showed that host genes in profile 1 are mainly related to organelle organization, cellular macromolecule metabolic process, and anatomical structure morphogenesis biological processes (Supplementary Figure S5A, Supplementary Table S5). The genes of Profile 11 are mainly enriched in enzyme binding, and cellular macromolecule metabolic process ( Figure S5B, Table S6). In addition, KEGG enrichment analysis using GESA showed that profile 1 and profile 11 are mainly enriched in the MAPK signaling pathway ( Figure 4B). Moreover, Profile 1 mainly regulated follicular development in MF vs. PF and ML vs. PL through MAP3K5, SOS1, AKT3, and SOS2 in the MAPK signaling pathway. Profile 11 mainly functions through the MAPK signaling pathway involving IL1RL1, SOS2, AKT3, BRAF, MAPK1, and PAK2 to regulate follicle development in MF vs. PF and ML vs. PL ( Figure  4C). Figure S6 revealed that the MAPK signaling pathways had important regulatory effects on sheep with different fertility and different stages in different circRNAs host genes of profile 1 and profile 11.

Construction of WGCNA Network
In this study, 2212 high-expression genes were obtained by filtering out the lowexpression genes. The results shown in Figure 5A indicate that there was no outlier in the gene cluster tree among all samples and therefore no further processing was required. R 2 > 0.85, β = 14 was used as the soft threshold to construct the co-expression module ( Figure 5B). After determining the soft threshold β = 9, Highly correlated genes were allocated to the same module (Figure 5C), and a total of nine co-expression modules were obtained. The genes in a module were assigned to a cluster based on their expression levels. Also, the genes with higher clustering were assigned to a module. The turquoise module contained the largest number of host genes (271) while the pink module contained the smallest number of host genes (20). The grey module is a set of genes that were not assigned to other modules ( Figure S7). Furthermore, ME correlation analysis between two modules was carried out ( Figure 5D). It was found that the ME correlation between the Brown and Black modules was 0.74, while it was 0.73 between the Black and Green modules. The correlation between the eigenvectors of nine modules and the sample sources was calculated ( Figure 5E Table S7). Likewise, KEGG analysis revealed that the brown module genes are mainly enriched in the WNT signaling pathway. Green module genes are enriched in focal adhesion. Yellow module genes are enriched in the mTOR signaling pathway and PI3K-Akt signaling pathway. Turquoise module genes are enriched in the estrogen signaling pathway and insulin resistance (Table 5). In summary, host genes AKT3 and were found to play key regulatory roles in multiple modules. The circRNAs corresponding to the above host genes are shown in Table S1.  modules, which are color-coded. (D) ME correlation between different modules. Black dots represent circRNAs, and the higher the value, the higher the correlation. The more and larger the asterisk, the greater the pairwise correlation between the modules. (E) Association analysis of gene co-expression network modules with tissues. Each row corresponds to a module, whose name is displayed on the left, and each column corresponds to a particular sample. The color of each cell at the row-column intersection indicates the correlation coefficient between the module and the sample. High correlations between specific modules and samples are shown in red.  (D) ME correlation between different modules. Black dots represent circRNAs, and the higher the value, the higher the correlation. The more and larger the asterisk, the greater the pairwise correlation between the modules. (E) Association analysis of gene co-expression network modules with tissues. Each row corresponds to a module, whose name is displayed on the left, and each column corresponds to a particular sample. The color of each cell at the row-column intersection indicates the correlation coefficient between the module and the sample. High correlations between specific modules and samples are shown in red.

Discussion
Fetal lambing traits are directly related to the follicle development of sheep ovaries [8]. Sheep fecundity, a complicated biological process, is tightly regulated by coordinated genes, covering coding and non-coding RNA. Notably, in the mammalian genome, 70-90% of non-coding RNAs are universally transcribed [42], indicating that they constitute the majority of the transcriptome. Transcriptome sequencing can uncover the genes related to reproductive traits and elucidate the underlying molecular mechanisms. In this study, we identified the DE-circRNAs in the ovaries of Hanper sheep and explored their target genes which may regulate ovarian development. We studied the circRNA genome characteristics in the follicles and corpus luteum of monotocous and polytocous Hanper sheep and found that most circRNAs had three or four exons, which was in agreement with previous findings in sheep skeletal muscle and uterus [8,17]. Besides, a total of 4257 circRNAs were identified between MF vs. PF and ML vs. PL. Interestingly, we found a total of 183 and 34 differentially expressed circRNAs in MF vs. PF. This shows a large difference in the expression patterns of circRNAs in ML vs. PL, indicating cell type-dependent characteristics of circRNA expression and differential regulation of circularization events in different cells. A similar phenomenon was also observed in other studies [42]. It is evident that cumulus cells are transcriptionally highly active in the process of follicular development. Follicular development could produce much more transcripts relative to the luteal tissues, and circRNA biogenesis mainly depends on the back-splicing

Discussion
Fetal lambing traits are directly related to the follicle development of sheep ovaries [8]. Sheep fecundity, a complicated biological process, is tightly regulated by coordinated genes, covering coding and non-coding RNA. Notably, in the mammalian genome, 70-90% of non-coding RNAs are universally transcribed [42], indicating that they constitute the majority of the transcriptome. Transcriptome sequencing can uncover the genes related to reproductive traits and elucidate the underlying molecular mechanisms. In this study, we identified the DE-circRNAs in the ovaries of Hanper sheep and explored their target genes which may regulate ovarian development. We studied the circRNA genome characteristics in the follicles and corpus luteum of monotocous and polytocous Hanper sheep and found that most circRNAs had three or four exons, which was in agreement with previous findings in sheep skeletal muscle and uterus [8,17]. Besides, a total of 4257 circRNAs were identified between MF vs. PF and ML vs. PL. Interestingly, we found a total of 183 and 34 differentially expressed circRNAs in MF vs. PF. This shows a large difference in the expression patterns of circRNAs in ML vs. PL, indicating cell type-dependent characteristics of circRNA expression and differential regulation of circularization events in different cells. A similar phenomenon was also observed in other studies [42]. It is evident that cumulus cells are transcriptionally highly active in the process of follicular development. Follicular development could produce much more transcripts relative to the luteal tissues, and circRNA biogenesis mainly depends on the back-splicing of pre-mRNA transcribed in a cell [43]. We thus speculated that a larger amount of circRNAs produced in follicular tissues could be attributed to their higher transcriptional activity. About 20% of these circRNAs could produce more than one type of circRNA, which was lower than cows (40%) and honeybees (61%) [44]. This suggested that alternative splicing was indeed the main event in the formation of circRNAs [41]. We found that several known genes related to follicular development produced more than one circRNA. For instance, TGFBR3 and AKT3 had three predictive circRNAs, and LTBP1 could produce up to 15 circRNAs. It seems that a parental gene can produce several circRNA subtypes, even more than 10 types in some cases. However, only one or two subtypes are prominent in general. Nevertheless, most other subtypes with low expression may play a specific function in the follicular development of Hanper sheep.
circRNAs are known to enhance the transcription of host genes by interacting with U1 small ribonucleoproteins (snRNPs) and RNA polymerase II [45]. Therefore, identifying the potential function of host genes associated with circRNA may reveal their function. In addition to the common host genes that affect reproductive performance such as SMAD1, CYP11A1, INSR, and PGR, we also identified potential members of the TGF-β binding protein family, which encoded for transfer growth factor binding protein 1 (LTBP1), located on sheep chromosome 3. Sheep LTBP1, with 29 exons and 28 introns, serves as the target of TGF-β and is usually released as an inactive complex from the cell. Studies show that TGF-β releases from epithelial cells requiring the proteolytic processing of LTBP1, thereby a cell cannot produce TGF-β in the absence of LTBP1 [46]. This may affect cell maturation and differentiation [47]. Moreover, circRNAs also function as miRNA sponges and are involved in miRNA-mediated post-transcriptional regulation. However, only a limited number of circRNAs contain enough binding sites to bind specific miRNAs. For instance, circRNA ciRS-7 with more than 70 conserved binding sites acts as a sponge for miR-7 to regulate brain cell development [48]. On the contrary, recent evidence suggests that circRNA as a miRNA sponge does not require a large number of putative targeting sites. The circular RNA produced from exon 2 of the HIPK3 gene has only two binding sites for miR-124 which is enough for inhibiting its activity [49]. Liang et al. show that circSCIN can be sponged with miR-133 and miR-148a/b, and its regulation determines the reproductive performance of Meishan and large white pigs [50]. CirceGFR, being a competitive endogenous RNA of miR-125a-3p, regulates FYN expression in mouse ovarian granulocyte cells [51]. Here, based on miRNA target site prediction analysis, we found that oar-miR-27a could be a target of novel_circ_0005497. miR-27a is known to inhibit the proliferation and promote apoptosis of granulosa cells in mice. This is achieved by reducing the key enzyme of estrogen synthesis, named aromatase [52]. In conclusion, circLTBP1 is a potential competitive endogenous RNA (ceRNA) that may regulate gene transcription, and affect fertility traits as well as follicular development in sheep.
In addition to miRNA regulation, circRNA might have additional roles. Recent studies showed that circRNAs could also directly regulate protein synthesis through mRNA. Also, some of them can be converted into proteins for insertion via IRES (internal ribosome entry sites) [53]. Several factors and pathways, including TGFβ, PI3K-AKT, Wnt/β-catenin, and TGFβ-SMAD signaling pathways are related to follicular growth and development. Interestingly, Tao et al. showed that the host genes of ovarian follicles circRNAs participate in the generation of the ovarian corpus callosum and p53 signal transduction before ovulation in goats [54]. CircRNA host genes are also involved in the production of ovarian steroids and regulating their signals. These signals are essential for follicle growth, oocyte maturation, and ovulation. In this study, using STEM, WGCNA, and differential circRNA functional enrichment analysis, it was found that the signaling pathways of MAPK, Wnt, mTOR, PI3K-Akt, TGF-β, Rap1, and ovarian steroidogenesis were involved in circRNAs mediated regulation of sheep reproduction.

Conclusions
In this study, the circRNA expression profile of high and low litter size Hanper sheep was established. A total of 4752 circRNAs were identified, of which 183 circRNAs were differentially expressed in MF vs. PF and 34 circRNAs were differentially expressed in ML vs. PL. In total, 17 key circRNAs involved in the development of ovine follicles were identified, and they function through the signaling pathways of MAPK, WNT, mTOR, PI3K-AKT, TGF-β, RAP1, EGF-EGFR-RAS-JNK, and ovarian steroidogenesis to regulate the development of sheep follicles. Also, the target miRNA sites in circRNAs for oar-miR-27a, oar-miR-16b, oar-miR-200a/b/c, oar-miR-181a, oar-miR-10a/b, oar-miR-432 were revealed through competitive endogenous RNA network analysis. It seems that circRNAs mainly function as sponges of several reproductive-related miRNAs, which in turn affects the fetal litter traits of Hanper sheep.  Figure S4. The top 20 GO enrichments of differential expressed circRNAs (DE-circRNAs) in polytocous sheep in the follicular phase versus monotocous sheep in the luteal phase (ML vs. PL); Figure S5. The network plot showed the biological process GO terms produced by ClueGO in the profiles generated by STEM; Figure S6. MAPK signaling pathways at different stages in different circRNAs host genes of profile 1 and profile 11; Figure S7. Distribution of the number of genes in modules. Figure S8. Sequences information.   Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.