Genome wide identification and characterization of fertility associated novel CircRNAs as ceRNA reveal their regulatory roles in sheep fecundity
Journal of Ovarian Research volume 16, Article number: 115 (2023)
Reproductive traits play a vital role in determining the production efficiency of sheep. Maximizing the production is of paramount importance for breeders worldwide due to the growing population. Circular RNAs (circRNAs) act as miRNA sponges by absorbing miRNA activity through miRNA response elements (MREs) and participate in ceRNA regulatory networks (ceRNETs) to regulate mRNA expression. Despite of extensive research on role of circRNAs as miRNA sponges in various species, their specific regulatory roles and mechanism in sheep ovarian tissue are still not well understood. In this study, we performed whole genome sequencing of circRNAs, miRNA and mRNA employing bioinformatic techniques on ovine tissues of two contrasting sheep breeds "Small tail Han (X_LC) and Dolang sheep (D_LC)", which results into identification of 9,878 circRNAs with a total length of 23,522,667 nt and an average length of 2,381.32 nt. Among them, 44 differentially expressed circRNAs (DECs) were identified. Moreover, correlation between miRNA-mRNA and lncRNA-miRNA provided us with to prediction of miRNA binding sites on nine differentially expressed circRNAs and 165 differentially expressed mRNAs using miRanda. miRNA-mRNA and lncRNA-miRNA pairs with negative correlation were selected to determine the ceRNA score along with positively correlated pairs from lncRNA and mRNA network. Integration of ceRNA score and positively correlated pairs exhibit a significant ternary relationship among circRNAs-miRNA-mRNA demonestrated by ceRNA, comprising of 50 regulatory pairs sharring common nodes and predicted potential differentially expressed circRNAs-miRNAs-mRNAs regulatory axis. Based on functional enrichment analysis shortlisted key ceRNA regulatory pairs associated with reproduction including circRNA_3257-novel579_mature-EPHA3, circRNA_8396-novel130_mature-LOC101102473, circRNA_4140- novel34_mature > novel661_mature-KCNK9, and circRNA_8312-novel339_mature-LOC101110545. Furthermore, expression profiling, functional enrichments and qRT-PCR analysis of key target genes infer their implication in reproduction and metabolism. ceRNA target mRNAs evolutionary trajectories, expression profiling, functional enrichments, subcellular localizations following genomic organizations will provide new insights underlying molecular mechanisms of reproduction, and establish a solid foundation for future research.
Graphical abstract summarizing the scheme of study
Reproductive traits, such as prolificacy, multiple births, and fecundity, are threshold and crucial for the profitability of sheep production but have low heritability rate [1,2,3]. Therefore, developing reproductive traits has become a significant interest in sheep breeding , as it can increase the production of sheep-derived products, ranging from wool to lamb meat . The genetic background of a ewe primarily determines its fertility, which depends on various factors, including ovulation rate, estrus cycle, and litter size, with the fecundity of the ewe playing a central role, mainly concentrated on the ovary [5,6,7]. Of note, several studies have demonstrated that follicular development and ovulation can be controlled by the set of single or multiple genes with major and minor effects [8, 9]. BMPR1, BMP15, GDF9, and B4GALNT2 were well known for their contribution to high prolificacy [3, 10, 11]. However, lack of knowledge of the genetics behind the complicated physiological mechanisms governing fertility hinders selection accuracy.
Reproductive traits, such as prolificacy, multiple births, and fecundity, are crucial for sheep production profitability but have a low heritability rate [1,2,3]. Therefore, developing reproductive traits has become a significant interest in sheep breeding , as it can increase the production of sheep-derived products, ranging from wool to lamb meat . The genetic background of a ewe primarily determines its fertility, which depends on various factors, including ovulation rate, estrus cycle, and litter size, with the fecundity of the ewe playing a central role, mainly concentrated on the ovary [5, 7]. Several studies have shown that follicular development and ovulation can be controlled by the set of single or multiple genes with major and minor effects, including BMPR1, BMP15, GDF9, and B4GALNT2, which are known for their contribution to high prolificacy [3, 10, 11]. However, the lack of knowledge regarding the genetics behind the complicated physiological mechanisms governing fertility hinders selection accuracy.
Non-coding RNAs (ncRNAs), including microRNAs (miRNAs), long non-coding RNAs (lncRNAs), and circular RNAs (circRNAs), widely participate in post-transcriptional gene expression regulation and several biological processes [12,13,14]. Circular RNAs (circRNAs) are closed long non-coding RNAs that are formed by the back-splicing of exons from a single pre-mRNA, in which 5' and 3' termini are covalently bonded . CircRNAs are extensively conserved, widely expressed in mammalian cells and often exhibit cell type- or tissue-specific expression patterns [16, 17]. Increasing evidences suggested that circRNAs have been implicated in the regulation of a number of biological processes and participate in various pathophysiology processes . The biological role of circRNA is serving as a decoy for microRNA or competing ceRNA [19,20,21]. This is the most widely studied function of circRNAs that sponge and inhibit relevant miRNA by complementary base paring. They can release the inhibitory effect of miRNAs on their targeted mRNAs, thereby promoting the transcription and protein translation of target genes. Salmena et al. introduced a systematic "ceRNA-hypothesis", in which circRNAs would actively communicate with miRNA through a "ceRNA language" composing vast network of interactions and controlling their expression level . While, short-stranded RNA called microRNA (approximately 22nt) suppresses the expression of mRNA by either inhibiting protein translation or destroying the target mRNA .The ceRNA hypothesis holds that ceRNA regulates the expression of transcripts by competing with mRNA for the same microRNA response elements (MREs). Regardless of whether it encodes a protein, RNA transcripts can compete with each other for binding to microRNA and regulate each other, thus forming a huge ceRNA network (ceRNETs) .
Several circRNAs have been identified via high-throughput RNA sequencing in various tissues, such as the hypothalamus , pituitary , uterus , and oviduct . Many recent studies have shown that circRNAs may have unique and important functions during embryonic development, and their expression is the most abundant and complex in the cerebral cortex on the 60th day of pregnancy . Study of Zhang et al. determined key hypothalamic circRNAs in sheep and found oar_circ_0012110 can reduce the level of oar-miR-665-3p by acting as miRNA sponge, in turn, regulate the expression of DIAPH1 gene which function in mediating GnRH pulsatile secretion . Another example of novel_circ_0000417 that competitively bind with miR-136 and promote the apoptotic process in human granulosa cells . In ovine oviduct, novel_circ_0017815 determined as a molecular sponge for oar-miR-181a and regulate the estradiol synthesis and follicular apoptosis . ceRNA study on porcine follicles have highlighted the role of circINHA as a sponge for miR-10a-5p, thus upregulates the expression of connective tissue growth factor (CTGF) and prevent granulosa cells from apoptosis . Chi-circ_ 0,008,219 control the ovine follicles growth and development through competitively bind with miR-468-3p, miR-84c-5p, and miR-483a [32, 33]. CircDDXlO regulates the expression of sirtuin 3 (SIRT3) by inhibiting the activity of miR-1301-3p or miR-4660 and eventually slow down the ovarian aging in humans . CircRNA 11,396 acts as a molecular sponge for miR-187 in bovine cumulus cells to upregulate the expression of BMPR2 (bone morphogenetic protein 2 receptor) antibody . Above discussion suggested that circRNAs have a specific effect on fetal development. Despite substantial study into the prediction, quantification, and annotation of circRNAs based on competing endogenous model, in-depth understanding of the molecular mechanisms of ovarian-related functions are important to study the reproductive traits of Small Tail Han sheep. It is valuable for us to investigate the critical role of ovarian circRNAs as a ceRNA in sheep fecundity. However, sheep-fecundity-associated circRNAs, miRNAs, mRNAs regulatory networks remain unknown.
Small Tail Han (STH) sheep, a well-known Chinese sheep breed, is endemic polytocous with year-round estrus and has an average litter size of 2.29, making it a suitable model to investigate the molecular genetic mechanisms linked with high prolificacy [36, 37]. In comparison, Dolang sheep, which are widely bred in the south Xinjiang region of China , have strong adaptability but low fecundity compared to Small Tail Han sheep [39,40,41,42], with an average litter size of 1.4 [43, 44]. Multiple studies have been conducted by researchers used RNA sequencing to analyze the gene expression patterns in the ovaries of Small Tail Han sheep during the follicular and luteal phases of the estrous cycle [45,46,47]. These studies have identified several differentially expressed genes that are involved in follicular development, oocyte maturation, and luteinization, all of which are important for fertility and concluded that Small Tail Han sheep have a unique molecular mechanism for regulating ovarian function, which may contribute to their high fecundity. Given the difference in prolificacy between these two sheep breeds, this study aimed to gain an understanding of the molecular mechanisms related to prolificacy and explore the expression of circRNAs as competitive endogenous RNA (ceRNA) in reproduction through a systematic transcriptome analysis, co-expression networking, functional enrichment, and genomic organization to obtain in-depth insights into the structure–function relationship of reproductive trait-associated circRNAs, miRNAs, and their target mRNAs. Therefore, the current study aims to identify the sheep-fecundity-associated circRNAs, miRNAs, and mRNAs regulatory networks in Small Tail Han and Dolang sheep to explore the regulatory mechanisms of the ovarian circRNAs in sheep fecundity.
Material and methods
The experiment was carried out at the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences, with approval from the Animal Care and Use Committee. The Ministry of Agriculture of the People’s Republic of China authorized the conduct of experiment of this study.
Ovarian tissue samples
Six healthy experimental animals, including three Small Tail Han sheep (high fecundity group) and three Dolang Sheep (as a control- low fecundity group), were used in this study. All the animals were approximately 2 years old, weighed around 50 kg, and raised under the same environmental conditions. Ovarian tissues were excided during dissection and the samples were collected for RNA sequencing and expression profiling of coding (mRNA) and non-coding RNA (circular RNA, miRNA) to investigate their roles in the ovary. All collected ovary specimen were immediately placed in liquid nitrogen and stored at − 70 °C for total RNA extraction.
Construction of mRNA, miRNA and circRNA libraries and sequencing
Total RNA was isolated from the excised ovaries using Trizol (Invitrogen), and RNA integrity was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, USA). For library construction, approximately 10 μg of total RNA per sample were used to deplete rRNA using the TruSeq Stranded Total RNA with Ribo-Zero Gold kit, followed by TRIzol extraction. By using the RNA Library Prep Kit, the rRNA-depleted RNAs were then fragmented and reverse-transcribed to produce cDNA libraries, and these libraries were further sequenced on the Illumina sequencing technology (HiSeqTM 2500), generating 150 bp/125 bp paired-end reads. For the miRNA-seq library, RNA fragments of 18–35 nucleotides were isolated and purified using a 15% polyacrylamide gel electrophoresis, and 5ʹ- and 3ʹ-termini were ligated with proprietary indexed adapters. Reverse transcription and low-cycle PCR were carried out to produce enough products for Illumina sequencing.
To assess the quality of raw sequences, we used Fast QC 0.11.5 followed by Trimmomatic-0.36 software, which trimmed out low quality reads/bases in different ways from 3' and 5' ends. For miRNA quality control, cutadapt version. 1.14 software used to remove connector sequence from miRNA fastq files, and Q20 quality control was performed via fastx_toolkit version. 0.0.13 software, retaining sequences with q20 of 80% and above only. Reads containing N bases were filtered out using NGSQCToolkit version. 2.3.3, resulting in high-quality Clean Reads for subsequent analysis.
In order to classified small RNA in the sequencing results, we processed clean reads through several databases and libraries. Firstly, we used Rfam (version. 10.1) to identify non-coding RNAs such as rRNAs, tRNAs, snRNAs, and snoRNAs. Next, we aligned the filtered reads from miRNA-seq libraries using bowtie version. 1.1.1 software and conducted a BLAST search against (http://www.sanger.ac.uk/software/Rfam), GenBank (http://www.ncbi.nlm.nih.gov/genbank/) and Repbase databases. To identify known miRNAs, we aligned the reads against miRBase version. Twenty-twp database (http://www.mirbase.org/). For unannotated small RNAs, we used miRDeep2 package version. 126.96.36.199 with the Perl script ‘quantifer.pl’  to predict novel miRNAs and analyzed the expression of known and novel miRNAs. We identified comparable miRNA star sequences using the miRBase database and the hairpin structure of a pre-miRNA. Finally, we identified DE miRNAs using DESeq2 in R , and predicted putative target mRNA using miranda to inquire about the functional role of miRNAs .
To analyze mRNA expression, we first aligned clean reads from RNA-seq libraries to the Ovis Aries reference genome using HiSAT2 version. 188.8.131.52. We then used StringTie version. 1.3.3b to assemble the reads and splice them into new transcripts. The sequencing reads were aligned to known mRNA transcripts sequences using an mRNA database available at (ftp.ncbi.nlm.nih.gov/genomes/all/GCF_000298735.2_O), and measured expression levels, which were normalized as FPKM by using express version. 1.5.1. We standardize the counts and identified differentially expressed mRNA using Leverage DESeq . Finally, mRNA with |log2fold changes|≥ 1 and p-value < 0.05 were identified as differentially expressed genes (DEGs/DEmRNAs).
To predict circRNAs, we first aligned the sequencing reads of each sample with the reference genome using BWA  software and generated a SAM file. We then used CIRI version. 2.0.3  and circ software version. 1.2  to predict novel circular RNAs from the unmapped reads. Only circRNA candidates that were identified were further analyzed for subsequent analysis. The expression levels of circRNA candidates were calculated with back-splice junction reads in RPM algorithm to quantify and normalized the number of junction reads counts and fold change by DESeq2 . Finally, filtered the differentially expressed circRNAs according to the difference multiple and the difference significance test with a foldchange > 2, and a p-value < 0.5. Also analyzed the enrichment of differentially expressed circRNAs through the annotation information of circRNA source transcripts. Furthermore, analyzed the functional enrichment analysis (GO and KEGG) of target genes involved in ceRNA network.
Bioinformatics analysis for ceRNA regulatory networking
To investigate the ceRNA relationships between circRNA, miRNA, and mRNA associated with sheep reproduction, Pearson’s correlation coefficient (r) were calculated between miRNA—mRNA and miRNA—circRNA pairs based on the expression levels of mRNAs, miRNAs, and circRNAs. We selected negatively correlated pairs with a p value < 0.05 and Pearson’s r > 0.8 for further analysis in X_LC-vs-D_LC. Further, used miranda (version.3.3a) software to predict binding sites for circRNA-miRNA and miRNA-mRNA interactions. Subsequently, we built a co-expression network among differentially expressed circRNA-miRNA and miRNA-mRNA pairs using shared pairs from predicted pairs from binding sites and predicted pairs from the expression of mRNA, circRNA, and miRNA through metscape plugin in cytoscape (version. 3.5.1). Similarly, calculated Pearson’s r between differentially expressed circRNA-mRNA pairs and screened out positively correlated pairs based on the role of mRNA-circRNA in the ceRNA relationship. Based on the principle for ceRNA prediction, only used shared pairs of miRNA-mRNA and miRNA-circRNA to predict ceRNA score in conjunction with hypergeometric distribution calculations using the MUTAME  method. The smaller the p-value corresponding to the ceRNA relationship, the more significant the miRNAs shared between two ceRNAs (mRNA and circRNA). We considered the shared pairs from predicted pairs of circRNA-mRNA based on the expression of circRNA and mRNA Pearson's correlation coefficient, and the predicted pairs of circRNA-mRNA based on ceRNA score principle as the true ceRNAs.
ceRNA network construction, phylogenetic tree & enrichment analysis of mRNA in ceRNA
CircRNAs and mRNAs communicate through shared miRNAs, which suggests that they may have similar functions. Using the co-expression data obtained in the previous step, we mapped ceRNA regulatory networks via Cytoscape software (version. 3.5.1) for the obtained ceRNA pairs that shared miRNA mutually. The full-length amino acid sequences of target mRNAs involved in ceRNETs were retrieved from NCBI database (https://www.ncbi.nlm.nih.gov/). To align protein sequences, the ClustalW tool (version. 2.0) was used with the default settings, and then manually modified in MEGA version. 7.0. Then, in MEGA 7.0 with default parameter, the neighbor joining (NJ) tree was constructed with bootstrap replicates by using the neighbor joining method . We performed functional enrichment analysis (GO and KEGG) of mRNAs contributing to the ceRNA network using the R package cluster profiler , with a p-value < 0.05 used as the cut-off value for KEGG pathway enrichment analysis.
Genomic structure, expression profiling and subcellular localization of target genes
The CDS and genomic sequences of target genes were retrieved from NCBI database (https://www.ncbi.nlm.nih.gov/). The CDS and genomic sequences were used to draw picture of exon/intron organization at Gene Structure Display Server (GSDS) Program (http://gsds.cbi.pku.edu.cn/). To analyze the tissue specific expression profiles of target genes acquired from transcriptome analysis of ovarian tissues of two sheep breeds (X_LC and D_LC), we used FPKM values from RNA-seq data. The heat map along with Phylogenetic tree was generated using TBtools version. 1.108  for differential expression analysis. Subcellular localization prediction of ovarian tissues associated target genes of ceRNA, binding cites was carried out by using several websites, online servers and tools such as TargetP online web server (http://www.cbs.dtu.dk/services/TargetP/), WOLF-PSORTtool (https://wolfpsort.hgc.jp/), ProtComp (http://linux1.softberry.com/berry.phtml) and CELLO v.2.5 (http://cello.life.nctu.edu.tw/).
Quantitative (qRT) PCR verification
To validate the expression of differentially expressed circRNAs, miRNAs, and mRNAs involved in ceRNA networks of X_LC-vs-D_LC, we used quantitative real-time PCR (qRT-PCR). We randomly selected five circRNAs, miRNAs, and mRNAs for validation. The expression levels of the selected circRNAs, mRNAs, were normalized against the expression of a housekeeping gene, GAPDH and U6 was used as reference gene of miRNA respectively. Primers were designed and synthesized by Sangon Biotech Co., Ltd. (Shanghai, China). Total RNA was extracted using RNA extraction reagent RNAisoPlus (TakaRa, AA6620-1). For circRNAs, Lan Yi bio reverse transcription kit (Lan Yi Biologics, LY-0160) and miRNA first strand cDNA synthesis Kit (tail addition method) (Sheng Gong, B532451) for miRNA were used for reverse transcription. Quantitative real-time polymerase chain reaction (qRT-PCR) was performed in the lightcycler 96 system (Roche Applied Science, Mannheim, America) using SYBR green real-time PCR master mix (Chun Lei Jie Chong, MR2001) for circRNAs and miRNA fluorescent quantitative PCR Kit (Dye Method) (Sango, B532461) for miRNA following the manufacturer’s instructions. The qRT-PCR analysis was carried out in triplicate. The total 20 μL reaction mixture contained 5 μL 2 × iTaq™Universal SYBR@ Green Supermix, 1 μL cDNA, 10μL H2O, and 0.2 μL each of forward and reverse primers. The following programs were used; pre-incubation 95 °C for 3 min, amplification (40 cycles) of 95 °C for 3 s, 60 °C and for 20 s, then melting curves 95 °C for 15 s, 60 °C for 15 s, and 95 °C for 15 s. The relative change of gene expression of circRNAs, miRNAs, and mRNAs among the control and experimental groups was calculated by the 2 − ∆∆Ct method.
All the data presented here mean values ± standard deviations. Student’s T-test was performed for the comparison, and P-value < 0.05 was considered as statistically significant.
Identification of circRNAs, miRNAs, mRNAs in Ovary Tissue of X_LC-Vs-D_LC
We collected ovarian tissue from Small Tail Han sheep (X_LC) and Dolang sheep (D_LC) (n = 3each group) for RNA-seq analysis. We constructed and sequenced 6 ribosomal RNA (rRNA)-depleted libraries, named X_LC (X_LC_1, X_LC_2, X_LC_3) and D_LC (D_LC_1, D_LC_2, D_LC_3). We obtained a total of 91.04 G of clean data, with the effective data volume of each sample distributed between 14.79 and 15.82G. The Q30 base distribution was 91.51 to 92.0%, and the average GC content was 49.23%. We compared the reads of each sample to the reference genome and obtained a comparison rate of 92.44 to 93.29%. Number of circular RNA predicted in each sample. Specifically, in D_LC, we predicted 3,282, 2,232, and 3,416 circRNAs with non-unique numbers and 1,053, 492, and 1,198 with unique numbers in samples D_LC_1, D_LC_2, and D_LC_3, respectively. D_LC_1, D_LC_2, D_LC_3 with non-unique number;3282, 2232, 3416 and unique number; 1053, 492, 1198. In X_LC, we predicted 2,980, 3,211, and 2,351 circRNAs with non-unique numbers and 804, 923, and 575 with unique numbers in samples X_LC_1, X_LC_2, and X_LC_3, respectively (Fig. 1A).
Based on the annotation information of protein-coding genes and transcripts released in the database, the transcript with the largest overlap in position with circRNA obtained through position information comparison, and the sequence of circRNA is predicted based on the transcript information (Table 1). We determined 9878 new circRNAs with a total length of 23,522,667 nt and an average length of 2381.32 nt.
Among them, majority of circRNAs were characterized with a length of 200–600 nucleotides (nt) and more than 2,000nt (Fig. 1B). The circRNA_4335 predicted with highest circRNA length, whereas circRNA_2452 with the minimum circRNA length. The circRNAs distribution on chromosome identified and obtained top ten chromosomes with maximum number of circular RNAs including; NC_019458.2 (1182), NC_019459.2 (987), NC_019460.2 (946), NC_019461.2 (436), NC_019462.2 (390), NC_019463.2 (409), NC_019464.2 (547), NC_019465.2 (353), NC_019484.2 (338) and NC_019466.2 (301) (Fig. 1D). A total of 3,781 parental genes generated these circRNAs, with some parental genes producing single circRNAs and others generating multiple circRNAs. The identified circRNAs were categorized into five types based on their genomic location; including sense-overlapping-8970 (90.81%), intergenic-356 (3.60%), exonic-316 (3.20%), antisense-151 (1.53%), and intronic-85 (0.86%) circRNAs (Fig. 1E). Notably, only a limited number of circRNAs consisted of protein-coding exons. CircRNAs that consisted of one exon were significantly longer than those that comprised multiple exons (Fig. 1B). Identified circRNAs were equipped with various number of exons such as 1854 circRNAs had three exons, 1628 circRNAs had two exons, 1534 circRNAs had four exons, and 1134 circRNAs had one exon (Fig. 1C).
Differentially expressed circRNAs, miRNAs and mRNAs in X_LC-vs-D_LC ovarian tissue differentiation
By comparison of X_LC-vs-D_LC, a total of 44 DE-circRNAs were identified. The number of significantly upregulated circRNAs was higher than that of significantly downregulated circRNAs in comparison groups (Fig. 2A). For mRNAs, a total of 20,533 mRNAs, including 397 DEGs were screened out with the criteria of differantial-pvauel-0.05 and Foldchange-2. Among them160 DEGs were upregulated and 237 were downregulated respectively (Fig. 2D). In case of miRNAs, 1,186 miRNAs comprising of 146 known miRNAs and 1040 novel miRNAs were identified in X_LC-vs-D_LC. Of these, 35 DEmiRNAs were tested by diff-pval-0.05 and FC-2 in X_LC-vs-D_LC (Fig. 2G). The differences resulting from the comparison were reflected in the volcano map and heat map of differentially expressed circRNA (Fig. 2B-C), mRNA (Fig. 2E-F) and miRNA (Fig. 2H-I) in X_LC-vs-D_LC.
Co-expression and prediction of target genes of miRNA-mRNA and miRNA-circRNA
In X_LC-vs-D_LC, correlation (r) was calculated between 35 DEmiRNAs and 397 DE mRNAs (DEGs), resulting in the identification of 2510 miRNA-mRNA correlated pairs with a threshold of an absolute value of r ≥ 0.80 and a p-value ≤ 0.05. According to the principle of miRNA and mRNA action, 903 negatively regulated correlated pairs were screened out, and their target binding sites (MREs) on mRNA were further predicted using miranda software version. 3.3a. Thus, we obtained 337 miRNA-mRNA interaction pairs that could be involved in transcriptional degradation/translational inhibition, contributed by 166 DEGs and 22 DEmiRNAs (Supplementary Table 1 and Fig. 3A). Similarly, a correlation coefficient was measured between 35 DEmiRNAs and 44 DECs, resulting in 205 miRNA-circRNA correlated pairs. Following the principle of miRNA and circRNA action, 45 negatively regulated pairs were screened out, and miRNA binding sites (MREs) on circRNAs were predicted using miranda v3.3a (Supplementary Fig. 1). This yielded 12 miRNA-circRNA pairs contributed by 9 DECs and 8 DEmiRNAs (Fig. 3B and Table 2).
ceRNA data analysis & co-expression analysis and filtering of mRNA and circRNA
Furthermore, the ceRNA score was calculated, resulting in to 88 regulatory relationships between circRNAs and their target genes using the MuTaME method. The smaller the p-value corresponding to the ceRNA relationship, the more significant the shared target miRNAs between the two ceRNAs (circRNAs and mRNA). Meanwhile, Pearson correlation was determined among 397 DEmRNAs and 44 DECs in X_LC-vs-D_LC, which led to 4274 pairs of mRNA-circRNA. The co-expression network of the top 200 circRNA-mRNA interactions was constructed (Fig. 4A). Based on the role of mRNA-circRNA in the ceRNA relationship, filtered 3005 positively correlated mRNA-circRNA pairs respectively (Supplementary Table 2). These pairs were combined with the result of ceRNA score, in turn revealing 50 common intersections (Fig. 4B and Supplementary Table 3).
ceRNA network & phylogenetic tree of ceRNA target genes
In accordance with the ceRNA theory, ceRNAs share miRNAs as a junction where upregulated miRNAs are associated with downregulated circRNAs and mRNAs, and downregulated miRNAs are associated with upregulated circRNAs and mRNAs. The ceRNETs were constructed using the obtained top 50 interactions based on ceRNA scores, comprising six miRNAs: novel339_mature, novel579_mature, novel62_star, novel14_mature, novel130_mature, and novel34_mature > novel661_mature; eight circRNAs: circRNA_9564, circRNA_9690, circRNA_3497, circRNA_8312, circRNA_3257, circRNA_8401, circRNA_8396, and circRNA_4140; and 40 mRNAs (Fig. 5A). To understand the evolutionary, structural, and functional associations of ceRNA target mRNAs, a phylogenetic tree was constructed. Clustering into similar clades was observed for genes with similar structures. The best representation was exhibited by novel579_mature, which was common in circRNA_8401, circRNAs_9564, and circRNAs_3257, and novel339_mature, which was common in circRNA_8312 and circRNA_9564 (Fig. 5A and B).
Functional enrichment analysis of ceRNA
Gene ontology (GO) analysis of mRNA in ceRNETs revealed that out of 40 mRNAs, 28 were successfully annotated, providing insight into biological processes (BP), cellular components (CC), and molecular functions (MF). Fisher’s exact test was used to calculate the enrichment significance of each term in BP, CC, and MF, with detailed results provided in Supplementary Table 4. Each term is arranged in ascending order according to p-value, indicating that a smaller p-value corresponds to more significant enrichment. A total of 65 terms had a p-value ≤ 0.05, with the first three terms having the smallest p-values being "midbody" (TermID: GO:0,030,496, p-value: 0.0054) and "transcription factor activity, protein binding" (TermID: GO:0,000,988, p-value: 0.0081). Among these, 43 terms had p-values ≤ 0.05 in BP, 10 terms in CC, and 12 terms in MF, with none having FDR ≤ 0.05 (Supplementary Table 4 and Fig. 6A-B).
In the KEGG enrichment analysis, only 12 out of 40 mRNAs in ceRNETs were annotated, participating in 40 different pathways. We selected the top 20 pathways with p-values ≤ 0.05. Among them, the pathways with the smallest p-values were "Protein processing in the endoplasmic reticulum" (TermID: path: oas04141, p-value: 0.023) followed by "Alpha-linolenic acid metabolism" (TermID: path: oas00592, p-value: 0.030), as shown in Fig. 6C. The majority of pathways were associated with genetic and metabolic pathways, such as galactose metabolism, metabolism of xenobiotics by cytochrome P450, drug metabolism-cytochrome P450, type I diabetes mellitus, antigen processing and presentation, progesterone-mediated oocyte maturation, cell cycle, cell adhesion molecules (CAMs), axon guidance, phagosome, Foxo signaling pathway, p53 signaling pathway, and Ras signaling pathway. Furthermore, we observed similarities between the functional enrichment analysis of DECs host genes and the ceRNA target genes (Fig. 6D-E).
Gene structure, expression profiling and subcellular localization of target DEGs
Most animal genes contain exons and introns, and their arrangement can reveal their evolutionary relationships and suggest their involvement in specific functions . To further elucidate the evolutionary, structural, and functional relationships of ceRNA target genes, we conducted a comprehensive comparative study by constructing a phylogenetic tree and analyzing gene structures, tissue-specific expression patterns, and predicted subcellular localizations. The compelling illustration in Fig. 7A, B, C, D demonstrates the interconnection between gene structures and subcellular localization. Previous studies have reported a relationship between exon and intron distribution patterns and their biological roles [60, 61]. In this study, we identified SCML4, IGDC33, MEGF10, KCNK9, and EPHA3 as the longest genes compared to LOC101103023, which has the shortest genomic length of 3 Kb. To gain further insight into the structural evolution of the aforementioned genes, we developed a new phylogenetic tree and examined their structural features using the online web portal GSDS (http://gsds.cbi.pku.edu.cn/). Comprehensive analysis of exon–intron structures aligned with the phylogenetic tree sequence provides a holistic illustration of their relative lengths (Fig. 7B) and demonstrates their random distribution throughout the tree. However, clustering into similar clades was observed for genes with similar structures. The best representation was exhibited by novel579_mature, which was common in cirRNA-8401, cirRNA_9564, and cirRNA_3257 (Figs. 5B and 7A). The number of exon regions varies among all the genes (from 1 to 18); however, surprisingly, all closely associated genes exhibited similar structural organization and divergent exon–intron lengths.
To elucidate the potential difference in fecundity between X_LC and D_LC, we compared the mRNA expression profiles of target genes involved in the ceRNA regulatory network (Fig. 7C). A prominent difference was found in gene expression profiles between Dolang sheep (D_LC) and Small Tail Han Sheep (X_LC). EPHA3, KLF9, SLA, SLC45A1, and SCML4, the target genes of novel579_mature, are highly expressed and upregulated in the ovary tissue of Small Tail Han sheep compared to Dolang sheep. Similarly, differential expressions of LOC101102473, PLA2G2D, and LOC101110545 were noted during tissue-specific expression analysis of ovary tissues from both groups. It is evident that our results from tissue-specific expression profiling are in agreement with the ceRNA network, GO analysis, and KEGG enrichment analysis, which confirms the reliability and accuracy of our transcriptome data analysis. These findings are further corroborated by qRT-PCR analysis in the next step.
In order to predict the subcellular localization of ceRNA target genes, we used several websites and online tools to ensure the accuracy of results. Collectively, data suggested that the majority of the aforementioned genes were localized to the cytoplasm, the nuclear envelope, and the plasma membrane. Although few genes were also predicted in mitochondria, cytoskeletal and extracellular membranes. This is consistent with their associated functions (Fig. 7D). It has been observed that most of the circRNA_3257, circRNA_8312, circRNA_8401, and circRNA_9564 based target genes were grouped in same clades and represented similar subcellular localization along with analogues expression profiles. However, novel579_mature was identified as a common miRNA binding site of circRNAs_3257, cirRNA_8401 and circRNAs_9564 and contains the similar target genes. Novel339_mature have a target binding site for circRNAs (circRNAs_8312, circRNAs_9564) circRNAs_8396, circRNAs_9690 and circRNA_4140 have a binding site for novel130_mature, novel62_star and novel34_mature > novel661_mature, indicating the precision and accuracy of our transcriptomic study, revealing strong collaboration of gene structures and their subcellular localization with tissue specific expression patterns along the phylogenetic tree (Figs. 5B and 7A-D).
Validation of candidate circRNAs
Based on our predicted ceRNA circRNA-miRNA-mRNA interaction network and ceRNA target mRNA KEGG enrichment analysis, we selected key circRNAs (circRNA_8396|NC_019477.2:26838353_26929260_ + , circRNA_8312|NC_019477.2:3206957_3240924_-, circRNA_9690|NC_019484.2:62838118_62839910_-, circRNA_4140|NC_019463.2:66446753_66449897_-), miRNAs (novel130_mature, novel339_mature, novel62_star, novel34_mature > novel661_mature, novel579_mature), and mRNAs (LOC101102473, LOC101110545, KCNK9, EPHA3, PLA2G2D) for expression validation in the ceRNA regulatory network. We performed qRT-PCR and observed dynamic expression profiles of the aforementioned circRNAs, miRNAs, and mRNAs. The results showed downregulation of circRNA_8396 and circRNA_8312 in X_LC-vs-D_LC, whereas circRNA_3257 and circRNA_9690 exhibited significant downregulation, and circRNA_4140 demonstrated remarkable upregulation in ovary tissue of X_LC-vs-D_LC (Fig. 8). In the case of miRNAs, novel130_mature, novel339_mature, and novel579_mature were significantly upregulated, but novel62_star was not significantly upregulated in X_LC ovary tissue. Conversely, novel34_mature > novel661_mature was downregulated in X_LC (Fig. 8). The qRT-PCR results for mRNA revealed that the expression of LOC101102473, LOC101110545, EPHA3, and PLA2G2D were remarkably downregulated in X_LC. However, KCNK9 was significantly upregulated in X_LC ovarian tissue (Fig. 8). These findings were consistent with our predicted results from the transcriptome data.
Differentially expressed circular RNAs and their functional enrichment analysis
A large body of knowledge has emphasized the significant role of noncoding RNAs in sheep fertility [32, 62,63,64]. The study of ceRNA consider landmark in understanding the mutual regulatory relationship and interactions of RNA-RNA [64, 65]. In this study, high-throughput sequencing was used to identify and characterize ovarian tissue circRNAs from Small Tail Han Sheep (X_LC) and Dolang sheep (D_LC) and determined ceRNA regulatory network. In this study, majority of circular RNAs were predicted with exon 1–5 which are consistent with the sheep uterus circRNAs  but inconsistent with sheep pituitary gland circRNAs  thereby exhibit complexity and functional diversity. Besides, studies have been reported about the expression and potential biological functions of circRNAs in reproductive organs in goats , mice  and humans . Several findings have been reported about the circRNAs expression and their biological functions as miRNA sponge [25, 27, 28]. In sheep uterus study, 147 and 364 circRNAs out of 32,687 were differentially expressed in polytocous and monotocous groups in the follicular phase and luteal phases and DECs host genes significantly enriched with estrogen signaling pathway, oxytocin signaling pathway . While, in the comparative study between the follicular phase/luteal phase of sheep, 15 DECs out of 3223 were predicted and DECs host genes were enriched in the Rap1 signaling pathway, PI3K–Akt signaling pathway and neuroactive ligand–receptor interactions . In this study, we found 9,878 new circRNAs with a total length of 23522667nt and an average length of 2381.32nt respectively. 44 circRNAs were differentially expressed and majority of enriched GO terms of DECs host genes were related to cell proliferation and reproductive process and KEGG enrichment analysis enriched in cell adhesion molecules, phagosome, PI3K-Akt signaling pathway, neuroactive ligand-receptor interaction, glutathione metabolism, metabolism of xenobiotics by cytochrome P450, MAPK, axon guidance, valine, leucine and isoleucine degradation, and endocytosis pathways. A similar phenomenon was observed in other studies [27, 28, 32]. Metabolic changes reported during transition period from ovulation to the estrous and biomolecules including vitamins, amino acids, lipid, benzoic acid, carbohydrates and other intermediate and secondary metabolites are at their highest levels at the time of ovulation [27, 69, 70]. Prior stud have been reported that MAPK pathway have a great effect on granulosa and cumulus cells which plays essential role in oocyte maturation . Therefore, based on above evidences we suggested that DECs might contribute in sheep prolificacy.
CeRNA analysis and networking/ functional enrichment analysis of ceRNA regulatory genes
Circular RNAs (circRNAs) as ceRNA contribute to a various signaling pathways that are crucial for the development process . The fact that mRNA expression negatively regulated by miRNA activity, whereas, circRNAs inhibit or relieve repression of miRNA for translation . Study determined the role of circRNA as ceRNA in follicular development in Hanper sheep . Here, circRNA and mRNA target binding sites were predicted for miRNA (Fig. 3A-B and Table 2) which benefited to determined regulatory relationship among differentially expressed circRNAs, miRNA, mRNA (ceRNETs) (Fig. 5A). In this study, numerous circRNAs holding binding site for common miRNA such as; circRNA_9564, circRNA_8312 shared novel339_mature miRNA, circRNA_9564, circRNA_8401, circRNA_3257 shared novel579_mature miRNA respectively. Whereas, circRNA_8396, circRNA_9690, circRNA_4140 have potential binding sites for only single miRNA including novel339_mature, novel62_star, and novel34_mature > novel661_mature respectively. novel579_mature dominantly target to the number of mRNAs including; CHMP4C, EPHA3, SCML4, MSH4, LOC106991756, LOC101120875, HOXA11, RGS5, SLC45A1, C6H4orf32, LOC101104325, and PDK4 which significantly involved in pathways including endocytosis and axon guidance pathway. Prior study determined the crucial role of endocytic pathways in the development of the reproductive organs . circRNA_8396 have potential binding sites for novel339_mature which target to the following mRNAs; CENPF, CLEC19A, AURKB, and LOC101102473.circRNAs_4140 have potential binding sites for novel34_mature > novel661_mature which target to KCNK9. To further investigate potential circRNAs mechanisms of action in ceRNA regulatory network, we applied KEGG analyses. It was reported that foxO signaling pathway, cell cycle, p53 signaling pathway, endocytosis, progesterone-mediated oocyte maturation participate in reproduction . The results we obtained suggest that multiple signaling pathways form a complex regulatory network involved in prolificacy.
Coherence of ceRNA regulatory network with differential expression datasets to regulate reproduction
CircRNAs biological function as a microRNA (miRNA) sponge and regulating the target mRNA expression by forming the circRNA-miRNA-mRNA regulatory axis . Pearson correlation based networking and expression profiling provided ceRNA regulatory pairs as shown in figure (Fig. 5A). KEGG enrichment analysis of mRNA involved in ceRNETs revealed significant pathways such as; Galactose metabolism, Glutathione metabolism, Metabolism of xenobiotics by cytochrome P450, Drug metabolism—cytochrome P450, Cell adhesion molecules (CAMs), Antigen processing and presentation, p53 signaling pathway, and Progesterone-mediated oocyte maturation pathways respectively. Of note, studies have demonstrated the involvement of above-mentioned pathways in ovarian physiology and play essential role in follicular development, oocyte maturation, development of reproductive organ, proliferation, immunity, antioxidant, and metabolic process [74,75,76,77]. Previous studies demonstrated that the exposure of xenobiotics may destroy the primordial follicles which is responsible for premature ovarian failure and reduce fertility [78, 79]. Glutathione metabolism, well known for their spectacular role as a free-radical scavenger, intervenient in xenobiotics metabolism, cell-cycle regulation, and a reservoir of cysteine , as well as play key roles in cellular redox homeostasis . Study reported the involvement of Progesterone-mediated oocyte maturation pathways in oocyte development in sheep breeds ,and regulating the uterine receptivity and maintenance of pregnancy .
Researchers established the role of galactose metabolism in energy delivery, and galactosylation of complex molecules . The metabolism of galactose to UDP-glucose involves three major enzymes, galactokinase, galactose-1-phosphate uridyltransferase (GALT), and UDP-galactose-4-epimerase. In case of any disturbance or deficiency of UDP-galactose are considered to be important in the pathogenesis of ovarian dysfunction and in GALT deficiency  due to interference with ovarian apoptosis and gonadotrophin signaling, thus effecting fertility. Wu et al. reported that deregulation of glucose metabolism in diabetic individuals , which might induce ovarian anomalies . It is believed that these pathways may play a critical role in regulating ovarian physiology. Previously, study have also examined miRNAs as important components of the p53 transcriptional network . MiR-25 and miR-30d, have been shown to negatively regulate the transcription of P53 gene. Several other miRNAs, including miR-16–1, miR-143, miR-145, miR-34, miR-194, miR-192, miR-215, and miR-29, have been identified as miRNAs that are involved in the transcription P53 network, either by being directly altered by p53 or through their associations with downstream genes targeted by p53.
Significance of ceRNA target gene in regulating reproduction and associated metabolic syndrome
Through ceRNA circRNA-miRNA-mRNA interaction analysis and based on functional enrichment analysis established the following key ceRNA regulatory axis; circRNA_3257-novel579_mature-EPHA3, circRNA_8396-novel130_mature-LOC101102473, and circRNA_4140- novel34_mature > novel661_mature-KCNK9. These shortlisted circRNAs were the result of sense-overlapping and the length of circRNA_3257, circRNA_4140, circRNA_8312, circRNA_8396, circRNA_9690, circRNA_9564 were 758, 3145, 33,662, 44,924, 1793, and 76,127 nt respectively. We found that downregulated circRNA_3257 had a potential binding site for upregulated novel579_mature miRNA and influences EPHA3 gene expression, which is regulated by the protein kinase A (PKA) pathway . The ephrin-Eph gene family has a known physiological role in regulating mammalian reproductive function, such as in granulosa cells of bovine ovarian follicles , mouse , and human luteinizing granulosa cells . Prior studies have suggested that EPHA3 plays a role in treating ovarian endometriosis, potentially promoting apoptosis and autophagy of macrophages via the inhibition of the mTOR signaling pathway and reducing oxidative stress [93, 94]. In this study, we found that downregulated EPHA3, involved in the axon guidance pathway, functions through cell-to-cell contact , and regulates the expression of guidance proteins Ras and Rho GTPases during embryonic development . Furthermore, EPHA3 has been shown to interact with presenilin genes sel-12 (PS1), regulating axon guidance and kinesin-mediated axonal transport of motor neurons. Loss of function mutations in the presenilin genes sel-12 results in abnormal axonal projections, an effect attributed to altered Notch signaling pathway , which is integral to maintaining fertility in the ovaries through developmental regulation and granulosa cell function . Based on the ceRNA hypothesis, we suggest that circRNA_3257 may act as a sponge for novel579_mature miRNA, thus favoring the expression of repressed EPHA3, and potentially play a role in the fecundity of Small Tail Han Sheep.
LOC101102473 is a target gene of novel130_mature that is downregulated in ovary tissue of Small Tail Han Sheep , and is involved in pathways related to reproduction, such as FoxO signaling pathway, Cell cycle, p53 signaling pathway, Progesterone-mediated oocyte maturation, consistent with the previous research . FOXO signaling is the central pathway controlling growth and metabolism in all cells , and closely related with ovarian function . Low levels of p53 expression maintain cell cycle homeostasis and cell death . Dysregulation of the cell cycle, particularly the G1-S-phase transition, is implicated in epithelial ovarian cancer . If DNA damage occurs, p53 accumulates in the cells and induces p21-mediated inhibition of cyclinD/CDK. The transition to S phase is triggered by the activation of the cyclinD/CDK complex, which phosphorylates the retinoblastoma protein pRb, a known cell proliferation regulator . The transcription factor p53 functions as a suppressor of cell growth, and alterations in p53 lead to loss of this negative growth regulation and more rapid cell proliferation. Previous studies have demonstrated that dysfunction of the p53 signaling pathway contributes to the development of ovarian cancers , hence affecting fertility. These outcomes suggest that circRNA_8396 might affect sheep fecundity by regulating the expression of LOC101102473 genes linked with the above signaling pathways via sponging novel130_mature miRNA.
The gene KCNK9, which is linked to potassium channels in X_LC-vs-D_LC and has been previously identified in the membrane of cow oocytes , is upregulated as a target gene of novel34_mature > novel661_mature. While previous studies have shown its potential role as a therapeutic target for adenomyosis which induce infertility . KCNK9 also helps to control progesterone production in the ovary, and intracellular potassium and calcium concentrations have a significant impact on fertility . Our functional enrichment analysis has revealed that KCNK9 is involved in the aldosterone synthesis and secretion pathway, which is critical for maintaining blood pressure, circulating blood volume, and uteroplacental perfusion during pregnancy . Several other studies have shown that aldosterone synthesis and secretion levels are high throughout pregnancy, indicating potential involvement in the regulation of placental and fetal development [108, 109]. However, studies have also reported that aldosterone is involved in gynecological diseases due to metabolic alterations induced by the usage of hormonal contraceptives, PCOS, uterine fibroids and endometriosis, inflammation, and hypertension, which can result in an increase in the synthesis of angiotensinogen, activating all the RAAS and inducing the onset of sodium and water retention . Therefore, we hypothesize that circRNA_4140 may act as a sponge for novel34_mature > novel661_mature and regulate the expression of KCNK9, which plays a crucial role in fertility.
In addition to the upregulated novel339_mature, the downregulated LOC101110545 (HLA class II histocompatibility antigen, DO alpha chain) gene is also involved in reproduction, as it is a target gene involved in the Cell adhesion molecules (CAMs)  and Phagosome [112, 113] pathways crucial for reproduction. CAMs are a family of glycoproteins that play a crucial role in various physiological processes, including cell migration, tissue development, and immune response. They also play key roles in inducing leukocyte infiltration in an inflammatory site during ovulation . Studies have shown that CAMs are implicated in the regulation of ovarian follicular development, oocyte maturation, and fertilization in sheep [115, 116], all of which are critical for successful reproduction and fecundity. One of the key CAMs involved in sheep reproduction is integrin, which facilitates cell–cell and cell–matrix interactions, and is regulated during folliculogenesis and ovulation. Integrin-mediated signaling is essential for the survival and growth of ovarian follicles, as well as for the adhesion and migration of oocytes and granulosa cells during follicle development . Another important CAM in sheep fecundity is selectin. Selectin is another important CAM in sheep fecundity. It is a family of carbohydrate-binding proteins that mediate leukocyte-endothelial cell adhesion during inflammation and immune responses. Studies have reported that selectins play a role in follicular development and ovulation by facilitating the adhesion and migration of leukocytes to the ovary and promoting the release of inflammatory cytokines, chemokines and growth factors [117, 118]. Changes in these molecules are associated with the cyclic changes in the estrous cycle, compounding their role in the ovulatory process . Further research is needed to elucidate the precise mechanisms underlying CAM-mediated signaling in sheep reproduction and to develop novel interventions to improve sheep breeding and production.
In livestock animals, including sheep, phagosomes play a crucial role in regulating reproductive processes such as ovulation, fertilization, and embryo development . As ovarian follicles develop, some follicles undergo atresia, a process by which they degenerate and are eliminated from the ovary [121, 122]. This process involves the death of granulosa cells, which are engulfed and cleared by phagocytes via phagocytosis. Moreover, in mammals, the engulfment of apoptotic cells by macrophages induces the production of anti-inflammatory cytokines to suppress inflammatory responses . Failure to properly clear apoptotic granulosa cells can lead to inflammation and oxidative stress, which can negatively impact follicular development and ovulation [122, 124]. Based on the evidence presented, we can infer that circRNA_8312 may influence LOC101110545 gene expression by suppressing novel339_mature activity. Although the role of the LOC101110545 gene in reproduction has yet to be investigated, our data suggest its possible involvement in reproductive performance through ceRNA networks that regulate ovarian physiology. Comparison between our results and those of previous studies suggests that these selected circRNAs may act as ceRNAs and contribute to signaling pathways regulating reproductive traits. However, further research is required to identify the exact mechanisms involved. This study provides an important source of information that various circular RNAs can act as competitive endogenous RNA, containing common miRNA binding sites, and may act as miRNA sponges to regulate the expression of valuable prolificacy genes, thereby improving sheep fecundity.
In this study, analyzed the expression of circRNAs, miRNAs, and mRNAs in the ovine tissue of Small Tail Han sheep and Dolang sheep via using RNA sequencing. A total of 44 differentially expressed circRNAs, 397 differentially expressed mRNAs, and 35 differentially expressed miRNAs were identified. We went through a comprehensive study on ceRNA regulatory network analysis of ovine tissues from two different sheep breeds with contrasting fecundity traits. All of the DE circRNAs and mRNAs (fold change > 1, p < 0.05), combined with all of the miRNAs (differently or non-DE). By calculating the Pearson correlation between circRNA-miRNA and miRNA-mRNA, we were able to retain the regulatory relationship of circRNA (upregulation)-miRNAs (downregulation or invariably)-mRNA (upregulation) or circRNA (downregulation)-miRNAs (upregulation or invariably)-mRNA (downregulation), and eventually calculated the ceRNA score to identify circRNA–miRNA–mRNA pairs. Based on the role of mRNA-circRNA in the ceRNA relationship, we found a co-expression positive relationship between the differentially expressed circRNA and mRNA, with a correlation coefficient greater than 0.8. A total of 50 significant ternary relationship (circRNAs-miRNA-mRNA) were identified involved in fertility-related pathways such as FoxO signaling pathway, Cell cycle, p53 signaling pathway, Progesterone-mediated oocyte maturation, Cell adhesion molecules, Phagosome, Antigen processing and presentation, and Axon guidance. We identified several key circRNAs, including circRNA_3257, circRNA_8396, circRNA_4140, and circRNA_8312, that act as molecular sponges for miRNAs, and target genes including EPHA3, LOC101102473, KCNK9, and LOC101110545. Our study provides new insights into the molecular mechanisms underlying fertility in sheep and highlights the importance of ceRNA network analysis for understanding complex regulatory networks.
Availability of data and materials
The materials and datasets used and analyzed during the present study are available from the corresponding author upon reasonable request.
Hernández-Montiel W, et al. RNA-seq transcriptome analysis in ovarian tissue of pelibuey breed to explore the regulation of prolificacy. Genes. 2019;10(5):358.
El-Halawany N, et al. Investigating the effect of GDF9, BMP15, BMP6 and BMPR1B polymorphisms on Egyptian sheep fecundity and their transcripts expression in ovarian cells. Small Rumin Res. 2018;165:34–40.
Ma H, et al. Identification of novel genes associated with litter size of indigenous sheep population in Xinjiang, China using specific-locus amplified fragment sequencing technology. PeerJ. 2019;7:e8079.
Miao X, et al. An integrated analysis of miRNAs and methylated genes encoding mRNAs and lncRNAs in sheep breeds with different fecundity. Front Physiol. 2017;8:1049.
Miao X, Luo Q, Qin X. Genome-wide transcriptome analysis in the ovaries of two goats identifies differentially expressed genes related to fecundity. Gene. 2016;582(1):69–76.
Miao X, et al. Ovarian transcriptomic study reveals the differential regulation of miRNAs and lncRNAs related to fecundity in different sheep. Sci Rep. 2016;6(1):1–11.
Miao X, et al. Comparative DNA methylome analysis of estrus ewes reveals the complex regulatory pathways of sheep fecundity. Reprod Biol Endocrinol. 2020;18(1):77.
McNatty KP, et al. Ovarian characteristics in sheep with multiple fecundity genes. Reproduction. 2017;153(233):e40.
Tian D, et al. Comparative transcriptome of reproductive axis in Chinese indigenous sheep with different FecB genotypes and prolificacies. Anim Reprod Sci. 2020;223:106624.
Hernández-Montiel W, et al. RNA-seq transcriptome analysis in ovarian tissue of Pelibuey breed to explore the regulation of prolificacy. Genes (Basel). 2019;10(5):358.
Chen HY, et al. Differential gene expression in ovaries of Qira black sheep and Hetian sheep using RNA-Seq technique. PLoS One. 2015;10(3):e0120170.
Skaftnesmo KO, et al. Integrative testis transcriptome analysis reveals differentially expressed miRNAs and their mRNA targets during early puberty in Atlantic salmon. BMC Genomics. 2017;18(1):1–12.
Yang H, et al. Comprehensive analysis of long noncoding RNA and mRNA expression patterns in sheep testicular maturation. Biol Reprod. 2018;99(3):650–61.
Santosh B, et al. Non-coding RNAs: biological functions and applications. Cell Biochem Funct. 2015;33(1):14–22.
Zhang Z-C, Guo X-L, Li X. The novel roles of circular RNAs in metabolic organs. Genes Dis. 2018;5(1):16–23.
Li HM, Ma XL, Li HG. Intriguing circles: conflicts and controversies in circular RNA research. Wiley Interdiscip Rev RNA. 2019;10(5):e1538.
Kulcheski FR, Christoff AP, Margis R. Circular RNAs are miRNA sponges and can be used as a new class of biomarker. J Biotechnol. 2016;238:42–51.
Huang Y, et al. Biological functions of circRNAs and their progress in livestock and poultry. Reprod Domest Anim. 2020;55(12):1667–77.
Yu C-Y, Kuo H-C. The emerging roles and functions of circular RNAs and their generation. J Biomed Sci. 2019;26(1):1–12.
Enuka Y, et al. Circular RNAs are long-lived and display only minimal early alterations in response to a growth factor. Nucleic Acids Res. 2016;44(3):1370–83.
Li Z, et al. Exon-intron circular RNAs regulate transcription in the nucleus. Nat Struct Mol Biol. 2015;22(3):256–64.
Salmena L, et al. A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language? Cell. 2011;146(3):353–8.
Fiannaca A, et al. miRTissue ce: extending miRTissue web service with the analysis of ceRNA-ceRNA interactions. BMC Bioinformatics. 2020;21(8):1–21.
Xu W, et al. CeRNA regulatory network-based analysis to study the roles of noncoding RNAs in the pathogenesis of intrahepatic cholangiocellular carcinoma. Aging (Albany NY). 2020;12(2):1047.
Zhang Z, et al. Comparative transcriptomics identify key hypothalamic circular RNAs that participate in sheep (Ovis aries) reproduction. Animals. 2019;9(8):557.
Li X, et al. Comprehensive expression profiling analysis of pituitary indicates that circRNA participates in the regulation of sheep estrus. Genes. 2019;10(2):90.
La Y, et al. Differential expression of circular RNAs in polytocous and monotocous uterus during the reproductive cycle of sheep. Animals. 2019;9(10):797.
Li Z, et al. Analysis of expression profiles of CircRNA and MiRNA in oviduct during the follicular and luteal phases of sheep with two fecundity (FecB Gene) genotypes. Animals. 2021;11(10):2826.
Wang D, et al. The role of circRNA-SETD2/miR-519a/PTEN axis in fetal birth weight through regulating trophoblast proliferation. BioMed Res Int. 2020;2020:9809632.
Jia C, et al. Loss of hsa_circ_0118530 inhibits human granulosa-like tumor cell line KGN cell injury by sponging miR-136. Gene. 2020;744:144591.
Kallen AN, et al. The imprinted H19 lncRNA antagonizes let-7 microRNAs. Mol Cell. 2013;52(1):101–12.
Liu A, et al. Differential expression and functional analysis of CircRNA in the ovaries of low and high fecundity hanper sheep. Animals. 2021;11(7):1863.
Dong R, et al. Genome-wide analysis of long noncoding RNA (lncRNA) expression in hepatoblastoma tissues. PLoS ONE. 2014;9(1):e85599.
Lee JT, Bartolomei MS. X-inactivation, imprinting, and long noncoding RNAs in health and disease. Cell. 2013;152(6):1308–23.
Thore S. Structural views on the steroid receptor RNA activator. In: Acta Crystallographica A-foundation and advances. Chester: Int Union Crystallography; 2013.
Miao X, et al. Co-expression analysis and identification of fecundity-related long non-coding RNAs in sheep ovaries. Sci Rep. 2016;6(1):1–10.
Miao X, et al. Ovarian proteomic study reveals the possible molecular mechanism for hyperprolificacy of Small Tail Han sheep. Sci Rep. 2016;6(1):1–10.
Liu T, et al. Differential regulation of mRNAs and lncRNAs related to lipid metabolism in Duolang and Small Tail Han sheep. Sci Rep. 2022;12(1):11157.
Chang W, Cui Z, Wang J. Identification of potential disease biomarkers in the ovaries of Dolang sheep from Xinjiang using transcriptomics and bioinformatics approaches. Indian J Anim Res. 2021;1:8.
Zhang Z, et al. Transcriptome sequencing-based mining of genes associated with pubertal initiation in Dolang sheep. Front Genet. 2022;13:818810.
Chang W, Cui Z, Wang J. Identification of potential disease biomarkers in the ovaries of Dolang sheep from xinjiang using transcriptomics and bioinformatics approaches. Indian J Anim Res. 2020;10:412–9.
Zhong F, et al. Study on the polymorphism of BMPR-IB gene associated with litter size in small-tailed Han sheep and Xinjiang Duolang sheep. China Herbivore. 2005;25(6):15–6.
Mo F, et al. Polymorphisms in BMPRIB gene affect litter size in Chinese indigenous sheep breed. Anim Biotechnol. 2021;34(3):538–45.
Wen Y-L, et al. The expression and mutation of BMPR1B<? xmltex\break?> and its association with litter size in<? xmltex\break?> small-tail Han sheep (Ovis aries). Arch Anim Breed. 2021;64(1):211–21.
Wang H, et al. Genome -wide transcriptome profiling in ovaries of small-tail Han sheep during the follicular and luteal phases of the oestrous cycle. Anim Reprod Sci. 2018;197:212–21.
Zhou M, et al. Single nucleotide polymorphisms in the HIRA gene affect litter size in small tail Han sheep. Animals. 2018;8(5):71.
Wang C, et al. Genome-wide identification of mRNAs, lncRNAs, and proteins, and their relationship with sheep fecundity. Front Genet. 2021;12:750947.
Euesden J, Lewis CM, O’Reilly PF. PRSice: polygenic risk score software. Bioinformatics. 2015;31(9):1466–8.
Anders S. Analysing RNA-Seq data with the DESeq package. Mol Biol. 2010;43(4):1–17.
Enright A, et al. MicroRNA targets in Drosophila. Genome Biol. 2003;4(11):1–27.
Anders S, Huber W. Differential expression of RNA-Seq data at the gene level–the DESeq package. 2012;10:f1000research. Heidelberg, Germany: European Molecular Biology Laboratory (EMBL).
Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv preprint arXiv:1303.3997. 2013.
Gao Y, Wang J, Zhao F. CIRI: an efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol. 2015;16(1):1–16.
Ledford H. Circular RNAs throw genetics for a loop. Nat. 2013;494(7438):415.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):1–21.
Tay Y, et al. Coding-independent regulation of the tumor suppressor PTEN by competing endogenous mRNAs. Cell. 2011;147(2):344–57.
Kumar S, Stecher G, Tamura K. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol Biol Evol. 2016;33(7):1870–4.
Chen C, et al. TBtools: an integrative toolkit developed for interactive analyses of big biological data. Mol Plant. 2020;13(8):1194–202.
Malik WA, et al. Genome-wide expression analysis suggests glutaredoxin genes response to various stresses in cotton. Int J Biol Macromol. 2020;153:470–91.
Yousuf S, et al. Genome-wide expression profiling and networking reveals an imperative role of IMF-Associated Novel CircRNAs as ceRNA in Pigs. Cells. 2022;11(17):2638.
Wang X, et al. Differentially expressed bZIP transcription factors confer multi-tolerances in Gossypium hirsutum L. Int J Biol Macromol. 2020;146:569–78.
Li C, et al. Genome-wide analysis of circular RNAs in prenatal and postnatal pituitary glands of sheep. Sci Rep. 2017;7(1):1–10.
Li C, et al. Genome-wide analysis of circular RNAs in prenatal and postnatal muscle of sheep. Oncotarget. 2017;8(57):97165.
Jin C, et al. Changes in circRNA expression profiles related to the antagonistic effects of Escherichia coli F17 in lamb spleens. Sci Rep. 2018;8(1):1–11.
Li Y, et al. Circular RNA is enriched and stable in exosomes: a promising biomarker for cancer diagnosis. Cell Res. 2015;25(8):981–4.
Liu Y, et al. Differentially expressed circular RNA profile signatures identified in prolificacy trait of Yunshang black goat ovary at estrus cycle. Front Physiol. 2022;13:576.
Lin X, et al. Expression dynamics, relationships, and transcriptional regulations of diverse transcripts in mouse spermatogenic cells. RNA Biol. 2016;13(10):1011–24.
Qian Y, et al. Potential significance of circular RNA in human placental tissue for patients with preeclampsia. Cell Physiol Biochem. 2016;39(4):1380–90.
Tríbulo P, et al. Changes in the uterine metabolome of the cow during the first 7 days after estrus. Mol Reprod Dev. 2019;86(1):75–87.
Romero JJ, et al. Pregnancy-induced changes in metabolome and proteome in ovine uterine flushings. Biol Reprod. 2017;97(2):273–87.
Zielak-Steciwko AE, et al. Expression of microRNAs and their target genes and pathways associated with ovarian follicle development in cattle. Physiol Genomics. 2014;46(19):735–45.
Li W, et al. Circular RNA TGFBR2 acts as a ceRNA to suppress nasopharyngeal carcinoma progression by sponging miR-107. Cancer Lett. 2021;499:301–13.
Gao Y, et al. Identification and characterization of circular RNAs in Qinchuan cattle testis. Royal Soc Open Sci. 2018;5(7):180413.
Hammes A, et al. Role of endocytosis in cellular uptake of sex steroids. Cell. 2005;122(5):751–62.
Zhang H, et al. Expression regulation and physiological role of transcription factor FOXO3a during ovarian follicular development. Front Physiol. 2020;11:1413.
Miao X, Qin QLX. Genome-wide transcriptome analysis of mRNAs and microRNAs in Dorset and small tail Han sheep to explore the regulation of fecundity. Mol Cell Endocrinol. 2015;402:32–42.
Wang H, et al. Expression profile analysis of sheep ovary after superovulation and estrus synchronisation treatment. Vet Med Sci. 2022;8(3):1276–87.
Hoyer PB, Sipes IG. Assessment of follicle destruction in chemical-induced ovarian toxicity. Annu Rev Pharmacol Toxicol. 1996;36:307–31.
Bhattacharya P, Keating AF. Ovarian metabolism of xenobiotics. Exp Biol Med. 2011;236(7):765–71.
Nunes SC, Serpa J. Glutathione in ovarian cancer: a double-edged sword. Int J Mol Sci. 2018;19(7):1882.
Meister A. Glutathione metabolism and its selective modification. J Biol Chem. 1988;263(33):17205–8.
Wang H, et al. Genome-wide specific selection in three domestic sheep breeds. PLoS ONE. 2015;10(6):e0128688.
Fair T, Lonergan P. The role of progesterone in oocyte acquisition of developmental competence. Reprod Domest Anim. 2012;47:142–7.
Coelho AI, Berry GT, Rubio-Gozalbo ME. Galactose metabolism and health. Curr Opin Clin Nutr Metab Care. 2015;18(4):422–7.
Gentao L, Georgina EH, Claude LH. Galactose metabolism and ovarian toxicity. Reprod Toxicol. 2000;14(5):377–84.
Wu T, et al. Application of metabolomics in traditional chinese medicine differentiation of deficiency and excess syndromes in patients with diabetes mellitus. Evid Based Complement Altern Med. 2012;2012:968083.
Ruan S, et al. Compound fuling granule suppresses ovarian cancer development and progression by disrupting mitochondrial function, galactose and fatty acid metabolism. J Cancer. 2018;9(18):3382.
Kumar M, et al. Negative regulation of the tumor suppressor p53 gene by microRNAs. Oncogene. 2011;30(7):843–53.
Martha LS, et al. The p53-signaling pathway and colorectal cancer: Interactions between downstream p53 target genes and miRNAs. Genomics. 2019;111(4):762–71.
Buensuceso AV, Deroo BJ. The ephrin signaling pathway regulates morphology and adhesion of mouse granulosa cells in vitro. Biol Reprod. 2013;88(1):25, 1–12.
Hatzirodos N, et al. Transcriptome profiling of granulosa cells of bovine ovarian follicles during growth from small to large antral sizes. BMC Genomics. 2014;15(1):1–19.
Egawa M, et al. Ephrin B1 is expressed on human luteinizing granulosa cells in corpora lutea of the early luteal phase: the possible involvement of the B class Eph-ephrin system during corpus luteum formation. J Clin Endocrinol Metab. 2003;88(9):4384–92.
Xu H, et al. EPHA3 enhances macrophage autophagy and apoptosis by disrupting the mTOR signaling pathway in mice with endometriosis. Biosci Rep. 2019;39(7):BSR20182274.
Liu H, et al. Hypoxia-inducible factor-1α promotes endometrial stromal cells migration and invasion by upregulating autophagy in endometriosis. Reproduction (Cambridge, England). 2017;153(6):809.
Nishikimi M, et al. Segregation and pathfinding of callosal axons through EphA3 signaling. J Neurosci. 2011;31(45):16251–60.
Noren NK, Pasquale EB. Eph receptor-ephrin bidirectional signals that target Ras and Rho proteins. Cell Signal. 2004;16(6):655–66.
Wittenburg N, et al. Presenilin is required for proper morphology and function of neurons in C. elegans. Nature. 2000;406(6793):306–9.
Genna EM, Lucio M, Asgerally TF. Notch signaling in reproduction. Trends Endocrinol Metab. 2021;32(12):1044–57.
Goldbraikh D, et al. USP 1 deubiquitinates Akt to inhibit PI 3K-Akt-FoxO signaling in muscle during prolonged starvation. EMBO Rep. 2020;21(4):e48791.
Li L, et al. The signaling pathways involved in ovarian follicle development. Front Physiol. 2021;12:730196.
Saldanha SN, Tollefsbol TO. Pathway modulations and epigenetic alterations in ovarian tumorbiogenesis. J Cell Physiol. 2014;229(4):393–406.
D’Andrilli G, Giordano A, Bovicelli A. Epithelial ovarian cancer: the role of cell cycle genes in the different histotypes. Open Clin Cancer J. 2008;2:7–12.
Buchynska LG, et al. Expression of p53, p21WAF1/CIP1, p16INK4A and Ki-67 proteins in serous ovarian tumors. Exp Oncol. 2007;29(1):49–53.
Hur C-G, et al. Expression and localization of two-pore domain K+ channels in bovine germ cells. Reproduction. 2009;137(2):237–44.
Chu L-H, et al. Epigenomic analysis reveals the KCNK9 potassium channel as a potential therapeutic target for adenomyosis. Int J Mol Sci. 2022;23(11):5973.
Kim J-M, et al. Role of potassium channels in female reproductive system. Obst Gynecol Sci. 2020;63(5):565–76.
Świątkowska-Stodulska R, et al. Renin-angiotensin-aldosterone system in the pathogenesis of pregnancy-induced hypertension. Exp Clin Endocrinol Diabetes. 2018;126(06):362–6.
Gennari-Moser C, et al. Regulation of placental growth by aldosterone and cortisol. Endocrinology. 2011;152(1):263–71.
Sabbadin C, et al. Aldosterone in gynecology and its involvement on the risk of hypertension in pregnancy. Front Endocrinol. 2019;10:575.
Hassan E, et al. Clinical Implications of the ovarian/endometrial renin-angiotensin-aldosterone system. Ann N Y Acad Sci. 2000;900(1):107–17.
D’Occhio MJ, et al. Adhesion molecules in gamete transport, fertilization, early embryonic development, and implantation—role in establishing a pregnancy in cattle: a review. Mol Reprod Dev. 2020;87(2):206–22.
Wu R, et al. Macrophage contributions to ovarian function. Hum Reprod Update. 2004;10(2):119–33.
Carlock C, et al. Ovarian phagocyte subsets and their distinct tissue distribution patterns. Reproduction. 2013;146(5):491–500.
Duffy DM, et al. Ovulation: Parallels With Inflammatory Processes. Endocr Rev. 2019;40(2):369–416.
Ochota M, Wojtasik B, Niżański W. Survival rate after vitrification of various stages of cat embryos and blastocyst with and without artificially collapsed blastocoel cavity. Reprod Domest Anim. 2017;52(Suppl 2):281–7.
Murdoch WJ, Wilken C, Young DA. Sequence of apoptosis and inflammatory necrosis within the formative ovulatory site of sheep follicles. J Reprod Fertil. 1999;117(2):325–9.
Oakley OR, et al. Periovulatory leukocyte infiltration in the rat ovary. Endocrinology. 2010;151(9):4551–9.
Caillaud M, Gérard N. In vivo and in vitro effects of interleukin-1beta on equine oocyte maturation and on steroidogenesis and prostaglandin synthesis in granulosa and cumulus cells. Reprod Fertil Dev. 2009;21(2):265–73.
Townson DH, Liptak AR. Chemokines in the corpus luteum: implications of leukocyte chemotaxis. Reprod Biol Endocrinol. 2003;1:94.
Zhao X, et al. Comprehensive analysis of differentially expressed profiles of mRNA, lncRNA, and miRNA of Yili geese ovary at different egg-laying stages. BMC Genomics. 2022;23(1):607.
Henson PM, Hume DA. Apoptotic cell removal in development and tissue homeostasis. Trends Immunol. 2006;27(5):244–50.
Regan SLP, et al. Granulosa cell apoptosis in the ovarian follicle-A changing view. Front Endocrinol (Lausanne). 2018;9:61.
Erwig LP, Henson PM. Immunological consequences of apoptotic cell phagocytosis. Am J Pathol. 2007;171(1):2–8.
Nagata S, Hanayama R, Kawane K. Autoimmunity and the clearance of dead cells. Cell. 2010;140(5):619–30.
This work was supported by a grant from the National Natural Science Foundation of China (No. 31970541), The Major Science and Technology Project of New Variety Breeding of Genetically Modified Organisms (Nos. 2009ZX08008-004B), the Agricultural Science and Technology Innovation Program (NO. ASTIP-IAS05), the Basic Research Fund for Central Public Research Institutes of CAAS (Y2016JC22, Y2018PT68), and (2013ywf-yb-5, 2013ywf-zd-2).
Ethics approval and consent to participate
All the procedures involving animals were approved by the animal care and use committee at the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (NO. IAS2019-82), where the study was conducted. All the experiments were performed in accordance with the relevant guidelines and regulations set by the Ministry of Agriculture of the People’s Republic of China.
Consent for publication
The authors declare no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Correlation and Target Found between miRNA and mRNA.
Correlation of circRNA-mRNA.
ceRNA score of ceRNA regulatory pairs.
Gene ontology of mRNA in ceRNA regulatory network.
Primer Sequence of circRNAs, miRNAs, and mRNAs for qRT-PCR.
About this article
Cite this article
Yousuf, S., Malik, W.A., Feng, H. et al. Genome wide identification and characterization of fertility associated novel CircRNAs as ceRNA reveal their regulatory roles in sheep fecundity. J Ovarian Res 16, 115 (2023). https://doi.org/10.1186/s13048-023-01178-2