Combined Multiomics Analysis Reveals the Role of ANXA1 Methylation and miRNA-Targeted LDLR in Polycystic Ovary Syndrome.

This study has been approved by the Human Research Ethics Committee of Yuncheng Central Hospital (Approval number: KYLL2019073). The animal experiments performed in the 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. All procedures followed were in accordance with t the Helsinki Declaration of 1964 and its later amendments. Animal Studies: The mice used in the study were purchased from the Institute of Cancer Research (ICR) Animal Experimental Center of Yangzhou University. The mice were housed in the Drum Tower Hospital Animal Experimental Center with constant temperature, humidity, and light, and had free access to water and food. All of the animal experiments in this study followed the guidelines of the Laboratory Animal Control Committee (Jiangsu Province, China).

Studies have shown that compared to other PCOS animal models, the mRNA levels of PNA mice are closest to those of PCOS patients [ 25 ]. Therefore, PNA mice were selected as PCOS animal samples in this study. To confirm the establishment of the PNA mouse model, we randomly selected six PNA mice and six control mice and measured their estrus time, serum testosterone levels, and number of secondary follicles at 8 weeks of age (Table  1 ), which in the PNA mouse group were significantly higher than those in the control group. It is evident that we have successfully constructed the PNA mouse model. Comparison of characteristics of PNA group and control group. Abbreviation: PNA, prenatally androgenized. A total of 17 patients with PCOS and 17 control women of normal childbearing potential were candidates for ovarian GC samples. Patients with PCOS were selected according to the Rotterdam criteria [ 7 ]. Table  2 lists the clinical data for all participants. Comparing the clinical data of the PCOS patient group with the control group, the body mass index (BMI) index and LH levels of PCOS patients were significantly higher than those of the control group ( p  < 0.05), and there was no significant difference in other characteristics. Comparison of characteristics of PCOS Clinical Patients and controls. Abbreviations: β‐HCG, beta‐human chorionic gonadotropin; BMI, body mass index; FBS, fasting blood sugar; E2, estrogen; LH, luteinizing hormone; P, progesterone; PRL, prolactin; T, testosterone. High‐throughput sequencing MBD‐seq data from ovarian tissue samples from three PNA mice and two control mice were preprocessed. After obtaining . bam files of methylation level in these ovarian tissue samples, genome‐wide differentially methylated region (DMR) analysis was performed with a screening threshold of fold ratio ≥ 2 and p ‐value < 0.05. A total of 26 174 DMRs were obtained, of which 5411 were hypermethylated and 20 763 were hypomethylated (Tables  S3 , S4 ). There were significantly more hypomethylated regions than hypermethylated regions, and the proportion of hypomethylated regions on each chromosome was higher than that of hypermethylated regions (Figure  1a ). General analysis of differentially methylated patterns. (a) Distribution of DMRs on chromosomes in mouse group; (b) Distribution of DMRs on chromosomes in human group. ( X ‐axis: chromosome number; Y ‐axis: Number of DMRs; Red bar: Hypermethylated DMRs; Green bar: hypomethylated DMRs.) Similarly, the high‐throughput sequencing MBD‐seq data from three ovarian GC samples from PCOS patients and two control ovarian GC samples were preprocessed. Finally, 1 264 800 DMRs were screened, of which 504 605 were hypermethylated and 760 195 were hypomethylated. Figure  1b shows the distribution of hypermethylated and hypomethylated regions on each chromosome, showing that the proportion of hypermethylated regions on chromosomes 17 and 19 is higher than that of hypomethylated regions, but the proportion of hypomethylated regions on other chromosomes is higher than that of hypermethylated regions. To further understand the distribution of DMRs in genomic characteristic regions, 26 174 DMRs were annotated and visualized statistically. As shown in Figure  2a,b , DMRs were mainly distributed in the distal intergenic region and intron region, including the first intron and other introns. A total of 2574 (9.83%) DMRs were located in the promoter region, including 456 (1.74%) hypermethylation regions and 2118 (8.09%) hypomethylation regions. There are much more hypomethylated regions than hypermethylated regions in the promoter region. For human, similarly, as shown in Figure  2c,d , DMRs were mainly distributed in the distal intergenic region and intron region. Compared with the mouse group, the proportion of DMRs located in the promoter region was significantly higher in the human group, and the number of hypermethylated regions was higher than hypomethylation regions in the promoter region. Analysis of differentially methylated patterns in different gene regions. (a) Hypermethylated region distribution in mouse group; (b) Hypomethylated region distribution in mouse group; (c) Hypermethylated region distribution in human group; (d) Hypomethylated region distribution in human group. (Color: Different gene regions.) We annotated DMRs and mapped 456 promoter regions with hypermethylation to 334 genes and 2118 promoter regions with hypomethylation to 1386 genes. In order to explore the function and significance of differentially promoter‐methylated genes (DMGs) in the process of PCOS, the promoter‐hypermethylation genes and the promoter‐hypomethylation genes were analyzed respectively (Figure  3 ). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed that the function of promoter‐differentially‐methylated genes is enriched in cell proliferation and migration, such as regulation of mitotic cell cycle (GO:0007346) and nuclear division (GO:0000280). It is mainly the promoter‐hypomethylation genes that are enriched in these biological pathways. The functions of these genes overlap with those of DEGs previously analyzed, suggesting that there is an interaction and co‐regulation of methylation levels and gene expression levels in the same pathway. Enrichment analysis of DMGs' feature in mouse group. (a) and (b) GO analysis and KEGG analysis of promoter‐hypermethylation genes, respectively; (c) and (d) GO analysis and KEGG analysis of promoter‐hypomethylation genes, respectively. (For GO analysis, X ‐axis: Gene counts in each GO term; Y ‐axis: GO terms; Blue bar: Biological process; Orange bar: Cellular component; Green bar: Molecular function. For KEGG analysis, X ‐axis: Gene ratio for each description; Y ‐axis: Description in KEGG; Circular diameter: Gene counts; Color temperature: p ‐value.) For human group, by annotating the DMRs, 121 946 promoter regions with hypermethylation were mapped to 17 843 genes, and 102 738 promoter regions with hypomethylation were mapped to 19 950 genes. The results of promoter‐hypermethylation gene function analysis are shown in (Figure  4a,b ). The functions of promoter‐hypermethylation genes are mainly enriched in cell growth and proliferation, such as chromosome segregation (GO:0007059) and DNA replication (GO:0006260). In addition, they were also enriched in some pathways highly related to PCOS, such as cellular carbohydrate metabolic process (GO:0044262), insulin signaling pathway and thyroid hormone signaling pathway. The promoter‐hypomethylation genes were also functionally analyzed shown in (Figure  4c,d ). They are mainly enriched in pathways closely related to PCOS, especially metabolic pathways, such as lipid catabolic process (GO:0016042) and triglyceride metabolic process (GO:0006650), as well as response to insulin (GO:0032868), and thyroid hormone signaling pathway. Enrichment analysis of DMGs' feature in human group. (a) and (b) GO analysis and KEGG analysis of promoter‐hypermethylation genes, respectively; (c) and (d) GO analysis and KEGG analysis; of promoter‐hypomethylation genes, respectively. (For GO analysis, X ‐axis: Gene counts in each GO term; Y ‐axis: GO terms; Blue bar: Biological process; Orange bars: Pathways related to cellular component; Green bar: Pathways related to molecular function; Orange box: Cell growth and proliferation; Blue box: Pathways related to metabolism. For KEGG analysis, X ‐axis: Gene ratio for each description; Y ‐axis: Description in KEGG; Circular diameter: Gene counts; Color temperature: −log 10 ( p ‐adjust).) The miRNA‐seq data of ovarian tissue samples from four PNA mice and three control mice were preprocessed. After obtaining the miRNA expression matrix from the samples, differential analysis was performed and 18 differentially expressed miRNAs (DemiRs) were obtained, of which nine were upregulated and nine were downregulated (Table  S5 ). Subsequently, hierarchical cluster analysis of DEmiRs was performed and the results were visualized using volcano plots and hierarchical clustering heat maps (Figure  5a,b ) The classification between the PCOS and control groups was good, and there were significant differences in miRNA expression patterns in ovarian tissues. Results of differential analysis of DEmiRs. (a) Volcano plot of DemiRs in the mouse group; (b) Hierarchical clustering heat plot of DemiRs in the mouse group; (c) Volcano plot of DemiRs in the human group; (d) Hierarchical clustering heat plot of DemiRs in the human group. (For volcano maps, X ‐axis: LogFC; Y ‐axis: −log 10 ( p ‐value); Red dot: Upregulated miRNAs; Blue dot: Downregulated miRNAs; Gray dot: MiRNAs with no significant expression differences. For heatmaps, X ‐axis: Sample type; Y ‐axis: DEmiRs; Color temperature: Normalized DEmiR expression level.) For human group, the miRNA‐seq data of ovarian GC samples from five PCOS patients and five control ovarian GC samples were preprocessed. Thirty‐eight DEmiRs were obtained, of which 23 were upregulated and 15 were downregulated (Table  S6 ). They will be used for subsequent multiomics joint analysis. Subsequently, a hierarchical cluster analysis of DEmiRs was performed, and the results were shown in Figure  5c,d , which indicated that the classification between the PCOS and control groups was as good as the mouse group. The high‐throughput sequencing RNA‐seq data of ovarian tissue samples from three PNA mice and two control mice were preprocessed, and the quality test results showed that the quality score for each base was above 30. Differential analysis revealed 3338 differentially expressed genes (DEGs), of which 1877 were upregulated and 1461 downregulated (Table  S7 ). Hierarchical cluster analysis of DEGs was performed, and the results were visualized in (Figure  6a,b ), which showed that the classification between the PNA group and the control group was good, and there were significant differences in the expression patterns of genes in ovarian tissues. Results of differential analysis of DEGs. (a) Volcano plot of DEGs in the mouse group; (b) Hierarchical clustering heat plot of DEGs in the mouse group; (c) Volcano plot of DEGs in the human group; (d) Hierarchical clustering heat plot of DEGs in the human group. (For volcano maps, X ‐axis: LogFC; Y ‐axis: −log 10 ( p ‐value); Red dot: Upregulated DNAs; Blue dot: Downregulated DNAs; Gray dot: DNAs with no significant expression differences. For heatmaps, X ‐axis: Sample type; Y ‐axis: DEGs; Color temperature: Normalized DEG expression level.) For human group, the high‐throughput sequencing RNA‐seq data of ovarian GC samples from three PCOS patients and three normal women were preprocessed. 1107 DEGs were obtained, of which 436 were upregulated and 671 were downregulated (Table  S8 ). The results of hierarchical cluster analysis of DEGs were visualized using a volcano plot and a hierarchical cluster heat map (Figure  6c,d ), showing that there were significant differences in the expression patterns of genes in ovarian tissues. In order to further study the regulatory effect of DEGs on PCOS, we conducted GO enrichment analysis and KEGG pathway analysis of DEGs to explore the biological processes and pathways related to PCOS. For the mouse group, we selected the top 10 p ‐value entries in each category from the GO enrichment analysis shown in Figure  7a . It can be seen that the function of DEGs is mainly enriched in two aspects: cell proliferation (e.g., nuclear division, meiosis, cell cycle) and immune regulation (such as T cell activation, lymphocyte differentiation, cell activation in immune response). There are also GO entries related to the proliferation of immune cells, such as regulation of monocyte proliferation (GO:0032944) and regulation of lymphocyte proliferation (GO:0050670). Figure  7b shows the top 20 p ‐value results for the KEGG pathway, indicating that DEGs are involved in pathways highly associated with PCOS, such as steroid biosynthesis and oocyte meiosis, as well as pathways related to inflammation and immune response, such as Th1 and Th2 cell differentiation, phagosome, cytokine‐cytokine receptor interaction. It is also involved in cell proliferation pathways such as cell cycle. Results of enrichment analysis of DEGs. (a) Results of GO enrichment analysis of DEGs in the mouse group; (b) Results of KEGG enrichment analysis of DEGs in the mouse group; (c) Results of GO enrichment analysis of DEGs in the human group. (d) Results of KEGG enrichment analysis of DEGs in the human group. (For GO analysis, X ‐axis: Gene counts in each GO term; Y ‐axis: GO terms; Blue bar: Biological process; Orange bar: Cellular component; Green bar: Molecular function; Orange box: Cell growth and proliferation; Blue box: Pathways related to metabolism; Orange box: Pathways related to immune modulation; Blue box: Pathways related to GCs. For KEGG analysis, X ‐axis: Gene ratio for each description; Y ‐axis: Description in KEGG; Circular diameter: Gene counts; Color temperature: −log 10 ( p ‐adjust); Orange box: Pathways related to inflammation and immune response.) Figure  7c shows the results of GO enrichment analysis for the human group. The function of DEGs is mainly enriched in immune regulation and GC‐related biological processes. There are many pathways related to immune modulation, such as T cell activation (GO:0042110) and leukocyte migration (GO:0050900). The biological processes related to GCs are also closely related to the occurrence and development of PCOS, mainly focusing on GC differentiation (GO:0030851), granulocyte migration (GO:0097530), and chemotaxis of GCs (GO:0071621). Figure  7d shows the results of KEGG pathway analysis, showing that DEGs are also involved in pathways related to inflammation and immune response, such as Th1 and Th2 cell differentiation, phagosome, and cytokine‐cytokine receptor interaction. To further investigate the regulatory mechanism of methylation in PCOS, we combined the results of RNA‐seq and MBD‐seq data analysis. First, we combined 3338 DEGs screened from the PNA mouse model group and 1107 DEGs screened from the clinical PCOS patient group. The R‐package homologene was used to find out the human homologous genes corresponding to the DEGs in the mouse group and then intersected with the DEGs in the human group. We finally obtained 236 DEGs, of which 152 were upregulated and 84 were downregulated. Hierarchical clustering analysis was performed, and the final results were visualized in Figure  8a,b . It can be seen that the PCOS group and the control group were well classified, and there were significant differences in the expression profile of DEGs. Mouse and human homologous DEG analysis. (a) Volcano plot ( X ‐axis: LogFC; Y ‐axis: −log 10 ( p ‐value); Red dot: Upregulated DNAs; Blue dot: Downregulated DNAs); (b) Hierarchical clustering heat map ( X ‐axis: Sample type; Y ‐axis: Homologous DEGs; Color temperature: Normalized homologous DEG expression level); (c) GO functional analysis ( X ‐axis: −log 10 ( p ‐value) in each GO description; Y ‐axis: GO descriptions; Orange box: Pathways related to immunity and inflammation); (d) KEGG functional analysis ( X ‐axis: −log 10 ( p ‐value) in each KEGG description; Y ‐axis: KEGG descriptions; Orange box: Pathways related to immunity and inflammation.) We performed GO enrichment analysis and KEGG pathway analysis to explore the role of DEGs in the occurrence and development of PCOS. Figure  8c,d showed the top 20 and top 10 entries with p ‐values in the GO analysis and KEGG analysis, respectively. The function of DEGs is mainly focused on immunity and inflammation, such as T cell activation (GO:0042110), leukocyte‐mediated immune regulation (GO:0002703), and immune effector regulation (GO:0002697). Subsequently, we combined 1720 DMGs from the PNA mouse model group with 120 581 DMGs from the PCOS patient group using similar methods. Finally, 1211 DMGs were obtained. We screened the corresponding DEGs and DMGs in the RNA‐seq group and MBD‐seq group, respectively, and then carried out transcriptome‐methylome multiomics analysis in order to further explore their interactions in PCOS regulation. After obtaining 236 DEGs and 1211 DMGs, the two groups of genes were intersected and screened according to the criteria of negative regulation between methylation level and expression level in the human group. 10 promoter‐hyper‐methylation & hypo‐expression genes and 16 promoter‐hypo‐methylation & hyper‐expression genes were obtained (Table  3 ). Finally, we screened genes with the same trend for differential expression in mouse and human groups. One gene, ANXA1 , with promoter‐hypomethylation and high expression was obtained. Most genes exhibit genetic or epigenetic differences between mice and humans. For instance, CD300A , KCNMA1 , and SORBS1 exhibited upregulated expression in both mouse and human PCOS phenotypes. However, their promoter methylation trends were opposite between the two species: hypermethylation was detected in PNA mice, while hypomethylation was detected in PCOS patients. Genes whose methylation levels are inversely correlated with gene expression levels. miRNAs negatively regulate the expression of target genes by inhibiting or degrading the translation of target mRNAs. To further investigate the regulatory mechanisms of miRNAs in PCOS, we performed a combined analysis of miRNA‐seq and RNA‐seq data. First, we aligned 18 DEmiRs screened from the PNA mouse model group to human homologous genes, then intersected with 38 DEmiRs screened from the PCOS patient group, and finally obtained 1 DEmiR: hsa‐miR‐106a‐5p. We predicted the target genes of has‐miR‐106a‐5p, intersected them with the final 236 DEGs, screened them according to the criteria of negative regulation between miRNA expression levels and gene expression levels, and finally obtained four differentially expressed target genes in PCOS patients: LDLR , RUNX1 , FICD , and STAC2 (Table  4 ). We screened genes with the same trend for differential expression in mouse and human groups and LDLR with low expression was obtained. Differential expression of hsa‐miR‐106a‐5p and its target genes. Abbreviations: DEG, differentially expressed gene; DEmiR, differentially expressed microRNA. The results of PCR validation are shown in Figure  9 . ANXA1 was significantly differentially expressed in PCOS patients (Figure  9a ), which was consistent with the above data analysis results with high reliability. For LDLR , there was also a significant differential expression (Figure  9b ), matching the results of the previous analysis. Therefore, both genes were included in the later discussion. We also conducted PCR validation in human PCOS patient samples for CD300A , KCNMA1 , and SORBS1 shown in Figure  9c–e . In PCOS patients, their expression is significantly upregulated. PCR validation results for human. (a) PCR results for ANXA1 ; (b) PCR results for LDLR ; (c) PCR results for CD300A ; (d) PCR results for KCAMN1 ; (e) PCR results for SORBS1 . ( X ‐axis: Sample type; Y ‐axis: Relative fold change; Circle dot: Normal samples; Square dot: PCOS samples; ** p ‐value < 0.01.) Combined with the p ‐value of the independent‐sample t‐test, the FC value of the fold analysis method and the VIP value of the OPLS‐DA model (Figures  [Link] , [Link] ), p ‐value  2 or  1 were used to identify candidate differential lipids. 34 significantly changed follicular fluid lipid metabolites (as shown in Table  5 ) were jointly screened and used as potential biomarkers for PCOS for subsequent discussion. In follicular fluid positive ion samples, 16 lipids, including triglyceride (TG) (16:0/17:0/18:1), TG (18:0/18:0/18:1), and TG (16:0/14:0/16:1), were significantly increased in the follicular fluid of PCOS patients. 10 lipids, including TG (8:0/8:0/8:0), TG (8:0/10:0/10:0), and TG (10:0/12:0/14:0), were significantly decreased. In follicular fluid negative ion samples, only phosphatidylinositol (PI) (16:0/16:1) was significantly upregulated in the PCOS group, while the contents of seven lipids, including phosphatidylcholine (PC) (16:0/22:5) and phosphatidylethanolamine (PE) (16:1e/22:6), were significantly downregulated. Follicular fluid metabolomics differential analysis results. Abbreviations: ESI, electrospray ionization; PC, phosphatidylcholine; PE, phosphatidylethanolamine; PI, phosphatidylinositol; TG, triglyceride.

In this study, the enrichment analysis of GO and KEGG was plotted in the data analysis of each omics. It can be seen that DMGs are mainly enriched in cell growth and proliferation and lipid metabolism pathways, while DEGs are mainly enriched in immune regulation‐related pathways. After that, we took the homologous genes of mice and humans. Then we looked for DMGs and DEGs common in both species and took the intersection of them. The number of genes obtained after screening was drastically reduced and mainly enriched in metabolism‐ and inflammation‐related pathways. Similar results were obtained for the analysis of miRNA expression, with a sharp decrease in the number after taking the intersection. It is evident that some genes play important roles in PCOS across species as well as across omics and are ideal for study. DNA methylation is closely related to PCOS. The changes in genomic methylation status can affect gene expression, which in turn plays a role in the occurrence and development of polycystic ovary syndrome. Approximately 50 genes with abnormal methylation status are associated with PCOS [ 26 ]. These genes are associated with symptoms of PCOS, such as insulin resistance [ 27 ], ovulation regulation, [ 28 ] hyperandrogenism [ 29 ], and a range of concurrent cancers [ 30 ]. ANXA1 plays an important role in the negative regulation of innate immunity. It can inhibit the migration and recruitment of neutrophils at the site of inflammation, thereby reducing their pro‐inflammatory activity. Otherwise, severe inflammatory lesions may occur in the tissue [ 31 ]. Glucocorticoids are also involved in this process, promoting the anti‐inflammatory effects of ANXA1 [ 32 ]. ANXA1 can also promote neutrophil apoptosis and macrophage cellular behavior to clear necrotic neutrophils. These cells or cell fragments activate the body's adaptive immune system, and ANXA1 can shorten their exposure time, limiting their stimulation of the immune system [ 33 ]. Meanwhile, ANXA1 is able to inhibit pro‐inflammatory Toll‐like receptors in dendritic cells, giving them a degree of tolerance [ 34 ]. The effect of ANXA1 on T cells is still inconclusive. Under normal conditions, the expression level of ANXA1 gene in T cells is low, but it is upregulated in inflammatory environments [ 35 ]. Studies have shown that when T cells are stimulated, the expression of ANXA1 receptors is also upregulated [ 36 ]. Inflammatory biomarkers play a crucial role in regulating ovarian function. Any disturbances related to inflammation may lead to ovarian dysfunction, thereby contributing to the development of PCOS [ 37 ]. Therefore, the upregulation of ANXA1 , which possesses anti‐inflammatory properties, could be one of the body's adaptive responses to the onset of PCOS. Aberrant hypermethylation of ANXA1 is associated with its aberrant expression and a variety of cancers. Compared with the control group, the methylation level of the ANXA1 gene in nasopharyngeal carcinoma tissues was extremely high, reaching 92%. This may be a direct cause of the decline in ANXA1 expression in nasopharyngeal carcinoma, accelerating lymph node metastasis in nasopharyngeal carcinoma. The expression level of ANXA1 in B‐cell lymphoma, prostate cancer, and esophageal cancer was also significantly decreased [ 38 ]. Current research focuses on the effects of hypermethylation and low expression of ANXA1 on cancer occurrence and progression. In this study, we found that ANXA1 showed hypomethylation in the promoter region and high expression in PCOS patients. ANXA1 antagonizes PPARγ signaling in the peroxisome proliferator‐activated receptor (PPAR) pathway. PPARγ inhibition in adipocytes can reduce insulin sensitivity [ 39 ] which may be linked to the symptoms of IR in PCOS. However, it is worthy to note that none of the samples of PCOS patients in this study met the criteria for obesity (BMI > 30). Considering the role of the gene ANXA1 in metabolism and inflammation, future studies need to include samples of obese PCOS patients. This may make the role of ANXA1 more visible. On the other hand, the selection of non‐obese women enables this study to avoid the possible interference caused by obesity factors, which has certain reference value. Meanwhile, studies have shown that ANXA1 is associated with intrauterine growth restriction (IUGR) in rat models. IUGR can increase the risk of developing type 2 diabetes, which aligns with the frequent occurrence of insulin resistance symptoms in PCOS patients [ 40 ]. However, many homologous genes do not exhibit entirely consistent patterns of gene expression and methylation changes, such as CD300A , KCNMA1 , and SORBS1 . These genes show upregulated expression in both mouse and human PCOS phenotypes. Yet, their promoter methylation trends are opposite between the two species: hypermethylation is observed in PNA mice, while hypomethylation is detected in PCOS patients. We believe this reflects species‐specific differences between mice and humans. Additionally, it suggests that the expression of these genes may be regulated by more complex mechanisms, where promoter methylation might serve as merely one contributing factor—or potentially play no role at all, which requires further research to elucidate. Altered miRNA expression has a significant impact on the development of various diseases. There has been sufficient evidence to support the combination between miRNA dysregulation in GCs and PCOS, which can be used to identify potential therapeutic targets [ 41 ]. Dysregulation of miRNAs in follicular fluid may be one of the key factors in ovulatory dysfunction in women with PCOS [ 42 ]. Also, miRNAs affect many biological processes including glucose and lipid metabolism; thus, they play an important role in insulin signaling pathway and IR in PCOS [ 43 ]. In this study, we found that one miRNA, hsa‐miR‐106a‐5p, had a higher expression in PCOS patients. Some studies have confirmed that hsa‐miR‐106a‐5p is significantly associated with renal cell carcinoma (RCC) progression and may serve as a candidate biomarker and a novel therapeutic target to inhibit RCC metastasis [ 44 ]. It has also been found that hsa‐miR‐106a‐5p is highly expressed in lung adenocarcinoma (LUAD), regulates the phosphorylation of the AMPK pathway, and is involved in the proliferation, migration, and autophagy of LUAD cells [ 45 ]. This may indicate that miR‐106a‐5p is associated with multiple cancer complications in PCOS. This present study found significantly reduced LDLR expression in patients with PCOS, which is consistent with existing studies. The specific reduced expression of adipose LDLR in PCOS may be associated with hyperandrogenism, which may be related to the tendency of PCOS patients with concomitant liver disease to progress to non‐alcoholic steatohepatitis (NASH) [ 46 ]. Meanwhile, the protein convertase PCSK9 promotes LDLR degradation in hepatocytes, and its expression is elevated in patients with PCOS. The high expression of PCSK9 may have an inhibitory effect on both mRNA and a protein expression level of LDLR [ 47 ]. LDLR deficiency in cells leads to reduced estrogen biosynthesis and secretion, affecting follicular cell proliferation and its steroid hormone synthesis [ 48 ]. LDLR is a key enzyme that regulates lipid metabolism, and abnormalities in metabolic status may negatively impact the female reproductive system and contribute to obesity (manifested as elevated BMI). Therefore, the downregulation of LDLR expression may be associated with the development and progression of PCOS. Studies have shown that the phosphoinositide‐3‐kinase–AKT (PI3K‐AKT) pathway is linked to common symptoms of PCOS, including insulin resistance (IR), obesity, and polycystic follicles [ 49 ]. The upregulation of PI may be associated with abnormalities in the PI3K‐AKT pathway in PCOS. During the human metabolic process, insulin binding to its receptor initiates the PI3K‐AKT pathway to regulate glucose metabolism. Specifically, activated PI3K phosphorylates phosphatidylinositol‐4,5‐bisphosphate (PIP2) to generate phosphatidylinositol‐3,4,5‐triphosphate (PIP3). PIP3 indirectly activates AKT, which mediates most key metabolic effects of insulin. As PIP2 is derived from phosphorylation of PI, elevated PI levels may suggest impaired PIP2 production, thus causing reduced PIP3 generation and then failure of AKT activation, finally leading to disrupted insulin‐mediated metabolic functions, aligning with the IR symptom commonly observed in PCOS [ 50 ]. The expression of PE in the follicular fluid of PCOS patients decreased. Studies have shown that PE levels are negatively correlated with the free androgen index (FAI). FAI is closely related to hyperandrogenemia, a classic symptom of PCOS. Additionally, FAI can serve as a potential marker for insulin resistance, further supporting the association between PE and PCOS [ 51 ]. In this study, we jointly analyzed the genome‐wide RNA‐seq, MBD‐seq and miRNA‐seq sequencing data of ovarian tissues from PNA mouse models and GCs from PCOS clinical patients. We explored the regulatory mechanism of DNA methylation and miRNA in PCOS, the potential miRNAs and genes as diagnostic markers for PCOS, and their biological functions. Two genes ( ANXA1 and LDLR ) were screened out, which play an important role in the PCOS‐related pathways. We screened out potential biomarkers for PCOS diagnosis and elucidated the regulatory pathways involved in them, which also had suggestive significance for the therapeutic targets of PCOS. Combined with the results of multiomics detection in clinical patients, the mechanism of PCOS occurrence and development can be further elaborated, so as to provide a more solid theoretical basis for the development of diagnostic kits and targeted treatment regimens for PCOS. However, the small number of samples used in this study may result in less reliable results. In future studies, we can further increase the number of samples and select PCOS patient samples that better reflect the heterogeneity of PCOS (including but not limited to obesity). In practical application, doctors can collect micro blood samples from subjects for high‐throughput sequencing. The sequencing results will be analyzed for differential expression. If ANXA1 has a significantly high expression, and the expression of LDLR is found low, the subject can be considered to be at a higher risk of developing PCOS, achieving the purpose of early diagnosis of PCOS.

Participants (17 PCOS patients and 17 normal subjects) were recruited in Yuncheng Central Hospital of Shaanxi Province. All of them agreed to participate in the study and provided written informed consent.

Polycystic ovary syndrome (PCOS) is a common hormonal imbalance syndrome in women of reproductive age [ 1 ]. It is estimated that 8%–13% of women of childbearing age are affected by PCOS, but up to 70% of them remain undiagnosed [ 2 ]. Polycystic ovaries, anovulation, and hyperandrogenemia are the three main features of PCOS [ 2 ]. Studies have shown that patients with PCOS have impaired endometrial function and some degree of weakened oocyte capacity compared to healthy women [ 3 , 4 ]. This is closely related to infertility (regardless of ovulatory status) and pregnancy complications in PCOS patients. In addition, PCOS can lead to insulin resistance (IR), menstrual disorders, hirsutism, obesity, acne, and oily skin. PCOS also comes with a host of complications such as endometrial cancer, type 2 diabetes, and cardiovascular disease [ 5 , 6 ]. Due to the extreme heterogeneity of PCOS, there is currently no definitive diagnostic method. The Rotterdam criteria recommended by the Endocrine Society are generally used for diagnosis, which require at least two of the following three manifestations to appear: hyperandrogenism, ovulatory dysfunction, and polycystic ovaries [ 7 ]. PCOS cannot be cured, but its symptoms can be relieved with treatment [ 2 ]. The clinical treatment of PCOS is mainly drug therapy such as using clomiphene, letrozole, metformin, statins, or gonadotropins [ 8 ]. Surgeries can be used for patients who have resistance or severe side effects of drug treatments. For instance, a minimally invasive surgery, transvaginal hydrolaparoscopy (THL), provides an accurate image of the female pelvis. Thus, surgeries can be performed in the patient's ovarian capsule drilling [ 9 ]. The exact cause of PCOS has not been identified, but women with a family history of type 2 diabetes or obesity are at higher risk [ 10 , 11 ]. Recent research advances in genetics and epigenetics have improved the understanding of the etiology and pathophysiology of PCOS, and certain loci associated with PCOS have stimulated the search for new therapeutic strategies [ 12 ]. Several studies have confirmed changes in the levels of histone acetylation, non‐coding RNA (ncRNA), and DNA methylation in PCOS [ 13 ]. RNA‐mediated DNA methylation (especially in promoter regions) leads to gene silencing at the transcriptional level. In active genes, promoters are usually hypomethylated [ 14 ]. A methylation histology study showed that PCOS patients have significant methylation level changes in several genes involved in a variety of physiological functions [ 15 ]. Abnormal histone modifications are associated with insulin‐related pathway dysfunction and may contribute to the development of PCOS [ 16 ]. Large amounts of aberrant long‐stranded ncRNAs were detected in various rodent models of PCOS and samples from different tissues of women with PCOS (serum, granulosa cells, follicular fluid, and ovarian cumulus cells) [ 17 ]. Also, PCOS is strongly associated with metabolomics. Abnormalities in pathways associated with steroidogenesis are of importance in the development of PCOS [ 18 ]. PCOS is usually characterized by increased levels of luteinizing hormone (LH) and gonadotropin‐releasing hormone (GnRH), resulting in ovulatory dysfunction and hyperandrogenemia [ 19 ]. Elevated androgen levels also impair the metabolic function of the pancreas, leading to IR symptoms [ 20 ]. Although the study of various omics has made great progress in the pathological process of PCOS, the combined multiomics analysis is still lacking. In this study, we conducted methyl‐CpG binding domain‐based capture and sequencing (MBD‐seq) and miRNA sequencing (miRNA‐seq) data from ovarian tissues of PNA mouse models and ovarian granulosa cells (GCs) of patients with PCOS to screen differentially expressed genes and differentially promoter‐methylated genes in PCOS. These data were jointly analyzed with metabolomics analysis and transcriptome sequencing (RNA‐seq) data to find differentially expressed miRNAs and differentially methylated genes in PCOS. Finally, through functional analysis, we aimed to identify and screen the biomarkers related to the diagnosis of PCOS and studied their biological functions, which are widely involved in metabolic functions, hormone synthesis, and immune‐inflammation in PCOS.

The authors declare no conflicts of interest.

Female mice were randomly mated with males. The discovery of the copulatory plug was identified as the first day of pregnancy in the female mice. After 16–18 days of gestation, the mice were randomly divided into two groups: The pregnant group was subcutaneously injected with 70 μL containing 350 μg dihydrotestosterone (DHT) and the control group with isopycnic sesame oil every day. The offspring of the DHT‐induced group were PNA mice, while the offspring of the control group were control mice. The mice were anesthetized and euthanized after 8 weeks, and the ovarian tissues were collected. All of the participants (17 PCOS patients and 17 normal subjects) were recruited in Yuncheng Central Hospital of Shaanxi Province. Participants, between the ages of 20 and 35 from the same geographical region, agreed to participate in the study and provided written informed consent. The diagnostic criteria for PCOS patients are at least two Rotterdam ESHRE/ASRM (2003) criteria [ 21 ]. The control group were normal women of childbearing age and had normal menstrual cycles, no clinical symptoms related to PCOS, and no endometriosis or other chronic diseases, etc. Clinical data were collected about 4 days after menstruation, including the number of follicles, preoperative serum hormone levels, and BMI. Participants' fasting whole blood was collected for 10 h. The Beckman UniCel DxI800 fully automatic chemiluminescence analyzer (Beckman Coulter, CA, USA) measures prolactin (PRL), LH, testosterone (T), thyroid stimulating hormone (TSH), and other hormone levels. Participants were injected subcutaneously with 5000 IU of gonadotropin‐releasing hormone for 15 days, 7 days after ovulation, and follicle size was measured with ultrasound. When more than three follicles were detected with a diameter of 16 mm, participants were then injected with 6000 IU of human chorionic gonadotropin. There was no ovarian stimulating injection such as recombinant FSH or human menopausal gonadotropin. After 1.5 days, the oocytes were collected under ultrasound guidance. Then, oocytes were removed, and follicular fluid and GCs were collected by centrifugation. The collected GCs were washed twice with Dulbecco's modified Eagle medium. The follicular fluid lipid and GCs were quickly transferred to −80°C for preservation. Mouse ovarian tissue samples and human GC samples were extracted by TRIzol (Invitrogen, Waltham, MA, USA) as recommended by the kit manufacturer's protocol, and the quality of RNA was assessed by 2% agarose denatured gel electrophoresis. Genomic DNA was extracted from mouse ovarian tissue samples and human GC samples using the QIAamp DNA Mini‐Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions, which use NanoDrop One (Thermo Scientific, Wilmington, USA) to identify the quality of genomic DNA. The sequencing data for MBD‐seq were downloaded from our group's previous studies ( https://www.ncbi.nlm.nih.gov/geo/ ; Accession Number: GSE156961 & GSE138575 ) [ 22 , 23 ]. After MBD‐seq sequencing, trim_galore (version 0.4.3) was used to trim and filter the raw data to obtain clean reads. FastQC (version 0.11.5) was then used for quality inspection. Clean reads were then aligned to the mouse reference genome mm10 or the human reference genome hg38 using bowtie2 (version 2.1.0). Then, we use samtools(version 1.5) to filter and sort the data with MAPQ30 as the standard, and finally use picard (version 2.25.4) to deduplicate and get the . bam file. Differential analysis of methylated regions across the genome was performed using the edgeR (version 3.38.4) package with thresholds of |log 2 FoldChange| > 1 and p  < 0.05. Subsequently, the distribution of genome‐wide differential methylation signals on chromosomes was visualized, and hierarchical cluster analysis of DMRs was performed to explore the differences in DNA methylation patterns between the two groups. Finally, the distribution of DMRs in genomic characteristic regions was analyzed by gene annotation, and DMGs were screened. GO enrichment analysis and KEGG pathway analysis were performed on the screened DMGs to explore the related biological pathways and functions of them in the occurrence and development of PCOS. The miRNA sequencing data used in this study for mouse ovarian and human GCs were downloaded from our previous studies ( https://www.ncbi.nlm.nih.gov/geo/ ; Accession Number: GSE168166 & GSE138575 ) [ 23 , 24 ]. Cutadapt (version 1.18) was used to trim and filter the raw miRNA‐seq data to obtain clean reads. This was followed by FastQC (version 0.11.5). The clean reads were then aligned to mouse or human mature miRNAs in miRBase (Release 22.1) using bowtie2 (version 2.1.0). Finally, samtools (version 1.5) was used to count and obtain the expression matrix of the sample miRNAs. Subsequently, the miRNA differential analysis was performed by edgeR (version 3.38.4) and DESeq2 (version 1.36.0) methods, and the threshold was set as |log 2 FoldChange| > 1 and p  < 0.05, and the results of the two methods were intersected to obtain DemiRs. Then, the final DEmiRs were visualized by using the volcano map. Finally, hierarchical clustering analysis was performed on DEmiRs to explore the differences in miRNA expression patterns between the two groups, and a hierarchical clustering heat map was drawn. The potential target genes corresponding to the DEmiRs were predicted using the miRTarBase (Release 9.0) database. Subsequently, GO enrichment analysis and KEGG pathway analysis were performed to explore the relevant biological pathways and functions of the target genes in the occurrence and development of PCOS. The transcriptome sequencing data for mouse ovarian tissue and human GCs were derived from our previous studies ( https://www.ncbi.nlm.nih.gov/geo/ ; Accession Number: GSE156961 & GSE138575 ) [ 22 , 23 ]. After RNA‐seq sequencing, the raw offline data were trimmed and filtered using trim_galore (version 0.4.3) to obtain clean reads. FastQC (version 0.11.5) was then used for quality inspection. The clean reads were then aligned to the mouse reference genome mm10 or the human reference genome hg38 using HISAT2 (version 2.2.1). Finally, featureCounts (version 2.0.3) were used to count each gene according to the Ensembl gene‐level annotation to obtain the mRNA expression matrix. This is followed by subsequent data processing and analysis in R. The difference analysis was carried out by using edgeR and DESeq2 packages, and the threshold was set as |log 2 FoldChange| > 1 and p  < 0.05. The analysis results of the two packages were intersected to obtain DEGs, and the final DEGs were visualized by using the volcano map. Finally, hierarchical clustering analysis was performed on DEGs to explore the differences in gene expression patterns between the two groups, and a heat map of sub‐clustering was drawn for visualization. The function of genes was enriched by the GO database from function and biological pathways involved. Pathway enrichment analysis was performed using the KEGG to explore the biological processes and pathways related to PCOS. Genomic DNA extracted from clinical samples was converted by bisulfite using the EZ DNA Methylation‐Gold Kit. In the presence of bisulfite, methylated cytosine does not undergo conversion, while unmethylated cytosine is converted to thymine. Primers designed for methylated and unmethylated promoter DNA sequences were used for PCR amplification. The primers used in MSP are shown in Table  S1 . The MSP result is shown in Figure  S1 . The gene expression level of ANXA1 was detected by RT‐qPCR. 500 ng total RNA was reverse‐transcribed into cDNA using PrimeScript RT Master Mix (RR037; TAKARA, Japan). The obtained cDNA was quantitatively detected by The Luna Universal qPCR Mix (M3003S; NEB, USA) in QuantStudio 3 (Applied Biosystems) according to the manufacturer's protocol. The primers used in RT‐qPCR are shown in Table  S2 . Follicular fluid stored was restored at −80°C to room temperature and transferred 50 μL of follicular fluid into a 1.5 mL centrifuge tube. 200 μL of methanol precooled at 4°C, 400 μL of chloroform, and 170 μL of distilled waterwere thoroughly added and vortexed. Next, we centrifuged at a low temperature of 12 000 rpm at 4°C for 10 min to collect the lower clear liquid. Then, the lower clear liquid was transferred to a vacuum drying concentrator for low‐temperature freeze‐drying. The freeze‐dried sample was dissolved in 100 μL L (methylene chloride:isopropyl alcohol:methanol = 2:1:1 ([vol:vol:vol])), vortex thoroughly and centrifuged at 12 000 rpm for 10 min. 50 μL of supernatant was taken for UPLC‐MS mass spectrometry. The mass spectrometry data were compared with the Human Metabolome Database. Orthogonal partial least squares discriminant analysis (OPLS‐DA) modeling analysis was performed to determine the difference between the samples. The variable importance in projection (VIP) of lipid metabolites in follicular fluid was greater than 1, |log 2 (FC)| ≥ 1, and p ‐value < 0.05 were identified as differential metabolites. The software used for statistical analysis of data in this study is Statistical Package for Social Sciences (SPSS) (version 24.0). The difference between the two groups was determined by the two‐sided Student's t ‐test; p ‐value < 0.05 was considered a statistically significant difference.

Figure S1. The result of MSP. Figure S2. Analysis of follicular fluid positive ion PCA results. Figure S3. Analysis of follicular fluid negative ion PCA results. Figure S4. Results of multivariate statistical analysis of metabolomics of follicular fluid positive ions. Figure S5. Results of multivariate statistical analysis of metabolomics of follicular fluid negative ions. Table S1. Primers used for MSP. Table S2. Primers used for PCR. Table S3. Hypermethylated DMRs in mouse group. Table S4. Hypomethylated DMRs in mouse group. Table S5. DemiRs in mouse group. Table S6. DemiRs in human group. Table S7. DEGs in mouse group. Table S8. DEGs in human group.