Characteristics of Circular RNA Expression Profiles of Porcine Granulosa Cells in Healthy and Atretic Antral Follicles

Circular RNAs (circRNAs) are thought to play essential roles in multiple biological processes, including apoptosis, an important process in antral follicle atresia. We aimed to investigate the potential involvement of circRNAs in granulosa cell apoptosis and thus antral follicle atresia. CircRNA expression profiles were generated from porcine granulosa cells isolated from healthy antral (HA) and atretic antral (AA) follicles. Over 9632 circRNAs were identified, of which 62 circRNAs were differentially expressed (DE-circRNAs). Back-splicing, RNase R resistance, and stability of DE-circRNAs were validated, and miRNA binding sites and related target genes were predicted. Two exonic circRNAs with low false discovery rate (FDR) high fold change, miRNA binding sites, and relevant biological functions—circ_CBFA2T2 and circ_KIF16B—were selected for further characterization. qRT-PCR and linear regression analysis confirmed expression and correlation of the targeted genes—the antioxidant gene GCLC (potential target of circ_CBFA2T2) and the apoptotic gene TP53 (potential target of circ_KIF16B). Increased mRNA content of TP53 in granulosa cells of AA follicles was further confirmed by strong immunostaining of both p53 and its downstream target pleckstrin homology like domain family a member 3 (PHLDA3) in AA follicles compared to negligible staining in granulosa cells of HA follicles. Therefore, we concluded that aberrantly expressed circRNAs presumably play a potential role in antral follicular atresia.


Introduction
Each estrous cycle, a cohort of small antral follicles in, is recruited to enter the pool of growing antral follicles. From this pool, dominant follicles are selected by a process highly dependent on follicle-stimulating hormone (FSH) [1,2]. Follicular development is not a very efficient process as, depending on the species, around 90-99.9% of ovarian follicles degenerate before reaching the ovulatory ovarian tissue, we chose to use the pig as an animal model to investigate the role of circRNAs in antral follicular atresia. For this purpose, the features of circRNAs in porcine granulosa cells derived from healthy and atretic antral follicles were studied using RNA-seq, followed by qRT-PCR in combination with circRNA-miRNA-mRNA network construction and validation.

Histological and Hormone Characteristics between Healthy Antral (HA) and Atretic Antral (AA) Follicles
In antral follicles with a well-vascularized follicular wall and clear follicular fluid, the cells from the granulosa layer were tightly adhered to each other, indicative of a healthy antral follicle. In follicles with a poorly vascularized follicular wall and a rather opaque appearance, the granulosa cell layer was disorganized, and numerous apoptotic cells were present; in some cases, apoptotic granulosa cells were found to be dispersed throughout the antrum, indicative of follicular degeneration. No cleaved caspase 3 (cCASP3) immunostaining was observed in granulosa cells of healthy antral follicles ( Figure S1A), while clear staining was observed in granulosa cells of atretic antral follicles ( Figure S1B). On the contrary, no morphological differences in the oocyte surrounding cumulus cells were observed between HA and AA follicles.
Anti-Müllerian hormone (AMH) and Insulin-like growth factor 1 (IGF1) concentrations ( Figure 1A,B) were both significantly decreased in follicular fluid of AA follicles compared to HA follicles. Concomitantly, estradiol concentrations were significantly lower in the follicular fluid of AA follicles ( Figure 1C). These results further supported the correctness of the macroscopical classification employed to discriminate between HA (well-vascularized) and AA (opaque colored) follicles.
healthy and atretic antral follicles were studied using RNA-seq, followed by qRT-PCR in combination with circRNA-miRNA-mRNA network construction and validation.

Histological and Hormone Characteristics between Healthy Antral (HA) and Atretic Antral (AA) Follicles
In antral follicles with a well-vascularized follicular wall and clear follicular fluid, the cells from the granulosa layer were tightly adhered to each other, indicative of a healthy antral follicle. In follicles with a poorly vascularized follicular wall and a rather opaque appearance, the granulosa cell layer was disorganized, and numerous apoptotic cells were present; in some cases, apoptotic granulosa cells were found to be dispersed throughout the antrum, indicative of follicular degeneration. No cleaved caspase 3 (cCASP3) immunostaining was observed in granulosa cells of healthy antral follicles ( Figure S1A), while clear staining was observed in granulosa cells of atretic antral follicles ( Figure S1B). On the contrary, no morphological differences in the oocyte surrounding cumulus cells were observed between HA and AA follicles.
Anti-Müllerian hormone (AMH) and Insulin-like growth factor 1 (IGF1) concentrations ( Figure  1A, B) were both significantly decreased in follicular fluid of AA follicles compared to HA follicles. Concomitantly, estradiol concentrations were significantly lower in the follicular fluid of AA follicles ( Figure 1C). These results further supported the correctness of the macroscopical classification employed to discriminate between HA (well-vascularized) and AA (opaque colored) follicles.
After scraping the cells from the follicular wall, as described in the Materials and Methods, 98.0 ±1.4% of the cells could be identified as granulosa cells based on positive staining with the follicle stimulating hormone receptor (FSHR) antibody. Furthermore, to ensure that the isolated granulosa cells were not contaminated with theca cells, the mRNA content of cytochrome P450 family 17 subfamily A member 1 (CYP17A1), a specific theca cell marker, was measured by qRT-PCR (see below) [22]. CYP17A1 levels were below the detection limit of the assay, implicating that theca cell contamination was negligible, further confirming the purity of the granulosa cell isolates, which were next used for RNA-seq.  After scraping the cells from the follicular wall, as described in the Materials and Methods, 98.0 ± 1.4% of the cells could be identified as granulosa cells based on positive staining with the follicle stimulating hormone receptor (FSHR) antibody. Furthermore, to ensure that the isolated granulosa cells were not contaminated with theca cells, the mRNA content of cytochrome P450 family 17 subfamily A member 1 (CYP17A1), a specific theca cell marker, was measured by qRT-PCR (see below) [22]. CYP17A1 levels were below the detection limit of the assay, implicating that theca cell contamination was negligible, further confirming the purity of the granulosa cell isolates, which were next used for RNA-seq.

Profiles and Characteristics of circRNAs in Granulosa Cells of Ovarian Antral Follicles
In total, 9362 distinct circRNAs containing at least two unique back-spliced reads ( Figure 2A) were identified based on the circRNAs identification standard [23]. Although 67.3% of circRNAs had average coverage of less than 10 unique back-spliced reads, there were 441 circRNAs highly expressed containing a high-average read count of more than 50, including circ_KIF16B, circ_IL1R1, circ_SLC30A7, and circ_ANGPT1. We identified 2515 overlapping circRNAs and 6847 novel circRNAs that had not been identified so far ( Figure 2B; [24,25]). The 54.36% (5090/9362) of circRNAs identified were derived from exons and generated from individual linear cognate genes. The remainder were annotated to sense, among others, overlapping (2243/9362) and intergenic regions (1195/9362). All exonic circRNAs were spliced from protein-coding exons ( Figure 2C). The length of most exonic circRNAs was less than 1300 nucleotides (nt), with a median length of~900 nt ( Figure S2A). CircRNAs were widely distributed across all chromosomes, with chromosome 1 containing the highest number of circular transcripts, namely, >1000 circRNAs, including 709 novel transcripts ( Figure S2B).

Profiles and Characteristics of circRNAs in Granulosa Cells of Ovarian Antral Follicles
In total, 9362 distinct circRNAs containing at least two unique back-spliced reads ( Figure 2A) were identified based on the circRNAs identification standard [23]. Although 67.3% of circRNAs had average coverage of less than 10 unique back-spliced reads, there were 441 circRNAs highly expressed containing a high-average read count of more than 50, including circ_KIF16B, circ_IL1R1, circ_SLC30A7, and circ_ANGPT1. We identified 2515 overlapping circRNAs and 6847 novel circRNAs that had not been identified so far ( Figure 2B; [24,25]). The 54.36% (5090/9362) of circRNAs identified were derived from exons and generated from individual linear cognate genes. The remainder were annotated to sense, among others, overlapping (2243/9362) and intergenic regions (1195/9362). All exonic circRNAs were spliced from protein-coding exons ( Figure 2C). The length of most exonic circRNAs was less than 1300 nucleotides (nt), with a median length of ~900 nt ( Figure  S2A). CircRNAs were widely distributed across all chromosomes, with chromosome 1 containing the highest number of circular transcripts, namely, >1000 circRNAs, including 709 novel transcripts ( Figure S2B).
Further analysis revealed that most host genes could only produce one single circRNA; nevertheless, there were a number of genes that produced multiple circRNAs ( Figure 2D, 9362 circRNAs from 3515 host genes). Seventy-two host genes generated more than 10 different circRNAs per gene. The most striking example was the 36 exons containing HEAT repeat containing 5A (HEATR5A) located on chromosome 7, which generated 54 distinct circRNAs (with at least two unique back-spliced reads). The abundance of these 54 isoforms was further investigated, and it appeared that one dominant isoform was significantly higher expressed than the other 53 isoforms. Similarly, 36% of the host genes (26/72) with more than 10 distinct circular isoforms also showed significantly higher expression compared to the other isoforms. This suggested that in gene loci with more circRNA isoforms, one isoform seemed to predominate over the other circRNA isoforms linked to this gene.  Further analysis revealed that most host genes could only produce one single circRNA; nevertheless, there were a number of genes that produced multiple circRNAs ( Figure 2D, 9362 circRNAs from 3515 host genes). Seventy-two host genes generated more than 10 different circRNAs per gene. The most striking example was the 36 exons containing HEAT repeat containing 5A (HEATR5A) located on chromosome 7, which generated 54 distinct circRNAs (with at least two unique back-spliced reads). The abundance of these 54 isoforms was further investigated, and it appeared that one dominant isoform was significantly higher expressed than the other 53 isoforms. Similarly, 36% of the host genes (26/72) with more than 10 distinct circular isoforms also showed significantly higher expression compared to the other isoforms. This suggested that in gene loci with more circRNA isoforms, one isoform seemed to predominate over the other circRNA isoforms linked to this gene.

Identification of Differentially Expressed circRNAs (DE-circRNAs)
The analysis of circRNAs expression showed that there were 62 circRNAs significantly differently expressed between HA and AA follicles ( Table 1, of which 49 circRNAs were significantly down-regulated, and 13 circRNAs were significantly up-regulated in AA follicles (fold change ≥ 2 and FDR ≤ 0.05, Figure 3A,B). All up-regulated and down-regulated circRNAs were further evaluated by unsupervised hierarchical clustering analysis ( Figure 3C), which showed a clear distinction in samples of HA and AA groups. Most DE-circRNAs were spliced from exons; over 84.37% of exonic circRNAs consisted of 2-6 exons, while a smaller fraction of DE-circRNAs was spliced from sense overlapping, intergenic, and antisense regions of the genome ( Figure S2C). The size distribution of the exonic DE-circRNAs was, in general, less than 1000 nt, while the length of the DE-circRNAs from sense overlapping regions was more than 3000 nt ( Figure S2D). In addition, DE-circRNAs were widely distributed across chromosomes 1-18 and the X chromosome ( Figure S2E). The distribution of DE-circRNAs was not uniform among different chromosomes; there was a general trend that numbers of DE-circRNAs per chromosome increased with absolute chromosome length ( Figure 3D).

GO and KEGG Analysis of Host Genes of DE-circRNAs
For the biological process ( Figure S3A, Table S1), five of the top 10 gene ontology (GO) terms were involved in gene expression-related processes, including gene expression, RNA metabolic process, transcription and DNA-templated synthesis, RNA biosynthetic process, and nucleic acid-templated transcription. Three others of 10 GO terms were related to homeostasis-related processes, including regulation of the homeostatic process, cellular homeostasis, and regulation of ion homeostasis. Regarding molecular function, the dominant categories were organic cyclic compound binding, heterocyclic compound binding, and metal ion binding, respectively ( Figure S3B). For cellular components, the dominant processes were centrosome, transcription elongation factor complex, and protein serine/threonine phosphatase complex ( Figure S3C). Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis of the host genes of these DE-circRNAs showed that apoptosis was one of the significantly enriched pathways ( Figure S3D, Table S2).
The underlying host genes related to GO processes and KEGG pathway analysis were similar. The enriched pathways included seven repetitive host genes of the DE-circRNAs (Table S2). Among these genes, three host genes appeared in most pathways, namely, PPP3CB, RB1, and SMAD2. In addition, PPP3CB enriched in the apoptosis-related pathway also appeared in nine out of the top 10 GO processes (Table S1), while RB1 and SMAD2 were present in six out of the top 10 GO processes. IL1R1 was present in three pathway maps, including apoptosis, CPEB2, NDUFB2, while NR5A2 appeared in only one pathway map. CPEB2 and NR5A2 were also present in the top 10 GO processes. These results indicated that the GO analysis was consistent with the pathway analysis.

Validation of DE-circRNAs Using qRT-PCR Combined with RNase R
Outward-facing divergent primers (Table S3) were designed against 10 randomly selected DE-circRNAs, and their expression levels were measured by qRT-PCR ( Figure 4A). The expression levels of circ_RB1, circ_CBFA2T2, circ_SLC30A7, circ_ANGPT1, circ_ANKHD1, and circ_FAM49B were significantly lower, while the expression levels of circ_WDR7, circ_KIF16B, circ_GMPS, circ_UNKL, and circ_IL1R1 were significantly higher in granulosa cells of AA follicles compared to HA follicles, confirming the RNA-seq data. In addition, we further confirmed the correctness of the amplified PCR products with specific circRNA junctions for circ_KIF16B, circ_SLC30A7, circ_ANGPT1, circ_ANKHD1, circ_FAM49B, and circ_CBFA2T2 by Sanger sequencing ( Figure 4B).
Four circRNAs (circ_ANGPT1, circ_SLC30A7, circ_KIF16B, and circ_CBFA2T2) and four corresponding linear mRNAs (ANGPT1, SLC30A7, KIF16B, and CBFA2T2) serving as controls were randomly selected for RNase R resistance analysis ( Figure 4C). The qRT-PCR results showed that the expression levels of these four circRNAs following RNase R digestion were similar to the untreated controls, while the expression levels of all four liner mRNAs were significantly decreased after RNase R treatment. These results further confirmed that these circRNAs were indeed circular in form. Representative examples of PCR products purified and sequenced to confirm circRNA junction sites by Sanger sequencing, n = 3. (C) qRT-PCR for the abundance of circRNAs (n = 6) and mRNAs (n = 6) treated with RNase R. The amount of circRNAs and mRNAs was calculated as fold change of RNase R over controls. The expression values in the controls were set to 1 for up-regulated genes and −1 for down-regulated genes. *, p < 0.05; **, p < 0.01; ***, p < 0.001. Follicles were isolated in batches at three different time points, each batch containing 15-20 ovaries.

circRNA-miRNA Interaction Network
A circRNAs-miRNA interaction network was constructed using Cytoscape based on our RNAseq data ( Figure S4). The top 10 up-regulated circRNAs and the top 20 down-regulated circRNAs were selected to identify circRNA-miRNA interactions based on TargetScan and miRanda analysis. The top five miRNAs for each circRNA were selected based on their context+ score.
When applying our stringent selection criteria for candidate circRNAs, two exonic circRNAs, with miRNA binding sites and fulfilling all criteria, namely, circ_CBFA2T2 (FC = 13; FDR = 0.007, being one of the down-regulated circRNAs) and circ_KIF16B (FC = 13.5; FDR = 0.004, being one of the most up-regulated circRNAs) (Table 1), attracted our interest. Following KEGG pathway analysis of the predicted targeted genes of circ_CBFA2T2 ( Figure S5A), we observed that the top enriched process referred to response to oxidative stress, a key trigger of apoptosis in granulosa cells [26]. Similarly, the pathway analysis of putatively targeted genes of circ_KIF16B showed that the main enriched processes directly associated with ovarian follicular development were the apoptosisrelated pathways, including the phosphatidylinositol 3-kinase (PI3K)-protein kinase B (AKT) athway Representative examples of PCR products purified and sequenced to confirm circRNA junction sites by Sanger sequencing, n = 3. (C) qRT-PCR for the abundance of circRNAs (n = 6) and mRNAs (n = 6) treated with RNase R. The amount of circRNAs and mRNAs was calculated as fold change of RNase R over controls. The expression values in the controls were set to 1 for up-regulated genes and −1 for down-regulated genes. *, p < 0.05; **, p < 0.01; ***, p < 0.001. Follicles were isolated in batches at three different time points, each batch containing 15-20 ovaries.

circRNA-miRNA Interaction Network
A circRNAs-miRNA interaction network was constructed using Cytoscape based on our RNA-seq data ( Figure S4). The top 10 up-regulated circRNAs and the top 20 down-regulated circRNAs were selected to identify circRNA-miRNA interactions based on TargetScan and miRanda analysis. The top five miRNAs for each circRNA were selected based on their context+ score.
When applying our stringent selection criteria for candidate circRNAs, two exonic circRNAs, with miRNA binding sites and fulfilling all criteria, namely, circ_CBFA2T2 (FC = 13; FDR = 0.007, being one of the down-regulated circRNAs) and circ_KIF16B (FC = 13.5; FDR = 0.004, being one of the most up-regulated circRNAs) (Table 1), attracted our interest. Following KEGG pathway analysis of the predicted targeted genes of circ_CBFA2T2 ( Figure S5A), we observed that the top enriched process referred to response to oxidative stress, a key trigger of apoptosis in granulosa cells [26]. Similarly, the pathway analysis of putatively targeted genes of circ_KIF16B showed that the main enriched processes directly associated with ovarian follicular development were the apoptosis-related pathways, including the phosphatidylinositol 3-kinase (PI3K)-protein kinase B (AKT) athway and apoptosis pathways ( Figure S5B). These two circRNAs were, therefore, selected for further functional analysis.

The Putative Role of Circ_CBFA2T2 in Apoptosis of Granulosa Cells in Antral Follicles
Circ_CBFA2T2 was derived from four exons, including exon 3-6 of the gene CBFA2T2 ( Figure 5A). qRT-PCR analysis of nuclear and cytoplasmic expression demonstrated that circ_CBFA2T2 was mainly localized in the cytoplasm of granulosa cells ( Figure 5B). The stability of circ_CBFA2T2 was investigated in granulosa cells treated with actinomycin D, an inhibitor of transcription. Total RNA was harvested at the indicated time points after treatment with actinomycin D ( Figure 5C). Analysis of circ_CBFA2T2 and CBFA2T2 mRNA showed that circular transcript half-life exceeded 24 h, whereas the associated linear transcript exhibited a half-life of less than 6 h, indicating that the circ_CBFA2T2 was highly stable.
Interestingly, these analyses showed that expression of circ_CBFA2T2 was significantly associated with the mRNA content of GCLC, a key antioxidant gene controlled by ssc-miR-9847-3p ( Figure 5F). The transcript content of circ_CBFA2T2 in the cumulus cells of AA follicles was similar to that of cumulus cells in HA follicles ( Figure 5G).

The Putative Role of Circ_KIF16B in Apoptosis of Granulosa Cells in Antral Follicles
Circ_KIF16B was derived from five exons, including exon14-18 of the KIF16B gene ( Figure 6A). Besides circ_KIF16B, there were four other circular transcripts observed that were derived from the KIF16B gene (named circ_KIF16B2, circ_KIF16B3, circ_KIF16B4, circ_KIF16B5, respectively). Circ_KIF16B was the predominant circular isoform as this isoform had the highest number of back-spliced unique reads (circ_KIF16B, 54; circ_KIF16B2, 26; circ_KIF16B3, 6; circ_KIF16B4, 9; circ_KIF16B5, 12). The qRT-PCR analysis of circ_KIF16B expression demonstrated that the circular form of KIF16B preferentially localized in the cytoplasm ( Figure 6B). Similar to circ_CBFA2T2, the stability of circ_KIF16B was investigated using actinomycin D ( Figure 6C). The analysis of circ_KIF16B and KIF16B mRNA showed that the circular transcript half-life exceeded 24 h, whereas the associated linear transcript exhibited a half-life of less than 6 h, indicating that the expression of circ_KIF16B was also highly stable.
To determine whether the up-regulation of circ_KIF16B expression was associated with the expression of the above 10 targeted mRNAs, linear regression analyses were performed. Interestingly, these analyses showed that the expression of circ_KIF16B was significantly associated with the mRNA content of TP53, a hub gene controlling apoptosis suppressed by ssc-miR-493-3p ( Figure 6F). In line with mRNA expression data, immunohistochemical staining for p53 in granulosa cells of HA follicles was weak to absent ( Figure 7A). In contrast, moderate to strong staining was observed in both the nuclei and cytoplasm of apoptotic granulosa cells in antral follicles ( Figure 7B), indicating the activation of p53 signaling during granulosa cell apoptosis.
Next, the expression of circ_KIF16B in cumulus granulosa cells of HA and AA follicles was investigated. qRT-PCR results showed that the transcript content of circ_KIF16B in cumulus granulosa cells of AA follicles did not differ significantly from that of HA follicles, although there was an increasing trend ( Figure 6G).
These results suggested that circ_KIF16B might promote apoptosis of granulosa cells, possibly by binding ssc-miR-493-3p, resulting in the up-regulation of the pro-apoptotic TP53 gene.

p53/PHLDA3 Pathway and Antral Follicular Atresia
To further confirm the role of the p53 signaling pathway as the downstream target of circ_KIF16B in granulosa cells, qRT-PCR and immunohistochemical staining experiments were performed. Firstly, downstream targets of TP53 were identified. qRT-PCR was used to measure the mRNA content of B-cell lymphoma 2 binding component 3 (PUMA), phorbol-12-myristate-13-acetate-induced protein 1 (BIM), and pleckstrin homology like domain family a member 3 (PHLDA3) in the granulosa cells of HA and AA follicles. These genes were selected based on the fact that they are reported p53 target genes in other cell types, including tumor cells and fibroblast cells [27,28]. Only the mRNA content of PHLDA3 was increased by more than five times in AA follicles compared to HA follicles, indicative of being a possible target of p53 in apoptotic granulosa cells. Immunohistochemical staining was used to detect the protein presence of p53, PHLDA3, and cCASP3 in the same HA and AA follicles. Consistent with the gene expression of PHLDA3, the immunostaining for PHLDA3 in granulosa cells of HA follicles was faint to absent ( Figure 7C), while moderate to strong staining was observed in the apoptotic granulosa cells of AA follicles ( Figure 7D). Similarly, cCASP3 staining was absent in the granulosa cells of HA follicles ( Figure 7E); however, many apoptotic cells with cCASP3 positive staining were present in the granulosa cell layer of AA follicles ( Figure 7F).

Discussion
This study was to our knowledge the first investigation that compared circRNAs profiles in granulosa cells of HA and advanced stage AA follicles in the porcine ovary. A large number of circRNAs from diverse genomic locations were identified; 62 of these circRNAs were DE-circRNAs. KEGG analysis of the host genes of the DE-circRNAs showed that apoptosis was one of the significantly enriched pathways. The characteristics of the DE-circRNAs, including back-splicing, RNase R resistance, and stability, were validated and suggested that these transcripts were real circular in form. We further showed that circ_CBFA2T2 and circ_KIF16B were two of the most relevant down-regulated and up-regulated exonic DE-circRNAs with miRNA binding sites. We further confirmed the expression of the targeted genes of these circRNAs, the antioxidant gene GCLC (potential target of circ_CBFA2T2), and the apoptotic gene TP53 (potential target of circ_KIF16B). Lastly, we showed that the p53/PHLDA3 pathway might participate in granulosa cell apoptosis. Based on these results, we hypothesized that these circRNAs might potentially play a role in granulosa cells apoptosis and thus antral follicular atresia.
Interestingly, despite relatively low expression levels per circRNA species, the diversity of the observed circRNA species is high, and many circRNA isoforms can be expressed by a single gene [29]. In line with this observation, our results showed that one parental gene could produce several circRNA isoforms, sometimes even more than 10. Only one or two isoforms were usually expressed at dominant levels, while most other isoforms exhibited low expression levels, suggesting that the dominant form might be purposefully produced to serve a specific function in granulosa cells.
CircRNAs have been suggested to function as miRNA sponges and participate in miRNA-mediated post-transcriptional regulation. Only a limited number of such circRNAs, however, contain multiple binding sites to sequester a particular miRNA. For example, circRNA-ciRS-7 harbors more than 70 conserved binding sites and acts as a sponge for miR-7 to regulate brain cell development [30]. Accumulating evidence suggests, however, that these large numbers of putative targeting sites are not necessary for circRNA to act as a miRNA sponge. An abundant circRNA derived from exon2 of the HIPK3 gene harbors only two binding sites for miR-124, but can despite this low number of targeting sites directly inhibit the activity of miR-124 [23]. In line with these observations, the present study showed that most of the DE-circRNAs harbored only one to two targeting sites for miRNAs. For instance, circ_CBFA2T2 and circ_KIF16B, located mainly in the cellular cytoplasm, had only two putative binding sites for ssc-miR-9847-3p and one putative binding site for ssc-miR-493-3p, respectively.
Several studies have reported that oxidative stress caused by the presence of the xenotoxic stressors can lead to ovarian follicular atresia [6,26]. In the present investigation, the mRNA concentrations of GCLC, the putative targeted gene of circ_CBFA2T2, and GCLM were significantly lower in the AA follicles compared to HA follicles. GCLC and GCLM are two subunits of glutamate-cysteine ligase, which acts as a rate-limiting enzyme in glutathione synthesis, an important antioxidant [31]. Hatzirodos et al. [22] observed a six-fold decrease in GCLC mRNA concentrations in small AA follicles in the bovine ovary. These observations strengthen the assumption of a disturbed glutathione antioxidant capacity in antral follicular atresia. Based on these observations, we speculated that circ_CBFA2T2 might play a role in antral follicular atresia via targeting GCLC.
miRNAs play an important role in the regulation of apoptotic processes. Kleemann et al. showed in ovarian carcinoma cell lines that miR-493-3p mimic transfection led to the disruption of the mitochondrial membrane potential and activation of caspases and apoptosis [8]. ssc-miR-493-3p is predicted to bind to circ_KIF16B, one of the identified DE-circRNAs in the present study. A predicted targeted gene of ssc-miR-493-3p is TP53, a gene in which the mRNA content was significantly up-regulated in granulosa cells of AA follicles in our study. The up-regulation of TP53 mRNA content coincided with intense p53 immunostaining in apoptotic granulosa cells of AA follicles, implicating that TP53 gene translation occurred, while p53 immunostaining was faint to absent in granulosa cells of HA follicles. This p53 immunostaining pattern coincided with previous observations in the rat ovary [32].
Recently, Haraguchi et al. (2019) found that in murine ovarian cumulus granulosa cells, the mouse double minute 2 homolog (Mdm2) -steroidogenic factor 1 (SF1) pathway played an important role in oocyte maturation by suppressing p53 activation; deletion of Mdm2 gene expression led to the activation of p53 [33]. We, however, did not observe differences in the expression of Mdm2 and SF1 in mural granulosa cells of healthy and atretic antral follicles (data not shown). Our data thus suggested that the activation of p53 in mural granulosa cells of atretic antral follicles in the porcine ovary did not seem to be due to a decrease in Mdm2 and SF1 gene expression.
It has been reported for other cell types, including fibroblasts and skin keratinocytes, that p53 preferably induces apoptosis via the B-cell lymphoma 2 (BCL-2)-regulated intrinsic apoptotic pathway by directly binding to its target genes, such as PUMA and BIM [27]. In our study, granulosa cells apoptosis, however, did not seem to favor these target genes as the mRNA content of BIM and PUMA was not altered in AA follicles compared to HA follicles. This suggested that another target gene must be involved in the p53-dependent apoptotic pathway in granulosa cells. Recently, a novel target gene of p53, named PHLDA3, has been shown to induce apoptosis by inhibiting the activity of the PI3K-Akt pathway in tumor cells [28,34]. The PI3K-Akt pathway plays an essential role in granulosa cell growth and apoptosis [5]. PHLDA3, a PH domain-only protein, can repress Akt activity by competitively binding to phosphatidylinositol 4,5-bisphosphate (PIP)2 and PIP3, leading to the activation of caspases and apoptosis. Ablation of endogenous PHLDA3 results in enhanced Akt activity and a decrease of p53-dependent apoptosis [28]. Our results showed that the mRNA content of PHLDA3 was significantly increased in AA follicles compared to HA follicles. Concomitantly, clear immunostaining of PHLDA3 and cCASP3 was confined to apoptotic granulosa cells in AA follicles, whereas staining in the granulosa cells of HA follicles was faint to absent. This suggested a potential role for PHLDA3 in granulosa cell apoptosis. We, therefore, hypothesized that circ_KIF16B might be involved in antral follicle atresia by serving as a miRNA sponge, which sequesters ssc-miR-493-3p, activating the p53/PHLDA3 signaling pathway and thus promoting granulosa cell apoptosis. The underlying molecular mechanisms of circ_CBFA2T2 and circ_KIF16B in granulosa cell apoptosis, however, needs to be further investigated.
Investigations of the role of circRNAs in ovarian follicular development so far have mainly focused on either the physiology of antral follicular growth or on ovarian diseases, e.g., PCOS. Not much is known about the role of circRNAs in follicle atresia. In a recent mouse study conducted by Jia and colleagues, it is reported that circRNAs might participate in ovarian granulosa cell apoptosis [12]. RNA-Seq analysis of pubertal (six-week-old) and neonatal (5-day-old) mouse ovaries showed up-regulation of circ_EGFR in the pubertal ovary. In vitro overexpression of circ_EGFR promoted the proliferation of granulosa cells, while knockdown of circ_EGFR significantly decreased the rate of proliferation, possibly via increased granulosa cell apoptosis. Our observations from the in vivo perspective extended these findings, as we identified by RNA-Seq 62 DE-circRNAs from a total of 9632 circRNAs in isolated porcine granulosa cells of HA and AA follicles. KEGG pathway analysis of the host genes of these DE-circRNAs showed that the apoptosis pathway was significantly enriched. We observed an increased mRNA content and strong immunostaining of p53/PHLDA3, a target of circ_KIF16B, in granulosa cells of AA follicles compared to HA follicles. These results offered further support for the role of circRNAs in the apoptosis of ovarian granulosa cells. Apart from our study, Guo et al. [35] showed that circ_INHA, generated from the intron of the inhibin subunit alpha gene, was down-regulated in porcine granulosa cells of early AA follicles compared to healthy follicles. circ_INHA inhibited granulosa cell apoptosis, possibly through sponging miR-10a-5p, thus modulating connective tissue growth factor (CTGF, CNN2) expression. CTGF, a member of the calponin (CNN) matricellular protein family, plays a role in many biological processes, including cell proliferation and apoptosis [35]. In contrast to the study by Guo and colleagues, we did not observe differential expression of circ_INHA via deep circRNA-Seq in granulosa cells from HA and advanced stage AA follicles. This discrepancy might be due to differences in the size of the selected antral follicles, as the present study analyzed DE-circRNAs in freshly isolated granulosa cells from 4-7 mm healthy and advanced atretic antral follicles collected during the follicular phase of the estrous cycle, while Guo and colleagues used somewhat smaller 3-5 mm antral follicles and analyzed DE-circRNAs in granulosa cells of healthy and early atretic antral follicles collected at an unknown stage of the estrous cycle.
In ovarian antral follicles, granulosa cells are subdivided into two distinct classes that are anatomically and functionally distinct: mural granulosa cells, which line the follicular wall, and cumulus granulosa cells, which surround and are in direct contact with the oocyte, together forming the cumulus-oocyte-complex (COC) [36]. There are several observations that have demonstrated that in antral follicular atresia, mural granulosa cells undergo apoptosis, while, at the same time, the cumulus cells remain healthy [37,38]. We confirmed these observations as we did not see morphological differences between COCs retrieved from AA or HA follicles. The absence of a difference in transcript content of circ_CBFA2T2 and circ_KIF16B between cumulus granulosa cells from HA and AA follicles further supported the morphological observation that during follicular atresia when increasing numbers of mural granulosa cells underwent apoptosis, cumulus cells remained healthy. It is not until late-stage atresia when most mural cells have undergone apoptosis that the cumulus cells also become apoptotic [36].

Animal and Follicles Collection
The whole experiment design was depicted as shown in the Figure 8. Sixty porcine ovaries from 30 gilts (nulliparous, around 180 days old with a bodyweight of approximately 120 kg) were used in this study. The gilts were housed at a local farm and monitored once daily to determine the stage of the estrous cycle in the presence of a boar. The first day of the estrous cycle onset was counted as day 0 of the estrous cycle. The animals were slaughtered in a local abattoir on day 16-19 of the estrous cycle, the beginning to mid follicular phase of the cycle [39]. In total, 60 ovaries were collected in our study, of which six ovaries from six different animals were snap-frozen immediately after the slaughter in the nitrogen liquid and stored at −80 • C until further immunohistochemical processing. The other 54 ovaries were collected immediately after slaughter and washed in 0.01 M phosphate-buffered saline pH 7.4 (PBS), transferred to a 30 • C PBS solution containing 1% penicillin-streptomycin, and transported to the laboratory within 1 h after slaughter. To ensure that follicle dissection was finished within a reasonably short time frame, the sows were slaughtered at three different days (three different batches). At each slaughter day, 15 to 20 ovaries were collected for antral follicle dissection (in total, 54 ovaries were used for follicle dissection). While dissecting one ovary, the remaining ovaries of that batch were stored in dulbecco's phosphate buffered saline (DPBS) on ice. For each batch, follicles were dissected within 1.5 h after arrival in the lab. Antral follicle size was estimated by measurement of two perpendicular diameters with a millimeter scale. Individual antral follicles, approximately 4-7 mm in diameter, were dissected from the ovaries. Healthy antral (HA) follicles and atretic antral (AA) follicles were classified morphologically under a stereomicroscope, as described previously [40][41][42]. Briefly, HA follicles were characterized by a vascular sheath on the follicular surface with a pinkish color and clear follicular fluid in the antrum. AA follicles were identified by an opaque color, the absence of clear vascularization on the surface of the follicle, and a large amount of debris floating in the follicular fluid [42,43]. Based on the criteria from Gioia et al. [43], these AA follicles were considered to be at an advanced stage of atresia.
In order to confirm the classification, 10 randomly selected follicles classified as either wellvascularized or non-vascularized from total 245 dissected follicles were fixed in diluted Bouin's fluid (0.9% picric acid, 4% formaldehyde, 5% glacial acetic acid) for 24 h at 4 °C and embedded in paraffin. The follicles were sectioned at a thickness of 5 µm. Six sections of each follicle were randomly selected, mounted on glass slides, stained with periodic acid Schiff's reagent (PAS) and Mayer's hematoxylin, and examined by light microscopy. Antral follicles were considered to be atretic when more than 5% of the granulosa cells showed signs of apoptosis, while the theca layer of these follicles showed signs of hypertrophy. In more advanced atretic follicles, apoptotic granulosa cells were not only observed in the mural cell layer but also throughout the follicular antrum [5,26]. The status of follicular atresia was further confirmed by immunostaining for the presence of cCASP3 (see below). The histological analysis confirmed the macroscopical identification of HA and AA follicles ( Figure  S1). Healthy antral (HA) follicles and atretic antral (AA) follicles were classified morphologically under a stereomicroscope, as described previously [40][41][42]. Briefly, HA follicles were characterized by a vascular sheath on the follicular surface with a pinkish color and clear follicular fluid in the antrum. AA follicles were identified by an opaque color, the absence of clear vascularization on the surface of the follicle, and a large amount of debris floating in the follicular fluid [42,43]. Based on the criteria from Gioia et al. [43], these AA follicles were considered to be at an advanced stage of atresia.
In order to confirm the classification, 10 randomly selected follicles classified as either well-vascularized or non-vascularized from total 245 dissected follicles were fixed in diluted Bouin's fluid (0.9% picric acid, 4% formaldehyde, 5% glacial acetic acid) for 24 h at 4 • C and embedded in paraffin. The follicles were sectioned at a thickness of 5 µm. Six sections of each follicle were randomly selected, mounted on glass slides, stained with periodic acid Schiff's reagent (PAS) and Mayer's hematoxylin, and examined by light microscopy. Antral follicles were considered to be atretic when more than 5% of the granulosa cells showed signs of apoptosis, while the theca layer of these follicles showed signs of hypertrophy. In more advanced atretic follicles, apoptotic granulosa cells were not only observed in the mural cell layer but also throughout the follicular antrum [5,26]. The status of follicular atresia was further confirmed by immunostaining for the presence of cCASP3 (see below). The histological analysis confirmed the macroscopical identification of HA and AA follicles ( Figure S1). Based on the criteria described above, 120 follicles (60 HA follicles and 60 advanced AA follicles) were selected from a pool of 235 follicles; the 115 follicles that were not included in the follicle selection did not meet our selection criteria for healthy and advanced atretic antral follicles based and were thus excluded from sampling. These two pools of 60 follicles were randomly divided into 6 groups of 10 follicles each. Follicular fluid was collected and pooled per 10 follicles. In the advanced stage atretic antral follicles, apoptotic granulosa cells had become lose from the follicle wall and were floating in the antrum, mixing with the follicular fluid. When dissecting these follicles, the follicular fluid and floating apoptotic granulosa cells were collected together in a Petri dish. The mural granulosa cells of these atretic follicles were scraped from the follicular wall and mixed with this follicular fluid to obtain a complete sample of granulosa cells (floating and mural cells) from the atretic follicles. To create comparable granulosa cell samples, the mural granulosa cells scraped from the walls of healthy antral follicles were also mixed with the follicular fluid. In brief, follicles were dissected in a Petri dish in PBS and cut in smaller segments. Follicular wall fragments (consisting of theca and mural granulosa layer) and the cumulus-oocyte-complex (COC) were separated. COCs were collected with a Pasteur pipette and transferred to another Petri dish, where cumulus granulosa cells were mechanically stripped from the oocyte. The collected cumulus cells were washed in PBS, centrifuged, and then stored at −80 • C until RNA extraction. Simultaneously, the mural granulosa cells were scraped off from the dissected follicle wall by gentle rubbing with a glass round-smooth Pasteur pipette [22]. The scraped granulosa cells were mixed with the follicular fluid of the corresponding follicles, diluted 1:3 in PBS, transferred to a 1.5 mL Eppendorf (EP) tube, and centrifuged at 800 g for 5 min. The supernatant was collected and stored at −80 • C for subsequent hormone analysis. The residual granulosa cells remaining in the precipitate were stored at −80 • C until RNA isolation.
To test the purity of the isolated granulosa cells, cells were labeled with an antibody against the follicle-stimulating hormone receptor (FSHR, ab150557, Abcam), a specific marker for granulosa cells [44]. On average, approximately 1 × 10 4 cells in three biological replicates (three different pools of follicles derived from three different batches of sows) were analyzed. After FSHR immunostaining, pictures were taken by a light microscope (Leica DM6000B, equipped with a DFC356FX digital camera and LasX imaging software; Leica Microsystems, Guangzhou, China) and analyzed using Image J software.

Granulosa Cells Cultured In Vitro
To test the stability of circRNAs, in vitro cultured granulosa cells were treated with actinomycin D, a transcription inhibitor. Isolated granulosa cells were seeded into 6-well or 12-well plates and cultured in Dulbecco's Modified Eagle Medium/Nutrient Mixture F-12 supplied with 10% FBS and 1% penicillin-streptomycin (Gibco, Thermo Fisher, Shanghai, China) for 24 h. Transcription was blocked by adding 2 µg/mL actinomycin D or dimethylsulfoxide (Sigma-Aldrich, St. Louis, MO, USA) as a control, with three replicates for each treatment under five time-points of 0, 6, 12, 18, and 24 h, respectively. At each time point, cells were collected for RNA isolation, and qRT-PCR was used to quantify the content of circRNAs and their linear mRNAs (see below).

Circular RNA Library Construction and Illumina Sequencing
Total RNA from granulosa cells of three randomly selected samples of pooled granulosa cells from either 10 HA or 10 AA follicles was isolated using TRIzol (Life Technologies, Shanghai, China) for the purpose of RNA-sequencing (RNA-seq). High throughput transcriptome sequencing was carried out by Cloud-Seq Biotech (Shanghai, China). Briefly, total RNA extraction and rRNA depletion were performed using the Ribo-Zero rRNA Removal Kit (Illumina, Shanghai, China), according to the manufacturer's instructions. RNA libraries were constructed from the rRNA-depleted RNA using the TruSeq Stranded Total RNA Library Prep Kit (Illumina), according to the manufacturer's instructions. The quality and quantity of the libraries were checked by the BioAnalyzer 2100 system (Agilent Technologies, Palo Alto, CA, USA). Ten pM libraries were denatured, captured on Illumina flow cells, amplified in situ, and finally sequenced for 150 cycles using the Illumina HiSeq Sequencer, according to the manufacturer's instructions. All Illumina sequencing data have been submitted to the Gene Expression Omnibus (GEO) under accession number (GSE136589).
The reads were aligned to the porcine genome (Sscrofa10.2/SusScr3) by employing bowtie2 software (v. 2.2.4, http://bowtie-bio.sourceforge.net/bowtie2/index.shtml) [47]. CircRNAs were detected and identified with find-circ software [30] downloaded from circBase [48]. The find-circ pipeline was run as suggested by the developers, except for an increased filtering stringency, requiring that both anchor segments mapped to the genome with mapping scores of 40. Only circRNAs, with two or more supporting reads within single samples, were used for further analysis.
Host genes, giving rise to individual circRNAs, were identified by matching the genomic location of circRNAs with the location of genes detected by TopHat/Cufflinks using BEDtools [49]. Published circRNAs were used to annotate the identified circRNAs. Back-splice junction reads for every circRNAs were normalized by mapped read numbers for each sample, and then log2 transformed. To detect differentially expressed circRNAs between granulosa cells from HA and AA follicles, a multiple testing correction method was used to estimate the false discovery rate (FDR). GO and pathway enrichment analysis were performed for the host genes of the differentially expressed circRNAs, using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) v6.8 [50]. Differently expressed circRNAs were visualized as heatmaps using the MultiExperiment Viewer (MeV, http://mev.tm4.org/) [51].

miRNA Target Prediction and circRNA-miRNA-mRNA Network Construction
The predicted miRNA binding sites for the differentially expressed circRNAs were identified using the miRanda [52] and Targetscan [53] algorithms. Only those circRNA-miRNA interactions predicted by both algorithms were used for the downstream network construction and analyses. miRNA-mRNA interactions that were common in both miRanda and TargetScan were then used to determine the gene targets of each filtered miRNA. Using these data, the outline of a circRNA-miRNA-mRNA regulatory network was visualized using Cytoscape (version 3.8, https://cytoscape.org/).

Candidate circRNAs Selection
To confirm that these circRNAs indeed played a functional role via miRNA sponging in antral follicular atresia, candidate circRNAs were selected from differentially expressed (DE)-circRNAs according to the following criteria. To exert miRNA sponging function, these circRNAs needed to be stably located in the cytoplasm and have miRNA binding sites [23]. Thus exon-derived circRNAs, which were mainly cytoplasmic in location, should be suitable candidates for further investigation (Table 1) (criterium 1). Only DE-circRNAs with a fold change (FC) of at least 10 were considered to be suitable for selection (criterium 2). The DE-circRNAs should be consistently expressed in all samples with FDR values < 0.01 (criterium 3). Lastly, the functions of the target genes of the candidate circRNAs needed to be described in the literature and preferably should be directly associated with antral follicular atresia (criterium 4).

RNA Preparation and qRT-PCR
The nuclear and cytoplasmic fractions of three samples, each containing pooled granulosa cells of 10 HA and 10 AA follicles, were extracted using the PARIS KIT (Thermo Fisher, Shanghai, China) according to the manufacturer's instructions. Briefly, approximately 10 7 granulosa cells were collected, washed in PBS once, and then placed on ice. The granulosa cells were subsequently resuspended in 500 µL ice-cold cell fractionation buffer and incubated on ice for 10 min to induce complete lysis. The samples were centrifuged at 500× g for 5 min at 4 • C, after which the cytoplasmic fraction was carefully aspirated from the nuclear pellet. The nuclear fraction was further lysed in cell disruption buffer. Following this, total RNA from whole-cell lysates of the nuclear and cytoplasmic fraction, respectively, was isolated using TRIzol (Life Technologies).
Next, 3 µg of total RNA of each sample was incubated for 15 min at 37 • C with or without RNase R (3 U/µg RNA, Epicenter Technologies, Shanghai, China), followed by inactivation for 3 min at 95 • C. The resulting RNA was subsequently purified using the RNeasy MinElute Cleaning kit (Qiagen, Guangzhou, China). In order to guarantee similarly effective RNase R treatment, all RNA samples were treated simultaneously in one run.
To quantify the amount of circRNA and mRNA, the qRT-PCR analysis was performed; 1 µg RNA of each sample was used for cDNA synthesis using the PrimeScript RT Master Mix (Takara, Guangzhou, China). Quantitative RT-PCR reactions were performed employing SYBR Premix Ex Taq II (Takara, Guangzhou, China). Individual samples were measured in duplicate. A standard curve using serial dilutions of the pooled sample (cDNA from all samples), negative control without cDNA template, and a negative control without reverse transcriptase (RT) were included in every assay. Only standard curves with efficiency between 90% and 110% and a correlation coefficient above 0.99 were accepted. Data were normalized against the reference gene β-actin. Divergent primers annealing at the distal ends of circRNAs were used to determine the abundance of circRNAs, while the convergent primers were used to measure the mRNA content. To further confirm the junction sequence of circRNAs, PCR products of divergent primers were gel purified and handed in for Sanger sequencing at Sangon Biotech Co. Ltd. (Shanghai, China). The primers used and PCR annealing temperatures for each gene are summarized in Table S3.

Immunohistochemistry
Frozen ovaries (n = 6 from six different animals) were serial sectioned (7 µm). Serial sections from each of the six ovaries were selected at random and mounted on Superfrost plus glass slides (Menzel-Gläser, Braunschweig, Germany). To determine the presence of proteins (p53, PHLDA3, and cCASP3) in porcine ovaries, immunohistochemistry was performed, according to Meng et al., with some modifications [5,26]. For each antibody tested, all ovarian sections were stained in one run, in order to be able to compare the immunohistochemical staining among the different animals. Briefly, sections were air-dried for 30 min and fixed in 4% phosphate-buffered paraformaldehyde for 10 min. Slides were subsequently washed in water, microwaved in sub-boiling 0.1 M sodium citrate buffer (pH = 6) for 10 min for epitope antigen retrieval. Slides were cooled down to room temperature and rinsed with Tris-buffered saline (TBS) pH 7.4. Endogenous peroxidase activity was blocked with 3% (v/v) hydrogen peroxide in methanol. Aldehyde residues were blocked with 0.3% glycine in TBS for 30 min. After rinsing with TBS, sections were incubated with 5% (wt/v) normal goat serum in TBS for 60 min at room temperature. Following the removal of goat serum, sections were incubated overnight at 4 • C in a humidified chamber with the respective primary antibodies (p53, 1:200; PHLDA3, 1:200; cCASP3, 1:1000) diluted in TBS-BSAc (Aurion, Wageningen, The Netherlands). Sections were rinsed and incubated at room temperature with the corresponding secondary biotin-labeled antibody diluted 1:200 in TBS-BSAc. After rinsing with TBS, the avidin-biotin complex (ABC kit elite, Vector Laboratories, Burlingame, CA, USA) was diluted 1:750 (v/v) in TBS-BSAc. Bound antibodies were visualized using 3-3 diaminobenzidine (Impact DAB kit, Vector Laboratories) diluted 1:400 (v/v). Sections were counterstained with Mayer's hematoxylin. Control sections were incubated with isotype IgG (Vector Laboratories, Burlingame, CA, USA), instead of the respective primary antibodies, according to the manufacturer's instructions. Background staining in these controls was negligible.

Statistical Analysis
To perform the statistical analysis, as described above, GraphPad Prism version 7.00 (Graphpad Software, San Diego, CA, USA) was used. Data were expressed as the mean ± standard error of the mean (SEM). Data were checked for normality, and when the normality was confirmed, the Student's t-test was used for data analysis. If normality could not be assumed, data were log10 transformed. If data was still not normally distributed after log transformation, a Mann-Whitney non-parametric test was used. p-values < 0.05 were considered to be significantly different.

Conclusions
Taken together, our study provided an overview of the presence of circRNAs in granulosa cells from HA and AA follicles. The association between the expression of two differentially expressed candidate circRNAs (circ_CBFA2T2 and circ_KIF16B) and their respective targeted genes GCLC and TP53 implicated a potential role for these circRNAs in granulosa cell apoptosis. The identification of RNA circularization may further expand our understanding of the molecular regulation atresia in antral follicles.
Supplementary Materials: Supplementary Materials can be found at http://www.mdpi.com/1422-0067/21/15/ 5217/s1. Figure S1. Representative immunostaining of cCASP3 (brown staining) as a marker of apoptosis in healthy (A) and atretic (B) antral follicles (n = 6 follicles from different animals). (A) Healthy antral follicles with no cCASP3 staining in the granulosa cells; (B) Atretic antral follicles with numerous granulosa cell-derived apoptotic bodies, adjacent to and within the antrum, that stain positively for the presence of cCASP3 (detail shown in insert). Granulosa cells are indicated by an asterisk, theca cells by an arrowhead. Scale bar represents 50 µm.  Figure S3. GO and KEGG analysis performed for the host genes of the DE-circRNAs. Two-D bar graph of biological process (A), molecular function (B), and cellular component (C) categories; bubble chart of potential signaling pathways (D). The length and color of each bar represent the number of genes and p-value, respectively. Similarly, the size and color of each bubble represent the number of genes in each pathway and p-value, respectively. DE-circRNAs, differentially expressed circRNAs. n = 3. Figure S4. CircRNAs-miRNAs interaction network. The V shape nods with red and green color represent circRNAs, and ellipse nodes (blue color) represent miRNAs. The red and green color denotes the up-regulation and down-regulation of circRNAs, respectively. The grey solid line represents the interaction between circRNAs and miRNAs, and the thickness of the line indicates the combined score. Figure S5. Representative GO terms of targeted genes of circ_CBFA2T2 (A) and circ_KIF16B (B). The size and color of each bubble represent the number of genes in each pathway and p-value, respectively. Table S1 All the significant GO terms; Table S2 KEGG