WGCNA Analysis Identifies Polycystic Ovary Syndrome-Associated Circular RNAs That Interact with RNA-Binding Proteins and Sponge miRNAs

Introduction

Polycystic ovary syndrome (PCOS) is a common endocrine metabolic disorder in women of childbearing age.¹ PCOS patients are typically characterized by androgen excess, polycystic ovaries and anovulation, and are often accompanied by insulin resistance and obesity.² Consequently, patients with this syndrome not only suffer from infertility,³ but also have an increased risk of developing diabetes, hypertension, cardiovascular diseases and cancer.⁴ However, management of this condition is still challenging, with anti-androgen and ovulation stimulating hormones as well as insulin-sensitizing drugs predominantly, which may cause serious side effects.² These findings indicate the necessity to deeply explore the molecular mechanisms of PCOS to develop more effective therapeutic strategies.

Ovarian follicular development and ovulation is a complex process that involves the communication between the oocytes and the surrounding granulosa cells (especially the cumulus subset).⁵ The cumulus granulosa cells were reported to promote oocyte growth by providing energy source (eg pyruvate, free fatty acids), steroid hormone (conversion of cholesterol to pregnenolone, progesterone, estradiol in turn under the binding of luteinizing hormone receptor and follicle-stimulating hormone with their corresponding ligands) and various growth regulators (eg epidermal growth factor) through gap junctions.^5–8 Dysfunction of cumulus granulosa cells (including increased apoptosis, decreased proliferation, estrogen deficiency and androgen accumulation) led to abnormal folliculogenesis and the development of PCOS.^9–11 Therefore, investigation of molecular mechanisms associated with dysfunction of cumulus cells may provide important insights into the pathogenesis of PCOS and potential therapeutic targets.

Circular RNAs (circRNAs) are a class of non-coding RNAs that form a covalently closed-loop structure from back-spliced exons, without a canonical 5′-end cap and 3′-end poly (A) tails.¹² Mounting evidence indicates that circRNAs participate in various biological processes by functioning as microRNA (miRNA) sponge, interacting with RNA-binding proteins (RBPs) and regulating the expression of their host genes.^13,14 These studies reveal circRNAs may represent attractive players for PCOS. This hypothesis has been demonstrated previously. Ma et al identified 286 circRNAs that were differentially expressed between PCOS and non-PCOS groups by microarray analysis and confirmed hsa_circ_0043533 was upregulated, while hsa_circ_0097636 was downregulated by quantitative reverse transcription polymerase chain reaction (qRT-PCR). The expression levels of hsa_circ_0043533 and hsa_circ_0097636 were positively correlated with serum testosterone in PCOS patients.¹⁵ Similarly, Che et al¹⁶ and Zhang et al¹⁷ used the microarray and qRT-PCR techniques to validate four (hsa_circ_0083952, hsa_circ_0082709, hsa_circ_0002425, hsa_circ_0015168) and two (hsa_circ_0001577, hsa_circ_0020093) circRNAs in the granulosa cells may be related to PCOS disease, respectively. Che et al¹⁶ also predicted miRNAs that interacted with circRNAs. In the study of Liu et al, overexpression of circPSMC3 was proved to decrease the insulin release and increase the number of granulosa cells in PCOS model mice by sponging miR-296-3p and then regulating PTEN expression.¹⁸ Additionally, in vitro studies implied ciRS-126¹⁹ and circASPH²⁰ regulated the proliferation and apoptosis of ovarian granulosa cells by targeting the miR-21-programmed cell death 4-reactive oxygen species and miR-375/mitogen-activated protein kinase kinase 6 axis, respectively. However, PCOS-associated circRNAs remain rarely reported.

The present study aimed to further screen PCOS-associated circRNAs by re-analysis of the microarray dataset of Ma et al¹⁵ based on weighted gene co-expression network analysis (WGCNA)²¹ that is a powerful systems biology algorithm for finding modules of highly correlated genes and associating the gene modules with clinical features. Furthermore, we also integrated the high-throughput datasets of miRNAs²² and mRNAs in cumulus cells of PCOS patients to construct the circRNA-miRNA-mRNA competing endogenous RNA (ceRNA) network and circRNA-RBP interaction network to predict the function mechanisms of circRNAs (Figure 1). Our study may highlight the importance of crucial circRNAs in the etiology of PCOS. Targeted regulation of these circRNAs may represent underlying treatment strategies for PCOS.

Figure 1 A flow chart displaying the analysis process in this study. Abbreviations: GEO, gene expression omnibus; PCOS, polycystic ovary syndrome; circRNAs, circular RNAs; miRNA, microRNA; RBPs, RNA-binding proteins; ceRNA, competing endogenous RNA; DECs, differentially expressed circRNAs; DEMs, differentially expressed miRNAs; DEGs, differentially expressed genes; DAVID, Database for Annotation, Visualization and Integrated Discovery; GS, gene significance; MM, module membership; PPI, protein–protein interaction.

Materials and Methods

Dataset Collection

PCOS-associated expression profile data of GSE145296,¹⁵ GSE72274²² and GSE155489²³ were obtained from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database by searching with the following keywords: (“circRNA” or “miRNA” or “transcriptomic”) and (“polycystic ovary syndrome”) and (“cumulus cells”). GSE145296 dataset analyzed the circRNAs in cumulus cells from six PCOS patients and six normal control individuals using Agilent-084217 CapitalBio Technology Human CircRNA Array v2. GSE72274 explored the altered miRNA profiles in cumulus cells from five women with PCOS and five without PCOS using Agilent-038166 cbc_human_miR18.0 array. GSE155489 investigated the transcriptome changes of cumulus cells in four women with PCOS compared with four women without PCOS using HiSeq X Ten whole-genome sequencing technology.

Differential Analysis for circRNAs, miRNAs and mRNAs

Differentially expressed circRNAs (DECs), miRNAs (DEMs) and protein-coding genes (DEGs) between PCOS patients and normal controls were identified by using the Bioconductor limma package of R software (v3.34.7; https://bioconductor.org/packages/release/bioc/html/limma.html).²⁴ P-value <0.05 and |log₂fold-change (FC)| ≥ 0.5 were set as the statistical thresholds for DECs, DEMs and DEGs. Heatmap was plotted using the “pheatmap” package (v1.0.8; https://cran.r-project.org/web/packages/pheatmap).

Function Enrichment Analysis

To understand the functions of DECs and DEGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) (v6.8; http://david.abcc.ncifcrf.gov). P-value <0.05 was considered to be statistically significant.

WGCNA Analysis for Identification of Hub circRNAs for PCOS

The WGCNA R package (v1.61; https://cran.r-project.org/web/packages/WGCNA/index.html)²¹ was used to construct a weighted correlation network using the expression data profile of DECs. First, Pearson coefficient was calculated to evaluate the correlation between the two genes. To ensure a scale-free topology, the soft threshold power was selected using the pickSoftThreshold function, by which the adjacency matrix was generated. Then, the adjacency matrix was transformed into a topological overlap matrix (TOM). A hierarchical clustering tree was built according to a TOM-based dissimilarity measure. Significant co-expression modules were identified using the DynamicTreeCut function, with the parameters set as minModuleSize ≥ 50 and cutHeight = 0.995. Next, module eigengene (ME, the representation of the overall expression level of each module), gene significance (GS, the correlation of individual genes with clinical trait) and module membership (MM, the correlation of individual genes with ME) were calculated. Pearson’s correlation analysis was performed to explore the correlation between ME and clinical trait as well as between GS and MM, in which p < 0.05 indicated the modules or genes in the modules were significantly associated with PCOS. The hub circRNAs in the significant modules were filtered with the threshold of |GS| > 0.8 and |MM| > 0.8.

Construction of ceRNA and circRNA–RBP Interaction Networks

The circinteractome (https://circinteractome.irp.nia.nih.gov) and starbase (v3.0, http://starbase.sysu.edu.cn) databases were used to predict the miRNAs or RBPs that interacted with DECs. The DEC-DEM pairs predicted by both the circinteractome and starbase databases were suggested to be important, while the DECs-DEGs (RBPs) interaction pairs were selected if they were predicted by any one database (in starbase, CLIP-Data was set as ≥3). The miRwalk database (v2.0; http://www.zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2) that contains 12 prediction programs was utilized to predict the target genes of DEMs. Those DEGs, which were predicted by at least 5 algorithms, were considered as the true candidate target mRNAs of DEMs. Considering that the network would be larger, the DEGs only with adjusted p-value < 0.05 and |log₂FC| ≥ 0.5 were used as the target genes of DECs or DEMs. The DEC-DEM and DEM-DEG interaction pairs were used to construct a ceRNA regulatory network. The DEC-RBP interaction pairs were used for establishing a circRNA-RBP network. These two networks were visualized by using Cytoscape (v3.6.1; www.cytoscape.org/).

Construction of a Protein–Protein Interaction (PPI) Network

To screen hub genes and explore their possible function mechanisms, the interactions between DEGs in the ceRNA and circRNA-RBP networks were predicted using the STRING (Search Tool for the Retrieval of Interacting Genes; v10.0; http://stringdb.org/) database. The DEGs with a confidence score ≥0.4 were included to construct the PPI network using Cytoscape. In addition, degree centrality was calculated for each protein in the PPI network using the CytoNCA plugin in Cytoscape software (http://apps.cytoscape.org/apps/cytonca). The hub genes were determined if they possessed a high degree.

Results

Identification of DECs, DEGs and DEMs

A total of 134,751 circRNAs were identified in the GSE145296 dataset. Of them, 3614 circRNAs were differentially expressed in cumulus cells of PCOS patients compared to controls, including 1823 upregulated and 1791 downregulated DECs (Figure 2A). The GSE155489 dataset comprised 22,064 mRNAs, from which 3544 displayed significantly differential expression in cumulus cells of PCOS patients compared to controls, including 1880 upregulated and 1663 downregulated DEGs (Figure 2B). Of 1887 miRNAs in the GSE72274 dataset, 1469 were screened as DEMs, including 35 upregulated and 1434 downregulated (Figure 2C). The crucial DECs, DEGs and DEMs are listed in Table 1.

Table 1 Crucial circRNAs, miRNAs and mRNAs in PCOS

Figure 2 Heatmaps showing hierarchical clustering analysis of differentially expressed circRNAs, protein-coding genes and miRNAs between PCOS and controls. (A) for differentially expressed circRNAs; (B) for differentially expressed protein-coding genes; (C) for differentially expressed miRNAs. PCOS, polycystic ovary syndrome. The color key is the Z-score.

Function Enrichment Analysis for DECs and DEGs

To reveal the potential functions of identified DECs, the parental genes of circRNAs were annotated, and function enrichment analysis was performed. A total of 215 GO biological process terms (Table S1) and 36 KEGG pathways (Table S2) were enriched, in which DECs were mainly involved in cell proliferation and immune regulation, such as GO:0051301~cell division (TPX2), GO:0000086~G2/M transition of mitotic cell cycle (TPX2), hsa04914:Progesterone-mediated oocyte maturation, hsa04110:Cell cycle and hsa04662:B cell receptor signaling pathway (Figure 3).

Figure 3 Function enrichment analysis for host genes of differentially expressed circRNAs.

Function enrichment analysis was also performed for all DEGs. The results showed that 80 GO biological process terms (Table S3) and 29 KEGG pathways (Table S4) were enriched. In line with the results of DECs, DEGs were also mainly associated with inflammation-mediated cell apoptosis, such as GO:0042981~regulation of apoptotic process (PEA15, proliferation and apoptosis adaptor protein 15; caspase 2, CASP2), GO:0002223~stimulatory C-type lectin receptor signaling pathway [TAB2, TGF-beta activated kinase 1 binding protein 2; PAK3, p21 (RAC1) activated kinase 3], GO:0045893~positive regulation of transcription, DNA-templated (RUNX2, RUNX family transcription factor 2), hsa05205:Proteoglycans in cancer (FZD4, frizzled class receptor 4), hsa04668:TNF signaling pathway (TAB2) and hsa04660:T cell receptor signaling pathway (PAK3). Furthermore, DEGs also participated in steroidogenesis, including GO:0006695~cholesterol biosynthetic process (HMGCR, 3-hydroxy-3-methylglutaryl-CoA reductase), hsa01130:Biosynthesis of antibiotics (HMGCR) (Figure 4). It was reported that the high androgen accumulation stimulated inflammation and promoted the apoptosis of cumulus cells.²⁵ Therefore, circRNAs that regulate inflammation-, apoptosis- and steroidogenesis-related DEGs may be important for PCOS.

Figure 4 Function enrichment analysis for differentially expressed protein-coding genes.

WGCNA Analysis Identifies Hub DECs Associated with PCOS

In WGCNA analysis, the soft threshold power was selected as 22 to ensure a scale-free network distribution (scale free index R² > 0.9; connectivity = 1) (Figure 5A). Five modules were identified by the DynamicTreeCut algorithm (Figure 5B), of which 853, 196, 902, 1547 and 116 DECs were enriched in the blue, brown, grey, turquoise and yellow module, respectively. All these modules were significantly associated with PCOS (Figure 5C). Also, the GS for PCOS and MM in all the modules showed a strongly significant correlation, suggesting genes in these modules were highly correlated with PCOS (Figure 6). Subsequently, based on the threshold of |GS| > 0.8 and |MM| > 0.8, a total of 39, 34, 29, 85 and 3 DECs were screened from the blue, brown, grey, turquoise and yellow module, respectively (Table 1), indicating they may be hub DECs.

Figure 5 Identification of polycystic ovary syndrome-related modules by WGCNA. (A) determination of the soft threshold power; (B) clustering dendrograms of circRNAs; (C) heatmap of the module–trait relationship. r, correlation coefficient; p was corresponding p-value.

Figure 6 Scatter plots of gene significance (GS) vs module membership (MM) in the five significant modules. (A) blue module; (B) brown module; (C) grey module; (D) turquoise module; (E) yellow module.

Identification of Hub DECs That Function as miRNA and RBP Sponges

To link the DECs and DEGs, we constructed the circRNA-miRNA-mRNA ceRNA and circRNA-RBP regulatory networks. Both the circinteractome and starbase databases predicted 17 hub DECs identified in WGCNA modules could interact with nine DEMs. miRwalk database predicted 315 DEGs were target genes of these nine DEMs. These interaction pairs were used to construct a ceRNA network (Figure 7).

Figure 7 Construction of a circRNA-miRNA-mRNA ceRNA network. Genes in rhombus represent mRNAs; genes in circle represent circRNAs; genes in triangle represent miRNAs. Different colors for circRNAs represent different modules. Abbreviation: ceRNA, competing endogenous RNA.

Among all DEGs, only one (RBPMS, RNA binding protein, mRNA processing factor) and four (MBNL2, muscleblind-like splicing regulator 2; YTHDF2, YTH N6-methyladenosine RNA-binding protein 2; DKC1, dyskerin pseudouridine synthase 1; RBM47,-RNA binding motif protein 47) were overlapped with the RBPs in the circinteractome and starbase database, respectively. Thus, these RBPs were selected for predicting the interacted circRNAs. As a result, one upregulated DEC (hsa_circ_0099109) was predicted to bind with RBPMS; 53 DECs were predicted to bind with MBNL2. Thus, these interaction pairs were used to construct a DEC-RBP regulatory network (Figure 8).

Figure 8 Construction of a circRNA-RBP network. Genes in rhombus represent mRNAs; genes in circle represent circRNAs. Different colors for circRNAs represent different modules. Abbreviation: RBPs, RNA-binding proteins.

Due to the fact that RBPMS and MBNL2 were not enriched in function results, we overlapped the DECs in the above two regulatory networks to indirectly explain their roles in PCOS. The results showed that five DECs were shared (hsa_circ_0018108, hsa_circ_0054275, hsa_circ_0056196, hsa_circ_0063309, hsa_circ_0070987) between two regulatory networks. These findings revealed their related ceRNA and DEC-RBP axes may be particularly pivotal for PCOS, such as hsa_circ_0063309/hsa_circ_0054275/hsa_circ_0056196/hsa_circ_0018108/hsa_circ_0070987-hsa-miR-139-5p-CASP2/ FZD4/HMGCR, hsa_circ_0063309/hsa_circ_0054275/hsa_circ_0056196-hsa-miR-217-FZD4/TAB2/PEA15/RUNX2, hsa_circ_0063309-hsa-miR-346-FZD4/PEA15/RUNX2, hsa_circ_0063309-hsa-miR-421-PAK3/CASP2/RUNX2, hsa_circ_0056196/hsa_circ_0018108-hsa-miR-496-HMGCR/FZD4 and hsa_circ_0018108/hsa_circ_0054275/hsa_circ_0056196/hsa_circ_0063309/hsa_circ_0070987-MBNL2. Also, DEGs in these regulatory axes were demonstrated to be hub genes after analysis of the degree of each protein in the PPI network (PAK3, degree = 13; RUNX2, degree = 9; HMGCR; TAB2, degree = 8; FZD4, degree = 6; CASP2, degree = 5; PEA15, degree = 2) (Figure 9).

Figure 9 Construction of a PPI network for the genes in the circRNA-miRNA-mRNA ceRNA network. Abbreviations: PPI, protein–protein interaction; ceRNA, competing endogenous RNA.

Discussion

Using WGCNA method, we identified 190 DECs (that were included in five significant modules) associated with the development of PCOS. Six of them (hsa_circ_0063309, hsa_circ_0054275, hsa_circ_0056196, hsa_circ_0018108, hsa_circ_0070987, hsa_circ_0099109) were considered important by integrating with DEMs and DEGs to construct the circRNA-miRNA-mRNA ceRNA and DEC-RBP regulatory networks. Function analysis revealed these circRNAs were mainly involved in inflammation-mediated cell apoptosis and steroidogenesis by influencing the downstream target genes.

Although all these DECs were not reported previously, two of them encoded the same host genes with other circRNAs (hsa_circ_0063309, DDX17; hsa_circ_0070987, ELF2); while these reported circRNAs were demonstrated to play prominent roles in cell proliferation and apoptosis: Li et al found circDDX17 (hsa_circ_0002211) was significantly downregulated in colorectal cancer tissues. Silencing of circDDX17 promoted colorectal cancer cell proliferation and inhibited apoptosis.^5,26 Peng et al observed circDDX17 was lowly expressed in breast cancer tissues and cell lines. Overexpression of circDDX17 inhibited cell proliferation and promoted cell apoptosis in breast cancer.²⁷ In vivo experiments by Ren et al showed elevated expression of circDDX17 blocked colorectal cancer cell growth and elevated chemosensitivity of 5-Fu.²⁸ circRNA microarray and qRT-PCR detected circ-ELF2 (hsa_circRNA_103741) were significantly upregulated in glioma tissues. The collapse of circ-ELF2 led to growth arrest and accelerated cell apoptosis.²⁹ The host genes of our other identified DECs were also revealed to mediate the regulation of hormone synthesis, inflammation or cell apoptosis. Ectopic expression of metastasis associated 1 family member 3 (MTA3, host for hsa_circ_0054275) was reported to promote the cAMP-stimulated progesterone secretion and testosterone production.^30,31 Treatment with high concentration of insulin significantly attenuated the level of MTA3 protein inLeydig tumor cells.³⁰ Oroxyloside treatment markedly induced apoptosis in glioma cells through upregulating MTA3.³² Actin-related protein 3 (ACTR3, host for hsa_circ_0056196) was proved to amplify B cell receptor signaling and is required for maximal B cell activation and proliferation.³³ Estrogen receptor α-engaged macrophages initiated chronic inflammation by upregulating the expression of kinesin family member 5B (KIF5B, host for hsa_circ_0018108).³⁴ CPSF6 was revealed to be upregulated after HIV-1 infection.³⁵ These findings indicate that our identified circRNAs may be involved in PCOS by changing inflammation, apoptosis and steroidogenesis of cumulus cells.

Many studies have found that circRNAs contain miRNA response element and they could work as endogenous miRNA ‘sponges’ to bind to miRNAs and then exert negative regulation on the expression of downstream genes. circDDX17 was implicated to regulate miR-605-cell cycle-related factors (cyclin dependent kinase 1 and p21),²⁷ miR-346/phospholysine phosphohistidine inorganic pyrophosphate phosphatase,³⁶ miR-31-5p-KN motif and ankyrin repeat domains1²⁸ axes. Mechanistically, circ-ELF2 acted as a ceRNA for miR-510-5p to positively modulate the expression of mucin 15 at posttranscriptional level.²⁹ However, their roles remain not well understood. In this study, we predicted five circRNAs (hsa_circ_0063309, hsa_circ_0054275, hsa_circ_0056196, hsa_circ_0018108, hsa_circ_0070987) could sponge hsa-miR-139-5p to upregulate the expression of CASP2, FZD4 and HMGCR; hsa_circ_0063309/hsa_circ_0054275/hsa_circ_0056196 could regulate the expression of FZD4/TAB2/PEA15/RUNX2 by sponging hsa-miR-217; hsa_circ_0063309 also functioned as a ceRNA for hsa-miR-346 and hsa-miR-421 to influence the expression of FZD4/PEA15/RUNX2 and PAK3/CASP2/RUNX2, respectively; hsa_circ_0056196/hsa_circ_0018108 acted as a ceRNA to regulate HMGCR/FZD4 function by sponging hsa-miR-496. In these interaction pairs, only the negative regulatory relationships between some miRNAs and target genes were confirmed by wet experiments, including miR-139-5p-FZD4³⁷ and miR-217-RUNX2,^38,39 while the circRNA-miRNA and other miRNA-target gene interactions were newly identified. Nevertheless, the target DEGs in the above interaction pairs were revealed to be directly or indirectly related with PCOS by involving in inflammation, proliferation, apoptosis or steroidogenesis of granulosa cells. Zhang et al found RUNX2 was upregulated in cumulus cells from PCOS patients compared with controls. Silencing of RUNX2 promoted estradiol production by enhancing the expression of cytochrome P450 family 11 subfamily A member 1 (CYP11A1) and CYP19A1 which were rate-limiting enzymes for estradiol biosynthesis.¹¹ Compared with women without RUNX2, women with RUNX2 expression had lower number of follicles and peak estradiol levels.⁴⁰ In granulosa cells that had an upregulated expression of RUNX2 due to knockdown of CREB3/Luman gene, estradiol and progesterone synthesis were decreased, anti-apoptotic protein Bcl-2 was downregulated.⁴¹ The FZD member gene FZD3 was upregulated in cumulus cells from PCOS patients. This upregulation showed a tight association with estrogen deficiency. Overexpression of FZD3 compromised estrogen production in human granulosa cells.⁴² Inhibition of PAK3 prevented the secretion of testosterone, while stimulated oestradiol release.⁴³ Addition of PAK3-specific siRNAs reversed the decreased proliferation rate of granulosa cells induced by slit guidance ligand 2.⁴⁴ PCOS women possessed higher amount of PED/PEA-15 protein levels than that in controls. This increase was positively associated with high total testosterone.⁴⁵ Kaur et al identified a 1.86-fold increase in the expression of TAB2 in granulosa cells of PCOS non-insulin resistant patients compared with matched controls and function analysis showed this gene was inflammation-related.⁴⁶ This pro-inflammatory role of TAB2 was further validated in a previous study that showed TAB2 could activate nuclear factor κ-B to promote the production of interleukin-1.⁴⁷ mRNA expression for CASP2 was shown to be increased with the onset of apoptosis of granulosa cells.⁴⁸ HMGCR encodes a rate-limiting enzyme (HMG-CoA reductase) in the cholesterol biosynthetic pathway; while cholesterol is a precursor of androgen production. The upregulation of HMGCR was induced by inflammatory agents.⁴⁹ The addition of simvastatin reduced ovarian androgen production in PCOS women by blocking the expression of HMGCR and inflammatory mediators.^50,51 Although miRNAs were not studied in granulosa cells, the expressions of some miRNAs (miR-421,⁵² and miR-496,⁵³ and miR-346⁵⁴) were confirmed to be significantly downregulated in patients or mice with diabetes, a comorbidity of PCOS. Overexpression of miR-139-5p was shown to reduce the blood glucose level and alleviate the apoptosis of liver cells in diabetic mice.⁵⁵ In vitro experiments observed downregulation of miR-217 could induce human podocyte cells apoptosis.⁵⁶ In accordance with these studies, we also found all related miRNAs were downregulated and mRNAs were upregulated in cumulus cells of PCOS patients findings. These findings implied the possible functions of circRNAs in PCOS.

CircRNAs can also interact with RBPs. In this study, we predicted that the above five circRNAs (hsa_circ_0063309, hsa_circ_0054275, hsa_circ_0056196, hsa_circ_0018108 and hsa_circ_0070987) could interact with MBNL2; hsa_circ_0099109 could interact with RBPMS to promote their upregulated expression in PCOS. In line with our prediction, Che et al verified that RBPMS was a downstream target of lncRNA HUPCOS which was significantly elevated in granulosa cells to promote the accumulation of testosterone in PCOS patients.⁵⁷ Downregulation of MBNL2 was shown to enhance the ability of dendritic cells to inhibit virus replication and exert anti-inflammatory roles.⁵⁸ A significant decrease in the percentage of dendritic cells was found in patients with PCOS compared with control.⁵⁹ These roles of RBPMS and MBNL2 may also indirectly explain the functions of related circRNAs.

There are some limitations in this study. First, hub circRNAs were only preliminarily identified by analysis of high-throughput datasets with small sample size. Their expression and associations with PCOS patients’ pathophysiological characteristics (estrogen, luteinizing hormone, follicle-stimulating hormone testosterone, insulin resistance index, metaphase II oocytes, maturation rate, embryos, and fertilization rate) should be evaluated using a large number of clinical data collected from our hospital. Second, overexpression and knockdown experiments should be conducted in vitro and in vivo to explore whether these hub circRNAs play roles in inflammation, apoptosis and steroidogenesis of cumulus cells in PCOS and whether targeted regulation of them effectively alleviates the symptoms of PCOS. Third, the downstream mechanisms of circRNAs by functioning as miRNA and RBP sponges should be validated by dual-luciferase reporter gene assay, co-immunoprecipitation, gain and loss of function experiments. Fourth, the DEMs and DEGs that were not linked to circRNAs should also be given special attention to demonstrate their importance in the future.

In conclusion, our present study suggests dysfunction of these six circRNAs may be the underlying mechanism to induce the development of PCOS. These circRNAs may mediate the hormone imbalance and apoptosis of cumulus granulosa cells by interaction with miRNAs (hsa-miR-139-5p, hsa-miR-217, hsa-miR-346, hsa-miR-421, hsa-miR-496) to influence their downstream target genes (CASP2, FZD4, HMGCR, TAB2, PEA15, RUNX2, PAK3) or interaction with RBPs (RBPMS, MBNL2).