Whole transcriptome analysis and construction of gene regulatory networks of granulosa cells from patients with polycystic ovary syndrome (PCOS).

We enrolled Mongolian female patients who received assisted reproductive conception treatment at the Reproductive Center of the Affiliated Hospital of Inner Mongolia Medical University, China. Data were collected with the consent of patients. The inclusion criteria for the PCOS group (3 cases) were based on the diagnostic criteria for PCOS as per the 2003 Rotterdam consensus [ 9 ]. The control group (3 cases) was selected from patients treated for male factor or tubal factor infertility during the same period. The control group was selected according to the criteria of normal ovulatory women, infertility caused by male factor or tubal factor, menstrual regularity (menstrual cycle 21–25 days), follicle-stimulating hormone (FSH) > 10 IU/L, and ultrasound confirmed normal ovarian morphology. There was no significant difference between the two groups in age (< 30 years) and body mass index (18.5 kg/m 2  ≤ BMI) ≤ 25.0 kg/m 2 ). The exclusion criteria for both groups included endometriosis, other endocrine diseases (such as abnormal thyroid function and hyperprolactinemia), and unexplained infertility. All patients underwent their first cycle using the standard long protocol for ovulation. The study was conducted in accordance with the principles outlined in the Declaration of Helsinki. All patients provided informed consent by signing a written consent form. We collected granulosa cells from women who underwent intracytoplasmic sperm injection (ICSI) and fulfilled the enrollment criteria. The cells were separated using hyaluronidase digestion, mechanically stripped using the Pasteur pipette to obtain the granulosa cells, and then washed with PBS. Six samples of RNA were extracted with TRIzol reagent (Invitrogen, USA), and the quality of RNA extracted was qualified. Ribosomal RNA was removed from each sample using the Epicenter Ribo-zTMo TM rRNA Removal Kit (Epicent, USA). cDNA libraries were prepared using the NEBNext ® UltraTM Directional RNA Library Prep Kit for Illumina ® (NEB, USA), as recommended by the manufacturer. The constructed library was sequenced on Illumina HiSeq 2500 platform to obtain raw reads. We performed a sequencing error rate distribution check using FastQC software to assess the quality of the sequencing database. We used the fragments per kilobase of exon model per million mapped reads (FPKM) value to measure the abundance of gene expression, taking into account the impact of gene length and sequencing volume differences on the calculation of gene expression. We used the R package Ballgown/edgeR (v3.14.0) for performing the differential expression analysis of mRNA and lncRNA. We estimated the expression of circRNAs from different transcripts. Screening conditions were statistically significant with |log 2 (fold change)|≥ 1, P -value ≤ 0.05, and Q-value ≤ 0.05. The differentially screened genes were analyzed using the Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) to gain a deeper understanding of the co-regulatory relationships involved in the pathogenic process. The results of GO analysis can indicate the cellular localization of genes biological processes and molecular functions. KEGG analysis results can clarify the basis as an important means of participation. The standard of P < 0.05 is used to define the target base. The threshold for significant correlation between its related functions and pathways. Using differential lncRNA, differential miRNA, and differential genes as the core, additional information on other differential RNAs with targeted relationships was added and compiled into a table. Based on ceRNA analysis, lncRNA, miRNA, and mRNA were extracted, all of which had differences. We used the Cytoscape software to analyze the correlation of DEGs and construct the co-expression network map. The threshold value for the correlation coefficient was set at > 0.999 or < − 0.999 with P  < 0.05.

In this study, we obtained a total of 502,354,062 raw reads through sequencing. After processing for quality control, each library yielded approximately 80,355,901 valid reads. The ratio of valid reads exceeded 95% in all cases. Furthermore, the Q20 values for each library were consistently above 99%, the Q30 values were all above 98%, and the guanine-cytosine content (GC content) ranged from 48 to 50% (Table  1 ). The above results indicate the reliability of the sequencing data. Table 1 Statistical results of clean reads and mapped reads Sample Raw data read Valid data Valid ratio (reads) Q20% Q30% GC content% N1 82651922 78748876 95.28 99.99 98.59 49 N2 88133018 84563631 95.95 99.99 98.62 50 N3 84819850 81267990 95.81 99.99 98.54 48 PCOS1 81404728 77993870 95.81 99.99 98.63 50 PCOS2 83227372 79205136 95.17 99.99 98.55 48 PCOS3 82117170 78643614 95.77 99.99 98.59 49 Statistical results of clean reads and mapped reads Based on the screening conditions of |log 2 (fold change)|≥ 1 and a P -value ≤ 0.05, we obtained 705 differentially expressed genes, of which 331 were upregulated and 374 were downregulated. We found a total of 204 differentially expressed lncRNAs, of which 83 were up-regulated and 121 were down-regulated. These results are shown in Figs.  1 2 . We generated heatmaps for the top 50 differentially expressed genes and lncRNAs, as depicted in Fig.  3 . Fig. 1 The mRNA and lncRNA number of granulosa cells were different between PCOS Group and control group Fig. 2 The volcano of mRNA and lncRNAs which were significantly differentially expressed between PCOS group and control group ( P  < 0.05) Fig. 3 Heatmaps of mRNA and lncRNAs which were significantly differentially expressed between PCOS group and control group ( P  < 0.05) The mRNA and lncRNA number of granulosa cells were different between PCOS Group and control group The volcano of mRNA and lncRNAs which were significantly differentially expressed between PCOS group and control group ( P  < 0.05) Heatmaps of mRNA and lncRNAs which were significantly differentially expressed between PCOS group and control group ( P  < 0.05) Differentially expressed genes (DEGs) were significantly enriched in biological processes such as immune system regulation and cell adhesion, as well as in cellular components including granule membranes and other membrane components. Additionally, they were enriched in cellular functions such as protein-binding, nucleotide-binding, ion-binding, and other cellular functions (Fig.  4 ). Fig. 4 GO function classification of differential genes The horizontal axis represents the pathway, and the vertical axis represents the number of enriched genes ( P  < 0.05) GO function classification of differential genes The horizontal axis represents the pathway, and the vertical axis represents the number of enriched genes ( P  < 0.05) GO annotation of differential lncRNAs target genes revealed a significant enrichment of biological process entries involved in signal transduction, transcriptional regulation, and ion transport. The cellular component entries of biofilm and cytoplasm were associated with these lncRNAs. Furthermore, the molecular function entries identified included protein-binding, ion-channel-binding, and DNA-binding in the granulosa cells of patients with PCOS (Fig.  5 ). Fig. 5 GO function classification of lncRNA target genes. The horizontal axis represents the pathway, and the vertical axis represents the number of enriched genes ( P  < 0.05) GO function classification of lncRNA target genes. The horizontal axis represents the pathway, and the vertical axis represents the number of enriched genes ( P  < 0.05) Genes that were differentially expressed in the granulosa cells of patients with PCOS as compared to normal women were primarily enriched in the Wnt, Toll-like receptor, tumor necrosis factor (TNF), type 1 T helper (Th1), type 2 T helper (Th2) cell differentiation, and transforming growth factor-beta (TGF-β) signaling pathways, as shown in Fig.  6 . KEGG enrichment analysis of target genes of DElnc RNAs, as shown in Fig.  6 , revealed that signaling pathways such as the immune metabolism-related signaling pathway, PPAR, PI3K–AKT, and P53 were involved in granulosa cells. Fig. 6 KEGG classification of differential expression target genes. The horizontal axis represents the enrichment index, while the vertical axis represents the enrichment pathway ( P  < 0.05) KEGG classification of differential expression target genes. The horizontal axis represents the enrichment index, while the vertical axis represents the enrichment pathway ( P  < 0.05) As shown in Table  2 , the mapped sequence ratio of each sample exceeded 96%, and the specificity mapped ratio exceeded 75%. We screened a total of 111 DEcirc RNAs based on the conditions |log 2 (fold change)|≥ 1 and P  ≤ 0.05. Out of these, we found that 20 were upregulated and 91 were downregulated. Table 2 circRNA sequencing data statistics Sample N1 N2 N3 PCOS1 PCOS2 PCOS3 Valid reads 78748876 83680954 81267990 77179262 79205136 77821470 Mapped 76659518 (97.35%) 81274605 (96.87%) 78720436 (97.30%) 75095597 (96.87%) 76831066 (97.30%) 75388892 (96.87%) Unique 59408883 (75.44%) 61136030 (73.06%) 63830814 (78.54%) 58327028 (75.57%) 64224394 (81.09%) 64994954 (83.52%) circRNA sequencing data statistics As seen from Table  3 , the analysis in this study mainly focused on three types of circRNAs, namely exon, intron, and intergenic. Among them, exons accounted for the majority, indicating that the formation sites of circRNAs are located at the boundaries of exons on one end and within the exons on the other end. Introns indicate that the formation sites of circRNAs are exclusively within the intronic region, while intergenics indicate that the formation sites of circRNAs are exclusively within the intergenic region. Table 3 The statistic of circRNA types Type N1 % N2 % N3 % PCOS1 % PCOS2 % PCOS3 % Exon 91.71 92.82 93.8 92.47 94.81 92.59 Intron 5.54 4.67 4.48 4.81 3.60 5.36 Intergenic 2.75 2.51 1.72 2.72 1.58 2.05 The statistic of circRNA types Based on the condition |log 2 (fold change) ≥ 1 and P  ≤ 0.05, there were 88 DEmiRNAs screened in the PCOS group vs. the control group. Out of these, 12 were upregulated and 76 were downregulated. We plotted heatmaps for miRNAs that were significantly differentially expressed, as shown in Fig.  7 . Fig. 7 Heatmap of the miRNA which were differentially expressed between PCOS group and control group ( P  < 0.05). The horizontal axis represents the sample, and the vertical axis represents the clustering pattern (unsupervised clustering) Heatmap of the miRNA which were differentially expressed between PCOS group and control group ( P  < 0.05). The horizontal axis represents the sample, and the vertical axis represents the clustering pattern (unsupervised clustering) We further analyzed the genes within the network by integrating the sequencing data and combining GO and KEGG pathways. We then sequentially placed the higher-expressed RNAs into the lncRNA–mRNA and miRNA–mRNA to screen for the corresponding target genes. This allowed us to co-construct the lncRNA–miRNA–mRNA co-expression network map. As shown in Fig.  8 , the screening yielded 14 pairs of targeting relationships. These relationships involve miRNA (miRNA-150)-lncRNA (lncRNA_FBXO30, lncRNA_SNHG14, and lncRNA_PVT1), as well as mRNA (VLDLR, PPP2R2B, and MYOCD) gene regulatory networks. These networks may have a significant impact on the development and progression of PCOS. Fig. 8 lncRNA–miRNA–mRNA network in granulosa cells. Connection represents the strength of the correlation coefficient. ( P  < 0.05) lncRNA–miRNA–mRNA network in granulosa cells. Connection represents the strength of the correlation coefficient. ( P  < 0.05)

We constructed a network map of the lncRNA–miRNA–mRNA interactions by integrating the sequencing data and conducting further analysis of the genes within the network in combination with the GO and KEGG pathways. In this network, we identified miRNA-150 as a key node, with ncRNAs (lncRNA_FBXO30, lncRNA_SNHG14, and lncRNA_PVT1) and mRNAs (VLDLR, PPP2R2B, and MYOCD) forming a targeting relationship. The target genes of miRNA-150 were found to be enriched in the mTOR signaling pathway, the insulin signaling pathway, and the adipocytokine signaling pathway. Diseases of the immune system and blood system, as well as disruptions in the process of embryonic development, are all linked to the aberrant expression of this miRNA.

PCOS is a common endocrine and metabolic disorder that affects women of reproductive age [ 10 ]. It is a complex polygenic disorder that is influenced by both epigenetic and environmental factors. The communication between granulosa cells and oocytes within follicular microenvironment is crucial for proper follicular development and oocyte quality. In this context, specific cytokines and signaling pathways play pivotal roles in regulating follicular dynamics. Therefore, understanding the pathogenesis of abnormal granulosa cell function in PCOS can provide valuable insights for identifying novel diagnostic and therapeutic biomarkers [ 11 , 12 ]. Recent studies have increasingly focused on elucidating the pathogenic mechanism of PCOS by analyzing the multi-gene set related to the disease. In our study, we employed high-throughput transcriptome sequencing to compare granulosa cells from PCOS patients and women with normal ovulation. We successfully constructed six transcriptomic libraries, covering mRNA, lncRNA, and circRNA, obtained an average of 77 million valid reads per sample, ensuring the quality of the subsequent analysis. An emerging body of research underscores the significance of lncRNA in regulating granulosa cell function and influencing pregnancy outcomes [ 13 ]. In this study, we identified 204 differentially expressed lncRNAs with 83 upregulated and 121 downregulated in PCOS. KEGG enrichment analysis revealed that immune-related pathways, as well as PPAR, PI3K–AKT, and P53 signaling pathways, were significantly involved in granulosa cells. Specifically, the p53 pathway is known for its role in regulating cell cycle arrest and apoptosis in response to DNA damage [ 14 ]. The PI3K/Akt pathway, which is essential for insulin sensitivity, is notably altered in PCOS, particularly in cases associated with insulin resistance (IR), which affects approximately 50% of PCOS patients. Dysregulation of this pathway may contribute to impaired granulosa cell proliferation and follicular development, leading to atresia [ 14 ]. In our study, we observed a significant upregulation of lncRNA PVT1 in the PCOS granulosa cells. PVT1, located on chromosome 8q24, is a well-known oncogenic lncRNA implicated in various malignancies. It has been shown to regulate miRNAs, contributing to cell proliferation, apoptosis, and inflammatory responses [ 15 ]. In the context of PCOS, PVT1 has been linked to the apoptosis of granulosa cells. Gao et al. [ 16 ] demonstrated that PVT1 inhibition, combined with PTEN overexpression, reduced apoptosis in ovarian granulosa cells of PCOS patients. Our findings confirm the elevated expression of PVT1 in PCOS granulosa cells, suggesting it as a potential target for therapeutic intervention. MicroRNAs (miRNAs) are recognized as critical regulators of gene expression, influencing follicular development and oocyte maturation. Altered miRNA expression has been linked to the pathophysiology of PCOS [ 17 , 18 ]. In our study, we identified 88 DEmiRNAs between the PCOS group and the control group, including miR-30b-5p, miR-29b-3p, miR-24-2-5p, miR-150, miR-19a-3p, miR-1839-3p, miR-10399-3p, miR-18b-5p. Notably, miR-30b-5p has been shown to regulates autophagy, and its reduced expression in ovarian tissue of PCOS model rats may contribute to disease progression [ 19 ] Overexpression of miR-30b-5p in granulosa cells from PCOS models significantly inhibited proliferation and induced apoptosis, highlighting its potential as a therapeutic target. CircRNAs are a unique class of RNAs formed by back-splicing, and they play crucial roles in various cellular processes. Recent bioinformatics advances have uncovered their functions in reproductive tissues such as granulosa cells and follicular fluid [ 20 – 22 ]. In our study, we found that more than 90% of the circRNAs were exon-derived, which are recognized for their role in regulating gene expression during cellular differentiation and early embryonic development [ 23 ]. Several differentially expressed circRNAs were enriched in glucose metabolic pathways, suggesting a possible link between circRNA expression and embryo quality, which could have implications for assisted reproductive technology [ 24 ]. Through co-expression network analysis, we identified several key candidate genes, including VLDLR, PPP2R2B, and MYOCD, as well as associated ncRNAs like lncRNA PVT1 and miRNA-150. Of particular interest, PPP2R2B, a regulatory subunit of protein phosphatase 2A (PP2A), was upregulated in PCOS samples. PP2A is involved in numerous signaling pathways, including autophagy regulation [ 25 , 26 ]. Abnormal activation of autophagy pathways has been observed in PCOS ovaries, and dysregulated autophagy may be related to the pathogenesis of PCOS. Our results suggest that the interaction between lncRNA PVT1 and PPP2R2B could be crucial in the abnormal autophagy observed in PCOS granulosa cells, offering new potential targets for future therapeutic strategies.

Polycystic ovary syndrome (PCOS) is a common gynecologic endocrine disease with a prevalence of 5% to 15% in women of reproductive age and representing a primary cause of ovulatory infertility worldwide [ 1 ]. The incidence of PCOS has increased by nearly 65% over the past decade in the China, significantly compromising the reproductive health [ 2 ]. Globally, the incidence of PCOS is becoming more prevalent among younger women, making is as a pressing public health issue. This complex condition arises from the interplay between genetic predisposition, environmental factors and lifestyle factors [ 3 ]. Despite considerable research efforts, the precise mechanisms underlying PCOS remain incompletely understood. Advances in the post-genome era have led to the emergence of various genomics technologies, enabling the exploration of gene expression profiles across a broad range of biological processes. High-throughput transcriptomic techniques are at the forefront of reproductive research, allowing for a comprehensive investigation of gene expression changes [ 4 ]. Granulosa cells (GCs) play a crucial role in mediating intercellular communication with follicular cells and oocytes. Apoptotic factors and signaling pathways expressed in GCs are integral to follicular development, maturation, and oocyte health [ 5 ]. Therefore, understanding the molecular signature of GCs is essential for advancing our knowledge of PCOS and improving therapeutic outcomes. Non-coding RNAs (ncRNAs), including microRNAs (miRNAs), long non-coding RNAs (lncRNAs), and circular RNAs (circRNAs), have emerged as pivotal regulators of gene expression. miRNAs modulate gene expression by targeting specific mRNAs, influencing key processes such as follicular development [ 6 ]. In GCs, dysregulated miRNAs have been implicated in the pathogenesis of PCOS [ 7 ]. lncRNAs, which exceed 200 nucleotides in length, are also implicated in PCOS through their involvement in cell proliferation, apoptosis, and hormone secretion [ 8 ]. Recent studies have shown altered lncRNA expression in PCOS, suggesting their role in disease progression [ 9 ]. Moreover, circRNAs, which function as miRNA sponges, form regulatory networks influencing various biological functions, including ovarian function [ 10 ]. Given the complex etiology of PCOS and the emerging importance of non-coding RNAs in reproductive health, this study aims to explore the expression of mRNAs and ncRNAs in GCs using RNA-seq high-throughput sequencing technology. By comparing the transcriptomic profiles of PCOS patients and normal controls, we aim to identify key molecular markers and pathways involved in PCOS pathogenesis, providing new insights for future therapeutic interventions.

Additional file1. Additional file1.