Global Transcriptional Profiling of Granulosa Cells from Polycystic Ovary Syndrome Patients: Comparative Analyses of Patients with or without History of Ovarian Hyperstimulation Syndrome Reveals Distinct Biomarkers and Pathways

Ovarian hyperstimulation syndrome (OHSS) is often a complication of polycystic ovarian syndrome (PCOS), the most frequent disorder of the endocrine system, which affects women in their reproductive years. The etiology of OHSS is multifactorial, though the factors involved are not apparent. In an attempt to unveil the molecular basis of OHSS, we conducted transcriptome analysis of total RNA extracted from granulosa cells from PCOS patients with a history of OHSS (n = 6) and compared them to those with no history of OHSS (n = 18). We identified 59 significantly dysregulated genes (48 down-regulated, 11 up-regulated) in the PCOS with OHSS group compared to the PCOS without OHSS group (p-value < 0.01, fold change >1.5). Functional, pathway and network analyses revealed genes involved in cellular development, inflammatory and immune response, cellular growth and proliferation (including DCN, VIM, LIFR, GRN, IL33, INSR, KLF2, FOXO1, VEGF, RDX, PLCL1, PAPPA, and ZFP36), and significant alterations in the PPAR, IL6, IL10, JAK/STAT and NF-κB signaling pathways. Array findings were validated using quantitative RT-PCR. To the best of our knowledge, this is the largest cohort of Saudi PCOS cases (with or without OHSS) to date that was analyzed using a transcriptomic approach. Our data demonstrate alterations in various gene networks and pathways that may be involved in the pathophysiology of OHSS. Further studies are warranted to confirm the findings.


Introduction
Ovarian hyperstimulation syndrome (OHSS) is a rare, iatrogenic complication of ovarian hyperstimulation by assisted reproduction treatment and is potentially life-threatening [1][2][3]. Women with polycystic ovary syndrome (PCOS) are at a higher risk for developing OHSS 2 of 12 due to having a large number of follicles on their ovaries and the tendency to over-respond to the hormones used for inducing fertility [2]. The causes of OHSS are generally unknown. Since complications can be severe and may even be life-threatening, it is critical to predict whether a woman going through gonadotropin therapy may develop OHSS.
Previous studies have reported several risk factors for the development of OHSS, such as young age, PCOS, low body mass index (BMI), and previous history of OHSS [4,5]. Granulosa cells (GC) are essential for ovarian folliculogenesis and are known to play a critical role in follicular development and oocyte maturation [6,7]. Several studies have indicated that GC dysfunction in women with PCOS may contribute to abnormal folliculogenesis [6,7]. Recent advances in omics technologies (genomics, transcriptomics, proteomics, and others) have enabled researchers to better understand the molecular characteristics of diseases and to identify disease biomarkers [7][8][9][10][11].
Genome-wide studies in various cells, including granulosa cells from women with PCOS and non-PCOS controls, have revealed several potential genes and pathways, such as the ERK/MAPK and VEGF signaling pathways [6,11]. However, OHSS (and PCOS) has a complex etiology and genetic basis and the disease pathogenesis and mechanism still largely remain unclear [5,6,12,13]. In this study, we aimed to identify the potentially essential genes and pathways in PCOS patients who would be prone to developing OHSS and to identify potential markers that may be linked to the OHSS and its pathobiology using a transcriptomic approach. To this purpose, we investigated if the PCOS patients who had developed OHSS have a distinct molecular profile in their granulosa cells as compared to those who did not in order to understand the biology of its development.

Patients
Twenty-four infertile women with PCOS undergoing IVF treatment were identified and examined by an IVF physician at Gynecology/IVF clinics at King Faisal Specialist Hospital and Research Center (KFSHRC). PCOS was defined according to the criteria established by the American Society of Reproductive Medicine and the European Society of Human Reproduction and Embryology (ASRM/ESHRE) [14]. The PCOS patients who are younger than 40 years, with FSH ≤ 12, body mass index (BMI) ≤ 35 kg/m 2 and not undergoing any immunosuppressive therapy, were asked to participate in the study. The study was approved by the institutional review board (IRB) of the King Faisal Specialist Hospital and Research Center (no.08-MED604-2; RAC#2100002; RAC#2110006) and it is performed in accordance with the current version of the Helsinki Declaration. Written informed consent was obtained from all patients before participation in the study.
The study included samples from 24 patients (n = 24) with PCOS. Pituitary downregulation was performed by administering gonadotrophin-releasing hormone (GnRH) agonist long or short protocols as described earlier [10]. Six patients had a positive history of OHSS (four patients had a history of OHSS and two had it in the current cycle) (referred to as "OHSS" group; n = 6) and the control group was without any history of OHSS (referred to as PCOS; n = 18). All OHSS cases were severe or moderate. Granulosa cells (3 mL) were collected into 15 mL regular falcon conical tubes from participating subjects undergoing control ovarian stimulation as described previously [10] for microarray and quantitative RT-PCR (qRT-PCR) experiments.

RNA Isolation
The total RNA was isolated from granulosa cells using Trizol or Total RNA Isolation Kit™ (ThermoFisher Scientific, Waltham, MA, USA) according to manufacturer's instructions. The quality and quantity of RNA were determined by measuring absorbance spectra on a UV/Vis spectrophotometer, NanoDrop ® ND-1000 (ThermoFisher Scientific, Waltham, MA, USA). Further quality check was carried out using the RNA 6000 Nano Assay and 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The high-quality RNA was either immediately used in the experiments (for the real-time RT-PCR and microarray experiments) or stored at −80 • C for further use.

Genome-Wide Gene Expression Profiling
Whole-genome gene expression profiling of samples from the OHSS (n = 6) and unrelated PCOS (n = 18) was performed using Affymetrix's Human Genome U133 Plus 2.0 Array (ThermoFisher Scientific, Waltham, MA, USA). High-quality RNA was converted into labeled cRNA, fragmented, hybridized onto the chip surface according to the manufacturer's protocols as described previously [15]. Briefly, the total RNA was converted into double-stranded cDNA, and then the cDNA was used for cRNA synthesis during which biotinylated UTPs and CTPs were incorporated into the cRNA. Target-labeled cRNA was fractionated and then hybridized onto the chip's surface. The experimental procedures and quality control procedures at each critical step (before hybridization as well as post-hybridization) were strictly followed according to the manufacturer's guidelines. Washing, staining, and scanning were performed according to the manufacturer's instructions and guidelines.
Microarray data normalization was performed using GC Robust Multi-array Average (GC-RMA) algorithm [16,17]. We performed independent two-sample t-test to identify genes whose expression significantly varied between OHSS and PCOS groups. Differentially expressed genes (DEGs) were defined as those with p-value < 0.01 and fold change (FC) ≥ 1.5. Two-dimensional hierarchical clustering is performed using Pearson's correlation with average linkage clustering. Functional annotation and biological term enrichment analysis were performed using DAVID Bioinformatics Resources [18]. Statistical analyses were performed using SPSS version 20 (SPSS, Inc., Chicago, IL, USA) and PARTEK Genomics Suite (Partek Inc., Chesterfield, MO, USA). All statistical tests were two-sided and p-value < 0.05 was considered statistically significant.

Functional, Pathway and Network Analyses
The functional, pathway, gene ontology enrichment, and network analyses were performed using Ingenuity Pathways Analysis (IPA) (QIAGEN Inc., Venlo, Netherlands; https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis (accessed on 18 September 2020). The DEGs were mapped to their corresponding gene object in the Ingenuity Pathways Knowledge Base and significant gene interaction networks are identified. Scores of ≥2 were considered significant after applying a 99% confidence level. A right-tailed Fisher's exact test was used to calculate a p-value determining the probability that the biological function (or pathway) assigned to that data set is explained by chance alone.

Quantitative RT-PCR (qRT-PCR)
To validate our microarray results, confirmatory qRT-PCR was performed using the ABI 7500 Sequence Detection System (ThermoFisher Scientific, Waltham, MA, USA). For this purpose, 50 ng total RNA procured from the same microarray study samples were transcribed into cDNA using Sensiscript Kit (QIAGEN Inc., Venlo, The Netherlands) under the following conditions: 25 • C for 10 min, 42 • C for 2 h, and 70 • C for 15 min in a total volume of 20 µL. Then, 2-5 µL of cDNA was amplified under the following conditions: Initial denaturation of 5 min at 95 • C followed by 34 cycles of "denaturation at 95 • C for 1 min, annealing at 60 • C for 1 min and extension at 72 • C for 1 min" and a final extension of 10 min at 72 • C.

Clinical Characteristics of the Patients
The study included 24 women with PCOS with or without a history of OHSS. Granulosa cells were isolated from six PCOS patients with a history of OHSS (referred to as OHSS) and 18 without OHSS (referred to as the PCOS group). The clinical features were similar between the two groups of patients in terms of age, weight, BMI, FSH, LH, total testosterone, and androstenedione (Mann-Whitney U-test p-value > 0.05) ( Table 1).

Identification of Differentially Expressed Genes
We analyzed the genome-wide mRNA expression profiling of 24 samples from OHSS and PCOS groups using Affymetrix's GeneChip ® Human Genome U133 Plus 2.0 Arrays, which includes over 47,000 transcripts and variants using more than 54,000 probe sets. The array technology is a well-established and reliable method to assess global gene expression profiling [15,21]. Comparison of transcriptomes of OHSS with those PCOS patients with no history of OHSS revealed significant dysregulation of 1520 probes, corresponding to 1188 genes (p-value < 0.01), of which 65 probes (corresponding to 59 genes) had a greater than 1.5-fold change (FC) between the two groups ( Table 2). The unsupervised principal component analysis (PCA) and hierarchical clustering in both dimensions (samples and genes) were performed using Pearson's correlation with average linkage clustering, where both analyses revealed clear discrimination of samples as PCOS and OHSS as well as the pattern of genes deregulation defining two main transcriptome clusters ( Figure 1A,B, respectively) that clearly demonstrated the significant differences in gene expression profiles of having a positive or negative history of OHSS among the PCOS patients.

Functional, Pathway and Gene Network Analysis of Dysregulated Genes
The gene ontology (GO) and functional analyses of differentially expressed genes in the OHSS group compared to PCOS were performed using IPA and DAVID Bioinformatics Resources [18,22]. The biological functions assigned to the data set were ranked according to the significance of the biological function to the dataset. As presented in Figure 1C, differentially expressed genes were enriched with functional categories including cellular development, connective tissue development and function, inflammatory and immune response, cellular growth and proliferation (including DCN, FOXP1, GRN, IL33, INSR, KLF2, FOXO1, FOS, VIM, PAPPA, and ZFP36). Significantly, altered canonical pathways in the OHSS included the PPAR, IL6, IL10, NF-κB and 14-3-3-mediated signaling pathways ( Figure 1D).
To gain an in-depth insight into the interactions of the dysregulated genes involved in the different pathways, genes that were significantly dysregulated in the OHSS were mapped to the gene networks using the Ingenuity Pathway Analysis [22,23]. The network analysis revealed potentially important hub genes, including DCN, IL33, VIM, VEGF, GPC4, KLF2, ELOVL5, KDM5A, Immunoglobulin, ZFP36, and the NF-κB, FOXO, and JAK/STAT signaling pathways that may be relevant to the pathophysiology of OHSS (Figure 2). pathways ( Figure 1D).
To gain an in-depth insight into the interactions of the dysregulated genes involved in the different pathways, genes that were significantly dysregulated in the OHSS were mapped to the gene networks using the Ingenuity Pathway Analysis [22,23]. The network analysis revealed potentially important hub genes, including DCN, IL33, VIM, VEGF, GPC4, KLF2, ELOVL5, KDM5A, Immunoglobulin, ZFP36, and the NF-κB, FOXO, and JAK/STAT signaling pathways that may be relevant to the pathophysiology of OHSS (Figure 2).

Confirmation of Gene Expression Using qRT-PCR
To validate the microarray results, confirmatory qRT-PCR was performed for randomly selected 21 significantly differentially expressed genes (p-value < 0.01), namely

Discussion
In this study, we performed a transcriptomic comparison of PCOS patients with a history of OHSS to those patients with no history of OHSS using granulosa cells in order to identify potential markers for OHSS and to understand its biology. To the best of our

Discussion
In this study, we performed a transcriptomic comparison of PCOS patients with a history of OHSS to those patients with no history of OHSS using granulosa cells in order to identify potential markers for OHSS and to understand its biology. To the best of our knowledge, this is the largest cohort of Saudi PCOS cases (with or without OHSS) to date that was analyzed using a transcriptomic approach.
Our results revealed the potentially important roles of genes related to cellular development, connective tissue development and function, inflammatory and immune response, cellular growth and proliferation, and reproductive system development and function, including genes such as FOXP1, FOXO1, DCN, IL33, INSR, KLF2, PAPPA, VIM, and ZFP36, that have significant gene expression changes in the OHSS group compared to PCOS patients [24][25][26]. Forkhead box (FOX) transcription factor family members have critical aspects in modulating the genes that are important for various cellular processes such as cell growth, differentiation, and longevity and some of which have a crucial role in embryonic development [27], particularly in female reproduction [28]. Interestingly, Forkhead box P1 (FOXP1) is known to cause estrogen-dependent endometrial cancers through the KRAS pathway. Moreover, altered FOXP1 expression and Wnt-related β-catenin acetylation were observed in endometriotic stromal cells from endometriosis patients [29]. Knockdown experiments pointed out the dysregulation of genes involved in Wnt signaling and recapitulated the endometriotic cellular activities such as reducing collagen gel contraction and inhibiting cell proliferation [29]. Several studies have shown that the deletion of the FOXO1 leads to embryonic cell death as a result of incomplete blood vessel development [30,31] and may have a critical role in placental morphogenesis in the developing embryo [32,33]. FOXO1 was implicated in putatively regulating genes involved in lipid and sterol biosynthesis, suggesting that it may play a role in follicular steroidogenesis (Liu et al., 2009). Furthermore, the expression of FOXO1 is shown to be associated with estrogen receptors alpha (ER-α) and beta (ER-β), which are produced primarily by the ovaries and the placenta during pregnancy. The follicle-stimulating hormone (FSH) stimulates the ovarian production of estrogens by the granulosa cells of the ovarian follicles and corpora lutea [34,35]. The pregnancy-associated plasma protein A (PAPPA) is significantly up-regulated in OHSS patients. PAPPA may have a role in female fertility by modulating ovarian function, preeclamptic placentae and steroidogenesis [24,36].
Zing finger protein (ZFP36), Kruppel-like factor 2 (KLF2), insulin receptor (INSR), elongation of very-long-chain fatty acid (ELOVL5) and Glypican-4 (GPC4) were found to be significantly altered in OHSS. GPC4, which is an adipokine that interacts with the INSR and influences insulin sensitivity [37], was down-regulated in OHSS, which may have a likely role in affecting the control of cell division, growth regulation, body fat distribution, insulin resistance, and arterial stiffness in OHSS [38]. Zing finger protein, ZFP36, was shown to be crucial for female fertility and early embryonic development [25] and influences ovulation and oocyte maturation [39,40]. Therefore, it was proposed as a promising candidate gene for obesity-associated metabolic complications [41]. Altered insulin functions have long been associated with abnormalities in female reproduction [42,43]. Conditions associated with insulin resistance, such as obesity and diabetes mellitus, are often accompanied by increased adiposity or hyperglycemia [44]. Obesity and diabetes are independently associated with an altered female reproductive function [13,45,46]. Interestingly, a recent study reported that ELOVL5 is involved in embryonic development and lipid metabolism [47]. KLF2 is expressed in endothelial cells and inhibited by the inflammatory cytokine interleukin-1, hence, implicated as a novel regulator response to proinflammatory stimuli [48]. KLF2 s inhibition for endothelial cell migration and angiogenesis is partly attributed to its ability to inhibit the VEGF receptor (VEGFR) 2/kinase insert domain-containing receptor (KDR) expression [49]. Moreover, G-protein signaling modulator 1 (GPSM1) and Myosin X (MYO10) also displayed dysregulation in OHSS patients. GPSM1 plays a critical role in regulating mitotic spindle orientation, cell polarity, and adenylyl cyclase activity [50] and MYO10 has key functions in filopodia. Experiments with fibroblast-like cells have revealed that MYO10 localizes to the tips of filopodia and undergoes intrafilopodial motility [51].
The pathway and gene network analyses revealed significant alterations in the PPAR, IL6, IL10, FOXO, JAK/STAT and NF-κB signaling pathways and potentially critical roles of IL33, VEGF, INSR, FOS, TGf-β, LIFR, and immunoglobulin that may be relevant to the pathophysiology of OHSS [52][53][54]. Indeed, previous studies also reported the involvement of the immune system, cytokines, and growth factors in the pathogenesis of OHSS [52,53,55,56]. For example, VEGF was implicated as having a significant role in the development of OHSS [26,52,54,57]. IL33 acts as both an extracellular cytokine and an intracellular nuclear factor with transcriptional regulatory properties [58]. IL33 was shown to have higher levels in PCOS patients compared to the controls [59]. Interestingly, dehydroepiandrosterone (DHEA)-induced PCOS in rats was shown to respond to omega-6 fatty acid (γ-linolenic acid (GLA)) treatment that led to a significant decrease in IL-33 levels in the rat's ovaries, hence, implicating GLA's potential use for likely human-related treatments for inflammatory responses in PCOS via the PPAR-γ pathway [60]. Previous studies have reported that the JAK/STAT pathway plays a critical role in the regulation of functions including immune regulation, growth, fertility, and embryogenesis [61][62][63][64].
A limitation of this study was that we had a relatively small sample size within the OHSS group. As the occurrence of OHSS is quite low (about 5% of treated women may encounter moderate and severe OHSS), we were able to recruit six patients with a history of OHSS (n = 6) (four patients had a history of OHSS and two had it in the current cycle) and 18 without any history of OHSS (n = 18) for our whole-genome expression analysis. This was, indeed, achieved after several years of sample collection. For whole-genome studies of such rare human diseases, including PCOS/OHSS, this number is considered a reasonable number to be able to identify differentially expressed genes (DEGs) for discovery-related investigations [6,54,65]. In addition, we used a stringent cutoff for selecting the DEGs; a p-value (of less than 0.01) as well as the fold change criteria. Furthermore, we performed an independent experimental method (qRT-PCR) to validate the microarray results and there was a strong correlation that existed between the microarray and qRT-PCR results (r = 0.9). Future studies may focus on investigating functional mechanisms and further validating the findings on a larger group of patient cohorts.

Conclusions
In conclusion, in this study, we identified markers and altered pathways that may have the potential to differentiate the patients who may be prone to developing OHSS and understanding the biology of its development, which opens new avenues for further research to confirm these findings.
Informed Consent Statement: Written informed consent was obtained from all patients before participation in the study.
Data Availability Statement: All data generated and analyzed in this study are included in this manuscript and its Supplementary Materials files.