Multiomics Analysis Reveals Molecular Abnormalities in Granulosa Cells of Women With Polycystic Ovary Syndrome

Polycystic ovary syndrome (PCOS) is the most common complex endocrine and metabolic disease in women of reproductive age. It is characterized by anovulatory infertility, hormone disorders, and polycystic ovarian morphology. Regarding the importance of granulosa cells (GCs) in the pathogenesis of PCOS, few studies have investigated the etiology at a single “omics” level, such as with an mRNA expression array or methylation profiling assay, but this can provide only limited insights into the biological mechanisms. Here, genome-wide DNA methylation together with lncRNA-miRNA-mRNA profiles were simultaneously detected in GCs of PCOS cases and controls. A total of 3579 lncRNAs, 49 miRNAs, 669 mRNAs, and 890 differentially methylated regions (DMR)-associated genes were differentially expressed between PCOS cases and controls. Pathway analysis indicated that these differentially expressed genes were commonly associated with steroid biosynthesis and metabolism-related signaling, such as glycolysis/gluconeogenesis. In addition, we constructed ceRNA networks and identified some known ceRNA axes, such as lncRNAs-miR-628-5p-CYP11A1/HSD17B7. We also identified many new ceRNA axes, such as lncRNAs-miR-483-5p-GOT2. Interestingly, most ceRNA axes were also closely related to steroid biosynthesis and metabolic pathways. These findings suggest that it is important to systematically consider the role of reproductive and metabolic genes in the pathogenesis of PCOS.

Introduction

Polycystic ovary syndrome (PCOS) is a life-long reproductive, neuroendocrine, and metabolic disorder that affects up to 6–15% of women of reproductive age (Risal et al., 2019). Its main clinical manifestations are ovulatory dysfunction, hyperandrogenemia, and polycystic ovaries, which can lead to infertility (Fauser et al., 2012). In addition to the above reproductive disorders, PCOS is often accompanied by metabolic abnormalities, such as insulin resistance (IR). IR can increase pituitary luteinizing hormone (LH) secretion, testosterone secretion in theca cells, and P450scc activity in granulosa cells (GCs), which interferes with follicle maturation and leads to the development of PCOS (Li et al., 2019). Studies have shown that the abnormal ovarian hormone production is mainly attributed to the hypertrophy of follicular theca cells and the altered expression of key steroid biosynthesis enzymes in GCs (Aste et al., 1998).

As the most abundant cells in the ovary, GCs are closely associated with the development of oocytes and play an essential role in both normal folliculogenesis and steroidogenesis (Hummitzsch et al., 2015). Previous studies have shown that cumulus and mural GCs contribute to the process of oocyte maturation by tight regulation and controlled changes in steroid hormones in the pathogenesis of PCOS (Holesh et al., 2020). Oocytes lack the capacity to carry out some metabolic processes, such as glycolysis and amino acid uptake. They rely on GCs to deliver nutrients and remove waste. In addition, the metabolic profile of GCs is associated with the fate of their accompanying oocyte (Gioacchini et al., 2018; Yilmaz et al., 2018; Fontana et al., 2020).

Previous studies have separately screened differentially expressed mRNAs, miRNAs (DEMs) and lncRNAs (DELs) in GCs to explore the regulatory mechanism of PCOS. Few studies have indicated that genome-wide DNA methylation changes may affect the expression of different genes in PCOS ovaries, as revealed by methylated DNA immunoprecipitation (MeDIP) experiments (Yu et al., 2015; Xu et al., 2016). However, to date, no studies have been performed to identify the whole transcriptome and methylome in same GCs of women with PCOS. In this study, the lncRNA-miRNA-mRNA expression profiles and DNA methylation of GCs in PCOS were comprehensively analyzed. Our goal was to integrate multiomics data to identify differentially expressed genes (DEGs) and differentially methylated regions (DMRs) in PCOS and to construct molecular networks that could help us to better understand the etiology of PCOS.

Materials and Methods

Sample Selection

Five women with PCOS and five age/body mass index (BMI)-matched control subjects from the Center for Reproductive Medicine, Shandong University, Jinan, China, were included in this study. All patients gave informed consent, and the study was approved by the institutional review board of the Reproductive Hospital Affiliated to Shandong University. The definition of PCOS was based on the 2003 Rotterdam criteria (Rotterdam ESHRE/ASRM-Sponsored PCOS Consensus Workshop Group, 2004). Women considered as suffering from PCOS had at least two of the following three characteristics: polycystic ovaries on ultrasound, irregularity/absence of menses, and hyperandrogenism. Cases with congenital adrenal hyperplasia, androgen-secreting tumors, Cushing’s syndrome, thyroid disease, and hyperprolactinemia were excluded. The control group was selected from healthy women who attended the center for IVF with their husbands due to a male factor. All relevant clinical information was obtained from the Electronic Medical Records System. BMI was calculated as weight (kg)/height² (m). Peripheral blood was collected on the 3rd to 5th days of the menstrual cycle to measure serum hormone levels. The levels of follicle-stimulating hormone (FSH), LH, estrogen (E2), prolactin (PRL), and testosterone (T) were measured with a chemiluminescence analyzer (Beckman Coulter, United States). Type B ultrasound was used to determine the antral follicle counts (AFC).

Retrieval of GCs

GCs were collected from follicular fluid obtained via ultrasound-guided transvaginal oocyte retrieval after informed consent had been given by the patients who received the long gonadotropin-releasing hormone agonist protocol. Oocyte retrieval was performed 36 h after human chorionic gonadotropin (hCG) injection by transvaginal ultrasound-guided needle puncture for follicles >15 mm in diameter. At the time of oocyte retrieval, follicular fluid aspirates were collected in sterile tubes and centrifuged. GCs were isolated and purified from the follicular fluid with Ficoll-Percoll (Solarbio, Beijing, China) as previously described (Iwase et al., 2009), and then immediately stored at –80°C for further analysis.

RNA−Seq Analysis and Quality Assessment

Total RNA was extracted using TRIzol Reagent (Invitrogen, CA, United States) and purified using an RNeasy Mini Kit (Qiagen, CA, United States). The quality of RNA was assessed using an Agilent 2100 Bioanalyzer (Agilent, Palo Alto, CA, United States). The rRNA-depleted RNA samples were further processed in accordance with the Illumina protocol (New England Biolabs, Massachusetts, United States). After cDNA synthesis, the samples were sequenced with an Illumina HiSeq 2500 using the paired-end (PE) sequencing strategy. The raw data were recorded. The overall quality of the RNA−seq data was evaluated by FastQC. Clean reads were aligned to the reference genome (Ensembl release 95, Homo sapiens) using TopHat2 (v2.0.14) (Kim et al., 2013) with the default parameters.

MeDIP-Seq Library Construction

MeDIP is a method for immunoprecipitating the methylated portion of the genome using an antibody capable of recognizing 5mC (Wilson and Beck, 2016). Following the manufacturer’s instructions, MeDIP was performed to analyze genome-wide methylation using the Zymo Research DNA Methylation IP Kit (Cat #D5101; Zymo Research, CA, United States). Immunoprecipitated DNA was PCR-amplified, purified, quantified, and sequenced on the Illumina HiSeq 2500 platform. MeDIP−seq reads were mapped to the human genome using BWA software (Li et al., 2012). MACS2 was used to call peaks. To study the DNA methylation differences between two groups, DMRs were identified using the Cummerbund (Trapnell et al., 2012) and ChIPpeakAnno (Zhu et al., 2010) packages in R. Briefly, DMRs were assigned to genomic regions based on gene annotations available from JGI and in-house repeat annotation in GFF3 format. The following gene regions were included: 3’UTR, 5’UTR, promoter, coding DNA sequence (CDS), intron, upstream 1 kb, and downstream 1 kb.

Screening and Clustering Analysis of Differentially Expressed mRNAs, DELs, and DEMs

Data preprocessing and follow-up analysis were performed in the R programming environment (version 3.6.1), and Bioconductor packages were applied for the analysis of DEGs. The lists of DEGs, DELs, and DEMs between controls and PCOS cases were generated using the edgeR package (version 3.32.0) (Robinson et al., 2010). To normalize the raw data, log-fold change | (logFC)| > 1.2 (mRNA and miRNA), | (logFC)| > 2.0 (lncRNA) and p-value <0.05 were considered to indicate statistically significant differences between the PCOS and control groups. To generate an overview of the lncRNA, miRNA and mRNA expression profiles and compare them between the two groups, hierarchical clustering analysis was performed based on the expression levels of all transcripts and significantly differentially expressed transcripts using the pheatmap R package based on Euclidean distance.

lncRNA, miRNA, and mRNA Prediction and Coexpression Network Construction

The miRNA target genes were predicted using the prediction results of the TargetScan and miRcode (Jeggari et al., 2012) databases. The potential target genes transcribed within a 10-kb region upstream or downstream of the lncRNAs were paired and predicted using the UCSC Genome Browser¹ (Song et al., 2019). The expression of differentially expressed mRNAs, DEMs and DELs was analyzed by Pearson’s correlation coefficient using the stats and pcaPP R packages. The miRNA-mRNA network and the lncRNA-mRNA coexpression network were constructed based on analysis of the correlations among the differentially expressed mRNAs, DEMs, and DELs. A p-value of <0.05 for the miRNA-mRNA network and one of <0.01 for the lncRNA-mRNA network were considered statistically significant. The target genes that overlapped with the DEGs were then identified and used to construct the miRNA-mRNA network and lncRNA-mRNA network using Cytoscape software (version 3.8.1) (Saito et al., 2012).

Functional Enrichment Analysis

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis in each module and network was conducted using the Database for Annotation, Visualization and Integrated Discover (DAVID). DEGs and enriched pathways were mapped using KEGG pathway annotation with KOBAS3.0². The top 10 KEGG pathways were selected and ranked by the enrichment factor. To perform literature-based functional analysis, a total of 370 follicle development- and 437 steroid metabolism-related genes were obtained from the Ovarian Kaleidoscope Database (OKDB)³. To identify key transcription factors (TFs), a total of 1496 human TFs were obtained from the Human Protein Atlas Database (HPA)⁴. Subsequently, the Venn diagram tool was used to help identify the common genes that were the focus of this work.

Statistical Analysis

Regarding the clinical characteristics of PCOS patients and controls, quantitative variables are expressed as the mean ± SD. P < 0.05 was considered significant. The clinical data analyses were performed with Statistical Package for Social Science (SPSS 25.0; IBM Corp, Armonk, NY, United States).

Results

Clinical Features

Table 1 presents the basic statistics of both PCOS and control subjects regarding the most important characteristics, such as FSH, LH, E₂, T, P, PRL, and AFC levels, as well as age and BMI. Significant differences between the two groups were found for LH, T, and AFC, all of which had higher levels in PCOS cases (P < 0.05).

TABLE 1

Table 1. Clinical characteristics of women with polycystic ovary syndrome and controls.

Differential Expression Analysis

To identify the DEGs, GCs from five healthy women and five women with PCOS were studied. As indicated in Figure 1A, a correlation plot was used to determine the correlation between samples and to verify the homogeneity between biological replicates. As presented in the histogram in Figure 1B, 669 mRNAs, 49 miRNAs and 3579 lncRNAs were differentially expressed between the PCOS and control groups. Among them, 546 and 123 mRNAs, 31 and 18 miRNAs, and 2226 and 1353 lncRNAs were upregulated and downregulated in PCOS, respectively. Hierarchical clustering heatmaps of the differentially expressed RNAs are shown in Figures 1C,E,G. All differentially expressed mRNAs, DEMs and DELs are listed in Supplementary Table 1. The volcano plots showed the differential expression of mRNAs, miRNAs and lncRNAs between the PCOS group and control group (Figures 1D,F,H). DNA methylation analysis of the MeDIP-seq data showed 890 CpG sites that were differentially methylated in PCOS GCs compared with control GCs (Supplementary Table 1). In terms of the gene structures associated with the CpG sites, the proportions of CDS, intron, downstream 1 kb, upstream 1 kb, 3’UTR, and 5’UTR were 126 (14.16%), 717 (80.56%), 16 (1.8%), 14 (1.57%), 16 (1.8%), and 1 (0.11%), respectively (Figure 1I). We identified 545 hypomethylated and upregulated genes, 19 hypermethylated and downregulated genes, 306 hypermethylated and upregulated genes, and 39 hypomethylated and downregulated genes by integrating the DNA methylation and gene expression data (Figure 1J). Moreover, the chromosomal locations of the DMRs were examined, and they were found to be present on all chromosomes except the Y chromosome (Supplementary Figure 1).

FIGURE 1

Figure 1. Global differentially expressed mRNAs, lncRNAs, miRNAs, and differentially methylated regions (DMRs) identified in PCOS and control granulosa cells. (A) Correlation heatmap between PCOS and control samples. (B) The numbers of differentially expressed mRNAs, differentially expressed lncRNAs (DELs) and differentially expressed miRNAs (DEMs). (C) Hierarchical clustering presentation of DEGs in the PCOS and control groups. (D) Volcano plot of DEGs in the PCOS and control groups. (E) Hierarchical clustering presentation of DELs in the PCOS and control groups. (F) Volcano plot of DELs in the PCOS and control groups. (G) Hierarchical clustering presentation of DEMs in the PCOS and control groups. (H) Volcano plot of DEMs in the PCOS and control groups. (I) The distribution of DMRs. (J) The correlation of DMR-associated genes and mRNA expression.

Functional Enrichment Analysis of DEGs, DELs, DEMs, and DMR-Associated Genes

To investigate the key pathways, the DEGs, DELs, DEMs and DMR-associated genes were evaluated and compared in terms of potential functional pathways in the KEGG database (Supplementary Table 2). As shown in Figure 2A, the results revealed that the DEGs were mainly involved in steroid biosynthesis and many metabolism-related pathways, such as type II diabetes mellitus, glycolysis/gluconeogenesis, carbon metabolism, biosynthesis of amino acids, HIF-1 signaling. Combining the known genes with human TFs, we found that the expression of FOXA1, HIF3A, and STMN1 was upregulated in PCOS GCs (Supplementary Figure 2A). The most significantly enriched biological functions of these TFs were the TGF-beta signaling pathway, IL-17 signaling pathway, endocrine resistance, human T-cell leukemia virus 1 infection, estrogen signaling pathway, Toll-like receptor signaling pathway, IR, GnRH secretion and TNF signaling pathway (Supplementary Figure 2B). For the known genes related to follicle development and steroid metabolism, AMH, FSHR, ESR2, DDX4, and SMAD9 expression was upregulated in the GCs of PCOS patients, while INHA, SOD2, and CYP11A1 expression was downregulated (Supplementary Figures 2C,D). All these functions and pathways have been proven to be closely correlated with the pathogenesis of PCOS.

FIGURE 2

Figure 2. KEGG pathway analysis of RNA-seq and methylation results in PCOS and control GCs. The top 10 affected biofunctions are grouped by disease. (A) Enriched KEGG pathways of DEGs. (B) Enriched KEGG pathways of DELs. (C) Enriched KEGG pathways of DEMs. (D) Enriched KEGG pathways of DMRs.

To further study the role and potential mechanisms of DELs, we identified 124 of their target mRNAs. KEGG analysis identified a total of 40 significantly enriched pathways, including metabolic pathways and steroid biosynthesis (Figure 2B). Notably, the metabolic pathways included carbon metabolism, biosynthesis of amino acids, IR, HIF-1 signaling pathway, terpenoid backbone biosynthesis and glycolysis/gluconeogenesis. Analysis of DEM target genes also identified a number of pathways. Further investigation by KEGG revealed that these miRNAs participated in the regulation of metabolism and steroid synthesis, such as 2-oxocarboxylic acid metabolism, glycerolipid metabolism, arginine and proline metabolism, the FoxO signaling pathway, ovarian steroidogenesis and steroid hormone biosynthesis (Figure 2C). The DMR-associated genes were also involved in metabolic pathways and steroid pathways. Specifically, these genes were associated with glycosaminoglycan biosynthesis, collecting duct acid secretion, bacterial invasion of epithelial cells, adherens junction, steroid biosynthesis, biosynthesis of unsaturated fatty acids and the notch signaling pathway (Figure 2D). Notably, all four omics enrichment analyses identified the steroid biosynthesis pathway. In addition, all three RNA omics analyses revealed the enrichment of metabolic pathways, which are key players in steroidogenesis by acting as a source of energy and substrate for steroid production. These multiomics enrichment results suggest a common etiology of abnormal metabolism and abnormal ovarian steroid formation in GCs of women with PCOS.

Construction of a lncRNA-miRNA-mRNA ceRNA Network

According to the predicted correlations among lncRNAs, miRNAs and mRNAs, a competing endogenous RNA (ceRNA) network was constructed using ceRNA mechanism analysis. The miRNA-mRNA coexpression network was constructed based on the correlation analysis between the DEGs and DEMs. A total of 67 differentially expressed target genes were predicted for 13 DEMs, which were used to construct the miRNA-mRNA coexpression network (Figure 3A). This network included 10 interactions and was associated with metabolic pathways, as determined by searching the KEGG database. The interactions included hsa-miR-548i-SOD2/IDH1, hsa-miR-500a-5p-NSDHL, hsa-miR-483-5p-GOT2, and hsa-miR-214-5p-BRCA1/MKI67. Among all DEM interactions in this regulatory network, FOXO1-hsa-miR-324-5p to DGKA-hsa-miR-148b-5p, FAM160A1-hsa-miR-628-5p, and HOMER2-hsa-miR-130b-5p-PRLR were revealed to represent continuous network connections. Similar to the miRNA-mRNA network, the lncRNA-mRNA coexpression network was constructed based on analysis of the correlation between the DELs and DEGs. In total, 34 lncRNAs and 112 mRNAs involved in 326 interactions were selected to generate the network map (Figure 3B).

FIGURE 3

Figure 3. Construction of the competing endogenous (ceRNA) regulatory network. (A) miRNA-mRNA interaction network. (B) lncRNA-mRNA coexpression network. (C) Sankey diagram of integrative network analysis of multi-RNA-seq data. (D) ceRNA interaction network of miRNA-mRNA-lncRNA interactions. This plot shows the potential regulatory linkage of different RNAs and biological pathways. The four modules represent lncRNAs, miRNAs, mRNAs, and pathways.

There were 217 nodes in the ceRNA network, which consisted of 79 lncRNAs, 6 miRNAs and 11 mRNAs, forming 8 pathways (Supplementary Table 3 and Figure 3C). The top five pathways of lncRNAs, miRNAs and mRNAs are displayed in Figure 3D, indicating the important biological significance of these molecules. KEGG pathway analysis was also performed to determine the involvement of coexpressed genes in different biological pathways. Five pathways overlapped with the enriched genes in the integrated ceRNA network, namely, glycerolipid metabolism, metabolic pathways, biosynthesis of unsaturated fatty acids, steroid biosynthesis, and peroxisome. For example, the AY603498-hsa-miR-628-5p-CYP11A1 and BC036229-hsa-miR-628-5p-HSD17B7 ceRNA axes, which contribute to steroid hormone biosynthesis, were downregulated in PCOS. The AK097578-hsa-miR-548i-IDH1 and AK128202-hsa-miR-483-5p-GOT2 networks were also identified to be associated with metabolic pathways. These analyses identified coexpressed genes that were associated with PCOS development.

Discussion

In the present study, we systematically investigated the differences in the mRNA-miRNA-lncRNA transcriptome and methylation modifications in control and PCOS GCs. Previous studies have focused on only single methylation modifications (Xu et al., 2016; Sagvekar et al., 2019) or single transcriptomics (Jones et al., 2015; Lan et al., 2015) in GCs. In addition, there were some multiomics studies on whole blood (Li et al., 2017), follicular fluid (Naji et al., 2018) and adipose tissue (Kokosar et al., 2016; Pan et al., 2018). However, few studies have performed multiomics analyses of GCs in patients with PCOS. Although many factors have been proven to play important roles in PCOS development in recent decades, no multiomics study has been performed in ovarian GCs. In this study, we systematically investigated control and PCOS ovarian GC mRNA-miRNA-lncRNA-DNA profiles and their potential regulatory networks.

In fact, some studies have shown that in GCs of PCOS patients, there is notable disruption of the entrainment of hormones and metabolic rhythms during the menstrual cycle (Wawrzkiewicz-Jalowiecka et al., 2020). The interaction of glucose/lipid metabolism and steroid synthesis shows obvious effects on both the development and the clinical manifestations of PCOS, mainly by increasing androgen availability, changing the function of GCs and disrupting follicle development (Pasquali et al., 2006). Follicle development depends on the synchronization of oocyte maturation and GC proliferation and differentiation. At the same time, the maturation of oocytes relies on the steroids and nutrients provided by GCs (Sutton-McDowall et al., 2010). From the perspective of follicle development, some of the DEGs identified in this study were related to androgen excess, impaired ovulation, and oxidative stress. Examples of such DEGs include CYP11A1, HSD17B7, and FOXO1, which have been identified to participate in the occurrence and development of PCOS (Sagvekar et al., 2019). PCOS is also characterized by IR and hyperinsulinemia. In PCOS, we also identified the involvement of some metabolic genes, such as IRS1 and INSR, showing increased expression, and IDH1 and GOT2, showing decreased expression at the transcriptional level, which may confer a genetic predisposition to developing this condition (Feng et al., 2015). Of note, TF analysis showed a higher level of the ZBTB16 gene in PCOS GCs, which is consistent with the findings of recent PCOS susceptibility gene studies. A large-scale genome-wide meta-analysis of PCOS patients of European ancestry suggests that a variant at the ZBTB16 locus was strongly associated with ovulatory dysfunction and polycystic ovarian morphology (Day et al., 2018). The SNP rs1784692 in the ZBTB16 gene was associated with PCOS and BMI levels in Han Chinese women (Yang et al., 2020).

PCOS is a complex and heterogeneous condition that results from the interaction of diverse genetic and environmental factors (Rosenfield and Ehrmann, 2016). In our study, we identified some key common enriched pathways, including glycolysis/gluconeogenesis, steroid biosynthesis, and IR, in the KEGG analysis of lncRNA-miRNA-mRNA interactions and DMR-associated genes. By constructing a ceRNA network, we also observed that the most highly enriched ceRNA axes might play a role in regulating metabolic pathways and steroid biosynthesis in the development of PCOS. Among them, both lncRNA-miR-628-5p-CYP11A1 and lncRNA-miR-628-5p-HSD17B7 ceRNA regulatory axes are associated with steroid hormone biosynthesis and metabolic pathways. A differential expression study showed that an increase in miR-628-5p serum levels at 20 weeks of gestation was observed in women who developed severe preeclampsia (Martinez-Fierro et al., 2019). In this study, the miR-628-5p-CYP11A1/HSD17B7 network was downregulated in PCOS GCs. A randomized clinical trial showed that frozen embryo transfer resulted in an increased risk of preeclampsia in twin pregnancy in women with PCOS (Zhang et al., 2018). Other preliminary evidence described that women with PCOS and the highest maternal testosterone levels in the early second trimester had the highest risk of developing preeclampsia (Valdimarsdottir et al., 2020). A transcriptome study showed that HSD17B7 and CYP11A1 expression was repressed in DHT-treated ovaries and that the dysregulation of HSD17B7 and CYP11A1 expression was associated with the biosynthesis and metabolism of steroids and cholesterol and lipids (Salilew-Wondim et al., 2015). Therefore, miR-628-5p may be an important hub regulator of HSD17B7 and CYP11A1 and may increase the risk of pregnancy complications by affecting steroid hormone biosynthesis and metabolic pathways in PCOS patients.

We also observed a downregulated axis consisting of lncRNAs-miR-483-5p and GOT2 associated with metabolic pathways in the ceRNA network. A previous study reported that miR-483-5p expression was significantly decreased in cumulus cells of PCOS patients (Shi et al., 2015). Other miRNA expression profiles revealed that miR-483-5p can regulate Notch3/MAPK3 expression (Xu et al., 2015) and progesterone concentrations (Sang et al., 2013) in cumulus GCs and follicular fluid of PCOS patients. Although few studies have focused on the function of GOT2 in PCOS, an important conclusion was reached by Yang et al. (2015), who found that GOT2 participates in mitochondrial metabolism through acetylation (Borst, 2020). Moreover, miR-483-5p is associated with future onset of both diabetes and cardiovascular disease (Gallo et al., 2018) and increases hepatic LDL receptor levels by inhibiting PCSK9 production (Dong et al., 2020). According to these results, it can be speculated that miR-483-5p may regulate GOT2 to contribute to the IR of PCOS and is worthy of further investigation.

In summary, the results show that women with PCOS have multiple transcriptional and epigenetic changes in GCs that are related to steroid hormone synthesis and metabolic pathways. Several genes and pathways, such as lncRNAs-miR-628-5p-CYP11A1/HSD17B7 and lncRNAs-miR-483-5p-GOT2, play important roles in the etiology of PCOS and may be novel candidate biomarkers or treatment targets for PCOS.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI GEO GSE168404.

Ethics Statement

The studies involving human participants were reviewed and approved by Institutional review board of the Reproductive Hospital Affiliated to Shandong University. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

HZ contributed to conception and design of the study. SGZ revised the manuscript. RSZ performed bioinformatics analysis and wrote the manuscript. YHJ supervised the bioinformatics analysis and edited the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This research was supported by the National Key Research and Development Program of China (2018YFC1004303), the Basic Science Center Program (31988101) and the National Program (82071606, 31871509) of NSFC, and the Foundation for Distinguished Young Scholars of Shandong Province (JQ201816).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We sincerely thank the patients for their participation. We thank openbiox community and Hiplot team (https://hiplot.com.cn) for providing technical assistance and valuable tools for data analysis and visualization. We also thank Honghui Zhang and Yu Zhang from Shandong University, China for their proofreading of this manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.648701/full#supplementary-material

Supplementary Figure 1 | Chromosome graph of differentially methylated regions (DMRs) between PCOS and control GCs. Different colors represent different types of DMR loci.

Supplementary Figure 2 | Specific pathway analysis of differentially expressed genes (DEGs). (A) Heatmap of DEGs associated with transcription factor regulation. (B) KEGG pathways of DEGs associated with transcription factor regulation. (C) Heatmap of DEGs associated with follicle development. (D) Heatmap of DEGs associated with steroid metabolism.

Supplementary Table 1 | Differentially expressed miRNA-mRNA-lncRNA-CpG-associated genes.

Supplementary Table 2 | KEGG pathways of miRNA-mRNA-lncRNA-CpG-associated genes.

Supplementary Table 3 | Integrated analysis of identified signature genes in the miRNA-mRNA-lncRNA network.

Footnotes

1. ^ http://genome.ucsc.edu/
2. ^ http://kobas.cbi.pku.edu.cn/kobas3/
3. ^ http://okdb.appliedbioinfo.net/
4. ^ https://www.proteinatlas.org/

References

Aste, N., Pau, M., and Biggio, P. (1998). Kerion celsi: a clinical epidemiological study. Mycoses 41, 169–173. doi: 10.1111/j.1439-0507.1998.tb00319.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Borst, P. (2020). The malate-aspartate shuttle (Borst cycle): how it started and developed into a major metabolic pathway. IUBMB Life 72, 2241–2259. doi: 10.1002/iub.2367

PubMed Abstract | CrossRef Full Text | Google Scholar

Day, F., Karaderi, T., Jones, M. R., Meun, C., He, C., Drong, A., et al. (2018). Large-scale genome-wide meta-analysis of polycystic ovary syndrome suggests shared genetic architecture for different diagnosis criteria. PLoS Genet. 14:e1007813. doi: 10.1371/journal.pgen.1007813

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, J., He, M., Li, J., Pessentheiner, A., Wang, C., Zhang, J., et al. (2020). microRNA-483 ameliorates hypercholesterolemia by inhibiting PCSK9 production. JCI Insight 5:e143812. doi: 10.1172/jci.insight.143812

PubMed Abstract | CrossRef Full Text | Google Scholar

Fauser, B. C., Tarlatzis, B. C., Rebar, R. W., Legro, R. S., Balen, A. H., Lobo, R., et al. (2012). Consensus on women’s health aspects of polycystic ovary syndrome (PCOS): the Amsterdam ESHRE/ASRM-Sponsored 3rd PCOS consensus workshop group. Fertil. Steril. 97, 28–38. doi: 10.1016/j.fertnstert.2011.09.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, C., Lv, P. P., Yu, T. T., Jin, M., Shen, J. M., Wang, X., et al. (2015). The association between polymorphism of INSR and polycystic ovary syndrome: a meta-analysis. Int. J. Mol. Sci. 16, 2403–2425. doi: 10.3390/ijms16022403

PubMed Abstract | CrossRef Full Text | Google Scholar

Fontana, J., Martinkova, S., Petr, J., Zalmanova, T., and Trnka, J. (2020). Metabolic cooperation in the ovarian follicle. Physiol. Res. 69, 33–48. doi: 10.33549/physiolres.934233

PubMed Abstract | CrossRef Full Text | Google Scholar

Gallo, W., Esguerra, J., Eliasson, L., and Melander, O. (2018). miR-483-5p associates with obesity and insulin resistance and independently associates with new onset diabetes mellitus and cardiovascular disease. PLoS One 13:e206974. doi: 10.1371/journal.pone.0206974

PubMed Abstract | CrossRef Full Text | Google Scholar

Gioacchini, G., Notarstefano, V., Sereni, E., Zaca, C., Coticchio, G., Giorgini, E., et al. (2018). Does the molecular and metabolic profile of human granulosa cells correlate with oocyte fate? New insights by Fourier transform infrared microspectroscopy analysis. Mol. Hum. Reprod. 24, 521–532. doi: 10.1093/molehr/gay035

PubMed Abstract | CrossRef Full Text | Google Scholar

Holesh, J. E., Bass, A. N., and Lord, M. (2020). Physiology, Ovulation. Treasure Island, FL: StatPearls Publishing.

Google Scholar

Hummitzsch, K., Anderson, R. A., Wilhelm, D., Wu, J., Telfer, E. E., Russell, D. L., et al. (2015). Stem cells, progenitor cells, and lineage decisions in the ovary. Endocr. Rev. 36, 65–91. doi: 10.1210/er.2014-1079

PubMed Abstract | CrossRef Full Text | Google Scholar

Iwase, A., Goto, M., Harata, T., Takigawa, S., Nakahara, T., Suzuki, K., et al. (2009). Insulin attenuates the insulin-like growth factor-I (IGF-I)-Akt pathway, not IGF-I-extracellularly regulated kinase pathway, in luteinized granulosa cells with an increase in PTEN. J. Clin. Endocrinol. Metab. 94, 2184–2191. doi: 10.1210/jc.2008-1948

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeggari, A., Marks, D. S., and Larsson, E. (2012). miRcode: a map of putative microRNA target sites in the long non-coding transcriptome. Bioinformatics 28, 2062–2063. doi: 10.1093/bioinformatics/bts344

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, M. R., Brower, M. A., Xu, N., Cui, J., Mengesha, E., Chen, Y. D., et al. (2015). Systems genetics reveals the functional context of PCOS loci and identifies genetic and molecular mechanisms of disease heterogeneity. PLoS Genet. 11:e1005455. doi: 10.1371/journal.pgen.1005455

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, D., Pertea, G., Trapnell, C., Pimentel, H., Kelley, R., and Salzberg, S. L. (2013). TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 14:R36. doi: 10.1186/gb-2013-14-4-r36

PubMed Abstract | CrossRef Full Text | Google Scholar

Kokosar, M., Benrick, A., Perfilyev, A., Fornes, R., Nilsson, E., Maliqueo, M., et al. (2016). Epigenetic and transcriptional alterations in human adipose tissue of polycystic ovary syndrome. Sci. Rep. 6:22883. doi: 10.1038/srep22883

PubMed Abstract | CrossRef Full Text | Google Scholar

Lan, C. W., Chen, M. J., Tai, K. Y., Yu, D. C., Yang, Y. C., Jan, P. S., et al. (2015). Functional microarray analysis of differentially expressed genes in granulosa cells from women with polycystic ovary syndrome related to MAPK/ERK signaling. Sci. Rep. 5:14994. doi: 10.1038/srep14994

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, M., Wu, H., Luo, Z., Xia, Y., Guan, J., Wang, T., et al. (2012). An atlas of DNA methylomes in porcine adipose and muscle tissues. Nat. Commun. 3:850. doi: 10.1038/ncomms1854

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, S., Zhu, D., Duan, H., Ren, A., Glintborg, D., Andersen, M., et al. (2017). Differential DNA methylation patterns of polycystic ovarian syndrome in whole blood of Chinese women. Oncotarget 8, 20656–20666. doi: 10.18632/oncotarget.9327

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Zhu, Q., Wang, W., Qi, J., He, Y., Wang, Y., et al. (2019). Elevated chemerin induces insulin resistance in human granulosa-lutein cells from polycystic ovary syndrome patients. FASEB J. 33, 11303–11313. doi: 10.1096/fj.201802829R

PubMed Abstract | CrossRef Full Text | Google Scholar

Martinez-Fierro, M. L., Carrillo-Arriaga, J. G., Luevano, M., Lugo-Trampe, A., Delgado-Enciso, I., Rodriguez-Sanchez, I. P., et al. (2019). Serum levels of miR-628-3p and miR-628-5p during the early pregnancy are increased in women who subsequently develop preeclampsia. Pregnancy Hypertens. 16, 120–125. doi: 10.1016/j.preghy.2019.03.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Naji, M., Nekoonam, S., Aleyasin, A., Arefian, E., Mahdian, R., Azizi, E., et al. (2018). Expression of miR-15a, miR-145, and miR-182 in granulosa-lutein cells, follicular fluid, and serum of women with polycystic ovary syndrome (PCOS). Arch. Gynecol. Obstet. 297, 221–231. doi: 10.1007/s00404-017-4570-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, J. X., Tan, Y. J., Wang, F. F., Hou, N. N., Xiang, Y. Q., Zhang, J. Y., et al. (2018). Aberrant expression and DNA methylation of lipid metabolism genes in PCOS: a new insight into its pathogenesis. Clin. Epigenetics 10:6. doi: 10.1186/s13148-018-0442-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Pasquali, R., Gambineri, A., and Pagotto, U. (2006). The impact of obesity on reproduction in women with polycystic ovary syndrome. BJOG 113, 1148–1159. doi: 10.1111/j.1471-0528.2006.00990.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Risal, S., Pei, Y., Lu, H., Manti, M., Fornes, R., Pui, H. P., et al. (2019). Prenatal androgen exposure and transgenerational susceptibility to polycystic ovary syndrome. Nat. Med. 25, 1894–1904. doi: 10.1038/s41591-019-0666-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenfield, R. L., and Ehrmann, D. A. (2016). The Pathogenesis of Polycystic Ovary Syndrome (PCOS): the hypothesis of PCOS as functional ovarian hyperandrogenism revisited. Endocr. Rev. 37, 467–520. doi: 10.1210/er.2015-1104

PubMed Abstract | CrossRef Full Text | Google Scholar

Rotterdam ESHRE/ASRM-Sponsored PCOS Consensus Workshop Group (2004). Revised 2003 consensus on diagnostic criteria and long-term health risks related to polycystic ovary syndrome (PCOS). Hum. Reprod. 19, 41–47. doi: 10.1093/humrep/deh098

PubMed Abstract | CrossRef Full Text | Google Scholar

Sagvekar, P., Kumar, P., Mangoli, V., Desai, S., and Mukherjee, S. (2019). DNA methylome profiling of granulosa cells reveals altered methylation in genes regulating vital ovarian functions in polycystic ovary syndrome. Clin. Epigenetics 11:61. doi: 10.1186/s13148-019-0657-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Saito, R., Smoot, M. E., Ono, K., Ruscheinski, J., Wang, P. L., Lotia, S., et al. (2012). A travel guide to Cytoscape plugins. Nat. Methods 9, 1069–1076. doi: 10.1038/nmeth.2212

PubMed Abstract | CrossRef Full Text | Google Scholar

Salilew-Wondim, D., Wang, Q., Tesfaye, D., Schellander, K., Hoelker, M., Hossain, M. M., et al. (2015). Polycystic ovarian syndrome is accompanied by repression of gene signatures associated with biosynthesis and metabolism of steroids, cholesterol and lipids. J. Ovarian Res. 8:24. doi: 10.1186/s13048-015-0151-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Sang, Q., Yao, Z., Wang, H., Feng, R., Wang, H., Zhao, X., et al. (2013). Identification of microRNAs in human follicular fluid: characterization of microRNAs that govern steroidogenesis in vitro and are associated with polycystic ovary syndrome in vivo. J. Clin. Endocrinol. Metab. 98, 3068–3079. doi: 10.1210/jc.2013-1715

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, L., Liu, S., Zhao, W., and Shi, J. (2015). miR-483-5p and miR-486-5p are down-regulated in cumulus cells of metaphase II oocytes from women with polycystic ovary syndrome. Reprod. Biomed. Online 31, 565–572. doi: 10.1016/j.rbmo.2015.06.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, J., Lu, Y., Sun, W., Han, M., Zhang, Y., and Zhang, J. (2019). Changing expression profiles of lncRNAs, circRNAs and mRNAs in esophageal squamous carcinoma. Oncol. Lett. 18, 5363–5373. doi: 10.3892/ol.2019.10880

PubMed Abstract | CrossRef Full Text | Google Scholar

Sutton-McDowall, M. L., Gilchrist, R. B., and Thompson, J. G. (2010). The pivotal role of glucose metabolism in determining oocyte developmental competence. Reproduction 139, 685–695. doi: 10.1530/REP-09-0345

PubMed Abstract | CrossRef Full Text | Google Scholar

Trapnell, C., Roberts, A., Goff, L., Pertea, G., Kim, D., Kelley, D. R., et al. (2012). Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat. Protoc. 7, 562–578. doi: 10.1038/nprot.2012.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Valdimarsdottir, R., Wikstrom, A. K., Kallak, T. K., Elenis, E., Axelsson, O., Preissl, H., et al. (2020). Pregnancy outcome in women with polycystic ovary syndrome in relation to second-trimester testosterone levels. Reprod. Biomed. Online 42, 217–225. doi: 10.1016/j.rbmo.2020.09.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Wawrzkiewicz-Jalowiecka, A., Kowalczyk, K., Trybek, P., Jarosz, T., Radosz, P., Setlak, M., et al. (2020). In search of new therapeutics-molecular aspects of the PCOS pathophysiology: genetics, hormones, metabolism and beyond. Int. J. Mol. Sci. 21:7054. doi: 10.3390/ijms21197054

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, G. A., and Beck, S. (2016). “Computational analysis and integration of MeDIP-seq methylome data,” in Next Generation Sequencing: Advances, Applications and Challenges, ed. J. K. Kulski (Rijeka: InTech), 159–169.

Google Scholar

Xu, B., Zhang, Y. W., Tong, X. H., and Liu, Y. S. (2015). Characterization of microRNA profile in human cumulus granulosa cells: identification of microRNAs that regulate Notch signaling and are associated with PCOS. Mol. Cell Endocrinol. 404, 26–36. doi: 10.1016/j.mce.2015.01.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, J., Bao, X., Peng, Z., Wang, L., Du, L., Niu, W., et al. (2016). Comprehensive analysis of genome-wide DNA methylation across human polycystic ovary syndrome ovary granulosa cell. Oncotarget 7, 27899–27909. doi: 10.18632/oncotarget.8544

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, H., Zhou, L., Shi, Q., Zhao, Y., Lin, H., Zhang, M., et al. (2015). SIRT3-dependent GOT2 acetylation status affects the malate-aspartate NADH shuttle activity and pancreatic tumor growth. EMBO J. 34, 1110–1125. doi: 10.15252/embj.201591041

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J., Zhao, R., Li, L., Li, G., Yang, P., Ma, J., et al. (2020). Verification of a ZBTB16 variant in polycystic ovary syndrome patients. Reprod. Biomed. Online 41, 724–728. doi: 10.1016/j.rbmo.2020.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Yilmaz, B., Vellanki, P., Ata, B., and Yildiz, B. O. (2018). Metabolic syndrome, hypertension, and hyperlipidemia in mothers, fathers, sisters, and brothers of women with polycystic ovary syndrome: a systematic review and meta-analysis. Fertil. Steril. 109, 356–364. doi: 10.1016/j.fertnstert.2017.10.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, Y. Y., Sun, C. X., Liu, Y. K., Li, Y., Wang, L., and Zhang, W. (2015). Genome-wide screen of ovary-specific DNA methylation in polycystic ovary syndrome. Fertil. Steril. 104, 145–153. doi: 10.1016/j.fertnstert.2015.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., Wei, D., Legro, R. S., Shi, Y., Li, J., Zhang, L., et al. (2018). Obstetric complications after frozen versus fresh embryo transfer in women with polycystic ovary syndrome: results from a randomized trial. Fertil. Steril. 109, 324–329. doi: 10.1016/j.fertnstert.2017.10.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, L. J., Gazin, C., Lawson, N. D., Pages, H., Lin, S. M., Lapointe, D. S., et al. (2010). ChIPpeakAnno: a bioconductor package to annotate ChIP-seq and ChIP-chip data. BMC Bioinformatics 11:237. doi: 10.1186/1471-2105-11-237

PubMed Abstract | CrossRef Full Text | Google Scholar