Skip to main content
Download PDF

Identification and expression profile analysis of circRNAs associated with goat uterus with different fecundity during estrous cycle

BMC Genomics volume 26, Article number: 349 (2025) Cite this article

Abstract

Background

The Yunshang Black Goat, a distinguished meat goat breed native to China, is renowned for its superior reproductive capabilities. Despite this, there is considerable phenotypic variability within the breed. During the reproductive cycle, the uterus plays a pivotal role, with its functions evolving in line with the different stages of the cycle. This study focuses on the uterine tissues, including both the endometrium and myometrium, of Yunshang Black Goats with high fecundity (HF) and low fecundity (LF) during the proliferative (FP) and secretory (LP) phases of the estrous cycle. By examining these tissues, we aim to elucidate the underlying molecular and physiological mechanisms of the observed differences in reproductive success.

Results

High-throughput sequencing was conducted, followed by bioinformatics analysis to identify the expression profiles of circRNAs. A total of 7,445 circRNAs were identified through the integration of findings from find_circ and CIRI2 software. Comparative analyses between the FPLF vs. FPHF and LPLF vs. LPHF revealed 149 differentially expressed (DE) circRNAs (94 up-regulated and 55 down-regulated) and 276 DE circRNAs (56 up-regulated and 220 down-regulated), respectively. The enrichment analysis indicated that the primary pathways involved were the Sphingolipid signaling pathway, MAPK signaling pathway, and GnRH signaling pathway, all of which are closely associated with cellular growth and development. Additionally, several key candidate genes were identified, such as FGF2 and MBTPS1. We also predicted a total of 281 miRNA-circRNA binding pairs, encompassing 263 circRNAs and 60 miRNAs, and simultaneously, 14 coding circRNAs were anticipated.

Conclusion

Based on the analysis, we have established the expression profiles of circRNAs during the follicular and luteal phases, respectively. Furthermore, using various analytical methods and data from high- and low-yield experimental control groups over different periods, we have identified multiple circRNAs that affect the high reproductive capacity of goats. Through enrichment analysis of the host genes of these circRNAs, we have discovered several key candidate genes. These findings provide fundamental data for the study of the molecular mechanisms underlying the fecundity of goats and pave the way for future genetic improvement strategies.

Peer Review reports

Introduction

The goat (Capra hircus) as one of the earliest domesticated livestock species, is renowned for its exceptional adaptability and versatility, making it a valuable livestock species worldwide. [1, 2]. It is believed to have originated from a domestication process involving various wild bezoar populations (Capra aegagrus) approximately 11,000 years ago [3]. Since then, complex intraspecific hybridization, especially the genetic contribution of domestic goats, has contributed to the development of diverse goat breeds [4,5,6]. The Yunshang black goat, a distinguished meat goat breed from Yunnan Province, China [7], is a crossbreed of the Nubi and Yunling goats. This breed is characterized by its superior meat production performance, robust adaptability, and notable fecundity [8]. Given that fertility directly impacts the economic efficiency of livestock, understanding the mechanisms behind the high reproductive output of the Yunshang black goat is of paramount importance to those in the livestock industry.

Numerous factors influence the number of offspring produced by goats, with much research focusing on the differential development of ovaries and the variance in ovulated follicle numbers [9]. However, in ruminants, pregnancy loss is a significant factor affecting female fertility [10, 11], which predominantly occurs during the implantation period. Embryo implantation is the process by which a blastocyst-stage embryo establishes structural and physiological connections with the endometrium after a period of free existence within the uterus [12]. The attachment phase of pregnancy is pivotal in the reproductive process of female animals, with the potential to enhance reproductive efficiency by increasing embryonic implantation rates and offspring numbers [13]. The determinants of attachment include embryo viability, endometrial receptivity, and the interaction between the embryo and the mother. It is reported that approximately 75% of natural conception pregnancy failures are due to an unsuitable uterine environment [14, 15]. Thus, the study of the female reproductive system, particularly the uterus, is crucial for improving offspring numbers.

The endometrial stromal and epithelial cells, influenced by estrogen and progesterone, proliferate and differentiate, thus becoming receptive to embryonic adhesion [13]. Numerous genes are involved in regulating this process, with Labarta et al. describing 140 endometrium-related genes, including GPX3, PAEP, and LIF [16]. Considering the dynamic expression of non-coding RNA (ncRNA) during the establishment of endometrial receptivity, it is essential to explore the regulatory role of ncRNA in this process. Circular RNAs (circRNAs) are a class of endogenous non-coding RNAs characterized by their closed loop structure, lacking free 5' or 3' ends, and are rich in mRNA regulatory elements (MREs), enabling them to perform various biological functions [13]. Study has indicated that circ-8073 acts as a molecular sponge for miR-449a, regulating centrosomal protein 55 (CEP55) and promoting the proliferation of endometrial epithelial cells via the PI3K/AKT/mTOR pathway [17]. CircRNA-9119 downregulates miR-26 expression, which in turn affects the endometrial epithelial cells of dairy goats in vitro by targeting the predicted sites and reducing the expression of prostaglandin-endoperoxide synthase 2 (PTGS2), thereby enhancing endometrial receptivity [18]. Additionally, Shen et al. confirmed that circ-BACH1 may influence endometrial receptivity as a competing endogenous RNA (ceRNA) for miRNAs [19]. Ma et al. discovered that circ-9110 could enhance the expression of the homeobox A1 gene (HOXA1) by sequestering miR-100-5p, subsequently regulating the PI3K/AKT/mTOR and ERK1/2 pathways to induce apoptosis in goat uterine stromal cells, improve uterine receptivity, and ultimately facilitate the process of embryo implantation [20]. In goats, circRNA has been identified as a key regulator of ovarian follicle development, and studies have shown that there are differences in circRNA expression in large and small follicles in Leizhou goats, indicating their potential role in follicle development [21]. In addition, circRNA is also involved in regulating the proliferation of goat granulosa cells. circCFAP20DC significantly promotes the transition of goat granulosa cells from the G1 phase to the S phase, thereby enhancing their proliferation [22]. These results provide valuable insights into the molecular regulation of mammalian endometrium and pregnancy in the future.

The aim of this manuscript is to investigate the expression profiles and regulatory roles of circRNAs in the uterus of Yunshang black goats with different fecundity levels during the estrous cycle. Specifically, we aim to identify differentially expressed circRNAs and elucidate their potential molecular mechanisms in regulating uterine receptivity and embryo implantation, thereby providing fundamental data for understanding the molecular basis of fecundity in goats.

Results

Results of goat uterine tissue morphology test

Collected uterine tissues were paraffin-embedded and stained, with the results depicted in Fig. 1. Morphological analysis revealed distinct differences in endometrial morphology between the follicular and luteal phases, evident in gland count, endometrial thickness, and uterine caruncle shape. The distance from the endometrial surface to the highest point of the caruncle was measured to determine uniform thickness. Statistical analysis revealed significant differences in gland count and endometrial thickness between the follicular and luteal phases, with the luteal phase having approximately twice the number of glands and endometrial thickness compared to the follicular phase (Fig. 1B, E). However, no significant differences were observed between high- and low-yield groups, although there was a trend towards higher gland counts and endometrial thickness in the high-yield group.

Fig. 1
figure 1

Histological sections of uterine tissue and statistical analysis of endometrial thickness and gland number. The figure shows H&E stained histological sections of uterine tissue from low and high yield samples at different stages. 'L' denotes the uterine lumen, 'Uc' represents the uterine caruncle, the black arrows indicate endometrial glands, 'En' stands for the endometrium, 'My' is the myometrium, 'Ua' signifies uterine arteries, and 'Pe' is the perimetrium. The circles and squares represent the uterine caruncle during the follicular and luteal phases, respectively. Scale bar: 500 Î¼m. 'NS' indicates no significant difference, while '***' indicates a highly significant difference (p < 0.001)

Full size image

Summary of sequencing quality control results in uterine tissue

Twenty cDNA libraries were constructed from uterine tissue samples of goats with high and low fecundity, collected during both the proliferative and secretory phases. High-throughput sequencing was performed using an Illumina Hiseq 2500 sequencer. The samples were divided into four groups, each with five replicates: the high-fecundity proliferative phase group (FPHF-1 through −5), the low-fecundity proliferative phase group (FPLF-1 through −5), the high-fecundity secretory phase group (LPHF-1 through −5), and the low-fecundity secretory phase group (LPLF-1 through −5). Before analyzing the sequencing data, the raw data were subjected to a stringent quality control and filtering process. The results of the filtering are shown in Fig. 1. The sequencing quality was assessed in three aspects: inspection of the base quality value distribution (Fig. 2A), assessment of the GC content distribution (Fig. 2B), and examination of the base content distribution (Fig. 2C) to ascertain the final clean reads. The results showed that the sequencing quality met the standards, and the clean reads were suitable for subsequent bioinformatics analysis.

Fig. 2
figure 2

Quality assessment of sequencing data. A Each base in the sequencing data has a corresponding quality value, which directly determines the accuracy of sequencing. B Distribution of GC content in the sequencing data. The GC content affects the efficiency of PCR amplification during library construction, so the sequencing process has a certain preference for fragments with different GC content. However, the overall GC content of the sequencing results should be consistent with the GC content of all expressed genes in the species. C Distribution of base content in the sequencing data. The distribution of ATGC can reflect the normality of the sequencing to some extent. D Schematic of the basic principle of find_circ software

Full size image

The comprehensive evaluation results for the clean reads are presented in Table 1. On average, the quality control process generated 121.1 million clean reads. Of these, 96.88% successfully aligned with the Capra hircus reference genome (ARS1). Additionally, the proportion of sequencing reads that aligned to multiple positions on the reference sequence averaged 5.39%, while those that aligned to unique positions accounted for 91.50%. Since the primary goal of this study is to identify circular RNAs (circRNAs), the number of junction reads, which are essential for circRNA identification, was also recorded, averaging 28.04%. Together, these clean reads meet the criteria for subsequent analysis. Further filtering criteria were applied to identify circRNA candidates, with details provided in Supplementary Table S1.

Table 1 Summary of mapping data from the uterine tissues
Full size table

Identification of circRNA in uterine tissues

A comprehensive analysis using both find_circ and CIRI2 software identified a total of 7,445 circRNAs after mapping to reference sequences and excluding linear RNA species (refer to Supplementary Table S2, Fig. 2D for details). All identified circRNAs were classified and annotated. The predominant category of circRNAs originated from exons, with each circRNA potentially containing one or more exons. Additionally, circRNAs were detected within intronic and intergenic regions, where multiple circRNAs could originate from the same genomic locations. Consistent with previous findings, this study revealed that over 95% of circRNAs were exon-derived on average (Fig. 3A).

Fig. 3
figure 3

Summary of circRNA identification. A Based on the identification results, the sources of circRNAs were mainly divided into exon, intron, and intergenic regions, and pie charts were generated to show their proportions. B Length distribution of circRNAs in four groups of uterine tissues. C Chromosome distribution of identified circRNAs in uterine tissues

Full size image

Length calculations indicated that the sequence lengths of circRNAs predominantly ranged from 200 base pairs (bp) to 10,000 bp (Fig. 3B). A chromosome-wise density analysis revealed that circRNAs were primarily distributed across chromosomes 1 through 11, with the highest density observed on chromosomes 1, 2, 3, and 11, which collectively accounted for approximately 22% of the circRNAs (Fig. 3C and Supplementary Table S2).

To provide a more nuanced view of the distribution, a detailed radar map of circRNA density was constructed, illustrating the distribution patterns of circRNAs among the 20 individuals across the ten chromosomes with the densest circRNA presence (Supplementary Fig. 1). This visualization aids in understanding the genomic landscape of circRNAs and their potential regulatory roles in the context of goat uterine tissue.

Results of expression analysis of DEcircRNA

Given the unique nature of circRNAs, their expression levels are conventionally estimated by counting the number of junction reads that span their splicing sites. For this study, circRNA expression levels were quantified for each sample and normalized to spliced reads per billion mapping (SRPBM). Violin plots depicting SRPBM for all transcripts were used to compare expression levels across various samples. The uniform widths of the violin plots in Fig. 4A indicate equivalent numbers of transcripts among the samples (Supplementary Table S3).

Fig. 4
figure 4

Analysis of circRNA expression levels: (A) Violin plot comparing SRPBM expression levels across different experimental groups. B Correlation matrix showing the expression level relationships between samples; a correlation coefficient closer to 1 indicates greater similarity in expression patterns. C Differential circRNA volcano plot for high- and low-yield groups during the proliferative phase; the topic shown here is host gene. D Differential circRNA volcano plot for high- and low-yield groups during the secretory phase. E Union of DEcircRNAs from each group

Full size image

Additionally, correlation analysis was performed on the samples based on their expression profiles to validate the scientific rationale of the sample grouping and the overall reliability of the experiment from a molecular standpoint. The clustering results revealed that samples FPLF-5, LPHF-4, and LPLF-2 did not conform to the grouping criteria, whereas the remaining samples were suitable for further analysis (Fig. 4B).

Differential expression analysis was used to identify circRNAs with significant expression differences between comparison groups. The genes targeted by these differentially expressed circRNAs may contribute to distinct phenotypes. Volcano plots were used to visually represent the number of differentially expressed circRNAs identified through various comparisons, primarily based on log2 fold change and p-value. A total of 149 differentially expressed (DE) circRNAs (94 up-regulated and 55 down-regulated) and 276 DE circRNAs (56 up-regulated and 220 down-regulated) were identified in the comparisons of FPLF versus FPHF and LPLF versus LPHF, respectively (Fig. 4C, D; Supplementary Table S4). The circRNAs derived from individual samples were consolidated to determine the differential expression profiles across various treatment cohorts. As depicted in Fig. 4E, the comparative analysis reveals distinct expression patterns between the disparate groups, which are both visually evident and statistically significant.

Enrichment analysis of the host gene for circRNA

Enrichment analysis of differentially expressed genes is crucial for deciphering gene functions and their broader biological implications. In this study, we performed enrichment analyses on differentially expressed genes across various comparison groups, identifying distinct counts of Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Enrichment was considered significant for GO terms with a Q-value ≤ 0.05 and for KEGG pathways with a P-value ≤ 0.05, based on our established screening criteria.

In the FPLF vs. FPHF comparison group, we identified 3,754 GO terms, of which 19 terms showed significant enrichment. These included categories such as biological process (n = 4), cellular components (n = 8), and molecular functions (n = 7). Similarly, in the LPLF vs. LPHF comparison, we identified 4,133 GO terms, of which 4 terms were significantly enriched. The top ten GO terms for biological process (BP), cellular components (CC), and molecular functions (MF) from both comparisons are depicted in Fig. 5 (upper panel). Notably, GO:0005622 and GO:0044424 were significantly enriched in the FPLF vs. FPHF comparison and also featured among the top ten terms in the LPLF vs. LPHF group. Although their Q-values exceeded our threshold, their P-values were remarkably low, as low as 0.0005. Ultimately, 26 terms were retained in the low-yield group, including biological processes (n = 14), cellular components (n = 5), and molecular functions (n = 7) (Supplementary Table S5).

Fig. 5
figure 5

GO term and KEGG enrichment analysis of the host gene for circRNA. Display content is consistent with the title of each figure

Full size image

For KEGG pathway enrichment analysis, a strategy analogous to GO enrichment was used. We identified 154 pathways in the FPLF vs. FPHF comparison and 158 pathways in the LPLF vs. LPHF comparison. The 20 most enriched pathways from each comparison were selected to create the KEGG bubble map (Fig. 5, lower panel). Among these, the Sphingolipid signaling pathway was significantly enriched in both comparisons. Ultimately, 12 pathways were retained in the high-yield group and 30 pathways in the low-yield group (Supplementary Table S6).

Gene set enrichment analysis for source mRNA of circRNA

Based on the KEGG enrichment results, we categorized the host genes of differentially expressed circRNAs and established distinct gene sets for further analysis. In the FPLF vs. FPHF comparison, 154 gene sets were formulated, excluding those with fewer than five genes. Ultimately, Gene Set Enrichment Analysis (GSEA) was performed on 39 gene sets. The analysis revealed that 10 gene sets were up-regulated in the FPLF group, but none met the significance threshold (P-value < 0.05, FDR < 0.25). Conversely, 10 gene sets were up-regulated in the FPHF group, with the gene set KO04141 significantly enriched (P-value = 0.03, FDR < 0.23; Fig. 6). The KEGG pathway associated with this gene set is "Protein processing in the endoplasmic reticulum," with genes MAPK8, MBTPS1, SIL1, and NPLOC4 identified as the core enrichment genes.

Fig. 6
figure 6

GSEA enrichment analysis for each pathway. Nominal p-value: represents the p-value and characterizes the credibility of the enrichment results. FDR stands for Q-value, which is the p-value corrected by multiple hypothesis tests. Note that GSEA uses p-value < 0.05 and q-value < 0.25 to filter the results. ES stands for Enrichment score, and Normalized ES represents the normalized Enrichment score

Full size image

Similarly, for the LPLF vs. LPHF comparison, the 158 KEGG-enriched pathways were transformed into gene sets, and 27 gene sets were selected for GSEA analysis. Eight gene sets were up-regulated in LPHF, with only one gene set, KO04360 (Axon guidance), crossing the significance threshold (P-value = 0.03, FDR = 0.24). The genes TRPC6, TRPC1, PAK3, and PTK2 were identified as the core enriched genes within this set.

Subsequently, 19 gene sets were up-regulated in the LPLF group, with only those surpassing the significance threshold depicted in Fig. 6. Genes such as FGF2, FLT1, PRKCB, RAP1B, VEGF, CREB5, PRKCE, and ARHGDIB recurred across various pathways, indicating their pivotal roles. The analysis concluded that these genes are central to the enrichment of the respective pathways and may influence the biological processes associated with fecundity in goats. These findings underscore the complexity of gene regulation and the multifaceted roles of non-coding RNAs in gene expression and cellular function (Supplementary GSEA).

Interaction regulation network about miRNA-circRNA

CircRNA molecules, rich in microRNA (miRNA) binding sites, act as miRNA sponges in cells, thereby alleviating the inhibitory effects of miRNAs on their target genes. Therefore, we conducted miRNA binding site analysis for the identified circRNAs to further investigate their functions. We predicted 281 miRNA-circRNA binding pairs, comprising 263 circRNAs and 60 miRNAs. These binding pairs include differentially expressed circRNAs from different comparison groups, such as FPLF vs. FPHF (3 up-regulated, 1 down-regulated) and LPLF vs. LPHF (1 up-regulated, 5 down-regulated), as shown in Fig. 7 and Supplementary Table S7.

Fig. 7
figure 7

Interaction regulation network of miRNA-circRNA. A Targeted miRNAs were predicted based on differentially expressed circRNAs. B The square represents miRNA, the blue outer circle represents downregulated circRNA, and the red outer circle represents upregulated circRNA. The thickness of the wire represents the absolute value of the predicted interaction correlation, while the color inside the circle indicates the p-value

Full size image

Identifying the coding potential of circRNAs

Although circRNA has been widely studied as a non-coding RNA, some studies have found that its coding ability is more stable than that of mRNA. Therefore, exploring the coding ability of circRNA is very important. The study of coding ability first requires IRES (internal ribosome entry site) prediction. We predicted IRES sites in all 7445 identified circRNAs, and 5014 circRNAs were found to have this site. Finally, we used three software tools to predict coding ability, and through CNCI analysis, 3027 circRNAs were found to have coding ability. CPC2 and CAPT identified relatively few circRNAs with coding potential, with 204 and 1001 circRNAs identified, respectively. Finally, the results of the four software tools were integrated, and the 14 condensable circRNAs jointly identified are shown in Fig. 8B.

Fig. 8
figure 8

Predicted results for the coding potential of circRNAs. A Results predicted by the four software tools. The second half shows the statistics of the total results predicted by each software. B Specific information of the 14 circRNAs identified by all four software tools after overlap

Full size image

Validation of RNA sequencing using RT-qPCR

To verify the accuracy and stability of the sequencing results, RT-qPCR was used to validate the RNA-seq data. The results showed that the quantitative results of circRNA were consistent with the sequencing results (Fig. 9), thereby confirming the reliability of the sequencing data. Relative gene expression was calculated using the 2−△△ct method.

Fig. 9
figure 9

Comparison of real-time quantitative PCR (RT–qPCR) and RNA-seq results about validation of differentially expressed circRNA, *p < 0.05,**p < 0.01

Full size image

Discussion

The reproductive traits of goats, which are economic traits, involve complex physiological processes and are controlled by various factors such as environment and genes [23, 24]. Currently, improving fecundity is the main breeding goal in goat breeding. With the rapid development of high-throughput sequencing technology, molecular marker-assisted breeding has been significantly enhanced. As a research hotspot in recent years, circRNA plays a significant role in various molecular processes in animals and plants [25]. It is an important regulator of gene expression, involved in transcriptional regulation, post-transcriptional translation, and affecting miRNA expression and other activities [26]. Several studies have analyzed the expression profiles of circRNAs in goat ovaries at different stages, identifying many differentially expressed circRNAs. Among them, key circRNAs can act as sponges for miRNAs, thereby regulating their expression and affecting gene transcription [27,28,29]. Additionally, the circRNA expression profile of the pituitary glands of sheep and goats was studied, and the differential circRNAs of sheep and goats with different fecundity were explored, providing a theoretical basis for high-yield research[30, 31].

The goat uterus is one of the most important organs of the reproductive axis. We fully excavated the circRNA expression profile of the goat uterus by high-throughput sequencing. A total of 7445 circRNAs were identified in this study, which was nearly 1400 and 3400 more than those identified in the ovary and fallopian tube by Liu et al., who previously studied Yunshang black goat [29, 32]. Additionally, Song et al. identified 21,813 circRNAs in the endometrial circRNA expression profile of Xinong Saanen dairy goats during early pregnancy [33]. In our study, 95% of the circRNAs in each comparison group originated from exon splicing, but in Song's study, the proportion of exon splicing was significantly lower, only nearly 50%. We infer that this difference may be due to different bioinformatics analysis methods. Song et al. mainly used TopHat-Fusion [34] for prediction and the FPKM normalization method, which is different from our software (find_circ, CIRI2) and SRPBM normalization method [25, 35, 36]. In summary, these differences, whether among different tissues of the same variety or between different varieties of the same tissue, suggest that circRNA plays a significant role in reproduction and exhibits species-tissue specificity.

A total of 149 DE circRNAs (94 up-regulated and 55 down-regulated) and 276 DE circRNAs (56 up-regulated and 220 down-regulated) were identified in the FPLF vs. FPHF and LPLF vs. LPHF comparisons, respectively. After comparing the overlap of up-regulated and down-regulated genes between the two periods, we found that 7 circRNAs were significantly differentially expressed in both the proliferative and secretory phases. They are derived from genes ROS1, VWA9, FGF2, MBTPS1, PRKD1, and PICALM, respectively. Orphan receptor tyrosine kinase ROS1 has a unique extracellular domain, consisting of fibronectin III repeat (FN-III), β propeller module (YWTD repeat) domain, and intracellular tyrosine kinase domain [37]. Study has shown that ROS1 plays an important role in epithelial cell differentiation and cell-to-cell signal communication [38]. The proteins encoded by the Vwa9 gene belong to the von Willebrand factor type A (vWA) domain family [39]. Proteins containing these domains are involved in various biological events, such as cell adhesion and migration [40]. We are all familiar with FGF2 genes, most of which are closely related to cell proliferation [41]. Through genome-wide association analysis of dairy goats, researchers found that the FGF2 gene contributes to animal ovary development and promotes corpus luteum angiogenesis and vascular endothelial cell differentiation [42]. There are special studies on the changes in FGF2 gene expression in the endometrium of goats during pregnancy. The data show that FGF2 and FGFR2 may play an important role in regulating uterine angiogenesis in the early stage of pregnancy [43, 44]. The differentially expressed circRNAs affect reproductive rates by influencing embryo implantation, growth, and differentiation. For example, circ-8073 acts as a molecular sponge for miR-449a, thereby regulating the proliferation of endometrial epithelial cells and affecting endometrial receptivity [17]. Additionally, circRNA-9119 downregulates miR-26 expression, thereby affecting endometrial epithelial cells by targeting PTGS2 and enhancing endometrial receptivity [18]. These findings underscore the role of circRNAs in modulating the uterine environment, thereby influencing embryo implantation rates and offspring numbers.

MBTPS1 is a member of a larger family of pro-protein convertases subtilisin kexins, responsible for the cleavage and activation of various protein substrates [45]. It is mainly located in the cis/medial Golgi apparatus, which is responsible for the cleavage of membrane-binding transcription factors and activation of proproteins [46]. Some studies have shown that it is expressed in the placental syncytiotrophoblast, which is the main component of the placenta and is mainly related to the mother [47]. Therefore, the differential expression of this gene in different reproductive groups is an important candidate gene to explore high fecundity. Serine/threonine-protein kinase D1 is an enzyme encoded by the PRKD1 gene, a member of the protein kinase D (PKD) family, and plays a role in many extracellular receptor-mediated signal transduction pathways [48, 49]. Several studies have found that PRKD1 is related to the first estrus in sows and plays an important role in affecting litter size [50, 51]. Steroids play a very important role in mammalian reproduction, controlling the transcription of most reproduction-related genes. Some studies have shown that PICALM, a receptor for clathrin-mediated endocytosis [52], is sensitive to estrogen in reproductive tissues, so whether it is related to the phenotype of our study needs to be further verified [53].

We then enriched and analyzed the host genes of differentially expressed circRNAs, including GO, KEGG, and GSEA. In the FPLF vs. FPHF comparison group, GO enrichment identified major cellular components such as cellular component organization, intracellular organelles, and organelle-related pathways, indicating that the main reason for the difference in high and low yield in the proliferative phase is the differential differentiation of uterine tissue cells that needs to be further explored. In the LPLF vs. LPHF comparison group, GO enrichment identified the main GTPase activator activity and GTPase regulator activity pathways. GTPases are a large family of hydrolase enzymes that bind to guanosine triphosphate (GTP) and hydrolyze it to guanosine diphosphate (GDP) [54]. They are involved in various cell processes and can respond to signal transduction activated by cell surface receptors and regulate cell differentiation and proliferation [55]. Among the overlapping GO pathways in our two comparison groups, only the intracellular components were significantly enriched. In the KEGG enrichment analysis, the pathways enriched in the proliferative phase included the Sphingolipid signaling pathway, MAPK signaling pathway, GnRH signaling pathway, and Apoptosis-fly, while the main enriched pathways in the secretory phase included Focal adhesion and Sphingolipid signaling pathway. Both analyses included the Sphingolipid signaling pathway. Sphingolipids are one of the important components of cell membranes, and the Sphingolipid signaling pathway is one of the most important secondary signal transduction systems in cells [56]. In particular, the balance regulation between ceramide and sphingosine 1-phosphate (S1P) is the rheostat of cell growth and programmed cell death [57]. Some studies have shown that S1P participates in Ca2+ movement and causes human airway smooth muscle contraction by phosphorylating myosin light chain kinase and forming stress fibers [58]. It is well known that the myometrium is a type of smooth muscle, and its role during pregnancy is self-evident. Thus, whether the movement ability of uterine smooth muscle in ewes can explain high fecundity is worth exploring. In the GSEA enrichment analysis, we defined the gene sets based on KEGG results to identify hub genes in the pathways. Finally, by overlapping different genes, we identified FGF2 and MBTPS1 as key candidate genes.

In addition to elucidating the functions and pathways of the host genes of differentially expressed circRNAs, we further explored whether these circRNAs possess miRNA binding sites and can function as miRNA sponges to modulate mRNA expression. A total of 281 miRNA-circRNA binding pairs were identified. Among them, circ-ARHGEF28 has been demonstrated to significantly drive bladder cancer (BC) progression by acting as a miRNA sponge. Specifically, circ-ARHGEF28 sequesters miR-548, thereby alleviating the inhibitory effects of miR-548 on its target gene KIF2C [59]. This interaction is particularly critical because miR-548 has been widely reported to function as a tumor suppressor in various cancers, including breast cancer, prostate cancer, and pancreatic cancer [60]. By sequestering miR-548, circ-ARHGEF28 effectively upregulates KIF2C expression, which in turn promotes cell proliferation, migration, and invasion in bladder cancer cells [59]. Similarly, circ-CSPP1 has been identified as a key player in both hepatocellular carcinoma (HCC) and renal cell carcinoma (RCC) by acting as a sponge for miR-493-5p. In HCC, circ-CSPP1 sequesters miR-493-5p, leading to the upregulation of HMGB1, a protein involved in inflammation and carcinogenesis [61]. The interaction between circ_CSPP1 and miR-493-5p is significant importance because miR-493-5p has been shown to suppress tumor progression by targeting oncogenic genes [62]. By sequestering miR-493-5p, circ-CSPP1 enhances HMGB1 expression, thereby promoting HCC development [61]. In RCC, circ-CSPP1 similarly sequesters miR-493-5p, resulting in the upregulation of RAC1, a protein involved in cell proliferation and metastasis [63].

Because circRNA has a ring structure composed of exons, some circRNAs contain open reading frames and IRES. Therefore, researchers speculate that under certain conditions, a large amount of circRNA can be converted into protein [64]. Fourteen circRNAs with coding potential were identified using IRESfinder, CPC2, CNCI, CPAT, and PLEK software. It is worth noting that the biological functions of these circRNAs need to be further explored from the perspective of protein translation.

Conclusions

This study identifies and characterizes differentially expressed circRNAs in the uterine tissues of Yunshang black goats with varying fecundity levels. Our findings highlight the involvement of circRNAs in key biological pathways related to cell growth and development, such as the Sphingolipid signaling pathway, MAPK signaling pathway, and GnRH signaling pathway. The identification of potential miRNA-circRNA interactions and coding circRNAs further underscores their regulatory roles in reproductive processes. These results provide valuable insights into the molecular mechanisms underlying fecundity in goats and pave the way for future genetic improvement strategies.

Materials and methods

Animals and ethics statement

Ethics approval (No. IAS2019-63) was granted by the Animal Ethics Committee of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (IAS-CAAS), Beijing, China. A cohort of 20 non-pregnant, adult Yunshang black doe-goats, aged 3 to 5 years and weighing 41 to 63 kg, was randomly selected from the Yixingheng Animal Husbandry Technology Co., Ltd's Tuanjie Township Base in Kunming City, Yunnan Province, China (24° 23′ North latitude).

These goats were divided into two groups based on their historical kidding records: the high-fecundity group (HF group, n = 10) with an average kidding number of 3.40 ± 0.42, and the low-fecundity group (LF group, n = 10) with an average kidding number of 1.80 ± 0.27. The difference in kidding numbers between the two groups was statistically significant (p < 0.01, T-test).

To synchronize estrus, all goats were administered a Controlled Internal Drug Releasing device (CIDR) containing 300 mg of progesterone (Inter Ag Co., Ltd., New Zealand). The CIDRs were removed after 16 days, marking the beginning of the experimental timeline at 0 h. Five goats from each group were humanely euthanized at 48 h (proliferative phase, FP groups), and the remaining ten goats were euthanized at 240 h (secretory phase, LP groups). For euthanasia, the goats were first anesthetized with 100 mg/kg sodium pentobarbital administered intravenously. Following anesthesia, the animals were exsanguinated to ensure a rapid and humane death. Blood samples were collected from these goats for hormone level determination using the competitive radioimmunoassay method, ensuring rigorous and scientifically valid data collection and analysis.

Sample collection, total RNA preparation and HE staining

Uterine tissues from the goats were immediately frozen in liquid nitrogen and stored at −80 Â°C for subsequent RNA extraction. Following established protocols from prior research [65]. total RNA for RNA sequencing was extracted from 20 uterine tissue samples, including both endometrium and stroma, using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's guidelines. The purity and concentration of the RNA samples were assessed using the NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA), and their integrity was confirmed by standard denaturing agarose gel electrophoresis to check for degradation and contamination. All samples exhibited an OD 260/280 ratio ranging from 1.8 to 2.2, indicating high purity. The integrity of the RNA samples was further validated using the Agilent Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, USA) with the RNA Nano 6000 Assay. All samples had an RNA Integrity Number (RIN) value above seven, ensuring their suitability for downstream analyses.

Thaw the frozen uterine tissues in pre-warmed saline containing antibiotics and antifungal drugs at 4 Â°C. Once thawed, trim the tissues to an appropriate size (including parts of the uterine horn and body) and fix them in 10% neutral buffered formalin for more than 48 h. Rinse the tissues in distilled water for 5 min. Dehydrate the uterus through a graded ethanol series (70%, 80%, 95%, and absolute ethanol), adjusting the time as needed based on the size of the uterus. After dehydration, clear the uterus in xylene until fully transparent, then immerse it in molten paraffin until saturated. Finally, embed the saturated tissue in a cassette. Use a rotary microtome to prepare 5-micron sections. Place the sections on a slide warmer until fully spread, then transfer the paraffin sections to microscope slides. Dry and stain them using the Hematoxylin and Eosin Staining Kit (Beyotime Co., Ltd., Beijing, China).

RNA library preparation and sequencing

The RNA library construction began with 3 Î¼g of high-quality RNA using the NEBNext® UltraTM Directional RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA), following the manufacturer's protocol closely. Ribosomal RNA (rRNA) was first removed using the Epicenter Ribo-Zeroâ„¢ Removal Kit (Epicenter, Madison, WI, USA), and the rRNA-depleted residue was purified through ethanol precipitation. The UNG enzyme was then used to degrade the second strand of the U-containing cDNA, followed by PCR amplification to generate the RNA library. The library products were purified using the AMPure XP system, and the library quality was assessed using the Agilent Bioanalyzer 2100 system. Finally, the RNA library was sequenced using a paired-end 150 bp (PE150) strategy on the Illumina Novaseq platform (Illumina, San Diego, CA, USA). All sequencing operations were conducted by Wuhan Frasergen Gene Information Co., Ltd. (Wuhan, China), ensuring professional data generation and management.

Data processing and transcriptome assembly

The imaging data from sequenced fragments generated by the high-throughput sequencer are converted into digital sequence data, known as reads, using CASAVA's base-calling algorithm. These reads are stored in fastq format, which includes both the sequence information and the corresponding quality metrics of the sequenced fragments. Prior to further analysis, the raw data undergo a critical filtration process:1) Reads containing adapter sequences were removed. 2) Reads with more than 1% of their length consisting of undetermined nucleotides (Ns) were removed, along with their paired-end counterparts. 3) Low-quality reads, where bases with a quality score of 20 or below constituted over 50% of the read length, were filtered out, along with their paired-end mates.

fter filtering, the quality of the clean reads was assessed, including calculations of quality score distributions (Q20 and Q30) and base composition, particularly GC content. The qualified clean reads were then aligned to the reference genome using HISat2 (version 2.1.0), and StringTie (version 1.3.5) was used for transcript assembly [66, 67]. HISat2 was executed with parameters tailored for strand-specific RNA-seq data, using options "-x –no-unal –un-conc", "-rna -strandness RF", and "-dta -t -p 4". StringTie was run with "-e -B", "-G ref.gtf -rf −1", and other parameters set to their default values, ensuring accurate and efficient mapping and assembly of the reads..

Identification and quantitation of uterine circRNA

Given the high rate of false positives in circRNA identification, we used a dual-software approach to rigorously screen and intersect the identification results based on the chromosomal positions of circRNAs. The find_circ software operates based on Bowtie2 alignment results, extracting a 20-nt anchor sequence from both ends of unaligned reads. These anchors are then re-examined against the reference sequence. A read is considered a candidate circRNA if the 5' end of the anchor aligns to the reference sequence with defined start and end sites (marked A3 and A4, respectively), the 3' end aligns upstream of this site (marked A1 and A2), and there is a canonical splice site (GT-AG) between A2 and A3 on the reference sequence. Only candidates with read counts of two or more are confirmed as identified circRNAs. The detailed parameters for the find_circ software are as follows: find_circ.py –genome = /path/to/ARS1.fa –prefix = chi_ –name = my_test_sample –stats = stats.txt –reads = splice_reads.fa > spliced_sites.bed.

Concurrently, CIRI2 detects PCC (paired-end mapping with a convergent signal) and PEM (pair-end mapping) signals, as well as junction reads indicative of the splicing event from the BWA alignment results. It then leverages dynamic programming alignment outcomes, along with circRNA read support numbers and genomic annotations, to refine the list of candidate circRNAs. The detailed parameters for the CIRI2 software are as follows: CIRI2.pl -T 20 -F ARS1.fa -A ARS1.gtf -I align.sam -O circRNA.xls.

For quantification, the spliced reads per billion mapped reads (SRPBM) for each circRNA were calculated, considering the number of junction reads spanning the circRNA's splice site [59]. The count of reads aligning to each transcript across samples was determined using StringTie software and subsequently normalized to ensure accurate representation and comparison of circRNA expression levels. This methodical approach ensures a more reliable and accurate identification of circRNAs, reducing the risk of false positives.

Enrichment analysis of host genes and DEcircRNA

Differential expression analysis was conducted using the DESeq2 R package [68]. Genes were considered differentially expressed if they had a p-value < 0.05 and an absolute log2 fold change > 1.5. Subsequently, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses for the parental genes associated with differentially expressed circRNAs. These analyses were performed using GOSeq (Release 2.12) and KOBAS (version 3.0) software [69]. The hypergeometric test was used to identify significantly enriched GO terms and KEGG pathways, with significance defined by adjusted Q-values < 0.05 or P-values < 0.01. Furthermore, Gene Set Enrichment Analysis (GSEA) was used to explore the differentially regulated pathways for distinct traits. The analysis was performed using the default parameters of the GSEA software, with significance determined by a p-value threshold < 0.05 and a false discovery rate (FDR) < 0.25. This comprehensive approach facilitates a deeper understanding of the biological functions and pathways influenced by differentially expressed circRNAs.

Prediction analysis of the targeting miRNA of DE circRNA

CircRNAs act as miRNA sponges, modulating gene expression by binding to miRNAs and thereby affecting their regulatory capacity. The miRanda algorithm is a valuable tool for predicting miRNA binding sites on circRNAs in animals, aiding in the exploration of these critical interactions. By analyzing circRNA-miRNA interactions, we can uncover the roles and mechanisms through which circRNAs regulate gene expression and their potential contributions to various biological processes. For the miRanda algorithm, we used the default sequence alignment threshold of 140 for the parameter -sc and the default free energy threshold of −10 for the parameter -en, ensuring that our predictions are based on established criteria.

Coding potential prediction of circRNA

To ascertain the protein-coding potential of circRNAs, we used four software tools: IRESfinder, CNCI, CPAT, and CPC2. Our approach began with IRESfinder, which predicts ribosome recruitment sites using its default parameters. We retained predictions with scores exceeding 0.5 for further analysis. Subsequently, we used CPC2 and CNCI for sequence-based predictions, following their default settings. For CPAT, we constructed predictive models using caprine coding sequences and non-coding RNA sequences, and then applied the software's default parameters for prediction.

The coding capacity of circRNAs was evaluated based on several criteria: the presence of an open reading frame (ORF) of sufficient length (three or more amino acids), ORF coverage (≥ 0.6), TESTCODE statistics (≥ 1.2), and hexamer usage bias (≥ 0.5). These metrics collectively indicate the likelihood of a circRNA encoding proteins. In conclusion, we determined the final set of circRNAs with coding potential by intersecting the predictions from all four software tools. This rigorous, multi-faceted approach ensures a high level of confidence in identifying circRNAs with potential protein-coding capabilities.

Verification of circRNA expression profiles with RT–qPCR

Several differentially expressed circRNAs were randomly selected to verify the reliability of the sequencing data and were confirmed by RT-qPCR, using RPL19 as the reference gene. For the RT–qPCR analysis, synthetic-detection templates were prepared using the PrimeScriptâ„¢ RT reagent kit (TaKaRa). RT–qPCR was performed using a Roche Light Cycler® 480 II system (Roche Applied Science, Mannheim, Germany). The primers were synthesized by Shanghai Sangon Biotech (Supplementary Table S8). RT–qPCR analysis of circRNA expression was conducted using the following procedure: 95 Â°C for 5 min, followed by 40 cycles of 95 Â°C for 5 s and 60 Â°C for 30 s. The data were analyzed using the 2−∆∆Ct method.

Statistical analysis

GraphPad Prism (version 5.0) software (San Diego, CA, United States) was used for statistical analyses of the RT-qPCR results and graphs. The statistical significance of the data was tested using paired t-tests. The results are presented as the means ± SEMs of three replicates, with statistical significance indicated by p values (*p < 0.05; **p < 0.01).

Data availability

The raw data for this study can be found in BioProject ID PRJNA731513 on NCBI.

References

  1. Zeder MA. Domestication and early agriculture in the Mediterranean Basin: Origins, diffusion, and impact. Proc Natl Acad Sci U S A. 2008;105(33):11597–604.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Zheng Z, Wang X, Li M, Li Y, Yang Z, Wang X, Pan X, Gong M, Zhang Y, Guo Y, et al. The origin of domestication genes in goats. Sci Adv. 2020;6(21):eaaz5216.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Daly KG, Delser PM, Mullin VE, Scheu A, Mattiangeli V, Teasdale MD, Hare AJ, Burger J, Verdugo MP, Collins MJ, et al. Ancient goat genomes reveal mosaic domestication in the Fertile Crescent. Science. 2018;361(6397):85–7.

    Article  CAS  PubMed  Google Scholar 

  4. Pidancier N, Jordan S, Luikart G, Taberlet P. Evolutionary history of the genus Capra (Mammalia, Artiodactyla): Discordance between mitochondrial DNA and Y-chromosome phylogenies. Mol Phylogenet Evol. 2006;40(3):739–49.

    Article  CAS  PubMed  Google Scholar 

  5. Hammer SE, Schwammer HM, Suchentrunk F. Evidence for introgressive hybridization of captive markhor (Capra falconeri) with domestic goat: cautions for reintroduction. Biochem Genet. 2008;46(3–4):216–26.

    Article  CAS  PubMed  Google Scholar 

  6. Grossen C, Keller L, Biebach I, Croll D, Consortium IGG. Introgression from domestic goat generated variation at the major histocompatibility complex of alpine ibex. PloS Genet. 2014;10(6): e1004438.

    Article  Google Scholar 

  7. Du X, Liu Y, He X, Tao L, Fang M, Chu M. Uterus proliferative period ceRNA network of Yunshang black goat reveals candidate genes on different kidding number trait. Front Endocrinol (Lausanne). 2023;14:1165409.

    Article  PubMed  Google Scholar 

  8. Tao L, He XY, Jiang YT, Lan R, Li M, Li ZM, Yang WF, Hong QH, Chu MX. Combined approaches to reveal genes associated with litter size in Yunshang black goats. Anim Genet. 2020;51(6):924–34.

    Article  CAS  PubMed  Google Scholar 

  9. Liu Y, Zhou Z, Guo S, Li K, Wang P, Fan Y, He X, Jiang Y, Lan R, Chen S, et al. Transcriptome analysis reveals key miRNA-mRNA pathways in ovarian tissues of yunshang black goats with different kidding numbers. Front Endocrinol (Lausanne). 2022;13: 883663.

    Article  PubMed  Google Scholar 

  10. Wang X, Guo X, He X, Liu Q, Di R, Hu W, Cao X, Zhang X, Zhang J, Chu M. Effects of FecB mutation on estrus, ovulation, and endocrine characteristics in small tail han sheep. Front Vet Sci. 2021;8: 709737.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Yang QY, Liu J, Wang Y, Zhao W, Wang WJ, Cui J, Yang JJ, Yue Y, Zhang S, Chu MQ, et al. A proteomic atlas of ligand-receptor interactions at the ovine maternal-fetal interface reveals the role of histone lactylation in uterine remodeling. J Biol Chem. 2022;298(1): 101456.

    Article  CAS  PubMed  Google Scholar 

  12. Dong GZ, Sun RD, Zhang R, Qin YF, Lu CC, Wang XR, Xia YK, Du GZ. Preimplantation triclosan exposure alters uterine receptivity through affecting tight junction protein(dagger). Biol Reprod. 2022;107(1):349–57.

    Article  PubMed  Google Scholar 

  13. Shekibi M, Heng S, Nie G. MicroRNAs in the regulation of endometrial receptivity for embryo implantation. Int J Mol Sci. 2022;23(11):6210.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Wang H, Dey SK. Roadmap to embryo implantation: clues from mouse models. Nat Rev Genet. 2006;7(3):185–99.

    Article  PubMed  Google Scholar 

  15. Norwitz ER, Schust DJ, Fisher SJ. Mechanisms of disease - Implantation and the survival of early pregnancy. New Engl J Med. 2001;345(19):1400–8.

    Article  CAS  PubMed  Google Scholar 

  16. Labarta E, Martinez-Conejero JA, Alama P, Horcajadas JA, Pellicer A, Simon C, Bosch E. Endometrial receptivity is affected in women with high circulating progesterone levels at the end of the follicular phase: a functional genomics analysis. Hum Reprod. 2011;26(7):1813–25.

    Article  CAS  PubMed  Google Scholar 

  17. Liu X, Zhang L, Liu Y, Cui J, Che S, An X, Song Y, Cao B. Circ-8073 regulates CEP55 by sponging miR-449a to promote caprine endometrial epithelial cells proliferation via the PI3K/AKT/mTOR pathway. Biochim Biophys Acta Mol Cell Res. 2018;1865(8):1130–47.

    Article  CAS  PubMed  Google Scholar 

  18. Zhang L, Liu X, Che S, Cui J, Liu Y, An X, Cao B, Song Y. CircRNA-9119 regulates the expression of prostaglandin-endoperoxide synthase 2 (PTGS2) by sponging miR-26a in the endometrial epithelial cells of dairy goat. Reprod Fertil Dev. 2018;30(12):1759–69.

    Article  CAS  PubMed  Google Scholar 

  19. Shen J, Chen L, Cheng J, Jin X, Mu Y, Li Q, Xia L, Gao Y, Xia Y. Circular RNA sequencing reveals the molecular mechanism of the effects of acupuncture and moxibustion on endometrial receptivity in patients undergoing infertility treatment. Mol Med Rep. 2019;20(2):1959–65.

    CAS  PubMed  Google Scholar 

  20. Ma L, Zhang M, Cao F, Han J, Han P, Wu Y, Deng R, Zhang G, An X, Zhang L, et al. Effect of MiR-100-5p on proliferation and apoptosis of goat endometrial stromal cell in vitro and embryo implantation in vivo. J Cell Mol Med. 2022;26(9):2543–56.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Liu J, Feng G, Guo C, Li Z, Liu D, Liu G, Zou X, Sun B, Guo Y, Deng M, et al. Identification of functional circRNAs regulating ovarian follicle development in goats. BMC Genomics. 2024;25(1):893.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Raza SHA, Wijayanti D, Pant SD, Abdelnour SA, Hashem NM, Amin A, Wani AK, Prakash A, Dawood MAO, Zan L. Exploring the physiological roles of circular RNAs in livestock animals. Res Vet Sci. 2022;152:726–35.

    Article  CAS  PubMed  Google Scholar 

  23. Luo J, Wang W, Sun S. Research advances in reproduction for dairy goats. Asian-Australas J Anim Sci. 2019;32(8):1284–95.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Ye J, Yao Z, Si W, Gao X, Yang C, Liu Y, Ding J, Huang W, Fang F, Zhou J. Identification and characterization of microRNAs in the pituitary of pubescent goats. Reprod Biol Endocrinol. 2018;16(1):51.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, Maier L, Mackowiak SD, Gregersen LH, Munschauer M, et al. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013;495(7441):333–8.

    Article  CAS  PubMed  Google Scholar 

  26. Patop IL, Wust S, Kadener S. Past, present, and future of circRNAs. EMBO J. 2019;38(16): e100836.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Tao H, Xiong Q, Zhang F, Zhang N, Liu Y, Suo X, Li X, Yang Q, Chen M. Circular RNA profiling reveals chi_circ_0008219 function as microRNA sponges in pre-ovulatory ovarian follicles of goats (Capra hircus). Genomics. 2017;S0888–7543(17):30129–35.

    Google Scholar 

  28. Liu X-r. Zhang L, Cui J, Yang L, Han J, Che S, Cao B, Li G, Song Y: circRNA landscape of non-pregnant endometrium during the estrus cycle in dairy goats. J Integr Agric. 2021;20:1346–58.

    Article  CAS  Google Scholar 

  29. Liu Y, Zhou Z, He X, Jiang Y, Ouyang Y, Hong Q, Chu M. Differentially Expressed Circular RNA Profile Signatures Identified in Prolificacy Trait of Yunshang Black Goat Ovary at Estrus Cycle. Front Physiol. 2022;13: 820459.

    Article  PubMed  PubMed Central  Google Scholar 

  30. Liu Y, Wang P, Zhou Z, He X, Tao L, Jiang Y, Lan R, Hong Q, Chu M. Expression profile analysis to identify circular RNA expression signatures in the prolificacy trait of yunshang black goat pituitary in the estrus cycle. Front Genet. 2021;12: 801357.

    Article  CAS  PubMed  Google Scholar 

  31. Yang J, Tang J, He X, Di R, Zhang X, Zhang J, Guo X, Chu M, Hu W. Comparative transcriptomics identify key pituitary circular RNAs that participate in sheep (Ovis aries) Reproduction. Animals. 2023;13(17):2711.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Sun Z, Hong Q, Liu Y, He X, Di R, Wang X, Ren C, Zhang Z, Chu M. Characterization of circular RNA profiles of oviduct reveal the potential mechanism in prolificacy trait of goat in the estrus cycle. Front Physiol. 2022;13: 990691.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Song Y, Zhang L, Liu X, Niu M, Cui J, Che S, Liu Y, An X, Cao B. Analyses of circRNA profiling during the development from pre-receptive to receptive phases in the goat endometrium. J Anim Sci Biotechnol. 2019;10:34.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Kim D, Salzberg SL. TopHat-Fusion: an algorithm for discovery of novel fusion transcripts. Genome Biol. 2011;12(8):R72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Gao Y, Zhang J, Zhao F. Circular RNA identification based on multiple seed matching. Brief Bioinform. 2018;19(5):803–10.

    Article  CAS  PubMed  Google Scholar 

  36. Zheng Q, Bao C, Guo W, Li S, Chen J, Chen B, Luo Y, Lyu D, Li Y, Shi G, et al. Circular RNA profiling reveals an abundant circHIPK3 that regulates cell growth by sponging multiple miRNAs. Nat Commun. 2016;7:11215.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Acquaviva J, Wong R, Charest A. The multifaceted roles of the receptor tyrosine kinase ROS in development and cancer. Biochim Biophys Acta. 2009;1795(1):37–52.

    CAS  PubMed  Google Scholar 

  38. Jun HJ, Roy J, Smith TB, Wood LB, Lane K, Woolfenden S, Punko D, Bronson RT, Haigis KM, Breton S, et al. ROS1 signaling regulates epithelial differentiation in the epididymis. Endocrinology. 2014;155(9):3661–73.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Bork P. Shuffled domains in extracellular proteins. FEBS Lett. 1991;286(1–2):47–54.

    Article  CAS  PubMed  Google Scholar 

  40. Colombatti A, Bonaldo P, Doliana R. Type A modules: interacting domains found in several non-fibrillar collagens and in other extracellular matrix proteins. Matrix. 1993;13(4):297–306.

    Article  CAS  PubMed  Google Scholar 

  41. Wen X, Hu G, Xiao X, Zhang X, Zhang Q, Guo H, Li X, Liu Q, Li H. FGF2 positively regulates osteoclastogenesis via activating the ERK-CREB pathway. Arch Biochem Biophys. 2022;727: 109348.

    Article  CAS  PubMed  Google Scholar 

  42. Wiles JR, Katchko RA, Benavides EA, O’Gorman CW, Escudero JM, Keisler DH, Stanko RL, Garcia MR. The effect of leptin on luteal angiogenic factors during the luteal phase of the estrous cycle in goats. Anim Reprod Sci. 2014;148(3–4):121–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Welter H, Wollenhaupt K, Einspanier R. Developmental and hormonal regulated gene expression of fibroblast growth factor 2 (FGF-2) and its receptors in porcine endometrium. J Steroid Biochem Mol Biol. 2004;88(3):295–304.

    Article  CAS  PubMed  Google Scholar 

  44. Presta M, Dell’Era P, Mitola S, Moroni E, Ronca R, Rusnati M. Fibroblast growth factor/fibroblast growth factor receptor system in angiogenesis. Cytokine Growth Factor Rev. 2005;16(2):159–78.

    Article  CAS  PubMed  Google Scholar 

  45. Nakagawa T, Suzuki-Nakagawa C, Watanabe A, Asami E, Matsumoto M, Nakano M, Ebihara A, Uddin MN, Suzuki F. Site-1 protease is required for the generation of soluble (pro)renin receptor. J Biochem. 2017;161(4):369–79.

    Article  CAS  PubMed  Google Scholar 

  46. Seidah NG, Prat A. The biology and therapeutic targeting of the proprotein convertases. Nat Rev Drug Discov. 2012;11(5):367–83.

    Article  CAS  PubMed  Google Scholar 

  47. Uhlen M, Fagerberg L, Hallstrom BM, Lindskog C, Oksvold P, Mardinoglu A, Sivertsson A, Kampf C, Sjostedt E, Asplund A, et al. Proteomics. Tissue-based map of the human proteome. Science. 2015;347(6220):1260419.

    Article  PubMed  Google Scholar 

  48. Owczarek CM, Portbury KJ, Kola I, Hertzog PJ. Assignment of protein kinase C mu (PRKCM) to human chromosome band 14q11 with somatic cell hybrids and radiation hybrids. Cytogenet Cell Genet. 2000;89(3–4):240–1.

    Article  CAS  PubMed  Google Scholar 

  49. Johannes FJ, Prestle J, Eis S, Oberhagemann P, Pfizenmaier K. PKCu is a novel, atypical member of the protein kinase C family. J Biol Chem. 1994;269(8):6140–8.

    Article  CAS  PubMed  Google Scholar 

  50. Sell-Kubiak E, Dobrzanski J, Derks MFL, Lopes MS, Szwaczkowski T. Meta-analysis of SNPs determining litter traits in pigs. Genes (Basel). 2022;13(10):1730.

    Article  CAS  PubMed  Google Scholar 

  51. Nonneman DJ, Schneider JF, Lents CA, Wiedmann RT, Vallet JL, Rohrer GA. Genome-wide association and identification of candidate genes for age at puberty in swine. BMC Genet. 2016;17:50.

    Article  PubMed  PubMed Central  Google Scholar 

  52. Tebar F, Bohlander SK, Sorkin A. Clathrin assembly lymphoid myeloid leukemia (CALM) protein: localization in endocytic-coated pits, interactions with clathrin, and the impact of overexpression on clathrin-mediated traffic. Mol Biol Cell. 1999;10(8):2687–702.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Kumar A, Dumasia K, Deshpande S, Balasinor NH. Direct regulation of genes involved in sperm release by estrogen and androgen through their receptors and coregulators. J Steroid Biochem Mol Biol. 2017;171:66–74.

    Article  CAS  PubMed  Google Scholar 

  54. Mizuno-Yamasaki E, Rivera-Molina F, Novick P. GTPase networks in membrane traffic. Annu Rev Biochem. 2012;81:637–59.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Shan SO. ATPase and GTPase tangos drive intracellular protein transport. Trends Biochem Sci. 2016;41(12):1050–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. van Kruining D, Luo Q, van Echten-Deckert G, Mielke MM, Bowman A, Ellis S, Oliveira TG, Martinez-Martinez P. Sphingolipids as prognostic biomarkers of neurodegeneration, neuroinflammation, and psychiatric diseases and their emerging role in lipidomic investigation methods. Adv Drug Deliv Rev. 2020;159:232–44.

    Article  PubMed  PubMed Central  Google Scholar 

  57. Laychock SG, Tian Y, Sessanna SM. Endothelial differentiation gene receptors in pancreatic islets and INS-1 cells. Diabetes. 2003;52(8):1986–93.

    Article  CAS  PubMed  Google Scholar 

  58. Kono Y, Nishiuma T, Nishimura Y, Kotani Y, Okada T, Nakamura S, Yokoyama M. Sphingosine kinase 1 regulates differentiation of human and mouse lung fibroblasts mediated by TGF-beta1. Am J Respir Cell Mol Biol. 2007;37(4):395–404.

    Article  CAS  PubMed  Google Scholar 

  59. Yang C, Li Q, Chen X, Zhang Z, Mou Z, Ye F, Jin S, Jun X, Tang F, Jiang H. Circular RNA circRGNEF promotes bladder cancer progression via miR-548/KIF2C axis regulation. Aging (Albany NY). 2020;12(8):6865–79.

    Article  CAS  PubMed  Google Scholar 

  60. Zhu S, He C, Deng S, Li X, Cui S, Zeng Z, Liu M, Zhao S, Chen J, Jin Y, et al. MiR-548an, transcriptionally downregulated by hif1alpha/hdac1, suppresses tumorigenesis of pancreatic cancer by targeting vimentin expression. Mol Cancer Ther. 2016;15(9):2209–19.

    Article  CAS  PubMed  Google Scholar 

  61. Yang G, Xu Q, Wan Y, Zhang L, Wang L, Meng F. Circ-CSPP1 knockdown suppresses hepatocellular carcinoma progression through miR-493-5p releasing-mediated HMGB1 downregulation. Cell Signal. 2021;86: 110065.

    Article  CAS  PubMed  Google Scholar 

  62. Wang X, Shi J, Niu Z, Wang J, Zhang W. MiR-216a-3p regulates the proliferation, apoptosis, migration, and invasion of lung cancer cells via targeting COPB2. Biosci Biotechnol Biochem. 2020;84(10):2014–27.

    Article  CAS  PubMed  Google Scholar 

  63. Zhang D, Yang X, Luo Q, Fu D, Li H, Zhang P, Tie C. Circular RNA CSPP1 motivates renal cell carcinoma carcinogenesis and the Warburg effect by targeting RAC1 through microRNA-493-5p. Acta Biochim Pol. 2023;70(3):693–701.

    CAS  PubMed  Google Scholar 

  64. Pamudurti NR, Bartok O, Jens M, Ashwal-Fluss R, Stottmeister C, Ruhe L, Hanan M, Wyler E, Perez-Hernandez D, Ramberger E, et al. Translation of CircRNAs. Mol Cell. 2017;66(1):9-21e27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Liang C, Han MC, Zhou ZY, Liu YF, He XY, Jiang YT, Ouyang YN, Hong QH, Chu MX. Hypothalamic transcriptome analysis reveals the crucial microRNAs and mRNAs affecting litter size in goats. Front Vet Sci. 2021;8:747100.

    Article  PubMed  PubMed Central  Google Scholar 

  66. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290-+.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT. StringTie and Ballgown Nat Protoc. 2016;11(9):1650–67.

    Article  CAS  PubMed  Google Scholar 

  68. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  Google Scholar 

  69. Xie C, Mao XZ, Huang JJ, Ding Y, Wu JM, Dong S, Kong L, Gao G, Li CY, Wei LP. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 2011;39:W316–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This work was financially supported by the National Natural Science Foundation of China (32102509), the Agricultural Science and Technology Innovation Program of China (CAAS-ZDRW202106 and ASTIP-IAS13), the China Agriculture Research System of MOF and MARA (CARS-38).

Author information

Authors and Affiliations

  1. State Key Laboratory of Animal Biotech Breeding, Institute of Animal Science, Chinese Academy of Agricultural Sciences (CAAS), Beijing, 100193, China

    Xiaolong Du, Yufang Liu, Xiaoyun He, Lin Tao & Mingxing Chu

  2. Department of Animal Genetics and Breeding, Laboratory of Animal Genetics and Breeding, College of Animal Science and Technology, National Engineering Laboratory for Animal Breeding, MARA, China Agricultural University, Beijing, 100193, China

    Xiaolong Du & Meiying Fang

Authors
  1. Xiaolong Du
    View author publications

    You can also search for this author inPubMed Google Scholar

  2. Yufang Liu
    View author publications

    You can also search for this author inPubMed Google Scholar

  3. Xiaoyun He
    View author publications

    You can also search for this author inPubMed Google Scholar

  4. Lin Tao
    View author publications

    You can also search for this author inPubMed Google Scholar

  5. Meiying Fang
    View author publications

    You can also search for this author inPubMed Google Scholar

  6. Mingxing Chu
    View author publications

    You can also search for this author inPubMed Google Scholar

Contributions

Conceptualization, MYF, and MXC; Methodology, XLD, YFL, MYF, and MXC; Validation: XLD, YFL, and XYH; Formal Analysis: XLD and MXC; Investigation: XLD, YFL, LT, XYH, MYF, and MXC; Resources, XLD, YFL, LT, and MXC; Data Curation, XLD, YFL, LT, MYF, and MXC; Writing-Original Draft Preparation, XLD, and MXC; Supervision, MXC; Project Administration, MXC; Funding Acquisition, MXC. All authors have read and agreed to the published version of the manuscript.

Corresponding authors

Correspondence to Meiying Fang or Mingxing Chu.

Ethics declarations

Ethics approval and consent to participate

All procedures performed in studies involving animals were in accordance with the ethical standards of the institutional and national research committees. Ethical approval was provided by the animal ethics committee of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (No. IAS2019-63). The manuscript adheres to the ARRIVE guidelines for the reporting of animal experiments, and the study was carried out in compliance with the ARRIVE guidelines.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Du, X., Liu, Y., He, X. et al. Identification and expression profile analysis of circRNAs associated with goat uterus with different fecundity during estrous cycle. BMC Genomics 26, 349 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-025-11489-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-025-11489-x

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us