Identification of Three Potential circRNA Biomarkers of Polycystic Ovary Syndrome by Bioinformatics Analysis and Validation

Introduction

PCOS is a common endocrine and metabolic disease in reproductive women, which is characterized by polycystic ovaries, hyperandrogenemia and chronic anovulation.¹ However, the health risks of PCOS go well beyond these commonalities, showing symptoms that extend into the reproductive period and even menopause.² In addition, patients with PCOS are at increased risk of developing metabolic diseases with increasing age, such as obesity,³ cardiovascular disease⁴ and so on. However, the treatment of PCOS still remains a challenge. Therefore, it is of great significance to study its pathogenesis and search for more effective biomarkers and treatment strategies. Follicular fluid (FF) has proved to provide an important microenvironment for follicular development and oocyte maturation. It is the medium for bidirectional communication between oocytes and the surrounding cells, such as cumulus cells.⁵ The components of FF are very complex including proteins, RNAs and metabolites,^6–10 which are secreted by theca cells, granulosa cells and oocytes.⁶ Recent studies have shown that circRNAs in FF play an important role in PCOS.¹¹ Huang et al reported that depletion of exosomal circLDLR in follicle fluid derepressed miR-1294 function and inhibits estradiol production via CYP19A1 in PCOS.¹ Otherwise, circDDX10 in granulosa cells (GCs) which derived from human follicular fluid may participate in the regulation of ovarian function by affecting the proliferation and apoptosis of GCs and steroid hormone synthesis.¹²

CircRNA is a kind of RNAs with loop structures produced by abnormal transcript splicing.¹³ The study of circRNAs has been applied to many important fields spans over 30 years. In recent years, most of the research work on circRNAs has revealed its molecular functions and mechanisms in more detail, and its biological functions include translation, the regulation of gene transcription and miRNA sponge activity.¹⁴ Among them, the functions of circRNAs - competing sponge to endogenous RNA to regulate RNA expression and affect RNA function - have attracted considerable scientific attention.¹⁵ All these functions suggested that circRNAs may be potential novel diagnostic and prognostic biomarkers of disease.¹⁶ So far, it has been reported that circRNAs are extensively involved in the processes of many human diseases, including a variety of cancers,^17–19 neurological disorders,^20,21 cardiovascular disease,^22–24 immune disease²⁵ and so on. Similarly, many circRNAs have been proved to be regulate the process of PCOS, such as circASPH,²⁶ circLDLR¹ and circPUM1.²⁷ However, the mechanism of the involvement of much more circRNAs in PCOS remains unclear. Otherwise, the structure of circRNA is stable, not easy to degrade,²⁸ and highly expressed in ovaries,²⁹ which made it to be a proper diagnosis strategy for PCOS.

In the present study, we selected three datasets of PCOS from the GEO datasets, aimed to find novel circRNAs as potential biomarkers for the early diagnosis of PCOS. We also performed a circRNA-miRNA-mRNA regulatory network to identify circRNA that are critical to the biological processes of PCOS. Finally, qRT-PCR was used to validate our bioinformatics analyses. The workflow of our study is shown in Figure 1.

Figure 1 The workflow of the study.

Materials and Methods

Clinical Samples

In our study, a total of 16 PCOS patients and 16 healthy people were collected from Fujian Maternity and Child Health Hospital, Affiliated Hospital of Fujian Medical University. Our study complies with the Declaration of Helsinki. All the research protocol for this study was approved by the Ethics Committee of Fujian Provincial Maternity and Children’s Hospital. PCOS patients were diagnosed according to the revised Rotterdam European Society of Human Reproduction and Embryology/American Society for Reproductive Medicine Criteria (2004). The PCOS patients who fulfilled at least two of the following criteria were included: chronic oligo-ovulation or anovulation, androgen excess and polycystic ovaries. The control group contained patients undergoing IVF or ICSI due to male factors or tubal factors of infertility, who had regular menstrual cycles, normal ovary sonographs and normal ovulating. All of the patients had no history of drugs affecting glucose and lipid metabolism and were without any known medical conditions or diseases, including Cushing’s syndrome, congenital adrenal hyperplasia, androgen-secreting tumors and endometriosis. The follicular fluid was separated from the subscribers. The FF samples were stored in a refrigerator at −80°C for RNA extraction. Informed consent was obtained from all the participants.

Data Sets Collection and Differential Expression Analysis

The RNA expression profiles of PCOS were downloaded from the National Center of Biotechnology Information Gene Expression Omnibus (http: //www.ncbi.nlm.nih.gov/gds/). The GSE145296³⁰ dataset is composed of 12 circRNA expression data of cumulus granulosa cells, including 6 PCOS patients and 6 normal control individuals. For dataset GSE155489, we chose 8 granulosa cells samples for differential mRNA analysis. GSE138572³¹ dataset contains the microRNA profiles of 10 samples, including 5 PCOS patients and 5 controls. The R software package was used to process the downloaded files, and the unqualified data are converted and removed. After the gene expression matrix was merged, the ensemble IDs were converted to gene names. The data were calibrated, normalized by using the edgeR package in R software, and then log₂ transformed. The differentially expressed circRNAs, mRNAs and miRNAs were screened by Limma package, with a threshold of |log₂(fold change) [FC]| > 0.58 and P-value < 0.05.

Function Enrichment Analyses

The differentially expressed circRNAs related genes were enriched by GO and KEGG analysis. GO contains three terms: biological processes (BP), cellular components (CC) and molecular functions (MF). The p value < 0.05 indicated significant GO term enrichment among differentially expressed genes. DAVID database (version 6.8; http://david.abcc.ncifcrf.gov) was used for pathway analysis.³² P value < 0.05 was considered as statistically significant.

ceRNA Network Construction

The circular RNA Interactome (CircInteractome) (https://circinteractome.nia.nih.gov/) database was used to predict the miRNA-binding sites (MREs). Next, we used miRTarBase³³ (http://mirtarbase.mbc.nctu.edu.tw) and TargetScan³⁴ (http://www.targetscan.org) databases to predict the interactions between miRNAs and mRNAs. Only the mRNAs recognized by both databases were considered as candidate mRNAs, and were intersected with DE-mRNAs to screen the DE-mRNAs targeted by DE-miRNAs. The circRNA-miRNA-mRNA regulatory network was constructed using a combination of circRNA-miRNA pairs and miRNA-mRNA pairs. Finally, the network was visualized and mapped using Cytoscape v3.6.1.

Validation with qRT- PCR Method

The potential expression of circRNAs was confirmed by qRT-PCR. Total RNA was extracted from the follicular fluid using TRIzol Reagent (Invitrogen) according to the manufacturer’s protocol. RNA purity was measured using a NanoDrop Spectrophotometer (Nanodrop, Thermo Fisher Scientific). The total PCR reaction of circRNA was 20μL, including 6 μL of cDNA, 10 μL of SYBR Green, and 0.5 μL of primers for forward (10 μmol/L) and reverse (10 μmol/L). GAPDH as an internal reference. Primers sequences are summarized in Table 1. All PCR reactions were conducted in triplicate. The relative expression of circRNAs was calculated by 2^−ΔΔCt method.

Table 1 Specific RNAs Primers for qRT-PCR Analysis

Statistical Analysis

Data were analyzed using GraphPad Prism 8.0 (GraphPad software, San Diego, USA). The data are expressed as means ± SD. Unpaired t-test was used to distinguish the differences between the two groups. P value < 0.05 was considered statistically significant with 2-tailed tests. Receiver operating characteristic (ROC) curve was used to analyze the diagnostic value, including the area under the curve (AUC), sensitivity, and specificity.

Results

Identification of Differentially Expressed Genes

Three datasets GSE145296, GSE155489 and GSE138572 were downloaded from the GEO database, respectively. The volcano and heatmap showed differential expression profiles of the PCOS and healthy samples (Figure 1). A total of 1591 DE-circRNAs were screened from GSE145296 dataset, including 945 upregulated and 646 downregulated circRNAs (Figure 2A and D). The DE-mRNAs and DE-miRNAs of PCOS patients and healthy controls, obtained from datasets GSE155489 and GSE138572, respectively. We totally identified 1688 DE-mRNAs (739 upregulated and 949 downregulated) (Figure 2B and E) and 38 DE-miRNAs (18 upregulated and 20 downregulated) (Figure 2C and F).

Figure 2 Differentially expressed circRNA, mRNAs and miRNAs in polycystic ovary syndrome patient samples compared with healthy controls. (A) Volcano plots of all circRNAs in dataset GSE145296. (B) Volcano plots of all mRNAs in dataset GSE155489. (C) Volcano plots of all microRNAs in dataset GSE138572. (D) Heatmap of differentially expressed circRNAs in dataset GSE145296. (E) Heatmap of differentially expressed mRNAs in dataset GSE155489. (F) Heatmap of differentially expressed microRNAs in dataset GSE138572. Red dots indicate significantly upregulated genes, and green dots indicate significantly downregulated genes in (A–C). The red diamonds indicate upregulated genes and blue diamonds indicated downregulated genes in (D–F).

Functional Enrichment Analysis of DE-circRNAs in PCOS

A total of 1591 DE-circRNAs were identified. Among these genes, hsa_circ_0017962, hsa_circ_0056003 and hsa_circ_0032879 were most significant. In order to investigate the potential function of the DE-circRNAs. GO and KEGG enrichment analysis was used to explore the functions of DE-circRNAs host genes. The top ten most significant enrichment BP, CC and MF terms are shown in Figure 3A–C. From our results, the top 3 enriched BP terms are response to oxygen levels, response to decreased oxygen levels and response to hypoxia. This may indicate that PCOS is in an anoxic environment. The top three enriched CC terms included collagen-containing extracellular matrix, mitochondrial protein complex, and mitochondrial matrix. The top three enriched MF terms included extracellular matrix structural constituent, peptide antigen binding and collagen binding. Otherwise, a total of 52 pathways were enriched. The top ten pathways were shown in Figure 3D. Among them, we found metabolism-related pathways were significantly enriched including Valine, leucine and isoleucine degradation, Propanoate metabolism, Pyruvate metabolism, Carbon metabolism, Glyoxylate and dicarboxylate metabolism, Citrate cycle (TCA cycle), which may suggest an active metabolism in PCOS.

Figure 3 Functional analysis of differentially expressed circRNAs host genes. (A) The top 10 ranked biological process. (B) The top 10 ranked cellular component. (C) The top 10 ranked molecular function. (D) The top 10 ranked KEGG pathway. The deeper the red, the more significantly the items are enriched.

A circRNA-miRNA-mRNA Regulatory Network of PCOS

Subsequently, we constructed a circRNA-miRNA-mRNA regulatory network to reveal the connections of circRNA, miRNA and mRNA in PCOS (Figure 4), which composed with 16 nodes and 14 edges. As shown in Figure 4, there were only three circRNAs (hsa_circ_0075691, hsa_circ_0075692 and hsa_circ_0085997) and three miRNAs (hsa-miR-320c, hsa-miR-455-5p and hsa-miR-346) in the network. The basic characteristics of the 3 circRNAs are listed in Table 2. Their basic structural pattern is shown in Figure 5.

Table 2 Basic Characteristics of the Candidate circRNAs

Figure 4 The ceRNA network of circRNA–miRNA–mRNA in polycystic ovary syndrome. The network composed with 16 nodes and 14 edges. Circular nodes indicate DE-mRNAs, triangular nodes indicate DE-circRNAs and diamond nodes indicate DE-miRNAs. Upregulated DE-miRNAs were labeled in red and downregulated were labeled in green.

Figure 5 The basic structural patterns of the three DE-circRNAs. Structural patterns of the 3 DE-circRNAs by circRNADb (http://reprod.njmu.edu.cn/circrnadb/circRNADb.php). (A) Structure of hsa_circ_0075691. (B) Structure of hsa_circ_0075692. (C) None.

Validation the Expression of the circRNAs

To further evaluate the expression of the three circRNAs, a total of 16 PCOS patients and 16 healthy controls were collected to validate the analyses results by qRT-PCR. Consistent with the omics data, hsa_circ_0075691 and hsa_circ_0075692 were significantly upregulated (Figure 6A and B), hsa_circ_0085997 was downregulated (Figure 6C) between PCOS and controls. To assess the potential diagnosis value of the three circRNAs for PCOS, we further implemented a ROC curve analysis. We found that ROC curve of hsa_circ_0075691 showed a significant distinguishing efficiency with an AUC value of 0.89 (95% CI: 0.65–0.96, **p <0.01) (Figure 6D), an AUC value of 0.80 for hsa_circ_0075692 (95% CI: 0.78–1.00, **p <0.01) (Figure 6E) and an AUC value of 0.75 for hsa_circ_0085997 (95% CI: 0.57–0.93, **p <0.01) (Figure 6F). These indicated that circRNAs could be candidate biomarkers for PCOS diagnosis.

Figure 6 RT-qPCR analysis of hsa_circ_0075691, hsa_circ_0075692 and hsa_circ_0085997. (A–C) Validation of hsa_circ_0075691, hsa_circ_0075692 and hsa_circ_0085997 by RT-PCR (**p< 0.01). (D–F) The ROC curve of hsa_circ_0075691, hsa_circ_0075692 and hsa_circ_0085997 in distinguishing PCOS and normal controls.

Discussion

PCOS, with its high incidence among women of reproductive age and multiple side effects, has become a global health burden and brings unpleasant life experiences to patients.³⁵ The etiology of PCOS is complex with many genetic mutations, so the clinical management is faced with difficulties.³⁶ Therefore, it is urgent to find effective and efficient targets for the treatment of PCOS. With the development of high throughput sequencing, studies applying this technique for assessing PCOS have revealed the differential expression of transcripts in PCOS. Recently, the expression profile of circRNAs in polycystic ovary syndrome has been partially reported.^29,37 Circular RNAs are spliced endogenous RNAs that can potentially be used as diagnostic and prognostic biomarkers for disease.^38–40 Nevertheless, there were only a few studies have reported the diagnostic value of circRNAs in PCOS. Recent studies showed that circASPH²⁶ and hsa_circ_0118530⁴¹ may have the potential to be used as diagnostic biomarkers of PCOS. These findings give us a direction for our research. In this study, the differentially expressed circRNAs hsa_circ_0075691, hsa_circ_0075692 and hsa_circ_0085997 have not been previously reported in other diseases, suggesting that the specificity and effectiveness of these circRNA could be potential biomarkers of PCOS.

As we know, circRNAs have a large number of binding sites to sponge miRNAs, which could competitively bind to the miRNAs, reducing the inhibition of target genes. In the study, we used three different datasets to construct a circRNA-miRNA-mRNA network in PCOS. Studying the function of these DE-circRNAs is a very interesting project. Through the GO analysis, we found that the biological process Response to hypoxia was significantly enriched. Otherwise, Wang and Wu has reported that hypoxia inducible factor-1 (HIF-1) decreased in granulosa cells of polycystic ovary syndrome, while the levels of reactive oxygen species increased.⁴² Similarly, Zhao et al found the protein and mRNA level of HIF-1α was significantly downregulated of PCOS endometrial tissues than that in healthy control.⁴³ These may indicated that the HIF pathway is inhibited in the ovary of PCOS and that PCOS causes the ovary to be in a peroxidation state. In addition, pathway analysis showed many metabolism-related pathways were significantly enriched, such as Valine, leucine and isoleucine degradation, Propanoate metabolism, Pyruvate metabolism, Carbon metabolism, Glyoxylate and dicarboxylate metabolism. As we know, PCOS is one of the most common endocrine diseases in women, so endocrine disruption obviously causes the changes in various metabolic pathways. And our findings are consistent with previous reports.^44–46

CeRNA network reveals a new mechanism of RNA interaction, greatly expanding the functional genetic information in the human genome and playing important roles in pathological conditions.⁴⁷ In the present study, we aimed to constructed a ceRNA network to reveal potential RNAs regulation relationship in PCOS. Through our bioinformatics analyses, three circRNAs was met to the ceRNA network. To validate the expression of the three genes, more samples were used to make a validation set. ROC analysis showed good specificity of the candidate genes. CircRNA-miRNA-mRNA network implied that circRNA may function by sponging miRNAs. MiR-346 has been reported to be upregulated under hypoxia in renal cell carcinoma.⁴⁸ Otherwise, miR-346 was upregulated in primary sclerosing cholangitis and associated with TNF-α signaling pathway.⁴⁹ So, we speculate that hsa_circ_0085997 may absorb miR-346 to involved in inflammatory regulation in PCOS. However, the function of hsa-mir-320c and hsa-mir-455-5p were not reported in detail. The functions and accurate regulatory mechanisms of the three circRNAs were need more validation. In a word, integrated bioinformatics and validation, we found three potential circRNA biomarkers to diagnose PCOS, but there remains further study of the complex regulatory mechanisms in PCOS.

Conclusion

In conclusion, through integrated bioinformatics, a circRNA related ceRNA network was constructed. Three circRNAs (hsa_circ_0075691, hsa_circ_0075692 and hsa_circ_0085997) were associated with the diagnosis and treatment of PCOS by RT-PCR and ROC curve analysis. However, there are still some limits about the molecular function and regulatory mechanism of the circRNAs in this study. We will focus on this in the following work.