Analysis of Methylome, Transcriptome, and Lipid Metabolites to Understand the Molecular Abnormalities in Polycystic Ovary Syndrome

Introduction

Polycystic ovary syndrome (PCOS) is one of the most common heterogeneous diseases of endocrine and metabolic origin, and it can seriously affect multiple aspects of the overall health of women.¹ The prevalence of PCOS in women of reproductive age ranges from 5% to 18% according to the 2018 International Evidence-based Guideline.² Patients with PCOS usually have several accompanying differential features such as hyperandrogenemia, insulin resistance, hypertension, overweight, type 2 diabetes, and infertility that have impacts across the patients’ lifespan.^3–5 Studies have shown that 70% women diagnosed with PCOS become infertile within the next few years.⁶ However, despite the severity of the condition, the pathogenesis of PCOS remains elusive.⁷

Granulosa cells (GCs) are a crucial component of the ovary, as they surround the oocyte and play a role in maintaining the homeostasis of the intracellular and extracellular environment.⁸ The gap junctions formed between GCs and oocytes allow the exchange of small molecule metabolites, which is vital for oocyte maturation.⁹ Insulin resistance in PCOS patients is closely associated with increased apoptosis of GCs.¹⁰ Dysfunction of GCs is likely to induce abnormal development of follicles, leading to the occurrence and development of polycystic ovaries. In addition, higher-than-normal levels of androgen in PCOS patients are associated with ovulation disorders, irregular menstruation, infertility, and hirsutism,^11,12 suggesting that androgen hyperactivity is also an important risk factor for PCOS.

Epigenetics is characterized by transmissible and reversible changes in gene function or activity while the DNA sequence remains unchanged.¹³ DNA methylation is involved in many diseases such as cancer, imprinting disorders, and repeat-instability disease and is one of the most widely studied and best-understood epigenetic phenomena.¹⁴ Researchers have observed that few genes such as PRDM1, PPARG, KCTD21, and LHCGR alter DNA methylation in the whole ovarian tissue, adipose tissue, peripheral blood, as well as GCs of PCOS patients.¹⁵ An abundance of differentially methylated CpG sites that were identified as being strongly associated with immunity, and inflammation was found in the ovaries and blood of patients with PCOS.^16,17 Recent studies have shown aberrant hypermethylation sites of gluteofemoral adipose-derived stem cell DNA in PCOS patients with abnormal lipids metabolism.¹⁸ Aberrant DNA methylation may promote or inhibit the dysregulation of pathways related to hormone synthesis, protein metabolism, and lipid metabolism.

In this study, we performed methyl-binding domain sequencing (MBD-seq) and RNA sequencing (RNA-seq) on GCs of PCOS patients and controls, respectively, to identify differentially methylated genes (DMGs) and differentially expressed genes (DEGs) to investigate new biomarkers for the diagnosis and treatment of PCOS. Subsequently, MBD-seq was also performed in prenatally androgenized (PNA) mice (PCOS mice model) and controls to integrate with human GCs MBD-seq sequencing data. The identified gene was verified by methylation-specific PCR (MSP) and real-time quantitative PCR (RT-qPCR). Finally, the serum of PCOS patients and control groups was analyzed by UPLC-MS to identify potential lipid metabolite markers. Our results revealed abnormal DNA methylation and lipid metabolism imbalance in PCOS patients, which provide valuable perspectives for the diagnosis and therapeutic targets of PCOS.

Materials and Methods

Participants

Women with PCOS (n=21) and healthy controls (n=21) were enrolled at Yuncheng Center Hospital in Shanxi Province, China. The participants were non-related and of Chinese Han ethnicity and aged between 20 and 35 years. All participants provided written informed consent. The study was approved by the Ethics Committee on Human Research of Yuncheng Center Hospital (Approval ID: KYLL2019073). All patients diagnosed with PCOS met at least two of the Rotterdam EHSRE/ASRM (2003) criteria.¹⁹ The control group comprised women with male-partner infertility or fallopian tube obstruction infertility. These women underwent in vitro fertilization treatment and had no insulin resistance, polycystic ovary morphology, or history of other chronic diseases, and they all had normal ovulation rhythm and androgen levels.

Clinical and Biochemical Measurements

Demographic and clinical data including height, weight, number of follicles, and preoperative serum hormone levels were collected from all patients, 2–5 days after menstruation. Ten-hour fasting blood samples were collected and centrifuged to separate the serum that was stored at −80°C until further use. The levels of follicle-stimulating hormone (FSH), thyroid-stimulating hormone (TSH), testosterone (T), prolactin (PRL), and luteinizing hormone (LH) were determined by the Beckman UniCel DxI800 fully automatic chemiluminescence analyzer (Beckman Coulter, CA, USA). The clinical parameters of all participants are shown in Table 1 and Supplementary Table S1.

Table 1 Clinical Parameters of the Polycystic Ovarian Syndrome (PCOS) Patients Compared with the Control Groups

Collection of GCs

All participants were subcutaneously injected with gonadotropin‐releasing hormone at the mid-luteal phase (7 days after ovulation) for about 15 days. Follicle size was monitored by ultrasound and serum estradiol. Participants were injected with 6000 IU of human chorionic gonadotropin when two or more follicles reached 15 mm in diameter. After 36 h, ultrasound-guided transvaginal oocyte aspiration was performed to obtain the GCs, then centrifuged and washed twice with Dulbecco’s improved Eagle’s medium. The separated GCs were rapidly preserved at −80°C.

Construction of PNA Mouse Model with PCOS

The animal experiments performed in this study were approved by the Experimental Animal Welfare & Ethics Committee of the School of Biomedical Engineering, Shanghai Jiao Tong University (No. 2019005), and conform to the animal care and use guidelines of Shanghai Jiao Tong University Experimental Animal Center. The Institute of Cancer Research mice were purchased from Shanghai Jiesijie Experimental Animal Co. Ltd., and housed in Shanghai Jiao Tong University Experimental Animal Center. The conditions of food, temperature (23±2°C), humidity (50%-60%), lighting (15–20 lux), and ventilation (0.15m/s) were suitable for survival. Females mated with male mice and the copulatory plug was checked daily. Discovery of the copulatory plug was considered day 1 of gestation. At 16–18 days of gestation, pregnant dams were subcutaneously injected with 70 μL sesame oil containing 350 μg dihydrotestosterone (DHT) (20649, MCE, China) or only 70 μL sesame oil vehicle. The offspring of the DHT-induced group were considered PNA mice, while the offspring of the control group were the control mice. The 8-week-old mice were anesthetized and then euthanized, and the ovary and serum materials were collected.

Assessment of Model Construction and Selection of Mice

PNA mice (n=8) and control mice (n=8) were randomly selected to verify the construction of the PCOS mouse model. We evaluated estrus time, secondary follicle number, and serum testosterone level at 8 weeks of age. The level of serum testosterone and the number of secondary follicles in the PNA group were significantly higher than in the control group, and the estrus time of the PNA mice was also significantly longer than that in the control group. The mice characteristics for model evaluation are shown in Table 2 and Supplementary Table S2.

Table 2 The Characteristics of Prenatal Androgenized (PNA) Mice Compared with the Control Group

Mouse Ovarian Tissue Processing and Histological Examination

After the mice were euthanized, the ovaries were dissected, and one portion was quickly stored at −80°C, while the other was fixed with 4% paraformaldehyde at 4°C overnight. The fixed ovaries were subsequently stored in 70% ethanol for histological examination. The fixed ovarian tissue was dehydrated with alcohol and cleared with xylene, then immersed in molten paraffin for embedding. Briefly, 5-µm-thick paraffin sections were mounted on 3-aminopropyltrisulfide silane-coated glass slides, which were then stained with hematoxylin-eosin (HE) to observe the growth of follicles; the number of follicles was counted.

Nucleic Acid Isolation

Total RNA from GCs samples was isolated using TRIzol (Invitrogen, Waltham, MA, USA), per the manufacturer’s instructions. RNA concentration and purity were determined via Nanodrop 2000 (Thermo Scientific, Waltham, MA, USA). Samples with acceptable absorbance ratios A260/A280 (1.9–2.2) and A260/A230 (2.0–2.4) were considered for further analysis. Next, 1% formaldehyde denaturing gel electrophoresis was employed to access the quality of RNA, and the gel was stained with 10,000×4S Red Plus Nucleic Acid Staining (Sango Biology, Shanghai, China).

Genomic DNA from mouse tissue samples and GCs samples were extracted using the QIAamp DNA Mini Kit (Qiagen, Hilden, Germany) according to the protocol described by the manufacturer. After extraction, the concentration and quality of DNA were quantitatively measured with the Nanodrop 2000 spectrophotometer.

MBD Sequencing and Bioinformatic Analysis

The sequencing data of MBD-seq was derived from our previous studies (https://www.ncbi.nlm.nih.gov/geo/; Accession Number: GSE138573 & GSE156961).^20,21 We performed MBD-seq to ascertain the differentially methylated regions (DMRs) in GCs and mouse ovaries to compare the methylation profiled regions. Briefly, 1 μg genomic DNA per sample was sheared using M220 Focused-ultrasonicator (Covaris, Woburn, MA, USA), and the sheared DNA with methylated regions was affinity-enriched using the MethylMiner^TM Kit (Invitrogen, Waltham, MA, USA), per the manufacturer’s protocols. The affinity-enriched fractions were used to construct an index library with the NEBNext Ultra DNA Library Prep Kit for Illumina (NEB, San Diego, CA, USA), and Agilent 2100 Bioanalyzer (Agilent, USA) was used to assess the quality of the libraries.

For the MBD-seq data of GCs, TrimGalore (version 0.6.7) was used to trim the original data, followed by using fastqc (version 0.11.8) to control the sequencing quality. Then, the clean reads obtained were mapped to the reference genome (Human hg38) using bowtie2 (version 2.3.5.1) with samtools (Release 1.5) to filter the data. DNA methylation datasets from PCOS and controls were analyzed and compared using the MEDIPS package (version 1.48.0). The methylated regions that reached the Adj.P<0.05 and |log₂ (fold change [FC])|≥1 were considered DMRs. For the MBD-seq data of mice, the data processing steps were similar, except that the filtered reads were mapped to the reference genome Mouse mm10. The quality control results of MBD-seq data are shown in Figures S1 and S2. Detailed MBD-seq data regarding the study were shown in Supplementary Table S3.

Transcriptome Sequencing and Data Analysis

The transcriptome sequencing was downloaded from our previous study (https://www.ncbi.nlm.nih.gov/geo/; Accession Number: GSE138518).²⁰ The protocol includes the construction and evaluation of RNA-seq libraries.

The RNA sequencing data was preprocessed and compared with UCSC Human Genome Browser hg38 assembly with HISAT2 (version 2.2.0). Subsequently, based on Ensembl gene-level annotations, Stringtie (version 2.2.0) was employed to ascertain the read counts of genes. The unnormalized counts were combined into a count matrix with R (version 4.2.1) and served as input to DESeq2 (version 1.36.0) and edgeR (version 3.38.4). The DEGs were identified by the intersections of the result of DESeq2 and edgeR analysis, and the expression changes with the value of |log₂ (FC)| ≥ 1. P<0.05 was considered to indicate statistical significance. The quality control results of RNA-seq data are shown in Figure S3. Detailed RNA-seq data regarding the study are presented in Supplementary Table S3.

Real-Time Quantitative PCR (RT-qPCR)

The TACSTD2 expression level was quantified by RT-qPCR. Total RNA (500 ng) was reverse-transcribed by PrimeScript RT Master Mix (RR037, TAKARA, Japan) to generate the complementary DNA (cDNA) for qPCR. The Luna^® Universal qPCR Mix (M3003S, New England Biolabs, USA) was used to carry out the qPCR in QuantStudio 3 (Applied Biosystems) following the manufacturer’s instructions, with glyceraldehyde-3-phosphate dehydrogenase (GAPDH) level as the loading control. The primers used to amplify TACSTD2 are shown in Supplementary Table S4.

Methylation-Specific PCR (MSP)

Genomic DNA of GCs was extracted for bisulfite conversion using EZ DNA Methylation‐Gold™ Kit (D5006, Zymo Research, USA). Unmethylated cytosines in the genome are converted to uracil and thymines after PCR amplification, while methylated cytosines are unchanged. Then, PCR was performed with primers specifically designed for methylation and methylated DNA, respectively. The primers are shown in Supplementary Table S4.

Lipid Metabolites Analysis

First, the collected serum samples were thawed on ice from −80°C, and a 50 μL serum was transferred to a 1.5 mL centrifuge tube. Then, 200 μL pre-cooled methanol, 400 μL chloroform, and 170 μL double-distillation water were added and blended thoroughly for 1 min. Next, the mixture was incubated at 4°C and centrifuged at 13,000 g for 10 min. The lower chloroform layer was collected. To the upper layer, 400 μL solvent mixture (CHCl₃:MeOH:H₂O=85:14:1 [vol:vol:vol]) was added, and after 10 min holding at 4°C, the solution was centrifuged at 13,000 g for 10 min, and the lower chloroform layer collected twice was combined and dried in a vacuum concentrator at a low temperature. Next, the dried samples were redissolved in 100 μL (methylene chloride:isopropyl alcohol:methanol=2:1:1 ([vol:vol:vol]), followed by vortex oscillation for 1 min and centrifugation at 13,000 g for 10 min. Finally, 50 μL of the supernatant was taken for UPLC-MS analysis.

The UPLC-MS analysis mass spectrometry data were compared with the human metabolomics database (Human Metabolome Database). The processed serum data matrix was imported into the SIMCA-P software for principal component analysis (PCA) for quality control of samples and orthogonal partial least squares discriminant analysis (OPLS-DA) modeling analysis to identify differences between groups. The lipid metabolites with variable importance in projection (VIP) score of more than 1, with the value of |log₂ (FC)|≥1 and P<0.05 were identified as differential metabolites. The differential metabolites were visualized via the volcano map and were clustered via the hierarchical clustering heatmap.

Statistical Analysis

Statistical analyses were performed using GraphPad Prism version 8.01 software (GraphPad Software, San Diego, CA, USA) and Statistical Package for Social Sciences (SPSS) version 24.0 software (IBM Corporation, Armonk, NY, USA). The differences between groups were determined by two-sided Student’s t-test, and P-value<0.05 was identified as statistically significant.

Results

Clinical Characteristics of the PCOS Patients and Control Groups

The clinical characteristics of 21 PCOS patients and 21 healthy women are shown in Table 1. The data showed that the level of testosterone (T) and luteinizing hormone (LH) was significantly increased in the PCOS group. The level of follicle-producing hormone (FSH) was decreased slightly, but there was no significant difference. Notably, the LH/FSH ratio of the control groups was 1.00, whereas this ratio was significantly elevated in the PCOS patients (2.11), which is consistent with the clinical features of PCOS.^22,23 There was no significant difference in other characteristics. Subsequently, GCs and the sera of PCOS patients and normal women were collected for further experiments. The complete workflow of this study is presented in Figure 1.

Figure 1 The complete workflow of the study.

Identification of PCOS Patients’ DMRs and Functional Analysis

To confirm the DMRs, genomic DNA of GCs was extracted from PCOS patients (n=3) and control groups (n=2) for standard MBD-seq. After processing and quality control of MBD-seq raw sequencing data of GCs, the methylation regions on the whole genome were analyzed using the MEDIPS package (version 1.48.0). Finally, 1,001,256 DMRs (Adj.P<0.05, |log₂ (FC)|>1) were obtained in PCOS patients, among which 599,396 (59.86%) were hypermethylated and 401,860 (40.14%) were hypomethylated (Supplementary Table S5). The DMRs have added gene annotations to ascertain the distribution in the genome (Figure 2). There was no significant difference in the overall distribution of the hypermethylated and hypomethylated regions, for the most portion in the distal intergenic region, followed by the intron region. However, the proportion of hypomethylated promoter regions (23.88%) was about 1.92-times that of the hypermethylated promoter regions (12.45%), which may have crucial implications for some genes in PCOS. Next, we accessed the distribution of hypermethylated regions’ (red) and hypomethylated regions’ (blue) DMRs in the whole genome (Figure 3A). The results showed that the hypermethylated regions of most chromosomes are higher than hypomethylated regions, with only a few chromosomes exhibiting the opposite result. In addition, we also analyzed the distribution of DMRs located in the CpG and genomic regions (Figure 3B). It is worth mentioning that the CpG in the DMRs of CpG islands, CpG shores, and CpG shelves reportedly have more hypomethylated than hypermethylated regions.

Figure 2 Gene annotation of the differentially methylated regions (DMRs) in GCs of PCOS patients vs control groups. (A) The distribution of hypermethylated DMRs; (B) The distribution of hypomethylated DMRs; (C) The distribution of overall DMRs. Abbreviations: GC, granulosa cell; PCOS, polycystic ovary syndrome.

Figure 3 Distribution and bioinformatic analysis of differentially methylated regions (DMRs) on the genome of PCOS vs control groups (A) The distribution of hypermethylated regions (red) and hypomethylated regions (blue) DMRs in the whole genome; (B) The distribution of DMRs located in CpG and genomic regions; (C) GO analysis of genes with DMRs promoters (D) KEGG analysis of genes with DMRs promoters. Abbreviations: GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

DNA methylation that leads to disease or metabolic disorders usually occurs at multiple CpG sites on CpG islands, many of which are found in the promoters of protein-coding genes.^24,25 Thus, we focused on DMRs in the promoter that may influence the transcription of some key mRNAs in PCOS. A total of 13,526 genes with DMRs in the promoter regions were identified, of which 9367 genes were hypermethylated and 4159 genes were hypomethylated (Supplementary Table S6). Then, we used the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases to perform ontology analysis on the selected DMGs. DMGs in the promoter region of GCs of PCOS patients were mainly enriched in the GO terms associated with catabolic processes and organelle, including small molecule catabolic process, organic acid catabolic process, establishment of organelle localization, and organelle fission (Figure 3C). Furthermore, the KEGG terms mainly included insulin signaling pathway, glycerophospholipid metabolism, and phosphatidylinositol signaling system (Figure 3D).

Ovary Morphology, and Serological Characterization of PNA Mice

The estrus duration, serum testosterone concentration, and number of secondary follicles were statistically monitored in the normal control groups and PNA mice as shown (Figure 4A). The duration of estrus (P<0.05) and the number of secondary follicles (P<0.01) were significantly higher in PNA mice than in normal control mice. Furthermore, we also observed that the serum testosterone level (P<0.001) of PNA mice was significantly increased, which is consistent with sex hormone dysregulation in PCOS patients. To explore the tissue morphology of the mouse ovary, we performed HE staining on the ovaries of the control groups and PNA mice (Figure 4B and C). The HE staining results showed that the ovaries of the control group showed more corpus luteum, indicating that the ovaries had recently undergone smooth ovulation. However, the ovaries of DHT-induced PNA mice showed a distinctive polycystic pattern with fewer corpus luteum, indicating no recent ovulation. Additionally, follicles with an atretic cyst-like appearance were observed in the ovaries of PNA mice, but not in control mice.

Figure 4 Serological and morphological identification of prenatal androgenized (PNA) mice. (A) The proportion of estrus duration, serum testosterone concentration, and number of secondary follicles in PNA mice and control groups (*P<0.05, **P<0.01, ****P<0.0001; multiply the proportion of estrus duration data by 100 to ensure the vertical height); (B) Control group mouse ovary HE staining section; (C) PNA mouse ovary HE staining section; the black arrows indicate the primordial follicles. Abbreviation: HE, hematoxylin-eosin.

Identification of PNA Mice DMRs and Functional Analysis

To ascertain the DMRs of mice, we also extracted genomic DNA from the ovarian tissues of 3 PNA mice and 2 control group subjects by MBD-seq. The MBD-seq data were processed similarly as the previous PCOS data. A total of 68,930 DMRs (Adj.P<0.05, |log₂ (FC)|>1) were obtained from PNA mice, of which 63,774 (88.07%) were hypermethylated and 8638 (11.93%) were hypomethylated (Supplementary Table S7). Similar to human GCs MBD-seq, the distal intergenic region, as well as the intron region, also accounted for the most DMGs. However, the proportion of hypomethylated promoter regions (8.02%) was only about 0.50-times that of hypermethylated promoter regions (16.11%) (Figure 5). We also evaluated the distribution of DMRs on each mouse chromosome (Figure 6A). Each hypermethylated DMRs ratio on chromosomes was greater than that of the hypomethylated DMRs ratio. Interestingly, the distribution of DMRs in CpG regions was much more different than that in human GCs (Figure 6B), indicating the complexity and heterogeneity of PCOS.

Figure 5 Gene annotation of the differentially methylated regions (DMRs) in PNA mice vs control mice. (A) The distribution of hypermethylated DMRs; (B) The distribution of hypomethylated DMRs; (C) The distribution of overall DMRs. Abbreviation: PNA, prenatal androgenized.

Figure 6 Distribution and bioinformatic analysis of differentially methylated regions (DMRs) on the genome of PCOS vs control groups. (A) The distribution of hypermethylated regions (red) and hypomethylated regions (blue) DMRs in the whole genome; (B) The distribution of DMRs located in the CpG and genomic regions; (C) GO analysis of genes with DMRs promoters; (D) KEGG analysis of genes with DMRs promoters. Abbreviations: GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

We also investigated the mouse DMRs located in the gene promoter region. A total of 5753 genes in the promoter region were found, including 5464 hypermethylated genes and 289 hypomethylated genes (Supplementary Table S8). The GO and KEGG databases were also used for the ontology analysis of DMGs selected from mice. The GO terms of DMGs in the promoter region of PCOS mice were mainly enriched in the endomembrane system organization, small GTPase-mediated signal transduction, establishment of organelle localization and cilium organization, and other biological processes (Figure 6C). The KEGG pathway of DMGs in the promoter region of PCOS mice was mainly enriched in tight junctions, AMPK signaling pathway, DNA replication, and fatty acids biosynthesis (Figure 6D). These data indicated that PCOS is an endocrine disease modulated by hormones and closely related to inflammation.

Identification of DEGs and Functional Analysis

We implemented RNA-seq of GCs from PCOS patients (n=3) and controls (n=2) to identify significant DEGs and performed functional enrichment analysis. The hierarchical clustering analysis revealed significant differences in expression patterns between PCOS patients and controls (Figure 7A). By setting the threshold parameters (|log₂ (FC)|>1, P<0.05), 2429 DEGs were identified in total via DESeq2 and edgeR (Figure 7B), of which 872 were significantly upregulated and 1557 were significantly downregulated (Supplementary Table S9). To evaluate the biological significance of DEGs, we performed GO ontology and KEGG pathway analysis. The GO ontology enrichment occurred in several biological and molecular processes, involved in the nuclear division, chromosome segregation, regulation of cell-cell adhesion, and positive regulation of T-cell proliferation (Figure 7C). Moreover, the KEGG terms were mainly enriched in some immune system and inflammation-related pathways, including signaling pathways regulating cytokine-cytokine receptor interaction, phagosome signaling pathway, and inflammatory bowel disease signaling pathway (Figure 7D).

Figure 7 Identification of differentially expressed genes (DEGs) in GCs of PCOS patients vs control groups. (A) The hierarchical cluster heatmaps of DEGs; (B) The Venn diagrams of up-regulated and down-regulated DEGs, respectively; (C) GO function analysis of DEGs; (D) KEGG function analysis of DEGs.

Identification of Hypermethylated Promoter Regions and DEG Downregulation

Integrating the data of differential methylation of promoter regions in ovarian tissue of the PCOS mice and GCs of the PCOS patients, we identified the following 8 genes: CDC42EP4, ERMN, EZR, PIK3R1, ARHGEF18, NECTIN2, TSC2, and TACSTD2 (Table 3). These genes may serve as epigenetic markers for the diagnosis of PCOS. Interestingly, several genes are reportedly involved in the pathogenesis of PCOS such as PIK3R1 and TSC2.^26,27 Analyzing 2429 DEGs with these genes, we found that the expression of TACSTD2 was negatively correlated with the methylation level of the promoter. This result was consistent with the commonly assumed hypermethylation and downregulation of a gene promoter.

Table 3 The Genes Were Identified by the Intersection of PNA Mice and PCOS Patients

Next, we used RT-qPCR and MSP to validate the TACSTD2 promoter methylation and gene expression levels using GCs from PCOS patients (n=5) and controls (n=5). The results of MSP are shown in Figure 8A. The promoter region of TACSTD2 was significantly hypermethylated in PCOS patients. Meanwhile, the representative results of RT-qPCR are shown in Figure 8B. The expression of TACSTD2 was significantly decreased in PCOS samples.

Figure 8 Validation of the promoter methylation and gene expression levels of TACSTD2. (A) MSP 2% agarose gel electrophoresis results; (B) RT-qPCR of TACSTD2 in PCOS patients GCs vs control groups GCs (**P<0.01). Abbreviations: MSP, methylation specific PCR; RT-qPCR, real-time quantitative PCR; M, methylated primer; U, unmethylated primer.

Serum Lipid Metabolites Profiling of PCOS

To compare lipidomics signatures in patients with PCOS, the UPLC-MS method was used on serum samples (n=12). The variable matrix and classification matrix were respectively generated in positive and negative ion modes to construct the OPLS-DA model so that the samples of the PCOS group and the control group could be separated (Supplementary Figure S4). Then, the VIP value of each lipid metabolite in the model was measured, which could be used to determine the ability to influence sample classification and identify differential lipids in PCOS patients’ sera. Based on VIP value>1 in the OPLS-DA model and candidate differential lipids (P<0.05, |log₂ (FC)|>1), 15 differential lipid metabolites were identified (Table 4).²⁸ Notably, some lipid metabolites such as lysophosphatidylcholines (LPC) and triacylglyceride (TG) have been reported as closely relevant to the pathogenesis of PCOS.^29,30 The screening results were visualized by the volcano map (Supplementary Figure S5), which can help to observe the lipid metabolites with significant differences. According to the VIP value in the OPLS-DA model, hierarchical clustering analysis was performed on differential lipids (Figure 9A). There were significant differences in the distribution patterns of differential lipids between the PCOS group and the control group, indicating that the lipid distribution of PCOS patients was similar. To investigate the biological functions of differential lipids in the human body and their impact on the occurrence and development of PCOS diseases, we used MetaboAnalyst 4.0 software to perform KEGG pathway analysis on the differential metabolites (Figure 9B), among which the glycerophospholipid pathway had the greatest impact (impact=0.11182, P=0.001). The results revealed that lipid metabolism was imbalanced in the sera of PCOS patients.

Figure 9 Analysis of serum differential lipids and metabolic pathways in PCOS patients vs control groups. (A) The heatmap of the hierarchical cluster of serum differential lipids; (B) Scatter plots of serum differential lipid KEGG metabolic pathway.

Discussion

DNA methylation, one of the most dominant epigenetic phenomena, often affects physiological processes such as proliferation, differentiation, and apoptosis of cells. As a common and highly heterogeneous reproductive endocrine disease, PCOS is considered influenced by various factors such as environment and inheritable genetics.³¹ Previous studies have shown significant differences in plasma lipidomics in PCOS patients.³² In addition, alteration of DNA methylation has been noticed in subcutaneous adipose tissues, ovarian tissues, GCs, plasma, and other samples from PCOS patients.^33–35 However, studies integrating methylome, transcriptome, and lipid metabolites in PCOS are rare.

In the present study, we screened 1,001,256 DMRs by MBD-seq from PCOS patients (Adj.P<0.05, |log₂ (FC)| > 1). The CpG in the DMRs of CpG islands, CpG shores, and CpG shelves had more hypomethylated than hypermethylated regions. However, except for introns and intergenic regions, which have a higher hypermethylated proportion, the promoter, coding sequence, exons, and other regions all have a higher hypomethylated proportion (Figure 3B). A total of 13,526 DMGs related to the promoter region were identified, which enriched physiological processes including the cell cycle, glycerophospholipid metabolism, and phosphatidylinositol signaling system (Figure 3D). Additionally, the constructed PNA mouse model was also sequenced by MBD-seq, and 68,930 DMRs were obtained. The distribution of DMRs in CpG regions was much more different from that in human GCs (Figure 6B). The ratio of hypermethylated genes in the CpG island, CpG shores, and CpG shelves was higher than that in the hypomethylated genes. In terms of DMRs gene distribution, except for introns and intergenic regions, the proportion of the hypermethylated region was higher than that of the hypomethylated region (Figure 6B). In addition, 5753 DMGs in the promoter region were related to biological functions such as insulin signaling pathways and fatty acid synthesis (Figure 6D). The lipid metabolic pathways are common features of DMGs in GCs and PNA mice, suggesting that DNA methylation is associated with lipid metabolism in PCOS.

Next, we investigated DNA methylation signatures in PCOS patients and PNA mice and identified 8 potential epigenetics markers—CDC42EP4, ERMN, EZR, PIK3R1, ARHGEF18, NECTIN2, TSC2, and TACSTD2 (Table 3). In particular, we noticed several key genes such as PIK3R1, NECTIN2, and TSC2. The PIK3R1 gene encodes the 85-kD phosphoinositide-3-kinase regulatory region, which was identified as a potential pathogenic variant in PCOS and accompanied by severe insulin resistance, lipodystrophy, and obesity.^36–38 NECTIN2 is a nectin cell adhesion molecule 2 associated with lipid-protein, was reported to be a novel marker in carotid atherosclerosis progression;³⁹ young patients with PCOS are likely to feature signs of atherosclerosis.⁴⁰ TSC2, which encodes for the growth protein tuberin, is activated by AMP-activated protein kinase to promote catabolic activities, and its physiological effect is consistent with insulin resistance in PCOS.⁴¹ This suggests that lipodystrophy, obesity, atherosclerosis, and insulin resistance observed in some of the PCOS patients might be induced by these genes.

Comparing the RNA profiles of GCs in PCOS patients and normal ovaries, 872 upregulated and 1557 downregulated RNAs were identified (Figure 7B). The enriched pathways are involved in steroid biosynthesis, ovarian steroidogenesis, and Type I diabetes mellitus, which also correlated with lipid metabolism in PCOS (Figure 7D). By integrating DNA methylation sequencing and transcriptome sequencing data, we found that the expression of TACSTD2 was negatively correlated with the methylation of the promoter region (Figure 8A and B), which was consistent with the results of a previous study.⁴² TACSTD2, an intronless cancer-associated antigen protein-encoding gene, was differentially expressed and identified as an independent marker in ovarian cancer.⁴³ Previous studies have reported that women with PCOS are at increased risk of ovarian cancer,⁴⁴ indicating that TACSTD2 is likely a potential biomarker of PCOS. Perhaps, such a unique gene has the potential to be diagnosed for PCOS, allowing highly sensitive and early detection. However, only few studies are currently focusing on TACSTD2 in PCOS, an attractive possibility that needs further investigation.

Considering that the signaling pathways and biological processes of DMGs and DEGs enrichment in the study are related to lipid metabolism, such as glycerophospholipid metabolism, fatty acid synthesis, and steroid biosynthesis, it is reasonable to postulate that there might be some connection between aberrant DNA methylation and disordered lipid metabolism in PCOS. Thus, combined with the UPLC-MS method of serum metabolic lipids in PCOS patients and normal groups, we found that CL (74:2), PC (16:0/16:1), TG (54:5), and MePC were significantly increased in PCOS, while CL (78:3), LPC (14:0e), SM (d35:3), and TG (52:3e) were significantly decreased (Table 4). Most PCOS patients have abnormal lipid metabolism and are accompanied by hyperandrogenism, chronic anovulation, insulin resistance (IR), and obesity.⁴⁵ CL, containing four fatty acid chains, is unique among phospholipids and plays an important role in energy metabolism. It also only exists in the inner mitochondrial membrane of eukaryotes.^46,47 Moreover, PC, LPC, and SM play important physiological functions in the human body, and the distribution patterns are associated with the developmental status of PCOS.⁴⁸ Furthermore, PC is one of the most stable and abundant (40%–50%) phospholipids in mammalian cell membranes and subcellular organelles, which plays a key intermediate role in glycerophospholipid biosynthesis.^49,50 Studies have shown that LPC (16:0), LPC (18:0), and LPC (18:1) can induce mitochondrial reactive oxygen species in aortic endothelial cells. Moreover, oxidative stress and endothelial cell activation are closely related to the occurrence of PCOS and IR.^51,52 In the process of steroid hormone synthesis, SM plays a role in regulating transcription, acting as a messenger, and performing nuclear and extracellular signal regulation.⁵³ In addition, TG has also been found to be a determinant of hyperandrogenemia in PCOS,²⁹ and hyperandrogenism is considered the defining feature of PCOS,⁵⁴ suggesting that TG is closely associated with the development of PCOS. These differential lipids are consistent with women with PCOS prone to metabolic diseases.

Table 4 Serum Differential Lipids in PCOS Patients Compared with Control Groups

The abnormal activation or inhibition of lipid metabolism pathways in women of reproductive age is closely associated with the occurrence of PCOS.^55–57 Interestingly, in our study, the serum differential metabolites were mainly enriched in the glycerophospholipid metabolism and glycerolipid metabolism pathways (Figure 9B). Furthermore, we noticed that glycerophospholipid, glycerolipid (Figure 3D), fatty acids biosynthesis, insulin signaling pathway (Figure 6D), and other lipid-related signaling pathways were also enriched in DMGs of GCs and PNA mice. In addition, several identified key genes of PIK3R1, NECTIN2, and TSC2 were also closely related to the pathogenesis of PCOS. Moreover, the DEGs of GCs are also enriched in lipid metabolism pathways such as ovarian steroidogenesis and steroid biosynthesis metabolic process (Figure 7C and D). These results indicated that the interaction between aberrant DNA methylation and abnormal lipid metabolites might contribute to the development of PCOS.

Despite many advances in the exploration of reproductive, metabolic, and psychological features, PCOS remains a poorly understood disease in terms of unified standards of diagnosis and therapy.⁵⁸ It is worth noting that there have been several prior multi-omics studies on PCOS that indeed provided important insights into understanding its diagnosis and treatment, but they were either on gut microbiome and serum metabolome, or methylome and competing endogenous axes profiles.^59,60 There are rare multi-omics studies on methylome, transcriptome, and serum lipid metabolites. In our study, combining the methylome and transcriptome data, the 8 DMGs identified are likely to contribute directly or indirectly to the progression of the disease, especially TACSTD2. While function enrichment analysis showed that multi-omics were enriched in lipid metabolism pathways, the 15 lipid metabolites may provide strategies for various therapies for PCOS. Perhaps, PCOS could be accurately diagnosed and targeted for treatment from the abnormal molecules identified in our study. This study also has some limitations. The sample size of the integration analysis is inadequate, and there is no include more comprehensive clinical data on PCOS patients. More sample size and comprehensive clinical data will be beneficial for the analysis of the results.

Conclusion

By integrating methylome and lipid metabolites profiling, we identified 8 potential epigenetic markers and 15 potential lipid metabolite markers for PCOS diagnosis and therapeutic targets. The function analysis indicated that these potential biomarkers were closely related to lipid metabolism. Combined with the transcriptome analysis, we showed that the expression of TACSTD2 was negatively correlated with the methylation of the promoter region. These findings may offer novel insights for an improved understanding of the etiology and diagnosis of PCOS.