Transcriptomic landscape of granulosa cells and peripheral blood mononuclear cells in women with PCOS compared to young poor responders and women with normal response.

STUDY QUESTION: Are transcriptomic profiles altered in ovarian granulosa cells (GCs) and peripheral blood mononuclear cells (PBMNCs) of women with polycystic ovary syndrome (PCOS) compared to young poor responders (YPR) and women with normal response to ovarian stimulation? SUMMARY ANSWER: RNA expression profiles in ovarian GCs and PBMNCs were significantly altered in patients with PCOS compared with normoresponder controls (CONT) and YPR. WHAT IS KNOWN ALREADY: PCOS is characterised by a higher number of follicles at all developmental stages. During controlled ovarian hyperstimulation, PCOS women develop a larger number of follicles as a result of an exacerbated response, with an increased risk of ovarian hyperstimulation syndrome. Despite the number of developing follicles, they are often heterogeneous in both size and maturation stage, with compromised quality and retrieval of immature oocytes. Women with PCOS appear to have a longer reproductive lifespan, with a slightly higher menopausal age than the general population, in addition to having a higher antral follicular count. As a result, the ovarian follicular dynamics appear to differ significantly from those observed in women with poor ovarian response (POR) or diminished ovarian reserve. STUDY DESIGN, SIZE, DURATION: Transcriptomic profiling with RNA-sequencing and validation using quantitative reverse transcription PCR (qRT-PCR). Women with PCOS (N = 20), YPR (N = 20) and CONT (N = 20). Five patients for each group were used for sequencing and 15 samples per group were used for validation. PARTICIPANTS/MATERIALS, SETTING,
METHODS: PCOS was defined using the revised Rotterdam diagnostic criteria for PCOS. The YPR group included women <35 years old with <4 mature follicles (at least 15 mm) on the day of the trigger. According to internal data, this group represented the bottom 15th percentile of patients' responses in this age group. It was consistent with Patient-Oriented Strategies Encompassing Individualize D Oocyte Number (POSEIDON) criteria for POR (Group 3). The young CONT group included women <35 years without PCOS or anovulation, who developed >14 mature follicles (at least 15 mm on transvaginal ultrasound). According to internal data, a threshold of >14 mature follicles was established to represent the top 25% of patients in this age group in this clinic.Overall, n = 60 GCs and PBMNCs samples were collected and processed for total RNA extraction. To define the transcriptomic cargo of GCs and PBMNCs, RNA-seq libraries were successfully prepared from samples and analysed by RNA-seq analysis. Differential gene expression analysis was used to compare RNA-seq results between different groups of samples. Ingenuity pathway analysis was used to perform Gene Ontology and pathways analyses. MAIN RESULTS AND THE ROLE OF CHANCE: In PBMNCs of PCOS, there were 65 differentially expressed genes (DEGs) compared to CONT, and 16 compared to YPR. In GCs of PCOS, 4 genes showed decreased expression compared to CONT, while 58 genes were differentially expressed compared to YPR. qRT-PCR analysis confirmed the findings of the RNA-seq. The functional enrichment analysis performed revealed that DEGs in GCs of PCOS compared to CONT and YPR were prevalently involved in protein ubiquitination, oxidative phosphorylation, mitochondrial dysfunction and sirtuin signaling pathways. LARGE SCALE DATA: The data used in this study is partially available at Gene Ontology database. LIMITATIONS, REASONS FOR CAUTION: The analysis in PBMNCs could be uninformative due to inter-individual variability among patients in the same study groups. Despite the fact that we considered this was the best approach for our study's novel, exploratory nature. WIDER IMPLICATIONS OF THE
FINDINGS: RNA expression profiles in ovarian GCs and PBMNCs were altered in patients with PCOS compared with CONT and YPR. GCs of PCOS patients showed altered expression of several genes involved in oxidative phosphorylation, mitochondrial function and sirtuin signaling pathways. This is the first study to show that the transcriptomic landscape in GCs is altered in PCOS compared to CONT and YPR. STUDY FUNDING/COMPETING INTEREST(S): This study was partially supported by grant PI18/00322 from Instituto de Salud Carlos III, and European Regional Development Fund (FEDER), 'A way to make Europe' awarded to S.H. M.C., S.H., S.T., L.R., M.R., I.R., A.P. and R.C. declare no conflict of interests concerning this research. E.S. is a consultant for and receives research funding from the Foundation for Embryonic Competence. TRIAL REGISTRATION NUMBER: N/A.

Introduction

After being recognized as a clinical entity more than 85 years ago by Stein and Leventhal (Azziz, 2006), there has been several attempts to refine the clinical diagnostic criteria for polycystic ovary syndrome (PCOS). Among the widely accepted definitions, the first one was the 1991 National Institute of Health (NIH) consensus definition which required oligo/amenorrhea and hypernadrogenemia/hyperandorenism (symptoms associated with response to androgens) in the absence of other causative factors. The Rotterdam criteria, which were later proposed in 2004 (Rotterdam ESHRE/ASRM-sponsored PCOS consensus workshop group, 2004a), added the presence of polycystic ovaries to the two previously described factors, and required two out of these three factors to be present for a PCOS diagnosis. The Androgen Excess and PCOS Society (PCOS-AE) definition recently proposed a revision of the Rotterdam criteria; accepting the two out of three-diagnostic paradigm, but requiring hyperandrogenism to be present (Azziz ). Polycystic ovary is defined as an ovary containing 12 or more follicles (or 25 or more follicles using new ultrasound technology) measuring 2 to 9 mm in diameter or alternatively as an ovary that has a volume of >10 ml on ultrasonography. A single ovary meeting either of these definitions is sufficient for diagnosis of polycystic ovaries (Rotterdam ESHRE/ASRM-sponsored PCOS consensus workshop group, 2004b). A high number of antral follicles is among the most important characteristics of women with this complex and potentially multi-factorial syndrome.

As expected, due to the increased antral follicle count (AFC) in women with PCOS (Webber ; Maciel ; Bhide ), these patients tend to respond to controlled ovarian stimulation (COS) with gonadotropins by producing a higher number of oocytes and achieving higher levels of maximum serum estradiol (Sun ). They are also at a higher risk of developing ovarian hyperstimulation syndrome (OHSS). In addition to having a higher AFC, women with PCOS seem to have a prolonged reproductive lifespan (Carroll ), with slightly higher menopausal age than the general population (Ramezani Tehrani ). Therefore, the ovarian follicular dynamics seem to show a significant contrast to that observed in women with poor ovarian response (POR) or diminished ovarian reserve (DOR).

POR/DOR represents the opposite end of a spectrum regarding follicular pool, and reaching a consensus on the definition of POR has also been challenging. A commonly used definition is the ESHRE (European Society for Human Reproduction and Embryology) Bologna criteria (Younis ). It requires two out of this three criteria to be present: a maternal age >40 years (or any other risk factor for low response); a previous poor response cycle with ≤3 oocytes retrieved; and an abnormal ovarian reserve test (AFC < 5–7 or anti-mullerian hormone (AMH) < 0.5–1.1 ng/ml). A more recent and detailed definition has come from the Patient-Oriented Strategies Encompassing Individualize D Oocyte Number (POSEIDON) group, considering age, AFC, serum AMH levels and past response to ovarian stimulation (Esteves ). POSEIDON Group 1 refers to young infertile women (<35 years old), with acceptable ovarian reserve markers (AFC ≥ 5; AMH ≥ 1.2 ng/ml), and unexpected poor (<4 oocytes retrieved) or suboptimal (4–9 oocytes retrieved) response to conventional ovarian stimulation, while POSEIDON Group 3 includes women <35 years old with low ovarian reserve markers (AFC < 5; AMH < 1.2 ng/ml). POSEIDON Groups 2 and 4 include women AMH ≥35 years old, with similar characteristics to Groups 1 and 3, respectively. Regardless of the diagnostic criteria used, it is widely accepted that women with POR/DOR tend to have lower AFC and respond poorly to COS, in contrast to women with PCOS.

Granulosa cells (GCs) surround the oocyte, support follicular development and are essential in the transition from the primordial into a mature follicle. The interactions between GCs and oocytes are critical for coordinated oocyte maturation (Thomas and Vanderhyden, 2006). Therefore, investigating the molecular pathways involved in the proliferation, differentiation and functional transformation of GCs is crucial for the understanding of mechanisms that regulate folliculogenesis. GCs not only play a critical role in normal folliculogenesis, but they are also involved in the abnormal folliculogenesis observed in disorders such as PCOS (Yilmaz ), and premature ovarian insufficiency (Collins ). Small antral follices in PCOS patients’ ovaries, failed to develop into larger dominant follicles. PCOS in the murine is characterised by an increase in the follicle number and GC proliferation (Cozzolino and Seli, 2020). Furthermore, increased GC proliferation in smaller follicles was observed in the ovaries of women with PCOS (Yildiz ; Liu ). In contrast, there is an increase in GC apoptosis in young poor responders (YPR). Therefore, GC dysfunction and/or proliferative dysregulation seems to be among the pathological features of women with PCOS and POR, and GCs constitute a useful cellular model to investigate these complex conditions.

In this study, we hypothesized that characterising the differences between women with PCOS and those with POR/DOR could be useful to identify clinically relevant molecular mechanisms that promote follicular growth and biomarkers of aberrant follicle development. To this end, we characterised the transcriptomic landscape of women with PCOS, comparing and contrasting it to women with POR and normoresponder controls (CONT). To achieve a detailed understanding of the differences between these opposite entities, we analysed both ovarian somatic cells (GCs) and non-reproductive somatic cells (peripheral blood mononuclear cells [PBMNCs]).

Materials and methods

Study patients

This was a prospective study and the samples were prospectively collected. All study procedures were approved by and conducted according to the Institutional Review Board of Hospital Universitario y Politécnico La Fe, Valencia, Spain (520/2018). Prior to collecting samples, all patients were informed, and written consent was obtained

PCOS was defined using the revised Rotterdam diagnostic criteria for PCOS (Rotterdam ESHRE/ASRM-sponsored PCOS consensus workshop group, 2004a). Other causes of oligomenorrhea or hyperandrogenism congenital adrenal hyperplasia, androgen-secreting tumors, Cushing’s syndrome were excluded on clinical grounds for the study. The YPR group included women <35 years old with <4 mature follicles (at least 15 mm) on the day of the trigger. According to internal data, this group represented the bottom 15th percentile of patients’ responses in this age group and was consistent with POSEIDON criteria for POR (Group 3). Women with ovarian cancer, endometriosis or severe systemic disease were excluded. The young normoresponder (CONT) group included women <35 years without PCOS or anovulation, who developed >14 mature follicles (at least 15 mm on transvaginal ultrasound). According to internal data a threshold of >14 mature follicles was established to represent the top 25% of patients in this age group. A total of 20 women with PCOS, 20 YPR and 20 CONT were recruited for the study at the Hospital La Fe (Valencia, Spain).

Controlled ovarian stimulation protocols, oocyte retrieval and granulosa cell collection

The patients included in the study received COS using the GnRH agonist standard long protocol or antagonist protocol. Ovarian stimulation was achieved using recombinant FSH (rFSH) and highly purified HMG (HP-HMG). COS was started after confirming the absence of dominant follicular development by ultrasound. Serum estradiol (E2) <30 mIU/ml. rFSH was administrated daily, with dosages ranging from 150 to 300 IU, alone or in combination with hMG. When the lead follicular diameter was greater than 18–19 mm, 5000–10 000 IU hCG and/or GnRH agonist (GnRHa) were used to achieve final follicular maturation. To reduce the risk of OHSS in PCOS patients, oocyte maturation was primarily achieved using GnRHa alone, whereas those in the YPR group used 10 000 IU hCG. Patients in the CONT group used one of these strategies based on their physician’s preference. Oocyte retrieval was scheduled 36 h after trigger administration.

Follicular fluid samples were carefully collected from the first aspirated follicle of each ovary. Samples containing blood contamination were excluded. Ovarian GCs were isolated and purified from each follicular fluid sample by centrifugation at 450 g for 5 min. Recovered GCs were then washed with 1 ml of phosphate buffer saline (PBS) and centrifuged at maximum speed for 5 min. The GCs were resuspended in 10 µl of PBS. After that, 50 µl of RNA were added and samples were stored at −20°C. Over 95% of purity was achieved for all GC samples before being used in subsequent experiments. The GCs were then processed for RNA extraction.

Collection of peripheral blood mononuclear cell samples

Blood was collected in EDTA BD Vacutainer® tubes (BD Diagnostics, Spain) to isolate PBMNCs by standard Ficoll-based centrifugation protocols. Briefly, 4 ml of blood was diluted to a final volume of 7 ml and added to 3 ml of Ficoll-Paque (GE Healthcare, Uppsala, Sweden) to undergo a 2500 rpm centrifugation for 30 min at RT. Buffy coat containing PBMNCs was then collected and PBMNCs were resuspended in PBS and washed twice. The pellet was resuspended in 100 µl of PBS, with 500 µl of RNA later and stored at −20°C until ready for RNA extraction.

RNA extraction and sequencing

Total RNA was extracted from PBMNCs using the RNeasy Micro kit (Qiagen). RNA integrity and concentration were measured using the Agilent High Sensitivity RNA screen tape system and Qubit2, respectively. All samples used for RNA-seq had an RNA integrity number (RIN) of 7 or higher. mRNA was purified from approximately 200 ng of total RNA with oligo-dT beads and sheared by incubation at 94 °C in the presence of Mg (Kapa mRNA Hyper Prep). Following first-strand synthesis with random primers, second-strand synthesis and A-tailing were performed with dUTP for generating strand-specific sequencing libraries. Adapter ligation with 3′ dTMP overhangs was ligated to library insert fragments. Library amplification amplified fragments carrying the appropriate adapter sequences at both ends. Strands marked with dUTP were not amplified. Indexed libraries that meet appropriate cutoffs for both were quantified by quantitative reverse transcription PCR (qRT-PCR) using a commercially available kit (KAPA Biosystems) and insert size distribution determined with the LabChip GX or Agilent Bioanalyzer. Samples with a yield of ≥0.5 ng/µl were used for sequencing.

RNA was extracted from purified GCs and cDNA was amplified using Smart-Seq v4 ultra-low input RNA kit for sequencing following the manufacturer’s protocol. cDNA was amplified and measured using dsDNA High Sensitivity Kit on Qubit 2.0 (Life Technologies, Carlsbad, CA, USA) for concentration and Tapestation 2200 (Agilent Technologies, Santa Clara, CA, USA) for size distribution. RNA-sequencing libraries were constructed using the Nextera® XT DNA Library Preparation Kit (Illumina, San Diego, CA, USA) and multiplexed using Nextera® XT Index Kit (Illumina). Libraries were quantified by Qubit2.0 and Tapestation 2200 (Agilent technologies). Indexed libraries were then pooled and sequenced on Illumina’s Hiseq 2500 platform with 75 bp pair-end reads. We analysed 15 samples in total, and approximately 50 million paired-end reads were achieved for each sample.

RNA-sequencing analysis

Image analysis, base calling and generation of sequence reads were produced using the HiSeq Control Software v2.0 (NCS) and Real-Time Analysis Software (RTA). Data were converted to FASTQ files using the bcl2fastq2 v1.8.4 software (Illumina Inc.). The reads were trimmed for quality and aligned with the reference human genome hg38 with gencode annotation (Frankish ). The normal annotation has ∼50K entries, but the gencode annotation has over 100K annotated regions on the genome. We used HiSAT2 for alignment, and StringTie and BallGown for transcript abundance estimation (transcript level expression analysis of RNA-seq experiments with HISAT, String Tie and Ball gown). We used DESEq2 for differential gene expression. Gene expression values were calculated as FPKM using Cufflinks 2.1.1. The genes were called as differentially expressed if the adjusted P-value was ≤0.05. We used R for downstream processing and visualization of the data.

Ingenuity pathway analysis

Ingenuity pathway analysis (IPA) Ingenuity Systems QIAGEN (Content version: 45865156, 2018, www.qiagen.com/ingenuity, Redwood City, CA, USA) software was used to perform Gene Ontology (GO) and pathways analyses. A differential expression Log2 fold change (FC) ≥0.584 false discovery rate (FDR) ≤0.05 was assigned as differentially expressed in different comparisons. Each gene symbol was mapped to its corresponding gene object in the Ingenuity Pathways Knowledge Base. Differentially expressed genes (DEGs) used in pathway analysis were determined between treated and control groups by using a filtering criteria FC >0.5 and Benjamini–Hochberg (B-H FDR) ≤0.05 (Benjamini and Hochberg,1995). IPA was used to identify pathways overrepresented in DEGs (Benjamini and Hochberg, 1995).

Quantitative reverse transcription PCR

Another set of PBMNC samples were used for the validation of selected genes that were seen to be differentially regulated in the RNA-sequencing studies. The qRT-PCR screening was performed on genes involved in the relevant clinical pathway. The primer sequences of the genes used for the validation are listed in Supplementary Table SI.

Total RNA was obtained from PBMNCs using RNAqueous Microkit (Thermo Fisher Scientific, Waltham, MA, USA) and was treated with DNase I (Thermo Fisher Scientific) for genomic DNA contamination. Reverse transcription was performed using the RETROscript kit (Thermo Fisher Scientific) in two steps: first, template RNA and random primers were incubated at 85°C for 3 min to eliminate any secondary structures, and then the buffer and enzyme were added, and the reaction was carried out at 42°C for 1 h. Gene expression analysis was done by real-time PCR using Syber Green on the ViiA7 real-time PCR machine from Applied Biosystems. To check reproducibility, each qPCR reaction was carried out in triplicate, and water was used as a negative control. The amplification was a 45-cycle run with annealing temperature at 60°C. GAPDH was used as the housekeeping gene. Moreover, at the end of the amplification, a melting curve for each sample was generated to assess the absence of primer dimers or contamination.

Real-time PCR (RT-PCR) was also performed to confirm differential expression of genes identified as altered in the cumulus cells. RNA was amplified, and cDNA was prepared using established protocol of SMART-Seq v4 Ultra Low Input RNA Kit for Sequencing. Gene expression analysis was done by RT-PCR using Syber Green on the ViiA7 real-time PCR machine by Applied Biosystems. The PCR reaction was prepared using 5 µl PowerUp™ SYBR™ Green Master Mix (2×), 1 µl forward and reverse primers, 1 µl of DNA template and 3 µl of nuclease-free water. The components were mixed thoroughly and briefly centrifuged. PCR cycling conditions were 95°C for 5 min followed by 45-cycle run with annealing temperature at 60°C. GAPDH was used for normalization. ΔΔCt method was used for the calculation of the difference in the expression of genes A melt curve analysis was subsequently performed to confirm the specificity of products.

Statistical analysis

All the characteristics of PCOS, YPR and healthy controls were given as mean ± SD. The expression data for GCs and PBMNCs obtained by RT-qPCR was analysed with the two-tailed t-tests using the GraphPad Prism version 5.00 software (GraphPad, San Diego, CA, USA; www.graphpad.com). Differences between groups were considered significant when the P-value was <0.05.

Results

RNAseq analyses were performed to compare mRNA expression levels in ovarian GCs and PBMNCs between the three experimental groups of women: (i) PCOS (n = 5) who were young (age 32.5 ± 0.7 years) with a robust response to ovarian stimulation during ART; (ii) YPR (n = 5; age 33.1 ± 1.4 years) with poor response to ovarian stimulation; (iii) CONT (n = 5; age 31.0 ± 2.1 years) who were women undergoing ART treatment. They did not have PCOS or anovulation, and showed a robust response to COS. Baseline characteristics and the number of oocytes retrieved are reported in Table I.

Table I

General characteristics of patients in the three groups.

	CONT	PCOS	YPR
Age (years)	31.0±2.1	32.5 ±0.7	33.1 ±1.4
Number of oocytes	22 ± 3.9	28 ± 6.4	5 ± 0.7
Number of metaphases II	21 ± 2.2	16 ± 3.8	4 ± 0.6
Infertility reason
Tubal	6	–	0
Endometriosis	1	–	4
Male infertility	9	–	3
Idiopathic	4	–	–
Diminished ovarian reserve	–	–	13
PCOS	–	20	–
Duration of infertility	3.3 ± 1.6	3.0 ± 1.3	3.1 ± 1.2
Tests for insulin resistance	1.7 ± 0.5	3.0 ± 2.1	0.9 ± 0.5
Antral follicle count	17.7 ± 9.8	26.3 ± 13.8	3.1 ± 1.8
AMH (pmol/l)	17.7 ± 10.3	37.7 ± 18.3	5.4 ± 2.7
Ovarian stimulation protocol	GnRH antagonist protocol	GnRH antagonist protocol	GnRH antagonist protocol
Duration of stimulation	10.0 ± 1.5	10.6 ± 1.8	10.9 ± 1.7
Total dose of gonadotropin	2183 ± 926	1798 ± 912	3911 ± 1526
Type of triggering
hCG	10	8	18
Agonist GnRH	9	9	1
Dual triggering	1	3	1
Estradiol trigger day	2238 ± 848	3335 ± 1863	1873 ± 696
Progesterone trigger day	1,26 ± 0.86	1.13 ± 0.57	0.83 ± 0.44
Rate/number OHSS*	2/20 (10%)	8/20 (40%)	0/20 (0%)
OTC therapy rate	6/20	14/20	0/20
PCOS sub-phenotype	0 (0%)		0 (0%)
Phenotype A with hyperandrogenism +oligomenorrhea +PCO		2(10%)
Phenotype B with hyperandrogenism +oligomenorrhea		4 (20%)
Phenotype C with hyperandrogenism+PCO		10 (50%)
Phenotype D with oligomenorrhea +PCO		4 (20%)

CONT, control group of normal responders; OHSS, ovarian hyperstimulation syndrome; OTC, over the counter; PCO, polycystic ovary; YPR, young poor responder.

Only cases of mild OHSS are reported in the study, no cases of moderate or severe OHSS.

To identify the genes and pathways differentially expressed in PCOS, YPR and CONT, an unbiased approach was adopted with comprehensive genome-wide transcriptomic analysis. The most significant genes for each comparison are reported in the Supplementary Table SII.

Differential gene expression in the PBMNCs of women with PCOS

Differentially expressed genes in PBMNCs of women with PCOS compared to CONT

Hierarchal clustering of the DEGs partitioned into two distinct clusters showed differential gene expression in PBMNCs of women with PCOS when compared with the CONT (Fig. 1a), suggesting high reproducibility of the sequencing data. A total of 65 genes showed changes in the expression between these two groups. 10 genes showed decreased gene expression in the PCOS samples when compared with controls and 55 genes showed increased expression (Fig. 1b). For genes that were over-expressed the FC ranged from 1.8 to 22.9. For under-expressed genes, the FC ranged from −1.8 to −6.49 (Fig. 1c). Zinc finger protein 57 homolog (ZFP57) and Immunoglobulin Lambda Constant (IGLC1) were significantly downregulated in PBMNCs of women with PCOS compared to the controls, FC = 22.9 and 9.31, respectively. Telomerase RNA component (TERC) and Small Nucleolar RNA, C/D Box 14D (SNORD14D) were overexpressed in PCOS compared to controls, FC = −4.08 and −3.23, respectively (Fig. 1c). The functional enrichment analysis based on data from IPA software were used to perform GO, and revealed that DEGs in PBMNCs of women with PCOS were mainly involved in the pathway of granulocyte adhesion and breast cancer regulation. All affected pathways are reported in Fig. 1d.

Figure 1.

Gene expression is altered in PCOS women's peripheral blood mononuclear cells (PBMNCs) compared to normoresponder controls (CONT). (a) The heat map illustration shows differentially expressed genes in CONT and PCOS in PBMNCs. The colour spectrum ranging from red to blue indicates normalized levels of gene expression from high to low. (b) Volcano plots for RNA-seq comparing CONT and PCOS. Red spots on the left upper box represent –log10 (P-value) ≥2; red spots on the right upper box represent the –log10 (P-value) <2. (c) Selected differentially expressed genes in CONT and PCOS, P < 0.05 for each. The transcripts per million (TPM) value represents the relative expression level comparable between samples. For the box plots, the bottom and top whiskers denote 5 and 95 percentile values, the bottom and top bounds of the rectangle denote the 25 and 75 percentile values, and the line in between denotes the median (50 percentile) value of the distribution. (d) Pathway analysis was evaluated using the Gene Ontology bioinformatics tool in PBMNCs. Log2 fold change (FC) ≥0.584 false discovery rate (FDR) ≤0.05.

Differentially expressed genes in PBMNCs of women with PCOS compared to YPR

Gene expression was found to be altered in PBMNCs of women with PCOS compared to YPR. The heat map showed the hierarchical cluster of the DEGs partitioned into two distinct clusters PCOS and YPR (Fig. 2a). Sixteen genes showed significant changes in the expression. Seven of them were upregulated in the PCOS when compared with the YPR samples, while nine genes showed decreased expression (Fig. 2b). FC ranged from 3.51 to 22.91 for upregulated and −0.6 to −22.94 for downregulated genes (Fig. 2c). The pathway analysis did not assign the genes to any specific pathways.

Figure 2.

Comparison in gene expression between PCOS and young poor responder (YPR) in peripheral blood mononuclear cells (PBMNCs). (a) The heat map illustration shows differentially expressed genes in PCOS and YPR in PBMNCs. The colour spectrum ranging from red to blue indicates normalized levels of gene expression from high to low. (b) Volcano plots for RNA-seq comparing PCOS and YPR. Red spots on the left upper box represent –log10 (P-value) ≥2; red spots on the right upper box represent the –log10 (P-value) <2. (c) Selected differentially expressed genes in PCOS and YPR, P < 0.05 for each. The transcripts per million (TPM) value represents the relative expression level comparable between samples. For the box plots, the bottom and top whiskers denote 5 and 95 percentile values, the bottom and top bounds of the rectangle denote the 25 and 75 percentile values, and the line in between denotes the median (50 percentile) value of the distribution.

Differential gene expression in the granulosa cells of women with PCOS

Differentially expressed genes in granulosa cells of women with PCOS compared to CONT

Hierarchal clustering of the DEGs partitioned into two distinct clusters in the heat map of significant genes. It showed differential GC gene expression between women with PCOS and control (Fig. 3a). Four genes showed decreased expression in the PCOS group (Fig. 3b and c). There were no genes found to be overexpressed in PCOS compared to controls. For the overexpressed genes, FC ranged from 1.77 to 5.1 as reported in Fig. 3c. The functional enrichment analysis performed GO revealed that DEGs in PCOS compared to controls were prevalently included in the protein ubiquitination pathway (Fig. 3d).

Figure 3.

Gene expression is altered in PCOS women's granulosa cumulus cells (GCs) compared to normoresponder controls (CONT). (a) The heat map illustration shows differentially expressed genes in CONT and PCOS in GCs. The colour spectrum ranging from red to blue indicates normalized levels of gene expression from high to low. (b) Volcano plots for RNA-seq comparing CONT and PCOS. Red spots on the left upper box represent –log10 (P-value) ≥2; red spots on the right upper box represent the –log10 (P-value) <2. (c) Selected differentially expressed genes in CONT and PCOS, P < 0.05 for each. The transcripts per million (TPM) value represents the relative expression level comparable between samples. For the box plots, the bottom and top whiskers denote 5 and 95 percentile values, the bottom and top bounds of the rectangle denote the 25 and 75 percentile values, and the line in between denotes the median (50 percentile) value of the distribution. (d) Pathway analysis was evaluated using the Gene Ontology bioinformatics tool in GCs. Log2 fold change (FC) ≥0.584 false discovery rate (FDR) ≤0.05.

Differential expression was seen in granulosa cells of PCOS patients compared to the YPR

Differential gene expression analysis in GCs showed a significant change in the expression of mRNAs in PCOS compared to YPR, as demonstrated by the heat map illustration (Fig. 4a). In total, 58 genes were differentially expressed between PCOS and YPR (Fig. 4b). Figure 4c shows the box plot representation of the top over-expressed genes with FC ranging from 2.07 to 2.33, and the under-expressed genes with FC ranging from −2.08 to −2.76. The IPA showed that the principal cluster of genes were involved in the pathways of oxidative phosphorylation, mitochondrial dysfunction and sirtuin signaling. The clusters of genes involved in oxidative phosphorylation were overexpressed in PCOS compared to YPR. In contrast, the genes involved in the pathway of sirtuin signaling were downregulated in PCOS compared to YPR (Fig. 4d).

Figure 4.

Gene expression is altered in granulosa cumulus cells (GCs) of women with PCOS and young poor responder (YPR). (a) The heat map illustration showing differentially expressed genes in PCOS and YPR in GCs. The colour spectrum ranging from red to blue indicates normalized levels of gene expression from high to low. (b) Volcano plots for RNA-seq comparing PCOS and YPR. Red spots on the left upper box represent –log10 (P-value) ≥2; red spots on the right upper box represent the –log10 (P-value) <2. (c) Selected differentially expressed genes in PCOS and YPR, P < 0.01 for each. The transcripts per million (TPM) value represents the relative expression level comparable between samples. For the box plots, the bottom and top whiskers denote 5 and 95 percentile values, the bottom and top bounds of the rectangle denote the 25 and 75 percentile values, and the line in between denotes the median (50 percentile) value of the distribution. (d) Pathway analysis was evaluated using the Gene Ontology bioinformatics tool in GCs. Log2 fold change (FC) ≥0.584 false discovery rate (FDR) ≤0.05.

Validation of differentially expressed selected genes by RT-qPCR

For validation of the RNA-sequencing data by qRT-PCR, two genes from each set were chosen for further analysis. qRT-PCRs were performed in independent samples (n = 15 for each group), and results showed a similar expression pattern to that observed in the RNA-sequencing.

In PBMNCs of women with PCOS compared to controls, there was a trend toward a decreased expression of RAB4 B (P = 0.30) and increased expression of TERC (P = 0.1) (Supplementary Fig. S1a). Similarly, a trend toward a decreased expression of peroxidasin (PXDN) (P = 0.3) and ZFP57 (P = 0.1) in the PBMNCs of women with PCOS compared to YPR was observed (Supplementary Fig. S1b).

When samples from women with PCOS were compared to controls, there was a decrease in the expression of membrane frizzled-related protein (MFRP) (P = 0.03) and a trend toward decreased expression of myotubularin-related protein 9 (MTMR9) (P = 0.3) in the GCs (Supplementary Fig. S2a). When GCs samples from women with PCOS were compared to those from YPR, there was an increase in the expression mitochondrially encoded 12S RRNA (MT-RNR1) (P = 0.03) and a trend toward increased expression of mitochondrially encoded ATP synthase membrane subunit 6 (MT-ATP6) (P = 0.1) (Supplementary Fig. S2b).

Discussion

This study describes the transcriptomic analysis of PBMNCs and GCs in women with PCOS compared to controls and young women with DOR, showing differences in gene expression landscapes and pathway alterations following ovarian stimulation.

In PBMNCs of women with PCOS, RNA-seq analysis found several genes to be altered compared to controls. Some of these genes were involved in the mechanisms of granulocyte and agranulocyte adhesion and diapedesis. In addition, TERC expression was increased in PBMNCs of patients with PCOS. There have been several studies suggesting that somatic cell telomere length may be affected in women with PCOS. Some suggested shortened telomeres in leukocytes in women with PCOS (Li ). Others did not find a difference in the leukocyte telomere length but described a negative correlation between inflammatory biomarkers and telomere length (Pedroso ) or detected altered telomere length in their GCs (Wei ).

The granulosa cumulus cells play an important role in follicular development, oocyte maturation, and ovulation. The crucial impact of follicular somatic cells on peri-ovulatory events is in large part mediated by the bidirectional communication between cumulus cells (CCs) and the oocyte through specialized gap junctions. Defective communication between the CCs and the oocyte may impair oocyte quality and result in poor embryo development and decreased fertility. Similarly, altered gene expression and thus function of GCs and CCs that affect extracellular matrix formation, cellular proliferation, differentiation or apoptosis could be associated with lower oocyte quality.

We found an increased expression of genes involved in oxidative phosphorylation and mitochondrial dysfunction in GCs of women with PCOS compared to YPR. This suggests that oxidative stress plays a more dominant role in the dysfunction of follicles from women with PCOS compared to similar age women with poor ovarian reserve. In contrast, the genes involved in the pathway of sirtuin signaling were downregulated in PCOS compared to YPR.

Sirtuins are complex proteins implicated in longevity. They take part in pathways that regulate oxidative stress, maintenance of metabolic homeostasis, DNA repair and mitochondrial function. The downregulation of sirtuin pathway in GCs of women with PCOS compared to YPR is quite interesting. It suggests that the metabolic dysfunction inherent to PCOS may negatively impact the sirtuin pathway and may downregulate longevity-promoting gene expression.

Sirtuin signaling responds to metabolic challenges, inflammatory signals or hypoxic/oxidative stress (Tatone ) and has emerged as a critical player in the regulation of key processes in oogenesis. In GCs, sirtuins may play a role in the regulation of proliferation and hormonal metabolism, as they are responsible for the activation of steroidogenesis associated with luteinization (Tatone ). Data suggest that sirtuin 1 (SIRT1) is involved in the maintenance of cellular redox balance with an important role in sensing and modulating oxidative stress in GCs (Tatone ). Furthermore, the mitochondrial sirtuins have been implicated as sensors of metabolic state in human GCs. Specific alterations to the mitochondrial sirtuins might affect mitochondrial proteins, resulting in metabolic alterations in the ovarian follicle (Pacella-Ince ), potentially contributing to follicular dysfunction in women with PCOS.

Several genes involved in the mechanism of oxidative phosphorylation were upregulated in the GCs of women with PCOS, showing a link to a possible increase of oxidative stress in their ovaries. Increased oxidative stress has been widely accepted as a pivotal pathological feature of PCOS. However, the regulatory mechanisms of oxidative stress in the ovaries have not been fully clarified. Oxidative phosphorylation followed by electron transport chain reaction represented the most important metabolic pathway related to energy production in cells. Nevertheless, it is associated with a high cost due to the production of reactive oxygen species (ROS) (Halliwell, 2007). Although ROS are pivotal in the biological processes and act as signaling molecules, ROS may interact with several biomolecules such as lipids and nucleic acids, causing damage in the CCs and GCs.

The normal redox homeostasis in the endoplasmic reticulum (ER), may be altered by the excessive ROS production. The imbalance in redox homeostasis induces ER stress, which causes the unfolded protein response by activating genes encoding factors involved in protein folding and antioxidative machinery to restore ER homeostasis (Zhang ). Altered ROS production in PCOS could be a causative factor reducing the functionality of GCs (Tarin, 1995) with consequent compromise follicle quality and aberrant reproductive endocrinology in the ovary (Qiao ). Our results showed that several genes associated with complex IV deficiency of the mitochondrial respiratory chain such as MT-RNR1, MT-RNR2 (Mitochondrially Encoded 12S RRNA) and COX2 (Mitochondrially Encoded Cytochrome C Oxidase II) were dysregulated in GCs from PCOS patients. These specific genes are considered essential components of the respiratory chain that catalyses the reduction of oxygen to water. In fact, decreased mitochondrial biogenesis and mtDNA have been described in GCs from PCOS women (Zhao ). Altogether, these results support the hypothesis that mitochondrial dysfunction of human GCs in women with PCOS may contribute to impaired steroidogenesis, fertilization, oocyte maturation and oocyte quality (Sreerangaraja Urs ).

The results of qRT-PCR analyses showed a similar expression pattern to that observed in the RNA-sequencing. In PBMNCs of women with PCOS compared to controls, there was a trend toward decreased expression of RAB4 B (P = 0.30) and increased expression of TERC (P = 0.1) (Supplementary Fig. S1a). Similarly, a trend toward a decreased expression of PXDN (P = 0.3) and ZFP57 (P = 0.1) in the PBMNCs of women with PCOS compared to YPR was observed (Supplementary Fig. S1b). In the GCs, a decrease in the expression of MFRP (P = 0.03) and a trend toward a decreased expression of MTMR9 (P = 0.3) was observed when samples of women with PCOS were compared to controls (Supplementary Fig. S2a). Comparing GC samples from women with PCOS to those from YPR, we also observed an increase in the expression MT-RNR1 (P = 0.03) and a trend toward increased expression of MT-ATP6 (P = 0.1) (Supplementary Fig. S2b).

In the current study, oocyte maturation was achieved using hCG, GnRHa or a combination of both. A study by Miller previously demonstrated that patients who receive hCG alone or as part of a double trigger (together with GnRHa) have higher levels of vascular endothelial growth factor and lower levels of pigment epithelium-derived growth factor. Others found increased epidermal growth factor-like proteins (amphiregulin and epiregulin) and decreased gap junction protein Connexin 43 in women triggered with a combination of hCG and GnRHa compared to hCG alone (Haas ). Therefore, in the current study, we cannot exclude the possibility that different trigger methods might contribute to differential gene expression between the groups. It is however noteworthy that none of the previous studies identified the genes implicated in our study to be affected by the medications used to induce oocyte maturation and ovulatory changes.

In many respects, follicular numbers and patterns of response to COS of women with PCOS and POR seem to be at the opposite ends of a spectrum. Therefore, it is worth investigating whether these two important fertility-related disorders develop as a consequence of a perturbation in a single pathway (i.e. up- vs downregulation of similar gene groups). Our study is a step in this direction and our preliminary findings suggest specific differences between women with PCOS and POR. Future studies with adequate sample size directed at well-defined subgroups of these extremely heterogeneous disorders are needed to delineate shared gene pathways that can be exploited to identify therapeutic targets.

Conclusion

RNA expression profiles in ovarian GCs and PBMNCs were significantly altered in patients with PCOS compared with CONT and YPR. GCs of PCOS patients showed altered expression of several genes involved in oxidative phosphorylation, mitochondrial function and sirtuin signaling pathways. These findings are consistent with prior proposed mechanisms involved in the known features of PCOS.

Data availability

The data underlying this article are available in the article and in its online supplementary material.

Authors’ roles

M.C. and E.S. take primary responsibility for the paper. S.T., S.H., I.P. and M.R. collected the data. M.C. and S.T. did the statistical analyses. M.C. and E.S. coordinated the research. S.H., L.R., R.S. and A.P. participated in the data analyses. M.C., S.T. and E.S. drafted the manuscript.

Funding

This study was partially supported by grant PI18/00322 from Instituto de Salud Carlos III and European Regional Development Fund (FEDER), ‘A way to make Europe’ awarded to S.H.

Conflict of interest

M.C., S.H., S.T., L.R., M.R., I.R., A.P. and R.S. declare no conflict of interests concerning this research. E.S. is a consultant for and receives research funding from the Foundation for Embryonic Competence.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.