Expression Profiles of Circular RNA in Granulosa Cells from Women with Biochemical Premature Ovarian Insufficiency

Abstract

Aim: To identify the expression profiles and potential functions of circular RNAs (circRNAs) in granulosa cells (GCs) from women with biochemical premature ovarian insufficiency (bPOI). Patients & methods: CircRNAs microarray analysis was performed to GCs from 8 patients with bPOI and 8 control women, followed by qRT-PCR in 15 paired samples. CircRNA–miRNA networks and the prediction of their enriched signaling pathways were conducted by bioinformatics analysis. Results: A total of 133 upregulated and 424 downregulated circRNAs was identified in women with bPOI. We constructed circRNA–miRNA networks and found that the most predominantly enriched signaling pathways were the FoxO signaling pathway and cellular senescence. Conclusion: CircRNAs are differentially expressed in bPOI, which might contribute to the pathogenesis of bPOI.

Keywords::

• circular RNA
• granulosa cell
• microarray
• miRNA
• premature ovarian insufficiency

Premature ovarian insufficiency (POI) is a pathological condition of the ovarian reserve exhausting before the age of 40 years, which manifests with amenorrhea or oligomenorrhea, hypoestrogenism and elevated serum follicle stimulating hormone (FSH) concentration (>25 mIU/ml) [1]. The prevalence rate of POI is about 1%, depending on ethnicity [2]. Patients with POI possess increased risk of cardiovascular diseases [3], osteoporosis [4] and cognitive impairment [5]. They are even found to have increased total and cancer-specific mortality [6].

Ovarian function is thought to decline progressively and irreversibly both physiologically and pathologically. Previous studies describe three clinical stages of POI, including the occult stage, in which patients manifest with subfertility and normal FSH levels, the biochemical stage, in which patients have raised FSH levels but spontaneous menstruation, and the overt stage, which is the typical POI condition [7]. The development of POI may last for a long period. However, most patients are diagnosed at the overt stage, when fertility is severely diminished [8]. It has been demonstrated that patients with POI usually experience a period of ovarian function decline manifested with irregular menstruation before the follicles are completely depleted [9]. Therefore, it is necessary to identify potential biomarkers for early detection and diagnosis in POI.

Although the etiology of POI remains obscure, a poor primordial follicle pool or accelerated follicle depletion are considered the major mechanisms of POI [10–12]. Granulosa cells (GCs) play a crucial role in folliculogenesis and follicles atresia [13]. Recently, different expression profiles of miRNAs, shortened telomere length and diminished telomerase activity in GCs from biochemical POI (bPOI) patients have been reported [14,15]. These observations point to the idea that the dysfunction of GCs may contribute to the pathogenesis and development of POI.

Circular RNAs (circRNAs) are a class of circular noncoding RNAs (ncRNAs) which have no free 3′- or 5′-end after reverse splicing [16]. Compared to traditional linear RNAs, circRNAs are characterized by structural stability, sequence conservation and cell type-specific and tissue-specific specificity that may function as miRNA sponges, protein sponges, nuclear transcriptional regulators and protein scaffolds to regulate gene expression [17]. An increasing number of studies have revealed that circRNAs are misregulated in multiple human diseases, including cancer [18], cardiovascular diseases [19] and neurological disorders [20]. Aberrant circRNA expression has also been implicated in some ovarian diseases such as polycystic ovary syndrome (PCOS) [21–23] and ovarian endometriosis [24–26]. However, little is known about the expression profiles and function of circRNAs in POI.

In this study, we identified circRNA expression profiles of GCs from patients with bPOI and control patients with normal ovarian reserves by performing circRNA microarray analysis and assessed their potential role in the etiology of POI by bioinformatic analysis. Our findings indicate that circRNAs may play a crucial role in the etiology of POI and serve as promising diagnostic biomarkers for the early and accurate diagnosis of POI.

Materials & methods

Ethical approval

This study was approved by the Ethics Committee of Nanfang Hospital of Southern Medical University (NFEC-2017-197). Written informed consent was obtained from all patients.

Participants & samples

GCs for this study were collected from 23 patients with bPOI and 23 age- and BMI-matched women with normal ovarian reserves. Both bPOI and control patients underwent in vitro fertilization/intracytoplasmic sperm injection treatment at the Center for Reproductive Medicine, Department of Gynecology and Obstetrics in Nanfang Hospital, Southern Medical University, Guangzhou, People’s Republic of China from July 2017 to July 2018. The inclusion criteria for bPOI patients included: <40 years of age, basal FSH (on days 2–4 of menstrual cycle) ≥10 mIU/ml, bilateral ovarian antral follicle count (AFC) ≤7, anti-Mullerian hormone (AMH) ≤1.0 ng/ml and having spontaneous menstruation. The inclusion criteria for control women were as follows: <40 years of age, basal FSH <10 IU/ml, bilateral AFC >7, AMH ≥1.26 ng/ml, regular menstrual cycles occurring every 25–35 days and entered the in vitro fertilization/intracytoplasmic sperm injection program for male factor or tubal factor. The exclusion criteria of the two groups were: an abnormal karyotype, a history of other endocrine diseases such as polycystic ovary syndrome, hyperprolactinemia, hyperthyroidism and a history of radiotherapy, chemotherapy and ovarian operation.

Ovarian stimulation

The flexible gonadotrophin-releasing hormone antagonist ovarian hyperstimulation treatment protocol that we used has been previously described [27]. In brief, on days 2–4 of their natural menstrual cycle, ovarian stimulation was started by daily injection of recombinant FSH (Gonal F, Merck Serono, Italy) or highly purified FSH (Lishenbao, Livzon, China), followed by co-treatment with cetrorelix (Cetrotide, Merck Serono, Germany) in a daily dose of 0.25 mg sc., when at least one of the following criteria was reached: the leading follicle was ≥14 mm in diameter, the serum estradiol level was ≥600 pg/ml or the serum luteinizing hormone (LH) level was ≥10 mIU/ml. Gonadotrophin and gonadotrophin-releasing hormone antagonist were continued until the day of trigger. When at least two follicles were ≥18 mm or three follicles were ≥17 mm in diameter (the follicle was ≥17 mm in diameter for patients with only one mature follicle), 0.2 mg triptorelin (Diphereline, IPSEN, France) combined with individualized dosages (5000–10,000 IU) of human chorionic gonadotropin (hCG, Livzon, China) were administered to trigger final oocyte maturation. Oocyte pick-up was performed 35–36 h later.

GCs collection

GCs were isolated from follicular fluid via density centrifugation as previously described [28]. Briefly, follicle fluid was centrifuged at 400 g for 10 min and cell pellets were resuspended in phosphate buffer solution (HyClone, UT, USA). The cell pellets were then slowly layered on 50% Percoll density gradient media (GE Healthcare, Uppsala, Sweden) and then centrifuged at 400 g for 20 min to remove erythrocytes. If erythrocytes were not completely removed, erythrocyte lysis buffer (Sigma, MO, USA) was used. The collected GCs layer was washed three-times with phosphate buffer solution. The GCs were then further stored at -80°C in 1 ml RNAiso Plus (Takara, Dalian, China) until RNA extraction.

Microarray analysis

Microarray analysis was performed by Kang Chen Bio-tech (Shanghai, China). Eight GC samples from bPOI patients and eight samples from control women were randomly selected for the Arraystar Human circRNA Array V2 analysis. Briefly, total RNA from GCs was isolated following the standard protocol of RNAiso Plus and quantified using the NanoDrop ND-1000 (Thermo Fisher Scientific, DE, USA). Total RNA was then digested with RNase R (Epicentre, CA, USA) to remove linear RNA and the enriched labeled circRNA was hybridized onto the Arraystar Human circRNA Array V2 (8 × 15K, Arraystar, MD, USA). The Agilent Scanner G2505C (Agilent Technologies, CA, USA) and Agilent Feature Extraction software (version 11.0.1.1, Agilent Technologies) was used to scan and analyze acquired array images. R software limma package was used for data processing.

Quantitative reverse transcription PCR

To validate the microarray, qRT-PCR was used to examine the relative expression levels of circRNAs on the LightCycler 480 (Roche Molecular Biochemicals, Germany) in the rest 15 GCs samples from bPOI patients and 15 samples from control women. Divergent primers were designed for distinctive backsplicing sites of each circRNA. A PrimeScript RT™ reagent Kit with gDNA Eraser (Takara) was used for reverse transcription. qRT-PCR was performed using TB Green^® Premix Ex Taq™ II kit (Takara) and GAPDH was used as an internal control for quantification of target circRNAs. Analysis of relative circRNA expression levels was calculated using the 2^-ΔΔCT method. Primer sequences are listed in Supplemental Table 1.

Construction of circRNA–miRNA network

Potential miRNA response elements (MREs) on circRNAs were predicted with Arraystar’s home-made miRNA target prediction software based on TargetScan and miRanda. The networks between circRNAs and the predicted top five targeted miRNAs were constructed by Cytoscape 3.4.0. TargetScan and Mirdb databases were used to predict the target genes of miRNAs. The overlapping results predicted by both TargetScan and Mirdb were included for further analysis.

Luciferase reporter assay

To verify miRNA-circRNA interactions, luciferase reporter assay was used to detect binding between hsa_circ_103903 and miRNAs (hsa-let-7c-3p, hsa-miR-18a-5p, hsa-miR-18b-3p, hsa-miR-323a-3p and hsa-miR-1301-3p) in a granulosa cell line, KGN cell. The full-length of hsa_circ_103903 3′UTR sequence was amplified and cloned into a pmiR-RB-ReportTM vector (Ribobio, Guangzhou, China), which contains renilla luciferase gene and firefly luciferase gene. MiRNA mimics for hsa-let-7c-3p, hsa-miR-18a-5p, hsa-miR-18b-3p, hsa-miR-323a-3p and hsa-miR-1301-3p were obtained from Ribobio. Each miRNA mimic or negative control oligonucleotide was cotransfected with pmiR-RB-Report vector with or without the sequence of hsa_circ_103903. The luciferase activity was measured after transfection by Dual-Glo^® Luciferase Assay System (Promega, MA, USA). Firefly luciferase activity was normalized against Renilla luciferase activity for each transfected well.

Kyoto Encyclopedia of Genes and Genomes pathway analysis

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed to annotate the function of miRNAs targeted genes for the top ten upregulated circRNAs and the enrichment score was shown as the -log (p-value).

Statistical analyses

Statistical evaluation was performed using SPSS version 20.0 (IBM, IL, USA). Data was presented as the means ± standard deviation, unless otherwise stated. Continuous variables were compared with the Student’s t-test or Mann–Whitney U test for independent samples, according to the normality of their distribution. Correlations between circRNAs and patients’ clinical characteristics were analyzed by Pearson’s rank correlation. A p-value of < 0.05 was considered statistically significant.

Results

Expression profiles of circRNAs in GCs of bPOI patients

The clinical characteristics of 23 patients with bPOI and 23 controls are summarized in Table 1. A total of 6522 upregulated and 5574 downregulated differentially expressed circRNAs were detected. Of these, 133 upregulated and 424 downregulated circRNAs were identified to be significantly differentially expressed by a fold change ≥2.0 and p < 0.05 in GCs of eight bPOI patients compared with eight control patients. Hierarchical clustering analysis provided an overview of circRNA expression (Figure 1A). A volcano plot and scatter plot also demonstrated the differential expression of circRNAs between two groups (Figure 1B & C). The characteristics of the top 20 upregulated and 20 downregulated circRNAs are shown in Table 2.

Table 1. Clinical characteristics of patients with biochemical premature ovarian insufficiency and controls.

Variable	bPOI (n = 23)	Control (n = 23)	p-value
Age (year)	32.22 ± 4.86	32.13 ± 4.20	0.949
BMI (kg/m^2)	20.29 ± 2.58	21.02 ± 2.06	0.307
AMH (ng/ml)	0.56 ± 0.31	2.99 ± 1.56	<0.001
AFC	4.29 ± 1.94	12.64 ± 2.53	<0.001
Basal FSH (mIU/ml)	15.93 ± 6.04	6.79 ± 0.96	<0.001
Basal LH (mIU/ml)	6.43 ± 3.47	5.49 ± 1.59	0.246
Basal FSH/LH	2.87 ± 1.34	1.33 ± 0.39	<0.001
Basal estradiol (pg/ml)	36.12 ± 30.31	34.90 ± 12.15	0.860
Basal progesterone (ng/ml)	0.55 ± 0.41	0.49 ± 0.22	0.495
Number of oocytes retrieved	3.04 ± 1.82	14.05 ± 6.72	<0.001

Data are presented as the mean ± SD; p < 0.05 is considered to be statistically significant.

AFC: Antral follicle count; bPOI: Biochemical premature ovarian insufficiency; AMH: Anti-Müllerian hormone; FSH: Follicle-stimulating hormone; LH: Luteinizing hormone; SD: Standard deviation.

Figure 1. Differential expression of circular RNAs in granulosa cells from eight patients with biochemical premature ovarian insufficiency and eight control patients.

(A) Hierarchical clustering analysis of differently expressed circular RNAs (circRNAs; fold change ≥ 2, p < 0.05). Red = upregulated, green = downregulated. (B) Volcano plot visualizing the differentially expressed circRNAs. The red point in the plot represents the differentially expressed circRNAs (fold change ≥ 2, p < 0.05). (C) Scatter plot visualizing the differentially expressed circRNAs. The circRNAs above the top green line and below the bottom green line indicate more than twofold change of circRNAs.

circRNA: Circular RNA.

Table 2. Characterization of the top 20 upregulated/downregulated circular RNAs in patients with biochemical premature ovarian insufficiency.

CircRNA	Regulation	Fold change	p-value	Gene symbol	Chromosome	Strand	Start	End
hsa_circRNA_104503	Up	4.3348856	0.031039861	TRIM24	chr7	+	138203933	138255748
hsa_circRNA_003785	Up	4.3191191	0.016443185	XLOC_012658	chr18	+	32902351	32907027
hsa_circRNA_104616	Up	4.1071068	0.0218736	SPIDR	chr8	+	48320422	48353104
hsa_circRNA_103903	Up	3.9442081	0.023707735	PAPD4	chr5	+	78952780	78964851
hsa_circRNA_101258	Up	3.9246349	0.020856298	VWA8	chr13	-	42385360	42393522
hsa_circRNA_036088	Up	3.8659703	0.029294692	RPLP1	chr15	+	69745158	69747884
hsa_circRNA_100236	Up	3.6513049	0.039709019	TMEM59	chr1	-	54506428	54509198
hsa_circRNA_402819	Up	3.4858104	0.030949297	GLB1	chr3	-	33114035	33118729
hsa_circRNA_001153	Up	3.4033414	0.028277148	MARCH1	chr4	+	164514761	164515093
hsa_circRNA_104299	Up	3.3711272	0.009393462	USP42	chr7	+	6150755	6155154
hsa_circRNA_103689	Up	3.339481	0.024133026	PTPN13	chr4	+	87693930	87696805
hsa_circRNA_405558	Up	3.3032197	0.040345806	SSH2	chr17	-	28120137	28120955
hsa_circRNA_403691	Up	3.296719	0.035614529	LOC101927768	chr6	+	114313056	114317224
hsa_circRNA_002131	Up	3.21306	0.035508723	BNIP3L	chr8	+	26248758	26265892
hsa_circRNA_001729	Up	3.2130321	0.012576279	ZNF646	chr16	-	31093273	31093358
hsa_circRNA_059571	Up	3.1824205	0.020124964	RALGAPA2	chr20	-	20553518	20620544
hsa_circRNA_407331	Up	3.1416441	0.043787336	ENOX2	chrX	-	129917520	129943235
hsa_circRNA_101115	Up	3.1134423	0.048630939	METAP2	chr12	+	95879654	95897937
hsa_circRNA_001065	Up	3.0765417	0.000857261	GYPC	chr2	-	127453811	127454077
hsa_circRNA_083948	Up	3.0694227	0.027726197	EIF4EBP1	chr8	+	37888019	37917883
hsa_circRNA_405965	Down	50.0742673	0.000388304	CCNT2-AS1	chr2	-	135625089	135626611
hsa_circRNA_000102	Down	34.546511	0.000170047	AKNAD1	chr1	-	109479800	109479932
hsa_circRNA_008389	Down	34.065194	0.000321687	DNAJC11	chr1	-	6727768	6741097
hsa_circRNA_402094	Down	18.2204676	0.00030518	HKR1	chr19	+	37838091	37838774
hsa_circRNA_030448	Down	16.1072968	0.000363936	LMO7	chr13	+	76301164	76415337
hsa_circRNA_004077	Down	13.7486624	0.000233947	VAT1L	chr16	+	77850817	77859358
hsa_circRNA_400027	Down	9.1762641	0.001987192	SNURF-SNRPN	chr15	+	25325262	25326442
hsa_circRNA_403482	Down	9.0520716	0.000680145	RBM22	chr5	-	150076352	150078193
hsa_circRNA_005198	Down	8.4638383	0.000280205	PARP4	chr13	-	25072253	25077915
hsa_circRNA_000542	Down	8.1853461	0.000344817	ARID4A	chr14	+	58785259	58796887
hsa_circRNA_000407	Down	8.1180835	0.000295619	SMARCC2	chr12	+	56558216	56558309
hsa_circRNA_000499	Down	7.964529	0.001099349	RBM26	chr13	+	79940776	79940857
hsa_circRNA_092523	Down	7.9593586	0.000980014	PANX2	chr22	+	50618381	50618469
hsa_circRNA_001405	Down	7.734526	0.000819786	PREX1	chr20	-	47316514	47316617
hsa_circRNA_001937	Down	7.4899568	0.002365583	CHD9	chr16	+	53155459	53155541
hsa_circRNA_001587	Down	7.0488002	0.000427757	XLOC_001374	chr2	+	19042277	19042456
hsa_circRNA_001481	Down	6.9676678	0.000368557	EMB	chr5	-	49694940	49707217
hsa_circRNA_103670	Down	6.9533279	0.037835328	CNOT6L	chr4	-	78694234	78697546
hsa_circRNA_100719	Down	6.8982817	0.000930526	DOCK1	chr10	+	128768965	128836080
hsa_circRNA_043602	Down	6.8486592	0.000239464	JUP	chr17	-	39679868	39925459

Validation of circRNAs by qRT-PCR

To verify the microarray analysis results, the expression levels of two upregulated circRNAs (hsa_circ_003785 and hsa_circ_103903) and two downregulated circRNAs (hsa_circ_008389 and hsa_circ_103670) were measured by qRT-PCR in GCs from another 15 patients with bPOI and 15 control patients. Consistent with the microarray data, the expression levels of hsa_circ_003785 (Figure 2A) and hsa_circ_103903 (Figure 2B) were significantly higher in the bPOI group and the expression of hsa_circ_008389 (Figure 2C) was significantly lower compared with the control group. However, no significant difference was observed in the expression of hsa_circ_103670 (Figure 2D) between the two groups.

Figure 2. Comparison of the expression of circular RNAs between patients with biochemical premature ovarian insufficiency and controls using quantitative real-time PCR.

(A) has_circ_003785 (B) has_circ_103903 (C) has_circ_008389 (D) has_circ_103670. *p < 0.05.

bPOI: Biochemical premature ovarian insufficiency.

Diagnostic value of novel circRNAs in bPOI

A receiver operating characteristic (ROC) curve analysis was used to analyze the expression levels of hsa_circ_003785, hsa_circ_103903 and hsa_circ_008389 in distinguishing patients with bPOI from control individuals (Figure 3). The highest area under the ROC curve (AUC) was found for hsa_circ_008389 (AUC: 0.784; 95% CI: 0.615–0.954; p = 0.008), followed by hsa_circ_103903 (AUC: 0.764; 95% CI: 0.593–0.936; p = 0.014) and hsa_circ_003785 (AUC: 0.729; 95% CI: 0.543–0.915; p = 0.033).

Figure 3. Receiver operating characteristic curve.

Receiver operating characteristic curve analysis of the expression levels of has_circ_003785 (red line), has_circ_103903 (green line) and has_circ_008389 (blue line) in biochemical premature ovarian insufficiency.

We further analyzed the associations between the expression levels of these three circRNAs and the clinical features of patients. As shown in Supplementary Table 2, significant correlations between the basal FSH/LH and the expression of both has_circ_003785 (r = 0.460; p = 0.012) and has_circ_103903 (r = 0.534; p = 0.003) were observed. Moreover, statistical analysis indicated correlations were observed between the expression of has_circ_008389 and the values of AMH (r = 0.461; p = 0.012) and AFC (r = 0.405; p = 0.032).

CircRNA–miRNA networks & pathway analysis

It has been widely proposed that some circRNAs can function as miRNA sponges by competing with miRNAs to modulate the translation of target genes [29]. Therefore, we built competing endogenous RNA networks between the circRNAs and the corresponding top five targeted miRNAs (Figure 4A). Interestingly, 19 miRNAs were matched with more than one circRNA in the networks. Among them, three downregulated circRNAs (has_circ_000542, has_circ_005189 and has_circ_030448) harbored miR-4778-3p MRE and three circRNAs (has_circ_001587, has_circ_103903 and has_circ_10115) harbored miR-1301-3p MRE.

Figure 4. The circular RNA–miRNA networks.

(A) The network comprising the top 20 upregulated (red nodes) and top 20 downregulated (deep blue nodes) circular RNAs and their target miRNAs (blue nodes) (miRNAs matched with more one circular RNAs are marked in yellow). (B) Binding sites of hsa-let-7c-3p, hsa-miR-18a-5p, hsa-miR-18b-3p, hsa-miR-323a-3p and hsa-miR-1301-3p in 3′UTR of has_circ_103903. (C) Hsa-let-7c-3p, hsa-miR-18a-5p, hsa-miR-18b-3p, hsa-miR-323a-3p and hsa-miR-1301-3p mimic or negative control oligonucleotide was co-transfected with pmiR-RB-Report vector with or without the 3′UTR sequence of has_circ_103903 to evaluate the effect of those miRNAs on the luciferase activities by luciferase reporter assay. *p < 0.05; **p < 0.01; ***p < 0.001.

circRNA: Circular RNA.

According to the MREs analysis, hsa_circ_103903 was observed to act as sponge for hsa-let-7c-3p, hsa-miR-18a-5p, hsa-miR-18b-3p, hsa-miR-323a-3p and hsa-miR-1301-3p. The positions of putative binding sites of these miRNAs in hsa_circ_103903 were showed (Figure 4B). The results of luciferase reporter assay revealed a significant decrease of luciferase intensity upon the co-transfection of the vector-hsa_circ_103903 and hsa-let-7c-3p, hsa-miR-18a-5p, hsa-miR-18b-3p, hsa-miR-323a-3p and hsa-miR-1301-3p mimics compared with co-transfected vector without the sequence of hsa_circ_103903 or miR-NC mimic (Figure 4C). The results showed that hsa-let-7c-3p, hsa-miR-18a-5p, hsa-miR-18b-3p, hsa-miR-323a-3p and hsa-miR-1301-3p directly bind and interact with hsa_circ_103903.

In order to explore potential pathways associated with the pathogenesis of POI, KEGG pathway analysis was conducted. Figure 5A showed the top ten relevant signaling pathways, including the FoxO signaling pathway, cellular senescence, PI3K-Akt signaling pathway and MAPK signaling pathway. Among these, we focused on the FoxO signaling pathway and cellular senescence, which had the highest enrichment scores (Figure 5B & C).

Figure 5. Prediction of miRNA-mediated signaling pathways.

(A) KEGG pathway analysis for the top ten upregulated circRNA-targeted genes. The vertical axis shows the annotated functions of the target genes. The horizontal axis shows the enrichment score (-log10 transformed p-value). (B) Mapping of the FoxO signaling pathway. (C) Mapping of cellular senescence pathway. Yellow marked nodes are associated with differentially expressed circRNAs.

circRNA: Circular RNA; KEGG: Kyoto encyclopedia of gene and genome.

Discussion

In this study, we investigated the circRNA expression profiles of GCs from patients with bPOI and analyzed the potential functions of misregulated circular RNA. We found a total of 133 upregulated and 424 downregulated circRNAs in women with bPOI and constructed circRNA–miRNA networks between the top 20 upregulated and top 20 downregulated circRNAs. The expression levels of hsa_circ_003785 and hsa_circ_103903 were upregulated and hsa_circ_008389 was downregulated in the bPOI group, consistent with qRT-PCR. We also explored further functions of the abnormal circRNAs. Our findings reveal new perspectives into the underlying mechanisms of POI.

POI is a highly heterogeneous condition with unclear pathology. The etiology studies of POI usually focus on factors, including genetic factors [30] and immune factors [31]. Currently, ncRNAs, including miRNAs and long noncoding RNAs (lncRNAs), are widely studied due to their epigenetic roles in various human diseases. The relationships between ncRNAs and POI have also been explored. For instance, miR-22-3p is a downregulated miRNA and is negatively associated with serum FSH in POI patients [32]. MiR-379-5p has been identified to be significantly upregulated in GCs from bPOI patients, which suppresses cell proliferation and impairs DNA repair function [14]. Previous studies also show that some misregulated lncRNAs influence cell proliferation of GCs [28,33]. CircRNAs are a type of ncRNAs discovered in eukaryotic cells as early as 1979, which were initially considered aberrant splicing byproducts with little functional potential [34]. New interest in circRNA research has sparked since the rapid development of microarray and RNA sequencing of ncRNAs. Several studies have proven that changes in circRNAs expression participate in ovarian development [35], ovarian aging [36,37] and aforementioned ovarian diseases [21–26]. In this study, we explored the alteration of circRNAs expression profiles in GCs isolated from bPOI patients compared with control patients. To our knowledge, this is the first study to identify the circRNAs expression profiles in the GCs of bPOI patients.

Due to their intrinsic circular characteristics, circRNAs can be more stable both inside cells and in extracellular plasma than linear RNAs. Thus, disease-associated circRNAs are considered promising diagnostic biomarkers [17]. In this study, we demonstrated that hsa_circ_003785 and hsa_circ_103903 were upregulated and hsa_circ_008389 was downregulated in the GCs of patients with bPOI by qRT-PCR (Figure 2). Spearman’s rank correlation indicates that both has_circ_003785 and has_circ_103903 were positively correlated with basal FSH/LH and has_circ_008389 is positively correlated with AMH values and AFC. The results of ROC curve analysis also showed that these circRNAs are potential diagnostic markers in POI.

The majority of circRNAs are derived from exons of protein-coding genes [38]. Thus, these circRNAs may alter gene expression outcomes by affecting the splicing of their precursor transcripts [39]. Hsa_circ_104616, an upregulated circRNA in GCs from bPOI patients, is derived from SPIDR gene, which has been reported to encode a protein involved in homologous DNA repair and play a role in autosomal recessive familial POI [40]. Unfortunately, we were unable to verify the altered expression of hsa_circ_104616 by qRT-PCR due to the failure to design a pair of specific primers. The expression of hsa_circ_008389 was confirmed by qRT-PCR and its precursor transcript is DNAJC11 gene, which encodes a mitochondrial protein. It has been reported that damaged mitochondrial function may be responsible for POI [41]. The hsa_circ_103903 is derived from PAPD4 gene, encoding a noncanonical poly(A) polymerase that is involved in gene expression regulation by post-transcriptional polyadenylation [42]. There is also some evidence supporting the idea that alterations in polyadenylation may contribute to POI [43]. Therefore, our results indicate that aberrant circRNAs may influence the etiology of POI via regulating their precursor transcripts.

CircRNAs have gained more attention because of their function as miRNA sponges by competitively binding to MREs. We constructed a circRNA–miRNA interaction networks in GCs from patients with bPOI. Among these miRNAs, has_miR-1301-3p, which matched with three circRNAs (hsa_circ_103903, hsa_circ_101115 and hsa_circ_001587), has been revealed to have a potential function on ovarian aging [36]. Another miRNA, hsa-miR-361-3p, has been found to match with two circRNAs (hsa_circ_001405 and hsa_circ_000499). Studies have also indicated hsa-miR-361-3p may participate in multiple reproductive processes, including follicle development [44], apoptosis of spermatogonia [45] and FSH synthesis regulation in pituitary cells [46]. Moreover, KEGG analysis based on miRNA-targeted genes demonstrates that the FoxO signaling pathway is the most predominantly enriched signaling pathway. Multiple components of the FoxO signaling pathway have been reported to be involved in the apoptosis of GCs and the pathogenesis of POI [47,48]. FOXO3A is regarded as a critical regulator and suppressor of ovarian follicle activation and female Foxo3a^–/– mice exhibit classic signs of hypergonadotropic hypogonadism secondary to POI [49]. Cellular senescence is another highly relevant signaling pathway shown by the results of KEGG analysis. As we know, POI is a condition of premature and pathological ovarian senescence. Previous studies have demonstrated accelerated telomere shortening and reduced telomerase activity, which are well-known senescence-linked factors, may contribute to follicular depletion in women with POI [15,50]. Therefore, we speculate that circRNAs may also play potential roles in POI via circRNA–miRNA-targeted gene regulation and its underlying mechanism requires further investigation.

Limitations of our study should be acknowledged. First, an important limitation of the study is the relatively small sample size, which is due to the difficulty in recruiting a large number of patients with bPOI and to collect enough GCs from them. Thus, our findings should be verified in a large-sample study. Second, GCs used in this study were stimulated by exogenous gonadotrophin; however, the circRNA expression profiles may be altered after ovarian hyperstimulation treatment. Finally, our study provides preliminary data on the mechanisms of circRNAs in POI and future studies should include the elucidation of how circRNA–miRNA networks participate in the etiology of POI.

Conclusion

In summary, this study reveals the circRNA expression profiles in GCs from bPOI patients by circRNA microarray analysis. The abnormal expression of circRNAs and their targeted miRNAs may be significantly related to the pathogenesis and development of POI. The findings of aberrantly expressed circRNAs in GCs may bring new perspectives on mechanisms underlying the etiology of POI.

Future perspective

CircRNA is a new research hotpot in the field of noncoding RNA, which has been reported to play an important role in multiple human diseases. The variety of biological mechanisms and special circular structure of circRNA give it talent to serve as a high-profile gene regulator and a promising diagnostic biomarker. However, the role of circRNA in reproductive endocrine disorders, especially in premature ovarian insufficiency, is still unclear. Further researches on how circRNA contribute to the pathogenesis and development of premature ovarian insufficiency should be conducted. The underlying mechanisms of circRNA in the proliferation, differentiation and apoptosis of ovarian granulosa cells are worth studying as well.

Summary points

• The circular RNA (circRNA) expression profiles in granulosa cells from biochemical premature ovarian insufficiency (bPOI) patients are different than those from women with normal ovarian reserves.

• A total of 133 upregulated and 424 downregulated circRNAs were identified as differentially expressed in bPOI.

• The expression levels of hsa_circ_003785 and hsa_circ_103903 were upregulated and hsa_circ_008389 was downregulated in the bPOI group, consistent with qRT-PCR.

• The results of receiver operating characteristic curve analysis showed that hsa_circ_003785, hsa_circ_103903 and hsa_circ_008389 are potential diagnostic markers in bPOI.

• Spearman’s rank correlation indicates that both has_circ_003785 and has_circ_103903 were positively correlated with basal follicle-stimulating hormone/luteinizing hormone and has_circ_008389 is positively correlated with anti-Mullerian hormone values and antral follicle count.

• Competing endogenous RNA networks were built between the circRNAs and the corresponding top five targeted miRNAs.

• CircRNA–miRNA interactions was confirmed by luciferase reporter assay.

• Kyoto Encyclopedia of Genes and Genomes analysis based on miRNA-targeted genes demonstrates that the FoxO signaling pathway and cellular senescence is the most predominantly enriched signaling pathway.

author-contributions

X-Y Zhou performed conception and design, acquisition of data, analysis and interpretation of data and drafting the article. S-L Chen performed conception and design of the study and revised the article. Y Li and Y-D Liu performed acquisition of data and drafting of the article. J Zhang, X Chen, J Zhe, Q-Y Zhang and Y-X Chen performed analysis and interpretation of data and revised the article. All authors gave their final approval of the version to be published.

Acknowledgments

We gratefully acknowledge all the patients who agreed to participate in the study and all the staff of the Center for Reproductive Medicine, Department of Gynecology and Obstetrics, Nanfang Hospital, Southern Medical University for their support and cooperation.

Ethical conduct of research

This study was approved by the Ethics Committee of Nanfang Hospital of Southern Medical University (NFEC-2017-197). Written informed consent was obtained from all patients.

Supplementary data

Financial & competing interests disclosure

This study was funded by the National Key Research & Developmental Program of China (2017YFC1001100), the National Natural Science Foundation of China (81671524), the Guangdong Science and Technology Department (2016A020218009), the Clinical Research Startup Program of Southern Medical University by High-level University Construction Funding of Guangdong Provincial Department of Education (LC2016ZD010), and the Clinical Research Program of Nanfang Hospital, Southern Medical University (2018CR016). The authors have no other relevant affiliations or financial involvement with any organization or entity with a financial interest in or financial conflict with the subject matter or materials discussed in the manuscript apart from those disclosed.

No writing assistance was utilized in the production of this manuscript.

Additional information

Funding

This study was funded by the National Key Research & Developmental Program of China (2017YFC1001100), the National Natural Science Foundation of China (81671524), the Guangdong Science and Technology Department (2016A020218009), the Clinical Research Startup Program of Southern Medical University by High-level University Construction Funding of Guangdong Provincial Department of Education (LC2016ZD010), and the Clinical Research Program of Nanfang Hospital, Southern Medical University (2018CR016). The authors have no other relevant affiliations or financial involvement with any organization or entity with a financial interest in or financial conflict with the subject matter or materials discussed in the manuscript apart from those disclosed. No writing assistance was utilized in the production of this manuscript.