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

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

Introduction

The kidding number of goats, an essential economic indicator of the goat industry, is mainly dependent on ovulation rate and ovarian growth (Notter, 2012). Many factors influence the kidding number and ovulation in goats. The reproductive performance of goats is different to that of sheep, and the genes conferring the major effects on prolificacy of goats are still unknown. Prolific and non-prolific breeds are characterized by differences in the estrus cycle. Determining the molecular mechanisms of follicular development has therefore been a focus of research in the genetic breeding of goats. In recent years, many studies of the mechanisms of prolificacy in goats have focused on candidate genes and their regulators, such as miRNAs, transcription factors, and varieties of regulated pathways. All these factors form a network of many interactions in vivo that regulate the goat reproduction trait (Lu et al., 2020; Zhao et al., 2020). The genes GDF9, Smad2, BMP15, and BMP4, and the transcription factor Sox9 are involved in kidding number and increased ovulation rate, and are proven to be crucial for the reproductive performance of goats (Hanrahan et al., 2004; Lai et al., 2016; Chinmoy et al., 2017; He et al., 2020; Wang et al., 2021). In addition, many factors affected the reproductive ability of goats including follicle recruitment, development, selection, and ovulation (Luo et al., 2019). A few signaling pathways, such as the Wnt/β-catenin, nodal mammalian target of rapamycin (mTOR), and bone morphogenetic protein (BMP)/Smad pathways, are involved in follicular development (Li et al., 2016).

CircRNAs are one type of long noncoding RNAs produced by pre-mRNA back-splicing that have been progressively identified and validated to form a covalently closed loop structure of single-strand, non-polyadenylated circular molecules (Memczak et al., 2013; Ashwal-Fluss et al., 2014), which show a previously neglected regulatory method in the form of the circRNA-miRNA-mRNA axes (Hansen et al., 2013). Usually, circular RNAs act as a sponge for related miRNAs, and compared with their corresponding colinear mRNA isoforms, circRNAs are more stable (Chuang et al., 2018). Recently, many studies on circRNA function have shown that regulation of circRNA expression can affect a variety of molecular and physiological phenotypes, mainly related to the nervous system, innate immunity, and diseases (Xiao et al., 2020). CircRNAs also act as a template for translation (Kristensen et al., 2019). Many circRNAs were shown to play important roles in the process of reproduction by regulating the expression of genes at the transcriptional and post-transcriptional levels. In the study of reproduction trait, circRalGPS2s were lowly expressed in follicle growth and predominant in the ovarian stroma, suggesting that circularization of circRNAs may not have a relationship with the development of granulosa cells (Shen et al., 2019). A few studies showed that the expression of circRNA_103827 and circRNA_104816 was negatively correlated with embryonic development (Cheng et al., 2017). Although circRNAs have various functions and act through mechanisms, the regulatory effect of circRNAs on the goat reproduction trait is uncharacterized. Hence, it is crucial to reveal the molecular mechanism of circRNAs regulating the prolificacy trait and the transition from the follicular phase to the luteal phase in goats.

The Yunshang black goat is an excellent breed of meat-producing goat, which has many favorable characteristics including tender meat and unique flavor, and live in Yunnan Province, China. Yunshang black goat is suitable for large-scale breeding in all of regions of China. Many studies have examined the molecular mechanisms of resource conservation, strain selection, and trait formation of this native variety. One excellent feature of Yunshang black goat is prolificacy, making it suitable for the study of prolific goats. The molecular mechanisms underlying prolificacy trait regulation in the Yunshang black goat have not yet been confirmed. In order to reveal the key molecular mechanism networks of the prolificacy trait and follicular development in this variety, in this study, we screened for differentially expressed (DE) circRNAs among high and low fecundity goats in two phases of the estrus cycle (four groups: high fecundity in the follicular phase, high fecundity in the luteal phase, low fecundity in the follicular phase, and low fecundity in the luteal phase, respectively) by RNA-seq. The crucial correlated miRNA-circRNA and circRNA-mRNA regulator networks and pathways in the prolificacy trait and the transition from the follicular phase to the luteal phase were identified. Furthermore, the function and mechanism of DE circRNAs were predicted by functional enrichment and coding potential analysis. These analyses will provide a basis for the molecular mechanism of follicular development and the transition from the follicular phase to the luteal phase in goats. Our purpose was to clarify the molecular mechanisms of follicular phase regulation of goat ovaries and to provide a research basis for further study on the interactions between circRNAs, miRNAs, host genes, and signaling pathways for prolific traits and the transition from the follicular phase to the luteal phase in goats. This study also provides basic data for Yunshang black goat prolificacy trait-related candidate gene screening and the conservation of native breed resources.

Materials and Methods

Sample Collection and Preparation

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

Preparation of RNA Library for Sequencing

According to the manufacturer’s protocol, total RNA was isolated from ground ovarian tissue powder using TRIzol reagent (Invitrogen, Carlsbad, CA, United States ). The degree of degradation and contamination of RNA samples was monitored by 1% agarose gels, and the NanoPhotometer^® spectrophotometer (IMPLEN, CA, United States) was used to check the purity of the RNA samples. The Qubit^® 2.0 Fluorometer (Life Technologies, CA, United States ) with the Qubit^® RNA Assay Kit was used to measure the concentration of RNA. The integrity of RNA was assessed by the RNA Nano 6000 Assay Kit of the Bioanalyzer 2,100 system (Agilent Technologies, CA, United States ).

The 3 μg of RNA per sample was used as the starting material for the RNA samples. Initially, the Epicentre Ribo-zero™ rRNA Removal Kit (Epicentre, United States ) was used to remove rRNA, and ethanol precipitation was used to purify the rRNA free samples. Secondly, the NEBNext^® Ultra™ Directional RNA Library Prep Kit for Illumina^® (NEB, United States ) was used to deplete RNA from rRNA and generate the sequencing libraries in accordance with the manufacturer. In brief, the divalent cations in NEB Next First Strand Synthesis Reaction Buffer (5X) under elevated temperature were used to perform fragmentation. The random hexamer primer and M-MuLV Reverse transcriptase (RNaseH-) were used to synthesize the first-strand cDNA. Subsequently, DNA polymerase I and rnase H were used to synthesize second-strand cDNA. The dUTP replaced the dNTPs with dTTP in the reaction buffer. The active exonuclease/polymerase converted the remaining overhangs into blunt ends. The 3′-ends of DNA fragment adenylation and the NEB Next Adaptor with hairpin loop structure prepared for hybridization were ligated. The AMPure XP system (Beckman Coulter, Beverly, United States ) was used to purify the library fragments and to select the cDNA fragments with different lengths. PCR was then performed by adding 3 μl of USER Enzyme (NEB, United States ) to size-selected, adaptor-ligated cDNA for 15 min at 37°C, followed by 5 min at 95°C. PCR was subsequently performed with Phusion High-Fidelity DNA polymerase, Index (X), and Primer universal PCR primers. Finally, the Agilent Bioanalyzer 2,100 system was used to assesses library quality, and the AMPure XP system was used to purify the products.

Identification and Analysis of Differentially Expressed circRNAs

CircRNAs were identified by find_circ and CIRI2 (Memczak et al., 2013; Gao et al., 2017). The following are the basic principles of find_circ: according to Bowtie2 alignment results, find_circ extracted the 20 nt anchor sequences from the ends of the reads that were not aligned to the reference sequence, and aligned each pair of anchor sequences with the reference sequence again. If the 5′-end of the anchor sequence aligned with the reference sequence (start and end sites A3 and A4, respectively), the 3′-end of the anchor sequence was aligned upstream of the site (start and end sites A1 and A2, respectively), and there was a splice over the site (GT-AG) between the A2 and A3 of the reference sequence, then this read could be considered a candidate circRNA. Finally, the candidate circRNAs with read counts greater than or equal to 2 were identified as circRNAs. CIRI2 searched for a PCC (paired end mapping) signal and PEM (pair-end mapping) signal and junction reads of the shear site of paired chiastic clipping according to the comparison of the results of BWA (Li and Durbin, 2009). Then, based on the results of dynamic programming alignment, circRNA read support numbers and genome annotation information was used to filter the candidate circRNAs.

Quantification was performed based on reads that mapped across the circularized junctions and were normalized using SRPBM (number of circular reads/number of mapped reads (units in billion)/read length) (Zheng et al., 2016). Differential circRNA expression was analyzed using DESeq2. CircRNAs with padj values of < 0.05 were considered to be significantly differentially expressed.

The GO and KEGG Enrichment of Host Genes With DE circRNAs

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

Prediction Analysis of the Targeting miRNA of DE circRNAs

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

Coding Potential Prediction of circRNAs

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

Verification of Sequencing Results by RT-qPCR

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

TABLE 1

TABLE 1. The primer information of DE circRNAs.

PPI Analysis of Host Genes

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

Results

Analysis of Transcriptome Sequencing of Ovaries in Yunshang Black Goats

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

TABLE 2

TABLE 2. Summary of data quality and genome alignment.

FIGURE 1

FIGURE 1. The identification of circRNAs in goat ovarian tissues. (A): Length distribution of circRNAs. (B): Function region statistics of circRNAs.

FIGURE 2

FIGURE 2. Density statistics of circRNAs in chromosomes in the four groups. HL: high fecundity goats in the luteal phase; LL: low fecundity goats in the luteal phase; HF: high fecundity goats in the follicular phase; LF: low fecundity goats in the follicular phase.

Analysis of the Expression of circRNAs

The expression of circRNAs per individual was calculated, and their levels were normalized to SRPBM (Figure 3A). Based on the normalized expression with |log2 (fold change) | > 1 and p-value < 0.05, 283 DE circRNAs (114 upregulated and 169 downregulated), 559 DE circRNAs (299 upregulated and 260 downregulated), 449 DE circRNAs (254 upregulated and 195 downregulated), and 314 DE circRNAs (210 upregulated and 104 downregulated) were identified in the HF vs LF, HF vs HL, HL vs LL, and LF vs LL comparisons, respectively (Figure 3B, Supplementary Table S1). The result showed that the host genes of circPRLR, circTGFBR2, and circFSHR were PRLR, TGFBR2, and FSHR, respectively. PRLR and FSHR are important genes in ovulation and these circRNAs probably related to follicular growth and development, or even to the goat prolificacy trait. To analyze the clustering model of the four groups with DE circRNAs, K-means and SOM clustering analyses were used to clustered the 1605 DE circRNAs (Figure 3C). For further validating the reliability of circRNA results, nine DE circRNAs were randomly selected for RT-qPCR across the four groups. The RT-qPCR results accorded with the RNA-seq data (Figure 4).

FIGURE 3

FIGURE 3. DE circRNA screening and clustering analysis. (A): The distribution of circRNA expression levels in each sample. The expression of circRNAs was normalized by SRPBM. (B): The volcano plot of DE circRNAs in the four comparisons. The X-axis represents the change in circRNA expression in different comparisons. The Y-axis represents the statistical significance of circRNAs. The blue dots represent the circRNAs that were not significantly different, the red dots represent circRNAs that were significantly upregulated, and the green dots represent circRNAs that were significantly downregulated. (C): The clustering diagram of DE circRNAs. The expression was gradually upregulated from blue to red.

FIGURE 4

FIGURE 4. Comparison of the results of DE circRNAs in RNA-seq and RT-qPCR analysis.

Enrichment Analysis of the Host Genes for DE circRNAs

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

FIGURE 5

FIGURE 5. GO and KEGG enrichment analysis of the host genes with DE circRNAs. (A): The results of GO enrichment. (B): The KEGG pathway enrichment.

Comparison of the host genes and KEGG pathway analysis showed that the host genes were annotated into 154, 175, 157, and 183 pathways for the HF vs LF, HL vs LL, LF vs LL, and HF vs HL comparisons, respectively (Figure 5B). Of these, progesterone-mediated oocyte maturation (ko04914), ovarian steroidogenesis (ko04913), oocyte meiosis (ko04114), the GnRH signaling pathway (ko04912), the estrogen signaling pathway (ko04915), and other pathways were involved in the development of ovarian and oocyte maturation.

Analysis of Targeting Relationship Between circRNAs and miRNAs

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

FIGURE 6

FIGURE 6. Predicted target relationships of circRNAs and miRNAs. Red triangles represent miRNAs. Green circles represent circRNAs.

Identifying the Coding Potential of circRNAs

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

TABLE 3

TABLE 3. The predicted results for the coding potential of circRNAs.

PPI Network for Goat Prolificacy Trait

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

FIGURE 7

FIGURE 7. The protein-protein interaction (PPI) network of host proteins with DE circRNAs in the comparisons.

Discussion

The process of reproduction is very complex. In particular, the study of the goat reproductive molecular mechanism is still in its infancy and there are still many unanswered questions. The reproduction traits in goats are governed by the precise regulation of genes, regulatory elements, and other factors, such as spatiotemporal expression, transcriptional regulation, and feedback mechanisms (Ye et al., 2018; Luo et al., 2019). Investigation of the ovary from goats in different estrus stages should provide information on potential follicle development and allow further improvement to kidding numbers. In a previous molecular biology study, hormone secretion and cell signaling transduction were found to play an important role in the mechanisms of ovarian and follicular development (Messina et al., 2016; Zou et al., 2020). Recently, circRNAs were reported to have a range of extremely important effects on many kinds of biological processes (Memczak et al., 2013). A lot of studies have showed that circRNAs are widespread in plants and mammals, however, little is known about circRNAs in goats, especially those related to the prolificacy trait and ovarian development during the different phases of the estrus cycle (Rybak-Wolf et al., 2015; Ye et al., 2015).

Here, the high-throughput sequencing of circRNAs in the ovaries from goats with different kidding number records was performed in different estrus phases. The circRNAs were distributed widely among chromosomes 1–11 in the four comparisons, and which showed that circRNAs had complexity and functional diversity (Figure 2, Additional file Supplementary Figure S1). The majority of circRNAs in this study were <10,000 nt in length, which was longer than those in previous studies in chicken granulosa cells (200–300 nt) (Shen et al., 2019). In total, 53.8% of the circRNAs were shorter than 10,000 nt and longer than 1,000 nt, the proportion was also higher than that in the studies of goat endometrium circRNAs (23% > 1,000 nt) (Song et al., 2019). Species and tissue differences may be the main reason for the variation in circRNA length. In our study, the main source of circRNAs in the ovary was exon splicing, which was supported by previous research (Lv et al., 2020). Given all of the above features, it was concluded that there were species- and tissue-specific differences in circRNA expression. The changes in circRNA abundance in the four comparisons might be associated with their specific roles during the different estrus phases of the ovary and kidding numbers. We identified 283, 559, 449, and 314 DE circRNAs in HF vs LF, HF vs HL, HL vs LL, and LF vs LL comparisons, respectively (Figure 3B, Supplementary Table S1). Then, these host genes of DE circRNAs were analyzed for GO terms and KEGG pathways. The majority of enriched GO terms were related to cell proliferation and reproductive process. There were 14 host genes of circRNAs related to the reproductive process; notably, mitogen-activated protein kinase 1 (MAPK1), which could mediate luteinizing hormone-induced communication disruption and oocyte maturation in ovarian follicles of rat (Gratao et al., 2008). MAPK1 was shown to play a basic role in regulating spindle assembly, microtubule organization, chromosome distribution, and meiosis resumption, and it could be hypothesized that activation of MAPK1 might improve meiosis resumption (Gordo et al., 2001). MAPK1 mRNA plays a crucial role in oocyte maturation in various species, including cattle and dogs, through its action on granulosa and cumulus cells (Quetglas et al., 2010; Lee et al., 2017). Moreover, some other reproduction-related pathways were also significantly enriched, including oocyte meiosis, the GnRH signaling pathway, and the estrogen signaling pathway.

In the HF vs LF comparison, circIGF2R, circMAPK8, circMYH10, and circPTEN were associated with the prolificacy trait. The IGF2R locus contains a highly abundant set of circIGF2R circRNAs, that showed multiple forms of different exon combinations in mice and plays significant roles in mammalian reproduction (Perera et al., 2018). MAPK8 might regulate the transcription of basal FSHb, and one study has shown that blocking its activation could inhibit the levels of major basal FSHb transcripts by half (Haisenleder et al., 2008). In the HL vs LL comparison, the host genes of circFSHR, circHMGCR, circSMAD3, and circFGFR2 were FSHR, HMGCR, SMAD3, and FGFR2, respectively. FSHR stimulates Gαs protein-dependent activation of cAMP/PKA, leading to phosphorylation of cAMP response element binding protein and steroidogenic protein, which is required for the production of estrogens, which in turn is known to be a growth stimulant (Casarini and Crépieux, 2019). In a study of common genetic variants affecting dizygotic (DZ) twins, a GWAS conducted in 1980 mothers of spontaneous DZ twins and 12,953 controls found that SNPs similar to SMAD family member 3 (SMAD3) were significant associated with DZ twinning (Mbarek et al., 2016). Collectively, these results suggested that those circRNAs were involved in ovarian and oocyte development and thus influence the number of kidding.

Among the DE circRNAs of the comparison between the estrus phases (HF vs HL and LF vs LL), we found that the host genes of circLDLR, circMAPK1, circSMAD1, and circSMAD2 were LDLR, MAPK1, SMAD1, and SMAD2, respectively. SMAD1 and SMAD2 are members of the SMAD family, which is related to signal transduction in the TGF-β superfamily and the inhibition of BMP activity, and they have been thought to control follicular selection by balancing luteinizing hormone receptor (LHR) transcription (Yu et al., 2017; Xu et al., 2018). In addition, LDLR has been shown to take part in the biosynthesis of steroid hormones by affecting the absorption of cholesterol substrates via circulating lipoproteins (Ying et al., 2013). The MAPK1 (p42ERK2) isoform of MAPK is found in mammalian oocytes and the kinase is shown to be activated during meiotic maturation (Fan and Sun, 2004). Together, these circRNAs might be associated with the difference in the estrus phases. In previous studies, similar results were found for non-coding RNAs, including long non-coding RNA and microRNA enrichment analysis, suggesting that these circRNAs might act as sponges of endogenous RNAs (ceRNAs) to regulate gene expression (Kutzler et al., 2001). The results obtained from this study and previous reports suggested that the prolificacy trait and estrus induction are complex processes that includes numerous events, and that the circRNAs distributed in the two estrus phases in high and low fecundity goats are mainly involved in follicle growth, reproduction events, and other biological processes.

Prediction of circRNA-miRNA binding relationships was completed by the miRNA target gene prediction method. In total, binding was predicted for 56 miRNAs, including the miR-133 family (such as miR-133a-3p and miR-133b), miR-129-3p, and miR-21. Several of these miRNAs have been shown to be closely related to ovarian development and the prolificacy trait. Differential expression of miR-133 occurs during meiotic maturation of oocytes (Song et al., 2014). MiR-129-3p serves as a tumor suppressor by targeting BZW1 in ovarian cancer cells and the restoration of miR-129-3p might be a novel therapeutic strategy for ovarian cancer (Liu et al., 2018). MiR-21 promotes follicle development in the ovary during ovulation in mice (Carletti et al., 2010). In the chicken reproductive process, circRalGPS2 can compete with miRNA to modulate mRNA expression that regulates downstream target genes (Chiu et al., 2017). In our study, miR-129-3p was predicted to target chi-circC2H2orf76, chi-circDUS2, chi-circGPR137B, chi-circMYSM1, and chi-circSVEP1, miR-21-3p was predicted to targeting chi-circUTP20, and miR-133a-3p and miR-133b were predicted to targets to chi-circAKAP10. Overall, these differentially expressed circRNAs may regulate the function of the ovary, for estrus and different kidding numbers, by regulating miRNA, although specific mechanisms require further study.

As circRNAs have a loop-like structure consisting of exons, some circRNAs consist of an open reading frame and IRES. Therefore, researchers hypothesized that numerous circRNAs could be translated into protein in certain conditions (Nagarjuna et al., 2017). IRESfinder, CPC2, CNCI, CPAT, and PLEK software identified 18 circRNAs with coding potential. It should be noted that the biological function of these circRNAs, from the perspective of protein translation, should be explored further. In addition, we demonstrated the interactions between host genes in the ovary of Yunshang black goats. A total of 69 PPI pairs for the comparisons and several host genes were associated with the prolificacy trait. The relationships between these interactions require further study.

Conclusion

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

Data Availability Statement

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

Ethics Statement

The animal study was reviewed and approved by All experimental procedures mentioned in the present study were approved by the Science Research Department (in charge of animal welfare issue) of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (IAS-CAAS) (Beijing, China). Ethical approval was provided by the animal ethics committee of IASCAAS (No. IAS 2019–63).

Author Contributions

Conceptualization: YL and MC; methodology: YL, ZZ, QH, and MC; validation: YL and ZZ; formal analysis: YL and MC; investigation: YL, ZZ, and MC; resources: YL, ZZ, XH, YO, and MC; data curation: YL, ZZ, and MC; writing-original draft preparation: YL and MC; supervision: MC; project administration: YL, QH, and MC; funding acquisition: MC and QH. All authors have read and agreed to the published version of the manuscript.

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), the Basic Research Foundation Key Project of Yunnan Province (202001AS070002), and the Major Science and Technology Project of Yunnan Province (202102AE090039).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2022.820459/full#supplementary-material

References

Ashwal-Fluss R., Meyer M., Pamudurti N. R., Ivanov A., Bartok O., Hanan M., et al. (2014). circRNA Biogenesis Competes with Pre-mRNA Splicing. Mol. Cel. 56, 55–66. doi:10.1016/j.molcel.2014.08.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Carletti M. Z., Fiedler S. D., Christenson L. K. (2010). MicroRNA 21 Blocks Apoptosis in Mouse Periovulatory Granulosa Cells1. Biol. Reprod. 83 (2), 286–295. doi:10.1095/biolreprod.109.081448

PubMed Abstract | CrossRef Full Text | Google Scholar

Casarini L., Crépieux P. (2019). Molecular Mechanisms of Action of FSH. Front. Endocrinol. 10, 305. doi:10.3389/fendo.2019.00305

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng J., Huang J., Yuan S., Zhou S., Yan W., Shen W., et al. (2017). Circular RNA Expression Profiling of Human Granulosa Cells during Maternal Aging Reveals Novel Transcripts Associated with Assisted Reproductive Technology Outcomes. PLoS One 12 (6), e0177888. doi:10.1371/journal.pone.0177888

PubMed Abstract | CrossRef Full Text | Google Scholar

Chinmoy M., Mangalika R., Mishra S., Sahoo S., Gangadhar N., Patra R. C. (2017). Genetic Polymorphism of Prolific Genes in Goat-A Brief Review. Explor Anim. Med. Res. 7, 132–141.

Google Scholar

Chiu H.-S., Martínez M. R., Bansal M., Subramanian A., Golub T. R., Yang X., et al. (2017). High-throughput Validation of ceRNA Regulatory Networks. BMC Genomics 18 (1), 418. doi:10.1186/s12864-017-3790-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Chuang T.-J., Chen Y.-J., Chen C.-Y., Mai T.-L., Wang Y.-D., Yeh C.-S., et al. (2018). Integrative Transcriptome Sequencing Reveals Extensive Alternative Trans-splicing and Cis-Backsplicing in Human Cells. Nucleic Acids Res. 46, 3671–3691. doi:10.1093/nar/gky032

PubMed Abstract | CrossRef Full Text | Google Scholar

Fahlgren N., Carrington J. C. (2010). miRNA Target Prediction in Plants. Methods Mol. Biol. 592, 51–57. doi:10.1007/978-1-60327-005-2_4

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan H.-Y., Sun Q.-Y. (2004). Involvement of Mitogen-Activated Protein Kinase Cascade during Oocyte Maturation and Fertilization in Mammals1. Biol. Reprod. 70 (3), 535–547. doi:10.1095/biolreprod.103.022830

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao Y., Zhang J. Y., Zhao F. Q. (2017). Circular RNA Identification Based on Multiple Seed Matching. Brief. Bioinform. 1, 8. doi:10.1093/bib/bbx014

CrossRef Full Text | Google Scholar

Gordo A. C., He C. L., Smith S., Fissore R. A. (2001). Mitogen Activated Protein Kinase Plays a Significant Role in Metaphase II Arrest, Spindle Morphology, and Maintenance of Maturation Promoting Factor Activity in Bovine Oocytes. Mol. Reprod. Dev. 59 (1), 106–114. doi:10.1002/mrd.1012

PubMed Abstract | CrossRef Full Text | Google Scholar

Gratao A. A., Dahlhoff M., Sinowatz F., Wolf E., Schneider M. R. (2008). Betacellulin Overexpression in the Mouse Ovary Leads to MAPK3/MAPK1 Hyperactivation and Reduces Litter Size by Impairing Fertilization1. Biol. Reprod. 78 (1), 43–52. doi:10.1095/biolreprod.107.062588

PubMed Abstract | CrossRef Full Text | Google Scholar

Haisenleder D. J., Burger L. L., Walsh H. E., Stevens J., Aylor K. W., Shupnik M. A., et al. (2008). Pulsatile Gonadotropin-Releasing Hormone Stimulation of Gonadotropin Subunit Transcription in Rat Pituitaries: Evidence for the Involvement of Jun N-Terminal Kinase but Not P38. Endocrinology 149 (1), 139–145. doi:10.1210/en.2007-1113

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanrahan J. P., Gregan S. M., Mulsant P., Mullen M., Davis G. H., Powell R., et al. (2004). Mutations in the Genes for Oocyte-Derived Growth Factors GDF9 and BMP15 Are Associated with Both Increased Ovulation Rate and Sterility in Cambridge and Belclare Sheep (Ovis aries). Biol. Reprod. 70 (4), 900–909. doi:10.1095/biolreprod.103.023093

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansen T. B., Jensen T. I., Clausen B. H., Bramsen J. B., Finsen B., Damgaard C. K., et al. (2013). Natural RNA Circles Function as Efficient microRNA Sponges. Nature 495, 384–388. doi:10.1038/nature11993

PubMed Abstract | CrossRef Full Text | Google Scholar

He L., Bi Y., Wang R., Pan C., Chen H., Lan X., et al. (2020). Detection of a 4 Bp Mutation in the 3′UTR Region of Goat Sox9 Gene and its Effect on the Growth Traits. Animals 10 (4), 672. doi:10.3390/ani10040672

PubMed Abstract | CrossRef Full Text | Google Scholar

Kristensen L. S., Andersen M. S., Stagsted L. V. W., Ebbesen K. K., Hansen T. B., Kjems J. (2019). The Biogenesis, Biology and Characterization of Circular RNAs. Nat. Rev. Genet. 20, 675–691. doi:10.1038/s41576-019-0158-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Kutzler M. A., Wheeler R., Volkmann D. H. (2001). “Canine Oestrus Induction Using the GnRH Agonist, Deslorelin,” in Proc. Ann. Symp. Eur. Vet. Soc. Small Anim. Reprod., 147–148.

Google Scholar

Lai F.-N., Zhai H.-L., Cheng M., Ma J.-Y., Cheng S.-F., Ge W., et al. (2016). Whole-genome Scanning for the Litter Size Trait Associated Genes and SNPs under Selection in Dairy Goat (Capra hircus). Sci. Rep. 6 (1), 38096. doi:10.1038/srep38096

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee S. H., Oh H. J., Kim M. J., Kim G. A., Choi Y. B., Jo Y. K., et al. (2017). Oocyte Maturation-Related Gene Expression in the Canine Oviduct, Cumulus Cells, and Oocytes and Effect of Co-culture with Oviduct Cells on In Vitro Maturation of Oocytes. J. Assist. Reprod. Genet. 34 (7), 929–938. doi:10.1007/s10815-017-0910-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Li H., Durbin R. (2009). Fast and Accurate Short Read Alignment with Burrows-Wheeler Transform. Bioinformatics 25 (14), 1754–1760. doi:10.1093/bioinformatics/btp324

PubMed Abstract | CrossRef Full Text | Google Scholar

Li X., Yue X., Pastor W. A., Lin L., Georges R., Chavez L., et al. (2016). Tet Proteins Influence the Balance Between Neuroectodermal and Mesodermal Fate Choice by Inhibiting Wnt Signaling. Proc. Natl. Acad. Sci. 113 (51), E8267–E8276. doi:10.1073/pnas.1617802113

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu F., Zhao H., Gong L., Yao L., Li Y., Zhang W. (2018). MicroRNA-129-3p Functions as a Tumor Suppressor in Serous Ovarian Cancer by Targeting BZW1. Int. J. Clin. Exp. Pathol. 11 (12), 5901–5908.

PubMed Abstract | Google Scholar

Liu Y., Zhou Z., He X., Tao L., Jiang Y., Lan R., et al. (2021). Integrated Analyses of miRNA-mRNA Expression Profiles of Ovaries Reveal the Crucial Interaction Networks that Regulate the Prolificacy of Goats in the Follicular Phase. BMC Genomics 22, 812. doi:10.1186/s12864-021-08156-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu T., Zou X., Liu G., Deng M., Sun B., Guo Y., et al. (2020). A Preliminary Study on the Characteristics of microRNAs in Ovarian Stroma and Follicles of Chuanzhong Black Goat during Estrus. Genes 11 (9), 970. doi:10.3390/genes11090970

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo J., Wang W., Sun S. (2019). Research Advances in Reproduction for Dairy Goats. Asian-Australas J. Anim. Sci. 32 (8), 1284–1295. doi:10.5713/ajas.19.0486

PubMed Abstract | CrossRef Full Text | Google Scholar

Lv X., Chen W., Sun W., Hussain Z., Chen L., Wang S., et al. (2020). Expression Profile Analysis to Identify Circular RNA Expression Signatures in Hair Follicle of Hu Sheep Lambskin. Genomics 112, 4454–4462. doi:10.1016/j.ygeno.2020.07.046

PubMed Abstract | CrossRef Full Text | Google Scholar

Mbarek H., Steinberg S., Nyholt D. R., Gordon S. D., Miller M. B., McRae A. F., et al. (2016). Identification of Common Genetic Variants Influencing Spontaneous Dizygotic Twinning and Female Fertility. Am. J. Hum. Genet. 98 (5), 898–908. doi:10.1016/j.ajhg.2016.03.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Memczak S., Jens M., Elefsinioti A., Torti F., Krueger J., Rybak A., et al. (2013). Circular RNAs Are a Large Class of Animal RNAs with Regulatory Potency. Nature 495, 333–338. doi:10.1038/nature11928

PubMed Abstract | CrossRef Full Text | Google Scholar

Messina A., Langlet F., Chachlaki K., Roa J., Rasika S., Jouy N., et al. (2016). Corrigendum: A microRNA Switch Regulates the Rise in Hypothalamic GnRH Production before Puberty. Nat. Neurosci. 19 (8), 1115. doi:10.1038/nn0816-1115a

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagarjuna R. P., Osnat B., Marvin J., Reut A. F., Christin S., Larissa R., et al. (2017). Translation of CircRNAs. Mol. Cel. 66 (1), 9–21. doi:10.1016/j.molcel.2017.02.021

CrossRef Full Text | Google Scholar

Notter D. R. (2012). Genetic Improvement of Reproductive Efficiency of Sheep and Goats. Anim. Reprod. Sci. 130, 147–151. doi:10.1016/j.anireprosci.2012.01.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Perera B. P. U., Ghimire S., Kim J. (2018). Circular RNA Identified from Peg3 and Igf2r. PLoS One 13 (9), e0203850. doi:10.1371/journal.pone.0203850

PubMed Abstract | CrossRef Full Text | Google Scholar

Quetglas M., Adona P., de Bem T., Pires P., Leal C. (2010). Effect of Cyclin-dependent Kinase (CDK) Inhibition on Expression, Localization and Activity of Maturation Promoting Factor (MPF) and Mitogen Activated Protein Kinase (MAPK) in Bovine Oocytes. Reprod. Domest. Anim. 45 (6), 1074–1081. doi:10.1111/j.1439-0531.2009.01498.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Rybak-Wolf A., Stottmeister C., Glažar P., Jens M., Pino N., Giusti S., et al. (2015). Circular RNAs in the Mammalian Brain Are Highly Abundant, Conserved, and Dynamically Expressed. Mol. Cel. 58 (5), 870–885. doi:10.1016/j.molcel.2015.03.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen M., Li T., Zhang G., Wu P., Chen F., Lou Q., et al. (2019). Dynamic Expression and Functional Analysis of circRNA in Granulosa Cells during Follicular Development in Chicken. BMC Genomics 20, 96. doi:10.1186/s12864-019-5462-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Song Y.-N., Shi L.-L., Liu Z.-Q., Qiu G.-F. (2014). Global Analysis of the Ovarian microRNA Transcriptome: Implication for miR-2 and miR-133 Regulation of Oocyte Meiosis in the Chinese Mitten Crab, Eriocheir Sinensis (Crustacea:Decapoda). BMC Genomics 15 (1), 547. doi:10.1186/1471-2164-15-547

PubMed Abstract | CrossRef Full Text | Google Scholar

Song Y., Zhang L., Liu X., Niu M., Cui J., Che S., et al. (2019). Analyses of circRNA Profiling during the Development from Pre-receptive to Receptive Phases in the Goat Endometrium. J. Anim. Sci Biotechnol. 10, 34. doi:10.1186/s40104-019-0339-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang F., Chu M., Pan L., Wang X., He X., Zhang R., et al. (2021). Polymorphism Detection of GDF9 Gene and its Association with Litter Size in Luzhong Mutton Sheep (Ovis aries). Animals 11 (2), 571. doi:10.3390/ani11020571

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao M.-S., Ai Y., Wilusz J. E. (2020). Biogenesis and Functions of Circular RNAs Come into Focus. Trends Cel Biol. 30, 226–240. doi:10.1016/j.tcb.2019.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie C., Mao X., Huang J., Ding Y., Wu J., Dong S., et al. (2011). KOBAS 2.0: a Web Server for Annotation and Identification of Enriched Pathways and Diseases. Nucleic Acids Res. 39, W316–W322. doi:10.1093/nar/gkr483

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu S.-S., Gao L., Xie X.-L., Ren Y.-L., Shen Z.-Q., Wang F., et al. (2018). Genome-Wide Association Analyses Highlight the Potential for Different Genetic Mechanisms for Litter Size Among Sheep Breeds. Front. Genet. 9, 118. doi:10.3389/fgene.2018.00118

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye C. Y., Chen L., Liu C., Zhu Q. H., Fan L. (2015). Widespread Noncoding Circular RNA S in Plants. New Phytol. 208 (1), 88–95. doi:10.1111/nph.13585

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye J., Yao Z., Si W., Gao X., Yang C., Liu Y., et al. (2018). Identification and Characterization of microRNAs in the Pituitary of Pubescent Goats. Reprod. Biol. Endocrinol. 16, 51. doi:10.1186/s12958-018-0370-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ying S. J., Xiao S. H., Wang C. L., Zhong B. S., Zhang G. M., Wang Z. Y., et al. (2013). Effect of Nutrition on Plasma Lipid Profile and mRNA Levels of Ovarian Genes Involved in Steroid Hormone Synthesis in Hu Sheep during Luteal Phase1. Anim. Sci. 91 (11), 5229–5239. doi:10.2527/jas.2013-6450

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu D., Chen F., Zhang L., Wang H., Chen J., Zhang Z., et al. (2017). Smad9 Is a Key Player of Follicular Selection in Goose via Keeping the Balance of LHR Transcription. bioRxiv. doi:10.1101/213546

CrossRef Full Text | Google Scholar

Zhao Z., Zou X., Lu T., Deng M., Li Y., Guo Y., et al. (2020). Identification of mRNAs and lncRNAs Involved in the Regulation of Follicle Development in Goat. Front. Genet. 11, 589076. doi:10.3389/fgene.2020.589076

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng Q., Bao C., Guo W., Li S., Chen J., Chen B., et al. (2016). Circular RNA Profiling Reveals an Abundant circHIPK3 that Regulates Cell Growth by Sponging Multiple miRNAs. Nat. Commun. 7, 11215. doi:10.1038/ncomms11215

PubMed Abstract | CrossRef Full Text | Google Scholar

Zou X., Lu T., Zhao Z., Liu G., Lian Z., Guo Y., et al. (2020). Comprehensive Analysis of mRNAs and miRNAs in the Ovarian Follicles of Uniparous and Multiple Goats at Estrus Phase. BMC Genomics 21, 267. doi:10.1186/s12864-020-6671-4

PubMed Abstract | CrossRef Full Text | Google Scholar