CircRNAs in Xiang pig ovaries among diestrus and estrus stages
Porcine Health Management volume 8, Article number: 29 (2022)
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.
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.
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.
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 . Ovary function directly influences the fecundity of female animals . 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 . 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 .
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 . CircRNAs are abundant in the cytoplasm and can be co-expressed with the linear transcripts from which they are derived . A great number of circRNAs have been identified in a variety of eukaryotic organisms by RNA-seq method . 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 , 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 . 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 . 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.
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.
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.
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.
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).
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 17, 18, 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).
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).
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 . 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 . 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 . 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 . 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 . 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 . 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 . 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 . 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 . And FOS, a critical downstream mediator of PGR and EGF signaling, plays an important role in ovulation . 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].
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
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).
Firstly, the reads were removed that contained low-quality, adaptor and more than 5% unknown nucleotides via Fastp software . 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 .
Expression analysis of circRNAs and mRNAs
The different expression patterns of circRNAs and mRNAs were calculated by edgeR package . 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 . The protein-coding gene expression level was counted by featureCounts software . 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 . The thresholds of DECs and DEGs were |log2 (fold_change)|≥ 1 and P < 0.05 .
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 .
Go and KEGG analysis
To analyze function of circRNAs, Go and KEGG analysis were performed by KOBAS online , 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 . 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.
Differentially expressed circRNA
Competitive endogenous RNA
Long noncoding RNA
- ASPH :
Chromodomain helicase DNA binding protein 2
Kyoto encyclopedia of genes and genomes
- FANCL :
FA complementation group L
- SUGCT :
- PAN3 :
Poly(A) specific ribonuclease subunit PAN3
E74 like ETS transcription factor 2
Centrosome and spindle pole associated protein 1
- IGF1R :
Insulin like growth factor 1 receptor
Cyclin dependent kinase 13
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
Dicer 1, ribonuclease III
F-box and leucine rich repeat protein 3
Differentially expressed gene
Reverse transcription-polymerase chain reaction
- HBEGF :
Heparin-binding epidermal growth factor-like growth factor
- PTK2B :
Protein tyrosine kinase 2 beta
Ribosomal Protein S6 Kinase A3
Phospholipase C Beta 1
- FOS :
Fos proto-oncogene, AP-1 transcription factor subunit
- PGR :
- EGF :
Epidermal growth factor
Knox RV. Impact of swine reproductive technologies on pig and global food production. Adv Exp Med Biol. 2014;752:131–60.
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.
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.
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.
Barrett SP, Salzman J. Circular RNAs: analysis, expression and potential functions. Development. 2016;143(11):1838–47.
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.
Granados-Riveron JT, Aquino-Jarquin G. The complexity of the translation ability of circRNAs. Bba-gene Regul Mech. 2016;1859(10):1245–51.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Marfella CG, Imbalzano AN. The CHD family of chromatin remodelers. Mutat Res. 2007;618(1–2):30–40.
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.
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.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Li X, Yang L, Chen LL. The biogenesis, functions, and challenges of circular RNAs. Mol Cell. 2018;71(3):428–42.
Bartsch D, Zirkel A, Kurian L. Characterization of circular RNAs (circRNA) associated with the translation machinery. Methods Mol Biol. 2018;1724:159–66.
Toms D, Pan B, Li JL. Endocrine regulation in the ovary by microRNA during the estrous cycle. Front Endocrinol (Lausanne). 2018;8:378.
Maalouf SW, Liu WS, Pate JL. MicroRNA in ovarian function. Cell Tissue Res. 2016;363(1):7–18.
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.
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.
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.
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.
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.
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.
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.
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.
Chen SF, Zhou YQ, Chen YR, Gu J. Fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):884–90.
Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9(4):357–9.
Gao Y, Wang JF, Zhao FQ. CIRI: an efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol. 2015;16(1):4.
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.
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.
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.
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.
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.
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.
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.
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.
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 (QKHPTRC5615, QKHRC4012, QKHZC2585, QKHZC2587, QKHPTRC 5781, and QJHKYZ152), the National High Technology Research and Development Program of China (863 Program) (No. 2013AA102503).
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
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
. Overview of RNA-Seq data.
. The circRNAs expressed in the ovaries of Xiang pig at diestrus and estrus phases.
. High confident circRNAs and DECs between diestrus and estrus ovaries.
. GO and KEGG for host genes of circRNAs in estrus and diestrus.
. The target miRNAs sponged by DECs.
. Target genes of miRNAs sponged by DECs.
. GO and KEGG for target genes of miRNAs sponged by DECs.
. Interaction of circRNA-miRNA-mRNA.
. CircRNAs and miRNAs in pathways involves in ovarian function.
. The DEGs between diestrus and estrus ovaries.
. GO and KEGG for DEGs.
. Interaction of DEC-miRNA-DEG.
. Interaction of DEC-miRNA-DEG associated with reproduction pathways.
. CircRNAs were derived from the gene CHD2.
. The information of verified 7 circRNAs.
Prediction of miRNA targets potentially sponged by DECs.
The circRNA-miRNA-mRNA interactive networks of circadian rhythm pathway.
The circRNA-miRNA-mRNA interactive networks of circadian entrainment pathway.
The circRNA-miRNA-mRNA interactive networks of estrogen signaling pathway.
The circRNA-miRNA-mRNA interactive networks of GnRH signaling pathway.
The circRNA-miRNA-mRNA interactive networks of progesterone-mediated oocyte maturation pathway.
The circRNA-miRNA-mRNA interactive networks of oocyte meiosis pathway.
The circRNA-miRNA-mRNA interactive networks of ovarian steroidogenesis pathway.
About this article
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