Ovary-derived circular RNAs profile analysis during the onset of puberty in gilts

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.


Background
Puberty is usually defined as the first estrus of mammals [1]. In human, the abnormal puberty has negative effects in diseases such as asthma [2], psychosocial disorder [3], hypogonadism [4] and reproductive system tumors [5,6]. It has been well recognized that the initiation of puberty is mainly driven by the hypothalamus-pituitaryovary (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 [17] and human [18]. 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 [23]. It is reported that most circRNAs are composed of only exonic sequences, while a few cir-cRNAs are composed of the exon-intronic or intronic sequences [24][25][26]. More recently, thousands of cir-cRNAs 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 [35], both of which are closely related to the abnormal puberty, and circRNAs are reported to involve in oocyte maturation and hormone synthesis [36]. For example, Cao et al. showed that circRNAs derived from oocytes exhibit with the characteristics of developmental stagespecific expression [37]. Xin et al. revealed that the depletion of circLDLR inhibits the expression of CYP19A1, thereby reducing the secretion of estrogen in polycystic ovary syndrome [38]. 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 [39]. 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.

AS of circRNAs in gilts' ovaries during puberty
The formation of circRNAs is dependent on AS [40]. In order to further explore the AS events involved in circRNAs, we identified the splicing events in cir-cRNAs. 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 twosample 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 [41] (Additional file 5). The isoforms spliced by A3SS events exist in pre-and in-puberty, Fig. 1 Overview of the identified circRNAs by RNA-seq analyses in ovaries of gilts. a CircRNAs were identified by two algorithms. b The Venn diagram shows the number of unique and common circRNAs in pre-, in-and post-puberty. c Circos plot of the circRNAs distribution in the whole genome of gilts. The outermost circle represents the distribution of the number of circRNAs, the blue circle represents the distribution of expression level of circRNAs, and the red inner circle represents the length of circRNAs. d Expression level of circRNAs in three stages, *p < 0.05. e Genomic dance of all detected circRNAs. f Transcript length of circRNAs. g Proportion of three types and exon number of the circRNAs in three stages. h The upset plot of three types of circRNAs and the corresponding parental genes 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 postpuberty) and "circ 1:100589850-100613174" (no uniquely expressed in any puberty) were derived from SMAD4 (Additional file 7). In order to further Fig. 2 The key signaling pathway of cirRNAs in pubertal transition. a KEGG analysis of all identified circRNAs (*P < 0.05). b Expression level of circRNAs involved in pubertal key pathways in three stages Fig. 3 The alternative splicing (AS) events of circRNAs and the presumed formation of cirRNAs in pubertal transition. a Number of four types of AS events of all detected circRNAs. b Differential IR events with the value of PSI value in three stages, *p < 0.05. c Two isoforms of circRNAs might were derived from PTCH2 by A3SS and IR splicing patterns 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 ovaryspecific circRNAs were exclusively expressed in the postpuberty (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 cir-cRNAs 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 Fig. 4 Analysis results of stage-specific and ovary-specific circRNAs. a Expression level of all circRNAs in three stages. b KEGG analysis of the parental genes of stage-specific circRNAs (P < 0.05). c Length of ovary-specific circRNAs and known circRNAs, *p < 0.05. d Expression level of ovary-specific circRNAs in three stages 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].

Discussion
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 [2]. We have previously demonstrated that pituitaryspecific circRNAs are related to reproduction-associated signaling pathways in pubertal gilts [48]. 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 [49]. Therefore, it is necessary to profile the expressions and changes of cir-cRNAs 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 [50]. Consistently, in this study, we demonstrated that exonic circRNAs were accounted for approximately 94 %. Sanger sequencing and RT-qPCR validation of circRNAs. a sanger sequencing of five circRNAs showed the back-splice junction. b-h seven circRNAs of differential expression and i one circRNA of insignificant difference was randomly selected for RT-qPCR. The primer information was listed in Additional file 11, *p < 0.05, **p < 0.01, *** p < 0.001 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 autophagyanimal 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 [54]. Therefore, cir-cRNA "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 cir-cRNAs can promote the transcription of their parental genes [24,26]. It is possible therefore that these 5 cir-cRNAs 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 [40]. For instance, compared with linear mRNAs, exons in Exon Skipping (ES) events were more likely included in circRNAs [26], and some introns are retained in circRNAs through alternative splicing [25]. 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 prepuberty, 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 [37]; Xia et al. comprehensively analysis circRNAs in human and mouse to establish a tissue-specific circRNAs database [27]; G. Maass et al. have shown that circRNAs exhibit a high degree of tissue-specificity in clinically relevant human tissues [61]. 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 [62]. There is evidence that FAF1, Fas-associated protein factor 1, is expressed in the cytoplasm of growing oocyte of the ovary [63]. Yang et al. have provided that LRP8 is the progestogenic-associated gene [64]. 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 [43]; AIG1 is induced by androgen and expressed at high levels in the ovaries [44]; SPTA1 might be related to oocyte maturation [65]; FSHR is expressed in granulosa cells that regulates proliferation of granulosa cell, maturation of follicular and production estrogen [46]; CERS3 is related to androgen production [66]. Studies have shown that circRNAs might act as a miRNA sponge to co-regulate the expression of related genes with miRNA [29]. 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 [67]. 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 [68]. 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 [69]. 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 downregulated 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.

Conclusions
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 Ribozero rRNA removal kits (Epicentre, Madison, WI, USA). We then used rRNA-depleted RNA to form a doublestranded 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 [70]. The clean data that after quality controlled was then mapped with two software, which were BWA and bowtie2 software [71], 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 [71], as well as using find_circ [72] 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 cir-cRNAs 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 postpuberty. 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) [73]. In addition, the alternative splicing events of circRNAs were determined by the CIRI-AS module [40], 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 [74], 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 [75] 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 mi-Randa, miRNAs with top 5 matching scores were eventually predicted. Furthermore, Bedtools [76] was used to extract the differentially up-regulated and downregulated mRNA sequences between any two pubertal pig groups (p.adj < 0.05, |log 2 FC| ≥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 [77].

Pathway analysis
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 [78]. 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 [29]. 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 Master-Mix (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, GAPD H 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.