Ovary-derived circular RNAs profile analysis during the onset of puberty in gilts
BMC Genomics volume 22, Article number: 445 (2021)
In mammals, the ovary is the essential system of female reproduction for the onset of puberty, and the abnormal puberty has negative outcomes on health. CircRNA is a non-coding RNA produced by non-canonical alternative splicing (AS). Several studies have reported that circRNA is involved in the gene regulation and plays an important role in some human diseases. However, the contribution of circRNA has received little known within the onset of puberty in ovary.
Here, the profiles of ovarian circRNAs across pre-, in- and post-pubertal stages were established by RNA-sEq. In total, 972 circRNAs were identified, including 631 stage-specific circRNAs and 8 tissue-specific circRNAs. The biological functions of parental genes of circRNAs were enriched in steroid biosynthesis, autophagy-animal, MAPK signaling pathway, progesterone-mediated oocyte maturation and ras signaling pathway. Moreover, 5 circRNAs derived from 4 puberty-related genes (ESR1, JAK2, NF1 and ARNT) were found in this study. The A3SS events were the most alternative splicing, but IR events were likely to be arose in post-pubertal ovaries. Besides, the circRNA-miRNA-gene networks were explored for 10 differentially expressed circRNAs. Furthermore, the head-to-tail exon as well as the expressions of 10 circRNAs were validated by the divergent RT-qPCR and sanger sequencing.
In summary, the profiles of ovarian circRNAs were provided during pubertal transition in gilts, and these results provided useful information for the investigation on the onset of puberty at the ovarian-circRNAs-level in mammals.
Puberty is usually defined as the first estrus of mammals . In human, the abnormal puberty has negative effects in diseases such as asthma , psychosocial disorder , hypogonadism  and reproductive system tumors [5, 6]. It has been well recognized that the initiation of puberty is mainly driven by the hypothalamus-pituitary-ovary (HPO) axis [7,8,9]. Several studies have shown that it is possible to treat abnormal puberty by regulating HPO axis [10,11,12,13] Furthermore, the ovary is the female reproductive system and the endocrine organ, which produces the steroids and peptide hormones necessary for the onset of puberty [14,15,16]. Previous studies have reported that the female puberty failure is presented with the decreased ovarian weight and denormal hormone levels in mice  and human . Accumulating studies support that the noncoding RNA (ncRNAs) play a vital role in the expression of essential genes that regulate the oocyte growth and ovarian endocrine function [19, 20] as well as the onset of puberty [21, 22].
Circular RNAs (circRNAs) are the category of ncRNA molecules in the cytoplasm of eukaryotes, which are produced by non-canonical alternative splicing (AS) named back-splicing . It is reported that most circRNAs are composed of only exonic sequences, while a few circRNAs are composed of the exon-intronic or intronic sequences [24,25,26]. More recently, thousands of circRNAs are identified through the high-throughput RNA sequencing (RNA-seq). In mammals, it was found that circRNAs are tissue-specific and stage -specific as well as evolutionarily conserved [27,28,29,30], and it has been shown that exonic circRNAs show miRNA sponge activity, and intronic circRNAs are likely to regulate the transcription of their host genes [31, 32]. These findings suggest that circRNAs may play a pivotal role in growth and development of mammals.
Recently, circRNAs have been found to involve in asthma [33, 34] and reproductive system tumors , both of which are closely related to the abnormal puberty, and circRNAs are reported to involve in oocyte maturation and hormone synthesis . For example, Cao et al. showed that circRNAs derived from oocytes exhibit with the characteristics of developmental stage-specific expression . Xin et al. revealed that the depletion of circLDLR inhibits the expression of CYP19A1, thereby reducing the secretion of estrogen in polycystic ovary syndrome . Moreover, Jia et al. found that overexpression of circEGFR increases the production of estradiol, while knockdown of circEGFR enhances the production of progesterone in mice . These observations suggested the potential importance and significance of circRNAs in the ovary, but the information of ovarian circRNAs remains little during the pubertal transition.
Collectively, in order to obtain more insights on the roles of circRNAs in ovaries during pubertal transition, the genome-wide analysis of circRNAs in ovaries across pre-, in- and post-puberty was performed by RNA-sEq. Then, the expression changes of circRNAs were explored as well as the stage-specific and tissue-specific circRNAs, and the potentially pubertal circRNAs were detected in this study. On the other hand, our study may provide a novel theoretical reference for the regulation of pubertal female by circRNAs.
Identification of ovary-derived circRNAs during the onset of puberty
To obtain a reliable result, CIRI2 and find_circ software were intersected to identify circRNAs. A total of 972 circRNAs candidates were identified in the pubertal transition of pubertal ovaries (Additional file 1) (Fig. 1a). Thereinto, 347, 293 and 827 circRNAs were identified in pre-, in- and post-puberty, respectively (Fig. 1b). These circRNAs were more distributed in Sus scrofa chromosome 1 (Fig. 1c). The expressions of circRNAs were the lowest in the post-puberty (Fig. 1d), compared to pre- and in-puberty. Meanwhile, the average genome distance of all circRNAs was 13,453 bp, and circRNAs shorter than 50,000 bp were accounted for 96.7 % (Fig. 1e). Besides, the average length of all circRNAs was 558 bp, and circRNAs with the length of 200–500 bp were accounted for 53.70 % (Fig. 1f). Notably, in the pre-puberty, exonic, intronic and intergenic circRNA occupied 93.66 %, 3.46 and 2.88 %, respectively; in the in-puberty, exonic, intronic and intergenic circRNA occupied 94.20 %, 4.10 and 1.70 %, respectively; in the post-puberty, exonic, intronic and intergenic circRNA occupied 95.41 %, 2.54 and 2.05 %, respectively (Fig. 1 g). Additionally, most circRNAs were made up of two and three exons. Specifically, circRNAs were made up of two exons occupied 26.80 %, 22.53 and 33.98 % in pre-, in- and post-puberty, respectively; circRNAs were made up of three exons occupied 34.87 %, 33.45 and 31.56 % in pre-, in- and post-puberty, respectively (Fig. 1 g). Notably, 722 genes were identified as the parental genes of these 972 circRNA. 922 exonic circRNAs were derived from 699 genes, and 23 intronic circRNAs were derived from 23 genes (Fig. 1 h). Taken together, 972 circRNAs were identified during onset of puberty in the ovaries of gilts.
Key pathways of cirRNAs in pubertal transition and circRNAs in pubertal genes
To further investigate the biological functions of circRNAs involved in pubertal ovary, the parental genes of all 972 circRNAs in puberty were used to the analysis of KEGG pathway (Additional file 2). Obviously, the functional pathways were significantly over-represented in pubertal ovaries, including steroid biosynthesis, autophagy-animal, MAPK signaling pathway, progesterone-mediated oocyte maturation and ras signaling pathway (Fig. 2a). In the steroid biosynthesis signaling pathway, “circ 14:14983402–14992789”, “14:14984501–14992789” and “14:14989270–15001827” derived from FDFT1 gene were uniquely expressed in the post-puberty (Fig. 2b). In the autophagy-animal signaling pathway, “circ 1:190648253–190654424” derived from HIF1A gene was uniquely expressed in the in-puberty. In the MAPK signaling pathway, “circ 15:25128581–25143159” and “circ 15:25140315–25143159” derived from MAP3K2 gene was uniquely expressed in the post-puberty. In the progesterone-mediated oocyte maturation, “circ 16:50437615–50478728” derived from CPEB4 gene was uniquely expressed in the post-puberty. In the ras signaling pathway, “circ 12:43673069–43691261” derived from NF1 gene was expressed in the in-puberty and post-puberty. Details of these circRNAs were shown in Additional file 3. Moreover, in order to further explore the circRNAs in pubertal genes, 10 puberty-related genes (ESR1, JAK2, NF1, ARNT, IGF1, KISS1, Gpr54, NKB, Mkrn3, GnRH) were selected and explored by manual reviewing the literature and databases, and 5 circRNAs derived from 4 pubertal genes (ESR1, JAK2, NF1, ARNT) were lastly found in the ovary of pubertal transition (Additional file 4). “circ 1:14373232–14374308” (uniquely expressed in the pre-puberty) and “circ 1:14416335–14457143” (expressed in the pre-, in- and post-puberty) were derived from ESR1; “circ 1:216914275–216951002” (uniquely expressed in the post-puberty) was derived from JAK2; “circ 12:43673069–43691261” (expressed in the in- and post-puberty) was derived from NF1; “circ 4:98369520–98372553” (uniquely expressed in the post-puberty) was derived from ARNT. Apart from “circ 1:14416335–14457143” and “circ 12:43673069–43691261”, other 3 circRNAs derived from 4 pubertal genes showed stage-specific expressions.
AS of circRNAs in gilts’ ovaries during puberty
The formation of circRNAs is dependent on AS . In order to further explore the AS events involved in circRNAs, we identified the splicing events in circRNAs. Compared with other events, A3SS events were the most splicing pattern in ovaries of gilts in puberty (Welch two-sample t-test, P < 0.05) (Fig. 3a). Strikingly, in the four types of AS events, IR showed more extreme post-pubertal tendency (Welch two-sample t-test, P < 0.05) (Fig. 3b). In other words, there were difference in IR event between pre-puberty and post-puberty. Furthermore, “circ 6:166505226–166505778” existed two isoforms, and its parental gene (PTCH2) was reported to regulate the follicle development  (Additional file 5). The isoforms spliced by A3SS events exist in pre- and in-puberty, but do not exist in post-puberty; the isoforms spliced by IR events exclusively exist in pre-puberty (Additional file 5) (Fig. 3c). The results above showed that the AS events might play a crucial role in formation of ovarian circRNAs during puberty.
Stage-specific and ovary-specific circRNAs in the pubertal transition
To further explore the stage-specific circRNAs during these pubertal stages, we investigated the expression of circRNAs in pre-, in- and post-puberty stages. Respectively, 72, 50 and 509 of circRNAs were uniquely expression in pre-, in- and post-puberty stages and considered to be stage-specific circRNAs (Fig. 4a). The parental genes of these stage-specific circRNAs were enriched in MAPK signaling pathway, progesterone-mediated oocyte maturation, oocyte meiosis and GnRH signaling pathway in post-puberty (Additional file 6) (Fig. 4b). Moreover, 154 circRNAs were expressed in all stages and considered as the co-existed circRNAs (Figs. 1b and 4a). Besides, some specific circRNAs and co-existed circRNAs were derived from the same gene. For instance, “circ 1:100589850–100603238” (uniquely expressed in the post-puberty) and “circ 1:100589850–100613174” (no uniquely expressed in any puberty) were derived from SMAD4 (Additional file 7). In order to further investigate the specific circRNAs in the ovaries, 964 known circRNAs that were found in nine tissues (brain, heart, kidney, liver, lung, skeletal muscle, spleen, testis, and retina) were excluded, leaving 8 circRNAs as being putative ovary-specific circRNAs which were only expressed in the ovary (Additional file 8). Subsequently, we found that the length of ovary-specific circRNAs were shorter than known circRNAs (Fig. 4c). Strikingly, apart from “circ 10:22806071–22812591”, other 7 putative ovary-specific circRNAs were exclusively expressed in the post-puberty (Fig. 4d). Besides, “circ 10:22806071–22812591” and “circ 10:22742781–22748221” were both derived from NR5A2 (Additional file 8). To sum up, the results showed that 64.92 % (631/972) of ovarian circRNAs expressed in the stage-specific means during pubertal transition, and 8 circRNAs were the ovary-specific circRNAs.
Potentially regulated network of differentially expressed circRNAs
Subsequently, 154 co-existed circRNAs were used to analysis differential expression between pair-wise comparison of three stages (Figs. 1b and 4a). In total, 10 circRNAs were identified as differentially expressed circRNAs (Additional file 9), of which 7 up-regulated circRNAs and 3 down-regulated circRNAs were identified in the pre- vs. post-puberty group; 2 up-regulated circRNA were identified in the in- vs. post-puberty group (Fig. 5a). To further explore the possible functions of the differentially expressed circRNAs, we tried to predict the miRNA binding sites of these circRNAs and explore the possible relationship between differentially expressed circRNAs and differentially expressed genes. The miRNAs with the top 5 highest score in miRanda-based circRNAs match were selected as potential miRNA targets and listed in Additional file 9 (see Methods for further details). We found that the differentially expressed circRNAs might interact with many of miRNAs, or might indirectly interact with differentially expressed genes (Fig. 5b-c). Noticeably, we also highlighted FSTL4, GAS2, AIG1, GNG2, FSHR and SPTA1 genes, which were associated with folliculogenesis or hormone [42,43,44,45,46,47]. For instance, “circ 9:131264261–131268491”, which was down-regulated in the pre- vs. post-puberty groups, might interact with FSTL4 via ssc-miR-320, might interact with GAS2 via ssc-miR-582-5p, and might interact with FSHR via ssc-miR-493-3p (Additional files 9 and 10). Taken together, the circRNA-miRNA-gene networks were explored for 10 differentially expressed circRNAs.
Validation of circRNAs
To verify the accuracy of our data, the divergent RT-PCR and sanger sequencing were utilized to identify the authenticity of circRNA, the head-to-tail splice junctions, as well as the expressions of circRNAs. The head-to-tail splice junctions of 5 circRNAs were determined by sanger sequencing, which proved that the circRNAs were circRNAs (Fig. 6a). Furthermore, 7 circRNAs of 10 differentially expressed circRNAs were selected for further investigation. The “circ 7:65198472–65198799” (Fig. 6b), “circ 16:44513270–44525240” (Fig. 6c), “circ 2:95657661–95677798” (Fig. 6d), “circ 9:131264261–131268491” (Fig. 6e), “circ 15:76166995|76177857” (Fig. 6f), “circ 1:121217326|121230821” (Fig. 6 g), and “circ 13:82256891|82274534” (Fig. 6 h) were significantly differentially expressed, which are in line with the RNA-sEq. Besides, the expression of “circ 2:151068704–151069641”, which was detected insignificantly differentially expressed, was insignificantly changed (Fig. 6i). Finally, results showed that the expressions of 8 selected circRNAs verified by divergent RT-qPCR were consistent with the trend of RNA-seq data (Additional file 9). Our results indicated that the existence of differentially expressed circRNAs, which further shows that our analysis is reliable.
In mammals, the onset of puberty is highly implicated in reproduction, and the abnormal puberty can cause various diseases. For instance, the precocious girls have twice the risk of asthma in adulthood than normally pubertal girls . We have previously demonstrated that pituitary-specific circRNAs are related to reproduction-associated signaling pathways in pubertal gilts . The ovary, as the important member of HPO axis, has been reported to guide female into puberty. Zhao et al. showed that circ_0023942 might inhibit the proliferation of human ovarian granulosa cells by regulating the expression of CDK-4 . Therefore, it is necessary to profile the expressions and changes of circRNAs in ovaries during pubertal transition. In this study, the circRNAs we obtained are widely distribution on 1 chromosome, which is consistent with previous literatures [48, 50]. It is worth noting that previous report has shown that circRNAs are divided into three categories, of which exonic circRNAs account for the majority . Consistently, in this study, we demonstrated that exonic circRNAs were accounted for approximately 94 %.
Moreover, we uncovered several circRNAs in some of the key pathways associated with the onset of puberty, including steroid biosynthesis, autophagy-animal, MAPK signaling pathway, progesterone-mediated oocyte maturation and ras signaling pathway. These pathways have all been reported to regulate the development of ovary and onset of puberty [41, 51,52,53]. Notably, for autophagy-animal pathway, “circ 1:190648253–190654424” was exclusively expressed in the in-puberty, and its parental gene HIF1A has been reported to have the highest expression in the early luteal phase . Therefore, circRNA “circ 1:190648253–190654424” might play a crucial role in autophagy in puberty. Furthermore, we found that 5 circRNAs derived from ESR1, JAK2, NF1 and ARNT, which are associated with the onset of puberty [55,56,57,58]. Previous studies have shown that circRNAs can promote the transcription of their parental genes [24, 26]. It is possible therefore that these 5 circRNAs might contribute to the onset of puberty by promoting transcription of ESR1, JAK2, NF1 and ARNT.
Interestingly, AS not only formed the linear mRNA, but also formed circRNAs, which caused the diversity of circRNAs . For instance, compared with linear mRNAs, exons in Exon Skipping (ES) events were more likely included in circRNAs , and some introns are retained in circRNAs through alternative splicing . In this study, we found that the AS events in circRNAs were consistent with previous reports, of which A3SS and ES events occurred in circRNAs in abundance. Furthermore, circRNA “6:166505226–166505778” from PTCH2 gene was spliced by IR and A3SS events in pre-puberty, while by A3SS in in-puberty (Fig. 3). It has been reported that PTCH2 is related to the functions of granulosa cells in ovaries of mammals [59, 60]. It is possible that AS events were differentially expressed in three pubertal stages.
Previous studies have shown that circRNAs have stage-specific and tissue-specific expression patterns ; Xia et al. comprehensively analysis circRNAs in human and mouse to establish a tissue-specific circRNAs database ; G. Maass et al. have shown that circRNAs exhibit a high degree of tissue-specificity in clinically relevant human tissues . According to our analysis of stage-specific circRNAs during the pubertal transition in gilts, we found that 509 of circRNAs specifically existed in post-puberty. Importantly, the parental genes of post-pubertal specific circRNAs were associated with the function of growth, development and reproduction. These results indicated that circRNAs might influence parental genes to regulate the puberty in female of mammals. Moreover, the parental genes of 8 circRNAs that were found uniquely expressed in ovaries might regulate the growth and development of ovary. For example, the parental gene of “circ 10:22806071–22812591” was NR5A2; the parental gene of “circ 6:161225244–161281707” was FAF1; the parental gene of “6:158955138–158958815” was LRP8. Guo et al. have shown that enhancing the expression of NR5A2 stimulates the progesterone synthesis in granulosa cell of porcine . There is evidence that FAF1, Fas-associated protein factor 1, is expressed in the cytoplasm of growing oocyte of the ovary . Yang et al. have provided that LRP8 is the progestogenic-associated gene . These results suggest that the tissue-specific circRNAs might regulate the transcriptions of their parental genes to involve in the onset of puberty.
Besides, many of target genes regulated by circRNAs-miRNAs network are the important puberty related genes and differentially expressed in the pubertal ovaries (Fig. 5b-c). For example, GAS2 is a regulator in the ovary during folliculogenesis and oocyte cyst breakdown ; AIG1 is induced by androgen and expressed at high levels in the ovaries ; SPTA1 might be related to oocyte maturation ; FSHR is expressed in granulosa cells that regulates proliferation of granulosa cell, maturation of follicular and production estrogen ; CERS3 is related to androgen production . Studies have shown that circRNAs might act as a miRNA sponge to co-regulate the expression of related genes with miRNA . For instance, circPUM1 acted as a molecular sponge of mir-760 and reduced the expression of mir-760, which in turn affected polycystic ovary syndrome . Moreover, Meng et al. characterized circRNA in healthy antral and atretic antral follicles, and found that abnormal circRNAs might play a crucial role in antral follicular atresia in pig . Another study showed that circINHA promoted the proliferation of granulosa cells and inhibited the apoptosis of granulosa cells by acting as a miR-10a-5p sponge in pig . In this study, we found that the differentially expressed genes interacted with differentially expressed circRNAs. For instance, in pre- vs. post-puberty group, the down-regulated “circ 9:131264261–131268491” interacted with the down-regulated FSHR, indicating that this circRNAs might as the miRNAs sponge to affect the expression of FSHR. These results suggested that circRNAs might affect the expression of genes by acting as sponge of miRNA, which in turn influence the onset of puberty in mammals. On the other hand, further investigations are needed to identify the functions of circRNAs.
In conclusion, we described the profiles of ovarian 972 circRNAs during pubertal transition in gilts, and these circRNAs were mostly enriched in steroid biosynthesis, autophagy-animal, MAPK signaling pathway, progesterone-mediated oocyte maturation and ras signaling pathway. 631 circRNAs were stage-specific, 8 circRNAs were tissue-specific, and 10 circRNAs were differentially expressed across pre-, in- and post-pubertal ovarian. Besides, 5 circRNAs were derived from four pubertal genes ESR1, JAK2, NF1 and ARNT. The isoforms of circRNAs spliced by IR were more likely to take place in post-puberty, and circRNA-miRNA-gene networks were explored for 10 differentially expressed circRNAs. Furthermore, several circRNAs were validated by the divergent RT-qPCR. These results suggest that circRNAs might play a crucial role in mammalian ovaries during onset of puberty, and further research should be undertaken to investigate the molecular mechanism of ovarian circRNAs in pubertal onset of mammals.
Preparation of samples
The gilts were purchased from Baishi Pig Farm, Zhongshan City, Guangdong Province, China. Three groups of Landrace × Yorkshire crossbred gilts were used: three gilts were served as pre-puberty gilts which were 160 days in age without any pubertal signs (weight = 81.38 ± 2.40 kg, age = 162 ± 3 d, no reddening, no swelling of the vulva, no standing reflex); three gilts were designated as the in-pubertal gilts which exhibited first pubertal signs (weight = 110.00 ± 2.00 kg, age = 212 ± 14 d, reddening, swelling of the vulva, standing reflex); three gilts were selected as the post-pubertal gilts which were 14 days beyond the pubertal phase (weight = 122.82 ± 9.11 kg, age = 216 ± 17 d). After euthanasia, the gilts’ ovaries (diameter ≥ 5 mm) were removed immediately to the liquid nitrogen, then stored at − 80 °C until further use.
RNA sequencing and data processing
Pre-, in-, and post-pubertal gilts’ ovaries were extracted the total RNA using the Trizol agent (Invitrogen, Carlsbad, CA, USA), the total RNA quality was then measured using the Agilent Bioanalyzer 2100 system (Agilent, Palo Alto, California, USA). Filter RNA samples with RNA integrity values greater than 7.0, and remove rRNA from qualified total RNA using Epicentre Ribo-zero rRNA removal kits (Epicentre, Madison, WI, USA). We then used rRNA-depleted RNA to form a double-stranded cDNA with the mRNA-Seq sample preparation kit (Illumina, SanDiego, USA). Later, the cDNA library of each sample was sequenced using the HiSeq 2500 sequencer based on the manufacturer’s instructions and further produced 150 bp paired-end reads. These raw reads were used the Cutadapt software to remove the low-quality reads and the 3’ adaptor-trimming for the quality control . The clean data that after quality controlled was then mapped with two software, which were BWA and bowtie2 software , and the reference genome used Sus scrofa11.1.
CircRNA identification and data analysis
Reports showed that CIRI2 has high sensitivity and low FDR, thus we used CIRI2 identified circRNA after BWA , as well as using find_circ  to identify circRNA after bowtie2 to minimize the number of false positives. The two programs look for potential circRNAs based on genomic comparisons. We screened circRNA with at least 2 unique junction reads as candidates, removed circRNAs with unclear break point, and filtered circRNA with a length greater than 100 kb (genome length, which defined as the distance from the first exon to the last exon in the circRNA). We eventually identified candidate circRNA in the gilts during pre-, in- and post-puberty. Thereinto, CIRI2 generated 1-base coordinates, but find_circ generated 0-base coordinates, thus we converted the two coordinates into a consistent 1-base for later analysis. Subsequently, we set the circRNA detected only in one pubertal stage as a stage-specific circRNA. In addition, the selection criteria for tissular specificity was as follows: the circRNAs identified in this study were matched with the known circRNAs in pigs by starting and ending the genome locations of circRNAs, and the new circRNAs were considered as the presumed tissue specific circRNAs. The known circRNAs were downloaded from circAtlas 2.0 (namely, the circRNAs database in vertebrates) which were included circRNAs of nine tissues (brain, retina, heart, kidney, liver, lung, skeletal muscle, spleen, and testis) . In addition, the alternative splicing events of circRNAs were determined by the CIRI-AS module , which classified the alternative splicing events into four forms: A3SS, A5SS, ES, and IR. The criteria for differential alternative splicing was as follows: PSI as the expression value, was subjected to the difference significance test (t-test) between any two pubertal pig groups. In this study, the EBSeq package was used to calculate the expression levels of circRNAs , which was quantified in RPM using the number of splicing junctions. The criteria for differentially expressed circRNAs was log2-fold_change- ≥ 1, adjusted p (p.adj) < 0.05. In addition, the value of any two pubertal pig groups was subjected to the difference significance test (Welch two-sample t-test) to analyze the significant differences.
Prediction of miRNA target and circRNA-miRNA-mRNA network construction
The interaction of circRNA-miRNA-gene was predicted by miRanda software  with a miRanda match score ≥ 175. The specific method is as follows: all the miRNAs sequence of Sus scrofa was obtained from miRBase database (http://www.mirbase.org/), all the circRNAs sequence was obtained using Bedtools, and the match score of miRNA and circRNA was scored utilizing miRanda, miRNAs with top 5 matching scores were eventually predicted. Furthermore, Bedtools  was used to extract the differentially up-regulated and down-regulated mRNA sequences between any two pubertal pig groups (p.adj < 0.05, |log2FC| ≥3 or ≤ − 3), respectively. Subsequently, miRanda software was used to predict the target genes of miRNA according to these sequences. Finally, the interoperability between circRNA-miRNA-gene was then described by the cytoscape software .
The parental gene of circRNA was analyzed by the KOBAS online software (http://kobas.cbi.pku.edu.cn/) with its GO function enrichment and KEGG pathway analyses . The hypergeometric test significance threshold P < 0.05 was considered for genes to indicate significant enrichment.
Validation of CircRNA
Back-spliced junction (BSJ) was a region consisting of canonical 5’ splice site sequence connected to upstream 3’ splice site sequence. The reliability of circRNA is usually verified by divergent primer flanking the BSJ utilizing RT and quantitative PCR (RT-qPCR) assays . The divergent primers of 10 circRNAs were designed to verify the accuracy of the RNA-sEq. Firstly, in accordance with the operator’s procedures of Taq PCR MasterMix (Tiangen, China), the cDNA template PCR were amplified by 35 cycles at 95 °C (30 s), 60 °C (30 s) and 72 °C (20 s), and the PCR product was visualized using a 2 % GelRed-stained agar glycogel. Furthermore, the sanger sequencing was further performed to directly examine the PCR product. In accordance with the manufacturer’s protocol, PrimeScript RT Reagent Kit (TaKaRa, Osaka, Japan) in a Mx3005P real-time PCR System (Stratagene, La Jolla, CA, USA) was used to qPCR, and the PCR standard procedure was denaturation 94 °C (5 min), 40 cycles at 94 °C (10 s), 52 to 62 °C (15 s), and 72 °C (30 s). As an internal reference, GAPDH normalized expression of circRNAs. We used the 2−ΔΔCt method to analyze the RT-qPCR data, and used the Student’s t test to assess differences in means of any two pubertal pig groups, the screening criteria for statistically significant were adjusted p < 0.05 (t-test). All of PCR primer were divergent design by the head-to-tail exon of circRNAs and listed in Additional file 11.
Availability of data and materials
The datasets used in this study have been submitted to the European Nucleotide Archive under accession number PRJEB39730 (https://www.ebi.ac.uk/ena/browser/view/PRJEB39730).
Alternative 3´ splice site
Alternative 5´ splice site
Polymerase chain reaction
Quantitative real-time PCR
Euling SY, Herman-Giddens ME, Lee PA, Selevan SG, Juul A, Sorensen TI, Dunkel L, Himes JH, Teilmann G, Swan SH. Examination of US puberty-timing data from 1940 to 1994 for secular trends: panel findings. Pediatrics. 2008;121(Suppl 3):S172–91.
Al-Sahab B, Hamadeh MJ, Ardern CI, Tamim H. Early Menarche Predicts Incidence of Asthma in Early Adulthood. Am J Epidemiol. 2011;173(1):64–70.
Hankin BL, Badanes LS, Abela JRZ, Watamura SE. Hypothalamic-Pituitary-Adrenal Axis Dysregulation in Dysphoric Children and Adolescents: Cortisol Reactivity to Psychosocial Stress from Preschool Through Middle Adolescence. Biol Psychiat. 2010;68(5):484–90.
Abreu AP, Kaiser UB. Pubertal development and regulation. Lancet Diabetes Endocrinol. 2016;4(3):254–64.
Bodicoat DH, Schoemaker MJ, Jones ME, McFadden E, Griffin J, Ashworth A, Swerdlow AJ: Timing of pubertal stages and breast cancer risk: the Breakthrough Generations Study (vol 16, R18, 2014). Breast Cancer Res 2020, 22(191):1–8.
Cesario SK, Hughes LA. Precocious puberty: a comprehensive review of literature. J Obstet Gynecol Neonatal Nurs. 2007;36(3):263–74.
Wildt L, Hausler A, Marshall G, Hutchison JS, Plant TM, Belchetz PE, Knobil E. Frequency and amplitude of gonadotropin-releasing hormone stimulation and gonadotropin secretion in the rhesus monkey. Endocrinology. 1981;109(2):376–85.
Watanabe G, Terasawa E. In vivo release of luteinizing hormone releasing hormone increases with puberty in the female rhesus monkey. Endocrinology. 1989;125(1):92–9.
Chongthammakun S, Terasawa E. Negative feedback effects of estrogen on luteinizing hormone-releasing hormone release occur in pubertal, but not prepubertal, ovariectomized female rhesus monkeys. Endocrinology. 1993;132(2):735–43.
Pepe GJ, Lynch TJ, Albrecht ED. Regulation of baboon fetal ovarian development by placental estrogen: onset of puberty is delayed in offspring deprived of estrogen in utero. Biol Reprod. 2013;89(6):132.
Hu K, Sun W, Li Y, Zhang B, Zhang M, Guo C, Chang H, Wang X. Study on the Mechanism of Sarsasapogenin in Treating Precocious Puberty by Regulating the HPG Axis. Evid Based Complement Alternat Med. 2020;2020:1978043.
Kendirci HN, Agladioglu SY, Onder A, Bas VN, Cetinkaya S, Aycan Z. Effects of GnRH analogue treatment on anterior pituitary hormones in children with central precocious puberty. J Pediatr Endocrinol Metab. 2015;28(9–10):1145–51.
Yang R, Wang Y, Zhang L, Zhao Z, Zhao J, Peng S. Prepubertal exposure to an oestrogenic mycotoxin zearalenone induces central precocious puberty in immature female rats through the mechanism of premature activation of hypothalamic kisspeptin-GPR54 signaling. Mol Cell Endocrinol. 2016;437(C):62–74.
Baudry M, Bi X, Aguirre C. Progesterone-estrogen interactions in synaptic plasticity and neuroprotection. Neuorscience. 2013;239:280–94.
Mlodawska W, Grzesiak M, Kochan J, Nowak A. Intrafollicular level of steroid hormones and the expression of androgen receptor in the equine ovary at puberty. Theriogenology. 2018;121:13–20.
Ernst E, Kjaersgaard M, Birkebaek NH, Clausen N, Andersen CY. Case report: stimulation of puberty in a girl with chemo- and radiation therapy induced ovarian failure by transplantation of a small part of her frozen/thawed ovarian tissue. Eur J Cancer. 2013;49(4):911–4.
Khristi V, Chakravarthi VP, Singh P, Ghosh S, Pramanik A, Ratri A, Borosha S, Roby KF, Wolfe MW, Rumi M. ESR2 regulates granulosa cell genes essential for follicle maturation and ovulation. Mol Cell Endocrinol. 2018;474:214–26.
Stagi S, di Tommaso M, Scalini P, Lapi E, Losi S, Bencini E, Masoni F, Dosa L, Becciani S, de Martino M. Triple X syndrome and puberty: focus on the hypothalamus-hypophysis-gonad axis. Fertil Steril. 2016;105(6):1547–53.
Ro S, Song R, Park C, Zheng H, Sanders KM, Yan W: Cloning and expression profiling of small RNAs expressed in the mouse ovary. RNA 2007, 13(12):2366–2380.
Pan Z, Zhang J, Li Q, Li Y, Shi F, Xie Z, Liu H. Current advances in epigenetic modification and alteration during mammalian ovarian folliculogenesis. J Genet Genomics. 2012;39(3):111–23.
Geach T. Neuroendocrinology: microRNAs regulate puberty timing. Nat Rev Endocrinol. 2016;12(7):372.
Stamou M, Ng SY, Brand H, Wang H, Plummer L, Best L, Havlicek S, Hibberd M, Khor CC, Gusella J et al: A Balanced Translocation in Kallmann Syndrome Implicates a Long Noncoding RNA, RMST, as a GnRH Neuronal Regulator. J Clin Endocrinol Metab. 2020;105(3):e231–44.
Zhang XO, Wang HB, Zhang Y, Lu X, Chen LL, Yang L: Complementary sequence-mediated exon circularization. Cell. 2014, 159(1):134–147.
Guo JU, Agarwal V, Guo H, Bartel DP: Expanded identification and characterization of mammalian circular RNAs. Genome Biol. 2014, 15 7):409.
Li Z, Huang C, Bao C, Chen L, Lin M, Wang X, Zhong G, Yu B, Hu W, Dai L et al: Exon-intron circular RNAs regulate transcription in the nucleus. Nat Struct Mol Biol. 2015; 22 (3):256–264.
Zhang Y, Zhang XO, Chen T, Xiang JF, Yin QF, Xing YH, Zhu S, Yang L, Chen LL: Circular intronic long noncoding RNAs. Mol Cell. 2013, 51(6):792–806.
Xia S, Feng J, Lei L, Hu J, Xia L, Wang J, Xiang Y, Liu L, Zhong S, Han L et al: Comprehensive characterization of tissue-specific circular RNAs in the human and mouse genomes. Brief Bioinform. 2017; 18 (6):984–992.
Salzman J, Chen RE, Olsen MN, Wang PL, Brown PO: Cell-type specific features of circular RNA expression. PLOS Genet. 2013, 9(9):e1003777.
Kristensen LS, Andersen MS, Stagsted L, Ebbesen KK, Hansen TB, Kjems J: The biogenesis, biology and characterization of circular RNAs. Nat Rev Genet. 2019, 20(11):675–691.
Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, Maier L, Mackowiak SD, Gregersen LH, Munschauer M et al: Circular RNAs are a large class of animal RNAs with regulatory potency. Nature 2013; 495 (7441):333–338.
Hansen TB, Jensen TI, Clausen BH, Bramsen JB, Finsen B, Damgaard CK, Kjems J: Natural RNA circles function as efficient microRNA sponges. Nature. 2013; 495 (7441):384–388.
Wilusz JE: Circular RNAs: Unexpected outputs of many protein-coding genes. RNA Biol. 2017; 14 (8):1007–1017.
Huang Z, Cao Y, Zhou M, Qi X, Fu B, Mou Y, Wu G, Xie J, Zhao J, Xiong W: Hsa_circ_0005519 increases IL-13/IL-6 by regulating hsa-let-7a-5p in CD4(+) T cells to affect asthma. Clin Exp Allergy. 2019; 49 (8):1116–1127.
Jia Y, Li X, Nan A, Zhang N, Chen L, Zhou H, Zhang H, Qiu M, Zhu J, Ling Y et al: Circular RNA 406961 interacts with ILF2 to regulate PM2.5-induced inflammatory responses in human bronchial epithelial cells via activation of STAT3/JNK pathways. Environ Int. 2020l 141:105755.
Xu H, Sun Y, You B, Huang CP, Ye D, Chang C: Androgen receptor reverses the oncometabolite R-2-hydroxyglutarate-induced prostate cancer cell invasion via suppressing the circRNA-51217/miRNA-646/TGFbeta1/p-Smad2/3 signaling. Cancer Lett. 2020; 472:151–164.
Xie S, Li M, Chen Y, Liu Y, Ma L, Sun X, Sun Y, Gao R, Huang T: Identification of circular RNAs in the ovarian follicles of Meishan and Duroc sows during the follicular phase. J Ovarian Res. 2020, 13 (1):104.
Cao Z, Gao D, Xu T, Zhang L, Tong X, Zhang D, Wang Y, Ning W, Qi X, Ma Y et al: Circular RNA profiling in the oocyte and cumulus cells reveals that circARMC4 is essential for porcine oocyte maturation. Aging (Albany NY). 2019; 11(18):8015–8034.
Huang X, Wu B, Chen M, Hong L, Kong P, Wei Z, Teng X: Depletion of exosomal circLDLR in follicle fluid derepresses miR-1294 function and inhibits estradiol production via CYP19A1 in polycystic ovary syndrome. Aging (Albany NY). 2020; 12 (15):15414–15435.
Jia W, Xu B, Wu J: Circular RNA expression profiles of mouse ovaries during postnatal development and the function of circular RNA epidermal growth factor receptor in granulosa cells. Metabolism. 2018, 85:192–204.
Gao Y, Wang J, Zheng Y, Zhang J, Chen S, Zhao F: Comprehensive identification of internal structure and alternative splicing events in circular RNAs. Nat Commun. 2016, 7:12060.
Lutz LB, Cole LM, Gupta MK, Kwist KW, Auchus RJ, Hammes SR: Evidence that androgens are the primary steroids produced by Xenopus laevis ovaries and may signal through the classical androgen receptor to promote oocyte maturation. Proc Natl Acad Sci U S A 2001, 98(24):13728–13733.
Wang Y, Cheng T, Lu M, Mu Y, Li B, Li X, Zhan X: TMT-based quantitative proteomics revealed follicle-stimulating hormone (FSH)-related molecular characterizations for potentially prognostic assessment and personalized treatment of FSH-positive non-functional pituitary adenomas. EPMA J 2019, 10(4):395–414.
York JP, Ren YA, Zeng J, Bin Z, Wang F, Chen R, Liu J, Xia X, Zhang P: Growth Arrest Specific 2 (GAS2) is a Critical Mediator of Germ Cell Cyst Breakdown and Folliculogenesis in Mice. Sci Rep 2016, 6:34956.
Seo J, Kim J, Kim M: Cloning of androgen-inducible gene 1 (AIG1) from human dermal papilla cells. Mol Cells 2001, 11(1):35–40.
Liang W, Ji L, Zhang Y, Zhen Y, Zhang Q, Xu X, Liu B: Transcriptome Differences in Porcine Alveolar Macrophages from Tongcheng and Large White Pigs in Response to Highly Pathogenic Porcine Reproductive and Respiratory Syndrome Virus (PRRSV) Infection. Int J Mol Sci 2017;18(7):1475.
Karakaya C, Guzeloglu-Kayisli O, Hobbs RJ, Gerasimova T, Uyar A, Erdem M, Oktem M, Erdem A, Gumuslu S, Ercan D et al: Follicle-stimulating hormone receptor (FSHR) alternative skipping of exon 2 or 3 affects ovarian response to FSH. Mol Hum Reprod 2014, 20(7):630–643.
Coyral-Castel S, Brisard D, Touze JL, Dupont M, Rame C, Uzbekova S, Dupont J: Analysis of in vivo oocyte maturation, in vitro embryo development and gene expression in cumulus cells of dairy cows and heifers selected for one fertility quantitative trait loci (QTL) located on BTA3. Theriogenology 2012, 77(9):1822–1833.
Chen Z, Pan X, Kong Y, Jiang Y, Zhong Y, Zhang H, Zhang Z, Yuan X, Li J: Pituitary-Derived Circular RNAs Expression and Regulatory Network Prediction During the Onset of Puberty in Landrace x Yorkshire Crossbred Pigs. Front Genet 2020, 11:135.
Zhao C, Zhou Y, Shen X, Gong M, Lu Y, Fang C, Chen J, Ju R: Circular RNA expression profiling in the fetal side of placenta from maternal polycystic ovary syndrome and circ_0023942 inhibits the proliferation of human ovarian granulosa cell. Arch Gynecol Obstet 2020, 301(4):963–971.
Tian J, Fu Y, Li Q, Xu Y, Xi X, Zheng Y, Yu L, Wang Z, Yu B, Tian J: Differential Expression and Bioinformatics Analysis of CircRNA in PDGF-BB-Induced Vascular Smooth Muscle Cells. Front Genet 2020, 11:530.
Leng X, Zhou H, Tan Q, Du H, Wu J, Liang X, He S, Wei Q: Integrated metabolomic and transcriptomic analyses suggest that high dietary lipid levels facilitate ovary development through the enhanced arachidonic acid metabolism, cholesterol biosynthesis and steroid hormone synthesis in Chinese sturgeon (Acipenser sinensis). Br J Nutr 2019, 122(11):1230–1241.
Quan C, Wang C, Duan P, Huang W, Chen W, Tang S, Yang K: Bisphenol a induces autophagy and apoptosis concurrently involving the Akt/mTOR pathway in testes of pubertal SD rats. Environ Toxicol 2017, 32(8):1977–1989.
Hiney JK, Srivastava VK, Vaden AD, Hartzoge NL, Dees WL: >Regulation of Kisspeptin Synthesis and Release in the Preoptic/Anterior Hypothalamic Region of Prepubertal Female Rats: Actions of IGF-1 and Alcohol. Alcohol Clin Exp Res 2018, 42(1):61–68.
Berisha B, Schams D, Rodler D, Sinowatz F, Pfaffl MW: Expression pattern of HIF1alpha and vasohibins during follicle maturation and corpus luteum function in the bovine ovary. Reprod Domest Anim 2017, 52(1):130–139.
Kugelberg E: Reproductive endocrinology: ESR1 mutation causes estrogen resistance and puberty delay in women. Nat Rev Endocrinol 2013, 9(10):565.
Wu S, Divall S, Hoffman GE, Le WW, Wagner KU, Wolfe A: Jak2 is necessary for neuroendocrine control of female reproduction. J NeurosciI 2011, 31(1):184–192.
Vossler MR, Yao H, York RD, Pan MG, Rim CS, Stork PJ: cAMP activates MAP kinase and Elk-1 through a B-Raf- and Rap1-dependent pathway. Cell 1997, 89(1):73–82.
Nocillado JN, Elizur A, Avitan A, Carrick F, Levavi-Sivan B: Cytochrome P450 aromatase in grey mullet: cDNA and promoter isolation; brain, pituitary and ovarian expression during puberty. Mol Cell Endocrinol 2007, 263(1–2):65–78.
Wijgerde M, Ooms M, Hoogerbrugge JW, Grootegoed JA: Hedgehog signaling in mouse ovary: Indian hedgehog and desert hedgehog from granulosa cells induce target gene expression in developing theca cells. Endocrinology 2005, 146(8):3558–3566.
Russell MC, Cowan RG, Harman RM, Walker AL, Quirk SM: The hedgehog signaling pathway in the mouse ovary. Biol Reprod 2007, 77(2):226–236.
Maass PG, Glazar P, Memczak S, Dittmar G, Hollfinger I, Schreyer L, Sauer AV, Toka O, Aiuti A, Luft FC et al: A map of human circular RNAs in clinically relevant tissues. J Mol Med (Berl) 2017, 95(11):1179–1189.
Guo R, Chen F, Shi Z: Suppression of Notch Signaling Stimulates Progesterone Synthesis by Enhancing the Expression of NR5A2 and NR2F2 in Porcine Granulosa Cells. Genes (Basel). 2020;11(2):120.
Peng H, Huo J, Gao Y, Chen J, Yu X, Xiao T: Fas-associated protein factor 1 is involved in meiotic resumption in mouse oocytes. J Reprod Dev 2018, 64(2):173–177.
Yang F, Wang M, Zhang B, Xiang W, Zhang K, Chu M, Wang P: Identification of new progestogen-associated networks in mammalian ovulation using bioinformatics. BMC Syst Biol 2018, 12(1):36.
Kranc W, Budna J, Chachula A, Borys S, Bryja A, Rybska M, Ciesiolka S, Sumelka E, Jeseta M, Brussow KP et al: Cell Migration” Is the Ontology Group Differentially Expressed in Porcine Oocytes Before and After In Vitro Maturation: A Microarray Approach. DNA Cell Biol 2017, 36(4):273–282.
Mizutani Y, Kihara A, Igarashi Y: LASS3 (longevity assurance homologue 3) is a mainly testis-specific (dihydro)ceramide synthase with relatively broad substrate specificity. Biochem J 2006, 398(3):531–538.
Deng L, Chen Q, Xie J, Wei W, Hui H: circPUM1 promotes polycystic ovary syndrome progression by sponging to miR-760. Gene. 2020, 754:144903.
Meng L, Teerds K, Tao J, Wei H, Jaklofsky M, Zhao Z, Liang Y, Li L, Wang CC, Zhang S: Characteristics of Circular RNA Expression Profiles of Porcine Granulosa Cells in Healthy and Atretic Antral Follicles. Int J Mol Sci. 2020;21(15):5217.
Guo T, Zhang J, Yao W, Du X, Li Q, Huang L, Ma M, Li Q, Liu H, Pan Z: CircINHA resists granulosa cell apoptosis by upregulating CTGF as a ceRNA of miR-10a-5p in pig ovarian follicles. Biochim Biophys Acta Gene Regul Mech 2019, 1862(10):194420.
Chen C, Khaleel SS, Huang H, Wu CH: Software for pre-processing Illumina next-generation sequencing short read sequences. Source Code Biol Med 2014, 9:8.
Gao Y, Wang J, Zhao F: CIRI: an efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol 2015, 16:4.
Hansen TB, Veno MT, Damgaard CK, Kjems J: Comparison of circular RNA prediction tools. Nucleic Acids Res 2016, 44(6):e58.
Ji P, Wu W, Chen S, Zheng Y, Zhou L, Zhang J, Cheng H, Yan J, Zhang S, Yang P et al: Expanded Expression Landscape and Prioritization of Circular RNAs in Mammals. Cell Rep 2019, 26(12):3444–3460.
Leng N, Dawson JA, Thomson JA, Ruotti V, Rissman AI, Smits BM, Haag JD, Gould MN, Stewart RM, Kendziorski C: EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments. Bioinformatics 2013, 29(8):1035–1043.
John B, Enright AJ, Aravin A, Tuschl T, Sander C, Marks DS: Human MicroRNA targets. >PLOS Biol 2004, 2(11):e363.
Quinlan AR, Hall IM: BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 2010, 26(6):841–842.
Su G, Morris JH, Demchak B, Bader GD: Biological network exploration with Cytoscape 3. Curr Protoc Bioinformatics 2014, 47:8–13.
Wu J, Mao X, Cai T, Luo J, Wei L: KOBAS server: a web-based platform for automated annotation and pathway identification. Nucleic Acids Res 2006, 34(Web Server issue):W720-W724.
We would like to thank the National Supercomputer Center in Guangzhou for its computing platform. Also, we are grateful to the editors and all the reviewers for their insightful comments and constructive suggestions that greatly improved our manuscript.
This research was supported by the China Agriculture Research System of MOF and MARA, the National Natural Science Foundation of China (Reference number: 31902131), the Special Fund for Science and Technology Innovation of Guangdong Province (Reference number: 2018B020203003), the National Natural Science Foundation of Guangdong Province (Reference number: 2019A1515010676), and the Science and Technology Program of Guangzhou (202002030071).
Ethics approval and consent to participate
All the animal experiments performed in this study were approved by South China Agriculture University’s Institutional Animal Care and Use Committee (approval number 2018B116), according to the regulations established by this committee and international standards for animal welfare.
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.
. List of the information of all identified circRNAs.
. List of the KEGG pathways enriched using parental genes of all CircRNAs.
. List of the key pathways related to the timing of puberty in parental genes of circRNAs.
. List of the circRNAs in pubertal genes.
. List of the alternative splicing in circRNAs.
. List of the KEGG pathways enriched using parental genes of stage-specific CircRNAs.
. List of the parental genes that are capable of producing stage-specific and non-specific circRNAs.
. List of the tissue-specific circRNAs.
. List of the differentially regulated circRNAs.
. List of the differentially expressed genes associated with puberty our development of ovary.
. List of primers used for validation.
About this article
Cite this article
Pan, X., Gong, W., He, Y. et al. Ovary-derived circular RNAs profile analysis during the onset of puberty in gilts. BMC Genomics 22, 445 (2021). https://doi.org/10.1186/s12864-021-07786-w
- Alternative splicing