Skip to main content

CircRNAs in Xiang pig ovaries among diestrus and estrus stages

Abstract

Background

The fecundity of sows is a trait of major economic in pig industry. The molecular regulation of estrus cycles can affect the fecundity of female animals. Compared with the other pig breeds, Xiang pig exhibits the special estrus behaviors. CircRNAs are thought to involve in regulation of multiple biological processes. However, the potential roles of circRNAs in ovary regulation on Xiang pig estrus are largely unknown.

Results

8,937 circRNAs were identified from eight libraries constructed from the ovarian samples of Xiang pig at estrus and diestrus stages by RNA sequencing method. Of which, 1,995 were high confidence circRNAs detected at least two junction reads in each ovary sample and seven circRNAs were validated by RT-PCR method. Furthermore, we identified 290 upregulated and 15 downregulated circRNAs in estrus ovaries. These differentially expressed circRNAs (DECs) derived from 273 host genes. And 207 miRNAs were identified to be targets sponged by 156 DECs with 432 binding sites, containing more than one miRNA binding site in each circRNA. Function enrichment analysis revealed that the host genes and the targets of miRNAs sponged by DECs were enriched in several reproduction-related signaling pathways, such as ovarian steroidogenesis, oocyte maturation, circadian rhythm, estrogen signaling pathway, GnRH signaling pathway, circadian entrainment, and oocyte meiosis. The circRNA-miRNA-mRNA networks revealed that 153 miRNAs interacting with 122 DECs and 86 miRNAs interacting with 84 DECs were involved in ovarian functions and ovarian circadian entrainment and circadian rhythm respectively. The DEC-miRNA-DEG (differentially expressed gene, DEG) networks associated with reproduction-related signaling pathways contained 22 DECs,18 miRNAs and 7 DEGs. 22 DECs were recognized as hub circRNAs during the estrus phase of Xiang pigs.

Conclusions

The circRNAs that function as miRNA sponges could play a key role in post-transcriptional regulation of gene expression during Xiang pig’s estrus cycle.

Background

Pork is the main meat consumed in human. Indigenous pigs contribute important part of pork production. Sow fecundity is pivotal in pig industry to ensure pork production amount, which is much diverse between pig breeds [1]. Ovary function directly influences the fecundity of female animals [2]. During each estrus cycle, the ovary undergoes a series of complex biological processes in morphological, hormonal, and biochemical changes. These biological processes in ovary are involved in the transcriptional and post-transcriptional regulation of many genes [3]. Previous studies have revealed that non-coding RNAs, including a variety of lncRNAs, miRNAs, and circRNAs are widely involved in post-transcriptional regulation of gene expression and various biological processes [4].

CircRNAs are a type of non-coding RNA with a covalently closed loop structure, which are derived from the back-spliced exonic or intronic sequences [5]. CircRNAs are abundant in the cytoplasm and can be co-expressed with the linear transcripts from which they are derived [6]. A great number of circRNAs have been identified in a variety of eukaryotic organisms by RNA-seq method [5]. CircRNAs contain miRNA binding sites and may function as miRNA sponges, or as transcriptional activators. Furthermore, circRNAs have been shown to segregate RNA binding proteins [7], and can even become translated into proteins through cap-independent translation initiation. Additionally, many evidences suggest that circRNAs are involved in a wide range of biological processes and function as ceRNA, which can influence the expression level of their parental genes [8, 9]. Recently, many studies using deep-sequencing approaches have reported that circRNAs are differentially expressed in ovary between pig breeds [8, 9]. Breed-specific circRNAs could be potentially associated to reproduction traits [10, 11].

Xiang pig is a miniature indigenous pig breed originated from the southeast in Guizhou province of China, and the meat were the dominant dietary intake sources for the mountain resident populations. It is featured by small size, early sexual maturity, lower litter size, excellent meat quality, and not clear exhibition of estrus behaviors [2]. The exhibition of estrus behaviors in Xiang gilts or sows was not very clearer than in Meishan or other European pigs [10, 11]. The difference on estrus expression between pig breeds could be affected by genes and could be improved by selection [12]. However, the potential roles of circRNAs in ovary regulation on Xiang pig estrus are largely unknown. The investigation of circRNAs in Xiang pig ovary may provide a valuable opportunity to understand the molecular basis of pig reproduction. To obtain more knowledge on the roles of circRNAs in Xiang pig ovary during estrus cycle, we performed a genome-wide analysis of transcripts of ovary tissues from Xiang pig sows at diestrus and estrus phases using RNA sequencing method. We investigated the expression profiles of circRNAs in ovarian tissues and identified differentially expressed circRNAs (DECs) (Fig. 1). Our results suggest that circRNAs probably participate in the regulation of gene expression in pig’s estrus cycle and reproduction.

Fig. 1
figure 1

The flow chart of present research

Results

Characterization of circRNAs expressed in Xiang pig ovaries at estrus and diestrus

To uncover the expression and function of circRNAs in porcine ovary, we performed total RNA sequencing in adult ovaries of Xiang pig during estrus and diestrus. The sequencing of the RNA libraries generated about 152 million 100-bp paired-end reads and an average of 124.3 million mapped reads per sample (Additional file 1: Table S1). We used find_circ and CIRI2 for circRNA detection and identified circRNAs with a junction reads greater than 2, which were detected in at least one sample of the biological replicates. We detected 4,404 circRNAs from the diestrus group, 7,680 circRNAs from the estrus group. In total, 8,937 unique circRNAs were detected from both software in two group samples (Additional file 2: Table S2), whereas 3,147 circRNAs were co-expressed in ovaries between two stages. The length of these circRNAs were shorter than 2 kb, which varied from 32 to 1,814 bases in estrus group, 29 to 1,146 bases in diestrus group, but mainly enriched between 100–400 bases (Fig. 2A). The composition and sources of these circRNAs between estrus and diestrus samples did not present significant difference. The percentages of the circRNAs that derived from exons, introns, and intergenic regions in estrus samples were 91.55%, 4.21%, and 4.6%, respectively, while they were 91.62%, 3.43%, and 5.4% in diestrus samples, respectively (Fig. 2B). According to their parental gene locations, these circRNAs were widely distributed on 1–18 autosomes, X chromosome, and mitochondria (Fig. 2C). It was found that the numbers of circRNAs from Chromosome 1, 6, and 13 were more than that from other chromosomes. Most (> 98%) circRNAs were derived from multiple exons, of which circRNAs with 2–3 exons accounted for > 87.3% within the same parental gene (Fig. 2D). Very few circRNAs were composed of more than 5 exons. The most exon composition (n = 11) was found in circ_2424 (host gene, ASPH). Additionally, 42.59% and 51.2% of the parental genes generated more than one circRNA in ovary samples at estrus and diestrus, respectively (Fig. 2E). We found that the gene DNA helicase (CHD2) had 28 predicted circRNAs.

Fig. 2
figure 2

Characterization of circRNAs expressed in Xiang pig ovaries at estrus and diestrus. A Amounts of circRNAs with different length. B Distribution of circRNAs in genic regions. C Distribution of circRNAs in chromosomes. D Exon numbers contained in circRNAs. E Characteristics of host genes producing circRNAs

Analysis of DECs in Xiang pig ovaries between estrus and diestrus

To evaluate the dynamic changes of circRNA expression in Xiang pig ovaries between estrus and diestrus, we focused on circRNAs that were detected at least two junction reads in each biological replicate of one specific tissue for high confidence. This analysis yielded 1,995 high confidence circRNAs. We performed differential expression analysis on the high confident circRNAs with edgeR. Compared with diestrus ovaries, 305 DECs were detected from estrus ovaries. Of these, 294 DECs were derived from 273 host genes (P < 0.05) with 290 upregulated and 15 downregulated DECs (Fig. 3A, Additional file 3: Table S3). In addition, we clustered the DECs in both diestrus and estrus ovaries. As shown in the heatmap (Fig. 3B), samples at the same stages were clustered together, and the expression levels of circRNAs exhibited dynamic changed during estrus cycle.

Fig. 3
figure 3

Volcano plot and heatmap analysis of DECs. A Volcano plot of DECs. B Heatmap of DECs Diestrus samples: 175,177,178, and 179; Estrus samples:193, 194,195, and 176

CircRNAs validation by RT-PCR method

To validate the reliability of predicted circRNAs, seven circRNAs were randomly chosen from the top 50 of CPM value in 1,995 high confidence circRNAs at estrus and diestrus phases and verified the region of spliced junction by RT-PCR method using a pair of divergent primers. These circRNAs consisted of 6 exons (circ_1952), 5 exons (circ_8664, circ_2414), 3 exons (circ_5597, circ_4508), 2 exons (circ_8670) and 1 exon (circ_0546) respectively. In Fig. 4A–G, the circularized states and the determined sequence at the junction of the 7 circRNAs were presented. The junction regions of circRNAs sequenced by Sanger method were consistent with those putative analysis of circRNAs by RNA-seq.

Fig. 4
figure 4

Diagram of circRNAs generation by a way of back-splicing and the sequence detected by Sanger method. A Circ_1952 was generated from exon 2, 3, 4, 5, 6, and 7 of gene FANCL. B Circ_8664 was generated from exon 2, 3, 4, 5, and 6 of gene SUGCT. C Circ_5597 was produced from exon 2, 3, and 4 of gene PAN3. D Circ_4508 was produced from exon 4, 5, and 6 of gene ELF2. E Circ_2414 was derived from exon 8, 9, 10, 11, and 12 of gene CSPP1. F Circ_0546 was derived from exon 2 of gene IGF1R. G Circ_8670 was generated at least from exon 2 and 5 of gene CDK13

Host gene function analysis of expressed circRNAs

To understand the biological functions of circRNAs, we performed GO and KEGG pathway analysis to predict the functions of circRNAs in Xiang pig ovaries during diestrus and estrus (Additional file 4: Table S4). GO analysis indicated that the host genes of circRNAs were significantly associated with cellular metabolic process, cellular component organization or biogenesis, regulation of catabolic process, biological regulation, binding and activity, and regulation of transcription by RNA polymerase II (P < 0.05) (Fig. 5A). Pathway analysis indicated that the host genes of circRNAs were involved in valine, leucine and isoleucine degradation, propanoate metabolism, human T-cell leukemia virus 1 infection, endocytosis, progesterone-mediated oocyte maturation, protein processing in endoplasmic reticulum, Fc epsilon RI signaling pathway, mitophagy–animal, and neurotrophin signaling pathway (P < 0.05) (Fig. 5B).

Fig. 5
figure 5

Go and KEGG analysis for host genes of circRNAs at estrus and diestrus phases. A Go functional analysis. B KEGG enrichment analysis

Prediction of miRNA targets potentially sponged by DECs

In this study, we used miRanda v3.3 and RNAhybrid v1.2 to predict the targeted miRNAs of DECs. We identified 432 binding sites in 156 DECs to bind with 207 miRNA molecules by all two prediction programs (Additional file 5: Table S5, Additional file 16: Fig. S1). We found that any one of circRNAs could contain 1–5 miRNA-binding sites. Accordingly, the same miRNA may bound with multiple circRNAs.

Many reports have proposed a ceRNA hypothesis that circRNAs can act as endogenous “miRNA sponge”, thereby modulating the expression of miRNA targets [7,8,9]. To acquire more knowledge concerning the biological role of DECs in ovary function, regulatory targets of those miRNAs sponged by DECs were predicted by miRanda v3.3 and PITA in the present study. A total of 9,307 potential targets (Additional file 6: Table S6) for the 207 miRNAs that bound by 156 DECs were predicted. Pathway analysis showed these predicted targets were involved in 292 possible KEGG pathways (P < 0.05) (Additional file 7: Table S7), including diseases (n = 85), signal transduction (n = 50), metabolism (n = 34), cell cycles, apoptosis, cell communication, and other pathways. The top 20 of KEGG enrichment were showed in Fig. 6A. Many pathways were associated with ovarian functions (Fig. 6B). Significantly, several pathways were associated with reproduction, such as oocyte meiosis, circadian entrainment, circadian rhythm, ovarian steroidogenesis, estrogen signaling pathway, progesterone-mediated oocyte maturation, and GnRH signaling pathway. Therefore, we further analyzed the interaction between circRNAs, miRNAs, and mRNAs of these pathways associated with reproduction in detail (Additional file 8: Table S8, Additional files 1718, 19, 20, 21, 22, 23: Figs. S2–S8). The circRNAs-miRNAs-mRNAs interactive networks involved in these KEGG pathways were constructed based on ceRNA mechanism. The resulting circRNA-miRNA-mRNA interactive networks provided nodes and linkages between circRNAs and their targets. Comparison of the circRNA-miRNA-mRNA networks between diestrus and estrus groups revealed that 153 miRNAs interacting with 122 DECs associated with ovarian functions and 86 miRNAs interacting with 84 DECs related with ovarian circadian entrainment and circadian rhythm (Additional file 9: Table S9).

Fig. 6
figure 6

KEGG enrichment analysis for protein-coding target genes of miRNAs sponged by DECs. A The top 20 of KEGG enrichment. B The KEGG pathways were related ovarian functions

To further understand the role of the key circRNAs in the regulation of gene expression during estrus stage, we performed the analysis of DEGs and constructed the DEC-miRNA-DEG interactive networks. At gene level, we identifieda total of 1,315 genes differentially expressed in ovaries between estrus and diestrus phases (Additional file 10: Table S10, Fig. 7A), of which 924 genes were upregulated and 391 genes were downregulated in estrus ovaries. Pathway analysis using KOBAS program showed DEGs participated in 51 possible KEGG pathways (P < 0.05) (Additional file 11: Table S11). Significantly, four pathways were associated with reproduction (Fig. 7B), such as steroid biosynthesis, cell cycle, ovarian steroidogenesis, and steroid hormone biosynthesis. And then we had predicted the 137 DECs-182 miRNAs-571 DEGs interactive networks (Additional file 12: Table S12). We found that 175 miRNAs which targeted 531 upregulated genes were sponged by 130 upregulated circRNAs and 12 miRNAs that also targeted 40 downregulated genes were competed by 7 downregulated circRNAs. At the same time, we further analyzed the interaction between DECs, miRNAs, and DEGs of above-mentioned pathways associated with reproduction (oocyte meiosis, circadian entrainment, circadian rhythm, ovarian steroidogenesis, estrogen signaling pathway, progesterone-mediated oocyte maturation, and GnRH signaling pathway). The DEC-miRNA-DEG networks associated with reproduction contained 22 DECs, 18 miRNAs, and 7 DEGs (Fig. 8, Additional file 13: Table S13).

Fig. 7
figure 7

Volcano plot and the KEGG pathways of ovarian functions for DEGs. A Volcano plot of DEGs. B The KEGG pathways of ovarian functions for DEGs

Fig. 8
figure 8

The DEC-miRNA-DEG networks associated with reproduction pathways

Discussion

Ovary is an important reproductive organ of female animal. They provide fertile oocytes, secrete reproductive hormones, and maintain the estrus cycle of female animals. In each estrus cycle, the ovary undergoes changes of proliferation, invasion, differentiation, and apoptosis, which involve in the transcriptional regulation of many genes. These physiological changes directly affect or determine the ovulation, fertilization rate and litter size of animals [13]. CircRNAs are a new class of endogenous non-coding RNAs that have been found to be widely expressed in human and animal cells and function in many biological processes [5, 14, 15]. Previous studies indicate that circRNAs are involved in regulation and may serve as novel regulators of ovarian follicle growth and development during porcine reproduction processes [10, 16,17,18]. In this study, we used RNA-seq to investigate the circRNA profiles in ovaries of Xiang pig sows during diestrus and estrus phases. We identified 8,937 circRNAs in the ovaries from diestrus and estrus groups. The circRNAs generated from estrus ovaries were more than that from diestrus ovaries. However, there were 3,147 circRNAs were co-expressed in ovaries between two stages. We found that over 91% of circRNAs were comprised of exonic sequences, containing one or more exons. Furthermore, about 5% of the circRNAs were generated from intergenic regions. The length of these circRNAs mainly enriched between 100 and 400 bases. These observations were similar to previous findings in pig ovaries [19]. It suggested that the circRNAs might play a key role in ovarian functions for Xiang pig.

Furthermore, 42.59 ~ 51.2% of the parental genes generated more than one circRNA in ovary samples at estrus and diestrus. The most predicted circRNAs (n = 28) were derived from the gene CHD2 (Table S14). Genetic studies demonstrate that Chromatin remodeling enzymes play critical roles in organizing genomic DNA within the native chromatin state. CHD2 is a member of the chromodomain helicase DNA-binding family of proteins and regulates gene expression through chromatin remodeling. Chd2-deficient mice have been demonstrated to exhibit a general growth delay and perinatal lethality [20]. These studies suggest that CHD2 plays an intrinsic role in normal mammalian development [21, 22]. In our study, we obtained 9 high confidence circRNAs from 28 predicted circRNAs derived from CHD2, of which 3 circRNAs were differentially expressed between estrus and diestrus ovaries. However, we didn’t predict any putative miRNA binding sites for these 3 DECs. These indicated these DECs derived from CHD2 do not act as miRNA sponges. Their function in ovaries should be further studied.

We identified 305 DECs (Additional file 3: Table S3). The identification of these DECs demonstrated that there was a larger difference in expression and regulation of genes between diestrus and estrus ovaries. At gene level, a total of 1,315 genes differentially expressed in ovaries between estrus and diestrus phases (Additional file 10: Table S10). By comparison with the expression in diestrus ovaries, 924 genes were upregulated and 391 genes were downregulated in estrus ovaries. Of the 273 host genes for 305 DECs, we found that only 20 host genes were differentially expressed. These observations suggested that gene expression during estrus cycle was regulated at different levels and the expression of circRNAs was highly regulated and controlled.

It is generally accepted that the functions of circRNAs are related to the functions of their parent genes [6]. Previous studies have indicated that many differentially expressed coding and noncoding RNAs are widely expressed in the diestrus or estrus stages and several genes could affect estrus of animals [23]. In our study, we found that eighteen genes, which have known to be involved in ovary functions or other reproductive processes, produced more than one circRNAs (Additional file 3: Table S3). For instance, CPEB1, MAPK9, PIK3CA, KDM1A, Dicer1, and FBXL3 contained two to three predicted circRNAs, respectively. These genes participate in several signaling pathways to regulate steroid hormone biosynthesis, ovarian circadian rhythm, oocyte development and maturation [24]. Several circRNAs (circ_3954, circ_1268, circ_6675, circ_6675, circ_4235, and circ_5758) were significantly upregulated in estrus ovary compared with diestrus ovary. The results illustrated that the circRNAs produced by these genes might involve in hormone biosynthesis, oocyte maturation and estrus cycle maintenance.

Even now the function of circRNAs is not well understood [25]. Many studies have revealed that some circRNAs can serve as a miRNA sponge and subsequently suppress its activity to regulate gene expression, or cross-talk with transcriptional machinery [7, 26].

We predicted 432 miRNA targeted by 156 of 305 DECs (Additional file 5: Table S5). Each circRNA contains one or more miRNA binding sites. For example, circ_1968 sponged ssc-let-7a, ssc-let-7d-5p, and ssc-let-7f-5p. Circ_2456 sponged ssc-miR-21-3p. Circ_2141 contained 9 potential binding sites to 4 different miRNAs, including ssc-miR-34a, ssc-miR-133a-5p, ssc-miR-138, and ssc-miR-7137-3p. Circ_1268 harbored two binding sites with ssc-miR-191. Circ_7558 functioned as sponges for ssc-miR-132 and ssc-miR-374b-5p. Previous studies have indicated that numerous miRNAs involved in regulating female reproductive hormone signaling during estrus [27,28,29,30]. MiR-191, miR-132, miR-370, and miR-181a were found to be associated with follicular development. Furthermore, miR-19b, miR-21, miR-31, miR-106a, and miR-224 were associated with follicular granule cell development. MiR-133 has been demonstrated to regulate oocyte meiosis and suppress ovarian cancer cell proliferation [31, 32].

To fully understand the biological role of miRNAs sponged by DECs in ovary function, we analyzed the targets of the miRNAs sponged by DECs. Pathway analysis showed that these predicted targets participate in a lot of KEGG pathways, including several pathways associated with reproduction, such as oocyte meiosis, circadian entrainment, circadian rhythm, ovarian steroidogenesis, estrogen signaling pathway, progesterone-mediated oocyte maturation, GnRH signaling pathway. Compare of the circRNA-miRNA-mRNA networks between diestrus and estrus groups revealed that 153 miRNAs interacting with 122 DECs associated with ovarian functions, and 86 miRNAs interacting with 84 DECs related with ovarian circadian entrainment and circadian rhythm processes (Additional file 9: Table S9).

To further illustrate the role of the circRNAs in the regulation of gene expression during estrus stage, we constructed the DEC-miRNA-DEG interactive networks (Additional file 12: Table S12). These networks revealed that so much circRNAs were interacted with miRNAs and may act as ceRNAs to medicate the expression of miRNA targets. Many of the target genes, such as PTK2B, RPS6KA3, CCNB3, PLCB1, FOS, and HBEGF, were associated with reproduction pathways. Thus, we hypothesized that these cirRNAs, miRNAs, and mRNAs in the DEC-miRNA-DEG networks could play critical roles in estrus regulation. For example, we found that thirteen circRNAs (circ_1429, circ_7842, circ_8618, circ_4891, circ_7904, circ_2367, circ_3536, circ_5076, circ_1401, circ_7986, circ_2720, circ_3338, and circ_1081) were interacted with twelve miRNAs (ssc-miR-10382, ssc-miR-1224, ssc-miR-149, ssc-miR-2320-5p, ssc-miR-326, ssc-miR-4339, ssc-miR-664-5p, ssc-miR-7138-5p, ssc-miR-7143-5p, ssc-miR-9802-3p, ssc-miR-9820-5p, and ssc-miR-9847-3p) as ceRNAs to medicate the expression of the related genes such as HBEGF, PTK2B, FOS, and ENSSSCG00000006719. Previous reports have shown that HBEGF may actively control the process of follicle growth and maturation in the zebrafish [33]. PTK2B can encode a cytoplasmic protein tyrosine kinase which was reported to participate in the development of the final stages of follicles and ovary development in Hu sheep [34]. And FOS, a critical downstream mediator of PGR and EGF signaling, plays an important role in ovulation [35]. In our study, the expression of HBEGF, PTK2B, and FOS gene were upregulated in estrus ovaries and the expression patterns were consistent with the thirteen circRNAs motioned above. The results showed that the ovary needs to express enough proteins, such as HBEGF, PTK2B, and FOS to meet the biological processes of follicular development and oocyte maturation during estrus cycle. Therefore, in order to remove the translation inhibition of target mRNA by the miRNAs, the cells would increase the expression level of related circRNAs as a miRNA sponge, thereby inhibiting the post transcriptional regulation of target mRNA by miRNAs.

Although some circRNAs function as miRNA sponges, 149 of 305 (48.85%) Xiang pig circRNAs identified in this study have no putative miRNA binding site. Several studies have suggested that most circRNAs do not act as miRNA sponges and they have functions including regulation of host gene transcription, protein binding, and translation [7,8,9].

Conclusions

In conclusion, our results demonstrated that ovaries generated abundant circRNAs during estrus, of which numerous circRNAs were differentially expressed between diestrus and estrus phases. We predicted 432 miRNA targets by 156 DECs. Each circRNA can contain one or more miRNA binding sites. Function enrichment analysis revealed that their host genes and the targets of miRNAs sponged by DECs were enriched in several reproduction-related signaling pathways, such as ovarian steroidogenesis, oocyte maturation, circadian rhythm, estrogen signaling pathway, GnRH signaling pathway. The DEC-miRNA-DEG networks associated with reproduction-related signaling pathways contained notes of 22 DECs, 18 miRNAs, and 7 DEGs. 22 DECs were recognized as hub circRNAs during the estrus phase of Xiang pigs. These results suggest that circRNAs probably participate in the regulation of gene expression in pig’s estrus cycle and reproduction.

Materials and methods

Sample preparation

Xiang pig sows after weaning were obtained from the Guizhou Dachang pig breeding company, Guizhou, China. The animal preparation and estrus detection were referred to the methods of Tang et al. (2018) and Ran et al. (2021) [2, 36]. Four animals from each group were regarded as biological replicates. The animals from post-weaning sows were monitored twice daily for behavioral estrus. On Day 10 and on Day 20 after estrus, the sows were considered in the diestrus and estrus phase, respectively. The ovarian samples were collected with surgery at 10 days before expected estrus and the day of the third estrus when the sows exhibited strong performance of reddening and swelling of the vulva. All samples were immediately frozen in liquid nitrogen and stored at − 80℃ until RNA extraction.

Library construction and RNA sequencing

The total RNA samples were isolated from ovarian tissues at diestrus and estrus stages using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. The quantity and integrity of the total RNAs were analyzed using a NanoDrop 2000 (Thermo Fisher Scientific Inc., Waltham, MA, USA) and an Agilent 2100 Size Bio-analyzer system (Agilent Technologies, CA, USA). The values of RNA integrity number (RIN) > 7.0 were used for RNA-seq analysis. About 5 μg total RNAs per sample were used for sequencing library preparation using a NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA). The ribosomal RNAs were removed from the total RNA using a Ribo-Zero™ GoldKits (Epicentre, Madison, WI, USA). The remaining RNA was fragmented and reverse-transcribed according to the description of TruSeq RNA LT/HT sample preparation kit (Illumina, USA). After the quality of the cDNA libraries were qualified by Bioanalyzer 2200 evaluation (Agilent, Santa Clara, CA), sequencing was conducted on an Illumina HiSeq 2500 instrument (Illumina, San Diego, CA, USA).

CircRNA identification

Firstly, the reads were removed that contained low-quality, adaptor and more than 5% unknown nucleotides via Fastp software [37]. The remaining high quality reads were then used for subsequent analysis. These reads were aligned to the reference genome of Sus scrofa (Sscrofa11.1) by employing Bowite2 or BWA [38, 39]. Bam files of unmapped reads from Bowite2 and sam files of mapped reads from BWA were input to find_circ and CIRI2, respectively. The candidates from two softwares were intersected to obtain the final circRNAs dataset based on chromosome location. The sequence splicing of circRNAs and the visualization of ring construction diagram were performed by using CIRI-full and CIRI-vis, respectively [40].

Expression analysis of circRNAs and mRNAs

The different expression patterns of circRNAs and mRNAs were calculated by edgeR package [41]. The reads numbers mapped to each circRNA were counted and the average of two softwares find_circ and CIRI2 was taken as expression level of each circRNA transcript. It was worth noting that the junction reads in a circRNA were greater than or equal to 2 in each sample. All CPM values of circRNA were added 0.1 for logarithm arithmetic. CPM = (circRNA read counts * 106) / the sum of circRNAs read counts [14]. The protein-coding gene expression level was counted by featureCounts software [42]. The expression level of mRNA was estimated by CPM value. CPM = (Count of reads mapped to a mRNA * 106) / Total count of mapped reads from the library [43]. The thresholds of DECs and DEGs were |log2 (fold_change)|≥ 1 and P < 0.05 [16].

Prediction of miRNA targets and circRNA-miRNA-mRNA network construction

Target miRNAs of DECs were predicted via MiRanda v3.3 (http://www.microrna.org/microrna/home.do) and RNAhybrid v1.2 (http://bibiserv.techfak.uni-bielefeld.de/rnahybrid). The intersections results from miRanda and RNAhybrid was identified as the target miRNAs with the minimal free energy of − 10 kcal/mol. MiRNAs sequences of pig were originated from miRBase database (http://www.mirbase.org/). The target mRNA by miRNA were determined as the shared mRNA between results from miRanda and PITA (http://genie.weizmann.ac.il/pubs/mir07/mir07_dat-a.html) with a set of minimal free energy to be − 10 kcal/mol. The diagram of circRNA-miRNA-mRNA regulatory interaction networks was depicted via Cytoscape software [44].

Go and KEGG analysis

To analyze function of circRNAs, Go and KEGG analysis were performed by KOBAS online [45], in which host genes of circRNAs and target genes of miRNAs sponged by DECs were taken as input. The established criteria (P < 0.05) was considered to indicate significant enrichment.

Validation of circRNAs by Sanger sequencing

Most of identified circRNAs were expressed at low levels with CPM value less than 1,000 especially those DECs between the estrus and diestrus phases. To confirm the reliability of the predicted circRNAs from RNA sequencing, seven circRNAs were randomly chosen from those circRNAs with CPM value (average CPM > 2500) ranked at the top 50. These seven circRNAs included circ_1952, circ_8664, circ_5597, circ_4508, circ_2414, circ_0546, and circ_8670. The information (eg: CPM, FDR, and so on) of those circRNAs was showed in Additional file 15: Table S15 in manuscript. The difference of these 7 circRNAs levels was not significant between estrus and diestrus stages. Moreover, we confirmed the reliability of the predicted circRNAs from RNA sequencing by verifying the region of spliced junction through RT-PCR method using a pair of divergent primers flanking the BSJ. The region of BSJ of circRNA is composed of a canonical 5′ splice site sequence joined to an upstream 3′ splice site sequence [46]. The PCR products were analyzed by 2% agarose gel electrophoresis and further proofed via Sanger sequencing (Sango Biotech, Shanghai, China).

Availability of data and materials

The data supporting the conclusions of this article are included within the article and its addition files. The raw sequencing data are available from SRA database in NCBI with accession number PRJNA835060.

Abbreviations

CircRNA:

Circular RNA

DEC:

Differentially expressed circRNA

MiRNA:

MicroRNA

MRNA:

Messenger RNA

ceRNA:

Competitive endogenous RNA

LncRNA:

Long noncoding RNA

ASPH :

Aspartate-β-hydroxylase

CHD2:

Chromodomain helicase DNA binding protein 2

GO:

Gene ontology

KEGG:

Kyoto encyclopedia of genes and genomes

GnRH:

Gonadotropin-releasing hormone

FANCL :

FA complementation group L

SUGCT :

Succinyl-CoA:glutarate-CoA transferase

PAN3 :

Poly(A) specific ribonuclease subunit PAN3

ELF2:

E74 like ETS transcription factor 2

CSPP1:

Centrosome and spindle pole associated protein 1

IGF1R :

Insulin like growth factor 1 receptor

CDK13:

Cyclin dependent kinase 13

CPEB1:

Cytoplasmic polyadenylation element binding protein 1

MAPK9 :

Mitogen-activated protein kinase 9

PIK3CA :

Phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit alpha

KDM1A :

Lysine-specific histone demethylase 1A

Dicer1:

Dicer 1, ribonuclease III

FBXL3:

F-box and leucine rich repeat protein 3

DEG:

Differentially expressed gene

RT-PCR:

Reverse transcription-polymerase chain reaction

BSJ:

Back-spliced junction

HBEGF :

Heparin-binding epidermal growth factor-like growth factor

PTK2B :

Protein tyrosine kinase 2 beta

RPS6KA3:

Ribosomal Protein S6 Kinase A3

CCNB3:

Cyclin B3

PLCB1:

Phospholipase C Beta 1

FOS :

Fos proto-oncogene, AP-1 transcription factor subunit

PGR :

Progesterone receptor

EGF :

Epidermal growth factor

PLC:

Phospholipase C

RNA-seq:

RNA sequencing

References

  1. Knox RV. Impact of swine reproductive technologies on pig and global food production. Adv Exp Med Biol. 2014;752:131–60.

    Article  PubMed  Google Scholar 

  2. Tang LT, Ran XQ, Mao N, Zhang FP, Niu X, Ruan YQ, et al. Analysis of alternative splicing events by RNA sequencing in the ovaries of Xiang pig at estrous and diestrous. Theriogenology. 2018;119:60–8.

    Article  CAS  PubMed  Google Scholar 

  3. Miao XY, Luo QM. Genome-wide transcriptome analysis between Small-tail Han sheep and the Surabaya fur sheep using high-throughput RNA sequencing. Reproduction. 2013;145(6):587–96.

    Article  CAS  PubMed  Google Scholar 

  4. Zhu LW, Li N, Sun LB, Zheng DW, Shao GF. Non-coding RNAs: The key detectors and regulators in cardiovascular disease. Genomics. 2021;113(1 Pt 2):1233–46.

    Article  CAS  PubMed  Google Scholar 

  5. Barrett SP, Salzman J. Circular RNAs: analysis, expression and potential functions. Development. 2016;143(11):1838–47.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, 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 

  7. Granados-Riveron JT, Aquino-Jarquin G. The complexity of the translation ability of circRNAs. Bba-gene Regul Mech. 2016;1859(10):1245–51.

    CAS  Google Scholar 

  8. Xu GX, Zhang JF, Li XJ, Hu H, Yang GS, Sun SD. Genome-wide differential expression profiling of ovarian circRNAs associated with litter size in pigs. Front Genet. 2019;10:1011.

    Google Scholar 

  9. Hu HY, Xi JZ, Zhou B, Zhang J, Li ZQ, Liu ZW, et al. Ovarian circular RNAs associated with high and low fertility in Large White sows during the follicular and luteal phases of the estrous cycle. Animals. 2020;10(4):696.

    Article  PubMed Central  Google Scholar 

  10. Xie S, Li MX, Chen YS, Liu YL, Ma P, Sun XM, et al. Identification of circular RNAs in the ovarian follicles of Meishan and Duroc sows during the follicular phase. J Ovarian Res. 2020;13(1):104.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Liang GM, Yan JY, Guo J, Tang ZL. Identification of ovarian circular RNAs and differential expression analysis between MeiShan and Large White Pigs. Animals. 2020;10(7):1114.

    Article  PubMed Central  Google Scholar 

  12. De Faria VR, Pinho RO, Camilo BS, Guimaraes JD, Silva FFE, Lopes PS, et al. Genes expression and phenotypic differences in corpus luteum and cumulus cells of commercial line and piau breed gilts. Theriogenology. 2019;136:111–7.

    Article  PubMed  CAS  Google Scholar 

  13. Lan DL, Xiong XR, Huang C, Mipam TD, Li J. Toward understanding the genetic basis of yak ovary reproduction: a characterization and comparative analyses of estrus ovary transcriptiome in yak and cattle. PLoS ONE. 2016;11(4): e0152675.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  14. Lu C, Shi X, Wang AY, Tao Y, Wang ZX, Huang CP, et al. RNA-Seq profiling of circular RNAs in human laryngeal squamous cell carcinomas. Mol Cancer. 2018;17(1):86.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  15. Cai HC, Li YM, Li HM, Niringiyumukiza JD, Zhang MD, Chen L, et al. Identification and characterization of human ovary-derived circular RNAs and their potential roles in ovarian aging. Aging-Us. 2018;10(9):2511–34.

    Article  CAS  Google Scholar 

  16. Cao ZB, Gao D, Xu TT, Zhang L, Tong X, Zhang DD, et al. Circular RNA profiling in the oocyte and cumulus cells reveals that circARMC4 is essential for porcine oocyte maturation. Aging-Us. 2019;11(18):8015–34.

    Article  CAS  Google Scholar 

  17. Guo TY, Huang L, Yao W, Du X, Li QQ, Ma ML, et al. The potential biological functions of circular RNAs during the initiation of atresia in pig follicles. Domest Anim Endocrinol. 2020;72: 106401.

    Article  CAS  PubMed  Google Scholar 

  18. Kim SH, Min KS, Kim NH, Yoon JT. Differential expression of programmed cell death on the follicular development in normal and miniature pig ovary. PLoS ONE. 2012;7(10): e46194.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Pan XC, Gong WT, He YT, Li NA, Zhang H, Zhang Z, et al. Ovary-derived circular RNAs profile analysis during the onset of puberty in gilts. BMC Genom. 2021;22(1):445.

    Article  CAS  Google Scholar 

  20. Marfella CG, Ohkawa Y, Coles AH, Garlick DS, Jones SN, Imbalzano AN. Mutation of the SNF2 family member Chd2 affects mouse development and survival. J Cell Physiol. 2006;209(1):162–71.

    Article  CAS  PubMed  Google Scholar 

  21. Marfella CG, Imbalzano AN. The CHD family of chromatin remodelers. Mutat Res. 2007;618(1–2):30–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Carvill GL, Heavin SB, Yendle SC, McMahon JM, O’Roak BJ, Cook J, et al. Targeted resequencing in epileptic encephalopathies identifies de novo mutations in CHD2 and SYNGAP1. Nat Genet. 2013;45(7):825–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. An XP, Zhang Y, Li F, Wang ZH, Yang SH, Cao BY. Whole transcriptome analysis: implication to estrous cycle regulation. Biol Basel. 2021;10(6):464.

    Article  CAS  Google Scholar 

  24. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Li X, Yang L, Chen LL. The biogenesis, functions, and challenges of circular RNAs. Mol Cell. 2018;71(3):428–42.

    Article  CAS  PubMed  Google Scholar 

  26. Bartsch D, Zirkel A, Kurian L. Characterization of circular RNAs (circRNA) associated with the translation machinery. Methods Mol Biol. 2018;1724:159–66.

    Article  CAS  PubMed  Google Scholar 

  27. Toms D, Pan B, Li JL. Endocrine regulation in the ovary by microRNA during the estrous cycle. Front Endocrinol (Lausanne). 2018;8:378.

    Article  Google Scholar 

  28. Maalouf SW, Liu WS, Pate JL. MicroRNA in ovarian function. Cell Tissue Res. 2016;363(1):7–18.

    Article  CAS  PubMed  Google Scholar 

  29. Tesfaye D, Gebremedhn S, Salilew-Wondim D, Hailay T, Hoelker M, Grosse-Brinkhaus C, et al. MicroRNAs: tiny molecules with a significant role in mammalian follicular and oocyte development. Reproduction. 2018;155(3):121–35.

    Article  Google Scholar 

  30. Cirillo F, Catellani C, Lazzeroni P, Sartori C, Nicoli A, Amarri S, et al. MiRNAs regulating insulin sensitivity are dysregulated in Polycystic Ovary Syndrome (PCOS) ovaries and are associated with markers of inflammation and insulin sensitivity. Front Endocrinol. 2019;10:879.

    Article  Google Scholar 

  31. Song YN, Shi LL, Liu ZQ, Qiu GF. 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 Genom. 2014;15(1):547.

    Article  CAS  Google Scholar 

  32. Guo JL, Xia BR, Meng FL, Lou G. MiR-133a suppresses ovarian cancer cell proliferation by directly targeting insulin-like growth factor 1 receptor. Tumour Biol. 2014;35(2):1557–64.

    Article  CAS  PubMed  Google Scholar 

  33. Liu KC, Ge W. Differential regulation of gonadotropin receptors (fshr and lhcgr) by epidermal growth factor (EGF) in the zebrafish ovary. Gen Comp Endocrinol. 2013;181:288–94.

    Article  CAS  PubMed  Google Scholar 

  34. Shabbir S, Boruah P, Xie L, Kulyar MF, Nawaz M, Yousuf S, et al. Genome-wide transcriptome profiling uncovers differential miRNAs and lncRNAs in ovaries of Hu sheep at different developmental stages. Sci Rep. 2021;11(1):5865.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Choi Y, Rosewell KL, Brännström M, Akin JW, Curry TE Jr, Jo M. FOS, a critical downstream mediator of PGR and EGF signaling necessary for ovulatory prostaglandins in the human ovary. J Clin Endocrinol Metab. 2018;103(11):4241–52.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Ran XQ, Hu FB, Mao N, Ruan YQ, Yi FL, Niu X, et al. Differences in gene expression and variable splicing events of ovaries between large and small litter size in Chinese Xiang pig. Porcine Health Manag. 2021;7(1):52.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Chen SF, Zhou YQ, Chen YR, Gu J. Fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):884–90.

    Article  CAS  Google Scholar 

  38. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9(4):357–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Gao Y, Wang JF, Zhao FQ. CIRI: an efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol. 2015;16(1):4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Zhong YY, Li LY, Chen ZT, Diao SQ, He YT, Zhang Z, et al. MiR143 inhibits steroidogenesis and induces apoptosis repressed by H3K27me3 in granulosa cells. Front Cell Dev Biol. 2020;8: 565261.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Robinson MD, McCarthy DJ, Smyth GK. EdgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    Article  CAS  PubMed  Google Scholar 

  42. Liao Y, Smyth GK, Shi W. FeatureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–30.

    Article  CAS  PubMed  Google Scholar 

  43. Corley SM, MacKenzie KL, Beverdam A, Roddam LF, Wilkins MR. Differentially expressed genes from RNA-Seq and functional enrichment results are affected by the choice of single-end versus paired-end reads and stranded versus non-stranded protocols. BMC Genom. 2017;18(1):399.

    Article  CAS  Google Scholar 

  44. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Xie C, Mao XZ, Huang JJ, Ding Y, Wu JM, Dong S, et al. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 2011;32:316–22.

    Article  CAS  Google Scholar 

  46. Jeck WR, Sorrentino JA, Wang K, Slevin MK, Burd CE, Liu JZ, et al. Circular RNAs are abundant, conserved, and associated with ALU repeats. RNA. 2013;19(2):141–57.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We would like to thank Guizhou Dachang pig breeding form, in Congjiang country, Guizhou province. We would like to thank TopEdit LLC (https://www.topeditsci.com) for the linguistic editing and proofreading during the preparation of this manuscript.

Funding

This research was supported and funded by the National Natural Science Foundation of China (No. 31960641 and No. 31672390), the Guizhou Provincial Science and Technology Projects (QKHPTRC[2019]5615, QKHRC[2016]4012, QKHZC[2017]2585, QKHZC[2017]2587, QKHPTRC [2018]5781, and QJHKYZ[2022]152), the National High Technology Research and Development Program of China (863 Program) (No. 2013AA102503).

Author information

Authors and Affiliations

Authors

Contributions

Conceptualization, XR JW; data curation, XN, YH, HL; investigation, SL, SH; writing-the draft, XN; writing-review and editing, XR, JW. All authors have read and approved to the final manuscript and agree to be personally accountable for author’s own contributions.

Corresponding authors

Correspondence to Xueqin Ran or Jiafu Wang.

Ethics declarations

Ethics approval and consent to participate

All animal procedures were approved by the guidelines of Guizhou University Subcommittee of Experimental Animal Ethics with no. of EAE-GZU-2020-P002.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have 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

Additional file 1. Table S1

. Overview of RNA-Seq data.

Additional file 2. Table S2

. The circRNAs expressed in the ovaries of Xiang pig at diestrus and estrus phases.

Additional file 3. Table S3

. High confident circRNAs and DECs between diestrus and estrus ovaries.

Additional file 4. Table S4

. GO and KEGG for host genes of circRNAs in estrus and diestrus.

Additional file 5. Table S5

. The target miRNAs sponged by DECs.

Additional file 6. Table S6

. Target genes of miRNAs sponged by DECs.

Additional file 7. Table S7

. GO and KEGG for target genes of miRNAs sponged by DECs.

Additional file 8. Table S8

. Interaction of circRNA-miRNA-mRNA.

Additional file 9. Table S9

. CircRNAs and miRNAs in pathways involves in ovarian function.

Additional file 10. Table S10

. The DEGs between diestrus and estrus ovaries.

Additional file 11. Table S11

. GO and KEGG for DEGs.

Additional file 12. Table S12

. Interaction of DEC-miRNA-DEG.

Additional file 13. Table S13

. Interaction of DEC-miRNA-DEG associated with reproduction pathways.

Additional file 14. Table S14

. CircRNAs were derived from the gene CHD2.

Additional file 15. Table S15

. The information of verified 7 circRNAs.

Additional file 16. Fig. S1.

Prediction of miRNA targets potentially sponged by DECs.

Additional file 17. Fig. S2.

The circRNA-miRNA-mRNA interactive networks of circadian rhythm pathway.

Additional file 18. Fig. S3.

The circRNA-miRNA-mRNA interactive networks of circadian entrainment pathway.

Additional file 19. Fig. S4.

The circRNA-miRNA-mRNA interactive networks of estrogen signaling pathway.

Additional file 20. Fig. S5.

The circRNA-miRNA-mRNA interactive networks of GnRH signaling pathway.

Additional file 21. Fig. S6.

The circRNA-miRNA-mRNA interactive networks of progesterone-mediated oocyte maturation pathway.

Additional file 22. Fig. S7.

The circRNA-miRNA-mRNA interactive networks of oocyte meiosis pathway.

Additional file 23. Fig. S8.

The circRNA-miRNA-mRNA interactive networks of ovarian steroidogenesis pathway.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, 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 changes were made. 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/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Niu, X., Huang, Y., Lu, H. et al. CircRNAs in Xiang pig ovaries among diestrus and estrus stages. Porc Health Manag 8, 29 (2022). https://doi.org/10.1186/s40813-022-00270-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40813-022-00270-1

Keywords