Skip to main content

N6-methyladenosine regulator-mediated methylation modification patterns and immune infiltration characterization in Polycystic Ovary Syndrome (PCOS)

Abstract

Background

Polycystic ovary syndrome (PCOS) is a multisystem-related disease whose pathophysiology is still unclear. Several regulators of N6-methyladenosine (m6A) modification were confirmed to play a regulatory role in PCOS. Nonetheless, the roles of m6A regulators in PCOS are not fully demonstrated.

Materials and methods

Four mRNA expression profiling microarrays were obtained from the Gene Expression Omnibus (GEO) database. Differentially expressed m6A regulators between PCOS and normal patients were identified by R software. A random forest modal and nomogram were developed to assess the relationship between m6A regulators and the occurrence risk of PCOS. A consensus clustering method was utilized to distinctly divide PCOS patients into two m6A subtypes (m6A cluster A/B). The patterns of differential expression and immune infiltration were explored between the two m6A clusters.

Results

In this study, 22 significant m6A regulators were identified between healthy controls and PCOS patients. The random forest model determined three optimal m6A regulators which are related to the occurrence risk of PCOS, including YTHDF1, RBM15 and METTL14. A nomogram was established based on these genes, and its predictive reliability was validated by decision curve analysis. The consensus clustering algorithm distinctly divided PCOS cases into two m6A subtypes. The ssGSEA algorithm found that the immune infiltration was markedly enriched in m6A cluster B than in cluster A. The m6A-pattern related differentially expressed genes (DEGs) of the two m6A subtypes were demonstrated by differential expression analysis. We found that they were enriched in immune-related genes and various infection pathways. Based on the m6A-pattern related DEGs, the PCOS patients were classified into two m6A-pattern related genomic subtypes (gene clusters A and B).

Conclusions

The present study provided evidence concerning the different modification patterns of m6A regulators in PCOS compared with normal patients. This study will help clarify the overall impact of m6A modification patterns and related immune infiltration on PCOS.

Introduction

Polycystic ovary syndrome (PCOS) is a multisystem-related disease of the reproductive system characterized by the pathological accumulation of follicular immaturity and atresia, ovarian and interstitial dysplasia, hyperandrogenemia (HA), hyperinsulinemia, insulin resistance (IR), metabolic abnormalities, endocrine disorders and polycystic ovary [1, 2]. Previous studies have shown significant advances in factors likely involved in PCOS pathogenesis, including the potential roles of androgen, insulin, AMH, TGFβ, oxidative stress, pro-inflammatory cytokines, advanced glycation end products (AGEs), and so on [3, 4], among which the HA seemed to be the most vital determinant of the PCOS pathophysiology and the associated metabolic dysfunction [5]. Due to the high heterogeneity and complexity of this multifactorial disease, the molecular basis of the etiology of PCOS is not fully understood yet. Therefore, a deeper investigation of the molecular mechanisms underlying the pathogenesis of PCOS is urgently needed to facilitate the development of therapeutic modalities and improve the prognosis of patients.

N6-methyladenosine (m6A), one of the most prominent and abundant epigenetic modification forms of mRNA and lncRNA in eukaryotic cells, is a dynamic, reversible, and highly conservative process under the regulation of methyltransferases, demethylases, and binding proteins [6,7,8]. These regulating proteins of m6A modification are also called "writers," "erasers," and "readers" of m6A [9,10,11]. Increasing pieces of evidence have shown that m6A modification plays crucial roles in physiological processes and disease progressions [12], such as embryonic development [13, 14], stress response, immunity [15], cell proliferation and self-renewal ability [16], tumorigenesis and metastasis [15, 17], and so on.

Previous studies have investigated and revealed the significant role of m6A RNA methylation regulators in the occurrence and progression of PCOS [18]. However, these researches concentrated predominantly on several of the 26 m6A regulators. Therefore, the function of m6A regulators in PCOS remains to be further investigated.

Since the m6A modification might play a part in the dysregulation of various systems, which provides an orientation for investigation in PCOS. Concurrently, the development of bioinformatics analysis and public databases, such as the Gene Expression Omnibus (GEO) [19], provides means to understand the molecular pattern of m6A regulators in PCOS. In this research, we comprehensively assessed the effects of m6A regulators on the occurrence risk and subtype categorization in PCOS. We constructed a nomogram to visually display the scoring coefficient and correlation between the occurrence risk of PCOS and the three optimal m6A regulators (YTHDF1, RBM15 and METTL14). As a result, PCOS patients could receive clinical benefits from this nomogram. Additionally, we identified two m6A subtypes closely associated with different immune infiltrations. This study uncovered two different patterns of m6A regulators in the granulosa cells of PCOS, which provides a theoretical foundation for further research into m6A modification and new therapeutic strategies for PCOS.

Materials and methods

Data extraction and processing

The RNA sequencing (RNA-seq) data of the GSE34526, GSE80432, GSE102293, and GSE137684 datasets was downloaded by using Perl (www.perl.org/) from the GEO database (http://www.ncbi.nlm.nih.gov/geo/) after searching for all expression microarrays that matched PCOS keywords with unnecessary information removed. The expression profiles of 19 controls and 25 PCOS patients from 4 datasets were combined, and the batch normalization of the read count values was completed to preserve data consistency by the "limma" and "sva" R packages.

Differential expression analysis of m6A RNA methylation regulators

A total of 26 m6A RNA methylation regulators were collected from previous studies, including 9 writers, 15 readers, and 2 erasers (Table 1). Differentially expressed analysis of these regulators based on the limma (Linear Models for Microarray Data) algorithm was performed between control and PCOS patients to identify the DEGs of the m6A regulators. A protein–protein interaction (PPI) analysis of the above DEGs was accomplished by the STRING website (https://string-db.org/), and the hub genes were extracted by Cytoscape's cytoHubba plug-in. The gene set variation analysis (GSVA) enrichment analysis was conducted to assess the variations in biological processes to match the biological function between control and PCOS patients.

Table 1 Three types of m6A modification regulators and their primary biological functions

Models comparison, selection, and establishment

Random forest (RF) and support vector machine (SVM) models were established respectively as training models to predict the prevalence of PCOS. Reverse cumulative distribution of residuals, boxplots of residuals, and receiver operating characteristic (ROC) curves were conducted to assess the accuracy of models. The m6A regulators with an importance value over 2 were considered as PCOS-specific genes.

Establishment of the nomogram

Based on the m6A regulators with an importance value over 2, a nomogram was built with the "rms" package of R to predict the occurrence risk of PCOS, as well as to demonstrate the correlation between the occurrence risk and selected m6A regulators. The calibration curve assessed the predictive accuracy and reliability of this nomogram. The decision curve analysis (DCA) was conducted to evaluate the benefits of decisions from the nomogram. A clinical impact curve was established to assess the rationality and benefit of decisions.

Identification of subtypes based on m6A regulators

Consensus clustering with K-means algorithms was applied to identify m6A regulator-related subtypes correlated with gene expression. The quantity and robustness of clusters were determined with a consensus clustering algorithm realized in the "ConsensuClusterPlus" R package. The principal component analysis (PCA) was used to confirm the m6A subtypes.

Identification of DEGs between m6A subtypes

The "limma" R package was applied to identify DEGs among different m6A subtypes, named m6A-related DEGs, with the criterion of adjusted p < 0.01 and |logFC|≥ 1.5. GO and KEGG enrichment analysis were utilized to investigate the potential function of the DEGs between m6A-subtypes responsible for PCOS via the "clusterProfiler" package in R. Based on the expression of m6A-related DEGs, the m6A-related gene subtypes were identified.

Establishment of the m6A gene signature

To quantify the m6A regulator patterns of PCOS, a scoring system—the m6A gene signature was constructed to evaluate the m6A regulator patterns of PCOS patients, which we termed as m6Ascore. The PCA was conducted to obtain m6Ascore of each PCOS patients. Both principal components 1 and 2 were selected to act as signature scores, and m6Ascore was acquired based on the following method: m6Ascore = ∑(PC1i + PC2i), of which PC1/2 refers to principal component 1/2 and i to the expression of m6A-related DEGs.

Estimation of immune cells infiltration

The single sample gene set enrichment analysis (ssGSEA) was applied to assess the infiltration of immune cells among different subtypes. Immune cells include immune-enhancing cells and immune-suppressive cells. The gene expression levels of the patients were sequenced with ssGSEA to acquire an infiltration grade. Then the expression data were summarized for immunological analysis. Differential expression analysis with moderated t-tests would be utilized to assess the enrichment of immune cells and the distribution of pro-inflammatory cytokines between different subgroups, and p < 0.05 was recognized as a significant result.

Statistics analysis

Linear regression analyses were applied to determine the relationship between m6A regulators. Kruskal–Wallis tests were utilized to identify a discrepancy between subtypes. All statistical analyses were carried out with two-tailed tests, and the significant value was considered p < 0.05. R × 64 (version 4.1.2) and packages in R were utilized to run and process data scripts, output results, and plot diagrams.

Results

The landscape of the m6A regulators in PCOS

Based on the selected datasets, the "limma" R package was applied to analyze differential expressions of 26 m6A regulators between control and PCOS patients. The complete expression data of 24 m6A regulators was successfully retrieved. The expression landscape of these m6A regulators was visualized in a histogram and a heatmap (Fig. 1A, B). Subsequently, 22 DEGs were observed in the differently expressed analysis of m6A regulators between control and PCOS samples. Most of the DEGs were overexpressed in PCOS compared to the control group, including METTL14, WTAP, VIRMA, ZC3H13, RBM15, CBLL1, YTHDC1, YTHDC2, YTHDF1, YTHDF3, HNRNPC, FMR1, IGFBP1, ELAVL1, and FTO, and several other DEGs such as METTL3, RBM15B, HNRNPA2B1, IGFBP2, IGFBP3, IGF2BP1, and ALKBH5 were downregulated in PCOS. The GSVA analysis demonstrated that the glyoxylate and dicarboxylate metabolism, glycosaminoglycan degradation, phenylalanine metabolism, fructose and mannose metabolism, and lysosome were obviously enriched in the control group than in PCOS patients (Fig. 1C). The PPI analysis demonstrated the interactions of 22 DEGs of m6A regulators, and Cytoscape illustrated three hub genes, METTL3, YTHDF1, and YTHDF3 (Fig. 1D, E). Additionally, the 22 DEGs of m6A regulators were mapped onto chromosomes by the "RCircos" package in R (Fig. 1F).

Fig. 1
figure 1

Landscape of m6A RNA methylation regulators in PCOS. A Differential expression analysis of m6A regulators between control and PCOS patients named “con” and “PCOS” respectively, and the gene expression has been converted by log2(counts + 1). B Expression heatmap of the identified 24 m6A regulators in two groups. C GSVA enrichment analysis between PCOS and control group. D PPI network analysis of the DEGs in m6A regulators. E Top 10 hub genes of m6A regulators in PPI network. F Chromosomal positions of the 22 DEGs of m6A regulators. *p < 0.05, **p < 0.01, and ***p < 0.001

Correlation analysis between m6A regulators in PCOS

The correlations between the m6A regulators are presented in Fig. 2A. The m6A regulators could display positive correlations, such as FMR1 and YTHDF3 (coefficient = 0.97), and negative correlations, such as YTHDC1 and IGF2BP1 (coefficient = -0.85). To investigate the possibility of co-expression between m6A regulators, we discovered a clear relationship between YTHDF3 and other regulators, with the greatest relevance for it with FMR1 as mentioned earlier, which was consistent with the PPI analysis. The correlation heatmap shows that the readers and erasers of m6A regulators in PCOS were highly related. Linear regression analyses were applied to further investigate the relationship between readers and erasers. Significant negative correlations could be observed between ELAVL1, FMR1, HNRNPC, IGFBP1, YTHDC2, YTHDF2, YTHDF3, and ALKBH5 in PCOS patients. Figure 2B-H illustrate that PCOS patients with low ALKBH5 expression levels tend to present high levels of above readers. And the expression of FTO demonstrated a negative association with IGFBP3 in PCOS (Fig. 2I).

Fig. 2
figure 2

Correlation analysis of m6A regulators in PCOS. A Correlation heatmap of 24 m6A regulators. B-I Correlation between m6A readers and erasers in PCOS. Reader genes: ELAVL1, FMR1, HNRNPC, IGFBP1, YTHDC2, YTHDF2, YTHDF3, IGFBP3; eraser genes: ALKBH5 and FTO

Evaluation of the RF model and SVM model

Based on 22 DEGs of m6A regulators, the RF model and SVM model were respectively constructed to identify optimal m6A regulators in order to evaluate the occurrence risk of PCOS. According to the "reverse cumulative distribution of residuals" and "boxplots of residuals" (Fig. 3A, B), the RF model retained smaller residuals, indicating that the RF model's predictive performance has more superiority to the SVM model. Moreover, the ROC curves were established to assess the accuracy of two models, and the AUC value also demonstrated the better performance of the RF model (Fig. 3C). Consequently, the RF model was applied to screen for the signature genes of PCOS. The RF model displayed the relationship between the error and the number of decision trees to find the point with the smallest cross-validation error for choosing the optimal number of decision trees (Fig. 3D). Based on the RF model, the DEGs of m6A regulators were given their respective gene importance, and gene importance rank demonstrated that YTHDF1, RBM15, and METTL14 had priorities in the model (Fig. 3E).

Fig. 3
figure 3

Establishment and comparison of the RF model and SVM model. A Reverse cumulative distribution of |residual| of RF and SVM model. B Boxplots of |residual| for displaying the residual distribution of RF and SVM model. C ROC curves demonstrated the accuracy of the RF and SVM mode. D The random forest model displayed the influence of the number of decision trees on the error rate. E The importance of the m6A regulators based on the RF model

Evaluation of a predictive nomogram

Three recommended m6A regulators with a gene importance score above 2 were selected to establish a nomogram for evaluating the occurrence risk of PCOS by the "rms" package of R (Fig. 4A). Based on the multi-factor regression analysis, the nomogram was built to integrate involved genes, score them by scaled line segments, and obtain the total points to predict the occurrence risk of PCOS. The expression levels of these selected genes were positively correlated with the risk score in the nomogram, indicating that they may be risk factors for PCOS patients. The nomogram was consistent with the above analysis based on the expression differences of m6A regulators in control and PCOS patients. The calibration curve displayed that the apparent line, the bias-corrected line, and the ideal line shared a close range and the same tendency, which indicated the superior predictive accuracy of the nomogram (Fig. 4B). The DCA curve demonstrated that the red line developed by the m6A regulator was always at the top of and away from the gray line, which stood for the high accuracy of this nomogram and patients' benefit from the decisions based on this nomogram (Fig. 4C). In addition, the clinical impact curve also proved the nomogram's high benefit and superior predictive performance (Fig. 4D).

Fig. 4
figure 4

Establishment of a nomogram model. A Establishment of a nomogram based on the three selected m6A regulators, YTHDF1, RBM15, and METTL14. B Calibration curves of the nomogram indicating predictive robustness of this nomogram. C DCA curve revealing the high accuracy and patients’ benefit of the decisions based on this nomogram. D Clinical impact curve of the nomogram

Analysis of subtypes based on m6A regulators

Based on the 22 differentially expressed m6A regulators, the consensus clustering algorithm was applied to identify different subtypes of PCOS patients. They were well categorized into two clusters when the cluster variable was 2 (Fig. 5A). The expressions of m6A regulators were compared in m6A cluster A and m6A cluster B. The expressions of METTL14, ZC3H13, CBLL1, YTHDC2, YTHDF3, HNRNPC, FMR1, ELAVL1, and FTO presented increased values in cluster A compared to those in cluster B, while the opposite performance was observed in RBM15B, IGFBP2, IGFBP3, and ALKBH5. Furthermore, METTL3, WTAP, VIRMA, RBM15, YTHDC1, YTHDF1, HNRNPA2B1, IGFBP1, and IGF2BP1 showed no significant differences between clusters (Fig. 5B, C). According to the PCA, we could distinctly differentiate the two m6A clusters (Fig. 5D).

Fig. 5
figure 5

Consensus clustering of 22 DEGs of m6A regulators. A Consensus matrix for k = 2. B Histogram of expressions of 22 DEGs in m6A regulators between two m6A clusters A and cluster B. C Heatmap of expressions of 22 DEGs in m6A regulators between two m6A clusters. D Principal component analysis of m6A clusters. *p < 0.05, **p < 0.01, and ***p < 0.001

Evaluation of the m6A related gene signature

To attest to the m6A subtypes, the DEGs were identified between the two m6A clusters with an adjusted p < 0.01 and a |logFC|≥ 1.5, which were termed the m6A-related DEGs. The GO and KEGG enrichment analyses were performed to explore the potential role of m6A-related DEGs in PCOS. The GO enrichment analysis revealed that the DEGs were especially affluent in immune-related biological processes (Fig. 6A, B). In KEGG enrichment analysis, DEGs were particularly abundant in a variety of infection-related pathways (Fig. 6C). The consensus clustering algorithm was employed to categorize the PCOS patients into two distinct gene subgroups based on these m6A-related DEGs. The two genomic subtypes were in accordance with m6A subtypes, and Fig. 6D displayed the differential expression of these m6A-related DEGs in two gene clusters. The results of differential expressions of m6A regulators between the two gene clusters were also similar to those between the m6A clusters (Fig. 6E).

Fig. 6
figure 6

Analysis of the m6A-related DEGs in PCOS. A Circle plot of the GO analysis revealing the potential gene functions of the m6A-related DEGs on the occurrence and development of PCOS. B Bubble plot of the GO analysis. C KEGG analysis of the m6A-related DEGs. D Expression heatmap of the m6A-related DEGs in gene cluster A and gene cluster B. E Differential expression analysis of the 22 identified m6A regulators in the two gene clusters. *p < 0.05, **p < 0.01, and ***p < 0.001

Evaluation of the immune cells infiltration patterns

The ssGSEA was conducted to assess the enrichment of immune cells in PCOS specimens and to demonstrate the relationship between the m6A regulators and immune cells (Fig. 7A). Prominent associations were observed between YTHDF1 and various immune cells. Therefore, further investigation was developed into the enrichment of immune cells in PCOS patients with high or low YTHDF1 (Fig. 7B). The boxplot demonstrated that patients with low YTHDF1 expression had obviously enriched immune cells. Ultimately, the differential immune cell enrichment between the m6A clusters was conducted to show the different immune responses of patients in different subtypes. The findings demonstrated that cluster B displayed higher infiltrating levels of immune cells than cluster A, particularly MDSC, activated B cells, and activated CD8 T cells (Fig. 7C), which indicated that patients in m6A cluster B might have a positive immune response for PCOS. The immune infiltration pattern of two gene clusters was in accordance with the pattern of m6A subtypes (Fig. 7D).

Fig. 7
figure 7

Immune infiltration analysis. A Correlation heatmap between infiltrating immune cells and the 22 DEGs of m6A regulators. B Differences in the abundance of immune infiltration between high and low YTHDF1 expression levels. C Differential expressions of immune cell infiltration between m6A cluster A and m6A cluster B. D Differential expressions of immune cell infiltration between gene cluster A and gene cluster B. *p < 0.05, **p < 0.01, and ***p < 0.001

Relationship between m6A and gene subtypes with cytokines

In the pathophysiological mechanism of PCOS, pro-inflammatory cytokines play an essential role and lead to inflammation and poor oocyte quality [4]. To further confirm the close correlation between the m6A subtypes and PCOS, a comprehensive exploration of the association between m6A subtypes or m6A-related gene subtypes and cytokines was conducted. As the results displayed, the cytokines presented significant discrepancies in both m6A clusters and gene clusters (Fig. 8B, C). Remarkably, IL1A, IL1RN, IL6, and TNF were overexpressed in m6Acluster B and genecluster B compared to m6Acluster A and genecluster A, which was consistent with existing reports [20,21,22]. These findings validated the close correlation between m6Acluster B or genecluster B in PCOS.

Fig. 8
figure 8

Relationships of m6A subtypes and m6A-related gene subtypes. A Sankey diagram of the m6A subtypes, gene subtypes, and m6A scores. B Differential expression analysis of cytokines between m6A cluster A and cluster B. C Differential expression analysis of cytokines between m6A-related gene clusters. *p < 0.05, **p < 0.01, and ***p < 0.001

Discussion

PCOS is a multifactorial disease caused by genetic, endogenous, endocrine regulation, metabolic disorders, and environmental factors, and the treatment strategies for PCOS involve a combination of lifestyle optimization, medical interventions, which include metformin, hormonal contraceptives, antiandrogens, and so on, and if needed, infertility management [1, 18, 23]. But the diagnosis and following management are often complex and unsatisfactory, highlighting the need for further investigation of the pathogenesis of PCOS. Existing researches have shown that m6A regulators play a role in the biological process of PCOS [18, 24]. As most studies focus on a single m6A regulator, the overall roles of multiple m6A regulators in PCOS are not comprehensively recognized. Identifying the m6A modification patterns in PCOS will contribute to enhancing our understanding of the pathogenesis of the disease and exploring new therapeutic strategies.

In this study, we comprehensively explored the m6A modification characters in PCOS patients, and the workflow of this study is demonstrated in Fig. 9. The expression profile data of 4 datasets were integrated and analyzed, and the expression levels of m6A regulators were obviously overexpressed in PCOS patients compared to healthy controls, which indicated that m6A modification might be closely related to the development of PCOS. Previous studies have demonstrated the m6A modification levels were increased and the m6A profile was altered in the luteinized GCs of PCOS patients. However, the knockdown analysis of METTL14 proved the absence of METTL14 / m6A / FOXO3 regulation in the luteinized GCs of PCOS patients. But the cell-specific effects of m6A modification are worth investigating and looking for other targets [24]. The elevated expression of FTO may conduce to the dysfunction of GCs by upregulating FLOT2 in PCOS [25]. The expression of IGFBP2 has been reduced in PCOS patients [26] and the levels of IGFBPs are involved in the selection of the dominant follicle [27, 28]. The PPI network was created by STRING to show the relationships between the differentially expressed m6A regulators and to discover hub genes. However, the identified hub genes, METTL3, YTHDF1, and YTHDF3, were obtained from protein interaction relationships of the m6A regulators, which indicated the internal connections of the regulators, and further exploration of the roles of m6A regulators in PCOS required a more comprehensive analysis of the expression levels of m6A regulators.

Fig. 9
figure 9

The flow chart of this study

As a result, RF and SVM models were constructed and compared based on the DEGs of m6A regulators to select a proper model to identify optimal regulators from the DEGs of m6A regulators for predicting the occurrence risk of PCOS. In addition, the RF model has become one of the most popular machine learning algorithms for its high accuracy, favorable robustness, and convenience. Based on the RF model's importance of the m6A regulators, three optimal m6A regulators, YTHDF1, RBM15, and METTL14, were selected to establish a nomogram by multivariate analysis rather than univariate analysis due to the possibility to ignore the multivariate association in the m6A regulators. However, the model could not integrate more disease features of PCOS yet due to the absence of sufficient clinical information and other features of PCOS in the GEO public database. Then, the calibration curve, the DCA curve, and the clinical impact curve further validated this nomogram model and demonstrated that PCOS patients may benefit from the decisions based on this nomogram. In this nomogram, YTHDF1, RBM15, and METTL14 were positively correlated with the risk for PCOS. Though there are few pieces of research on the correlation between three selected m6A regulators and PCOS, previous studies have demonstrated that the three optimal m6A regulators play an essential role in self-renewal and development of female germ stem cells, ovulation, placental insufficiency, and female reproductive system neoplasms such as cervical and endometrial cancer [29,30,31,32,33,34]. The functions of three selected m6A regulators and this nomogram in PCOS need to be confirmed in further research.

Furthermore, we identified two m6A-related subtypes in PCOS by consensus clustering, named m6A cluster A and m6A cluster B, and obtained m6A-related DEGs between two m6A clusters for identifying gene clusters. The expressions of m6A regulators in m6A cluster A and gene cluster A were significantly higher than those in m6A cluster B and gene cluster B, indicating that the distinction between clusters A and B is essential and warrants additional study. Studies have stated that hormone abnormalities in PCOS can overstimulate the immune system, and concurrently, inflammatory cell infiltration aggravates ovarian dysfunction and metabolic disorders [35, 36]. High infiltrating abundances of macrophages [37], neutrophils, activated dendritic cells [38], and eosinophils [39, 40] were observed in PCOS. Since chronic low-level inflammation plays a critical role in PCOS [41], we investigated and compared the immune infiltration and cytokine levels in different groups. The m6A cluster B presented higher immune cell infiltration behaviors consistent with existing researches, and a similar infiltration trend marked in the genomic subtypes. Additionally, the cytokine levels in m6A cluster B or gene cluster B were significantly higher than those in clusters A, consistent with existing studies [41,42,43] especially in the up-regulation of IL1A [20], IL1RN [21] and IL6 [22]. Therefore, these results indicated that m6A cluster B and gene cluster B with higher infiltration of immune cells could present more typical characteristics of PCOS, which manifested a higher risk and an inferior clinical performance. Consistent with the findings, PCOS patients in m6A cluster B or gene cluster B obtained higher m6A scores than clusters A.

Currently, the clinical importance and underlying mechanisms of m6A modification in the etiology of PCOS remain unknown, and our findings can help to better understand the roles of m6A involved in PCOS pathogenesis. Additionally, the development of precise gene-targeted therapy [44] possibly provides a new direction for m6A-targeted strategies specifically for PCOS patients. Nonetheless, there are several limitations in the present research. The relatively small number of PCOS patients included in the datasets and the absence of clinical information result in inadequate support to establish a more comprehensive assessment of the specific relationship between m6A regulators and PCOS. The findings are preliminary, and future collaborations on basic, clinical, and epidemiologic research will be critical in determining the roles of m6A in the occurrence and progression of PCOS.

Conclusion

In conclusion, the m6A regulators might play a role in the pathogenesis of PCOS. We characterized the PCOS with m6A regulators pattern and developed two different m6A subtypes of PCOS, one of which (clusterB) showed higher levels of immune infiltration and cytokines. This study will help clarify the overall impact of m6A modification patterns and related immune infiltration on PCOS. In addition, the developing precise therapy targeting m6A regulators is expected to exert a therapeutic effect on PCOS and the ensuing infertility. Further exploration of the roles of m6A modification in PCOS will help to investigate the pathogenesis and potential therapeutic method of PCOS.

Availability of data and materials

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Abbreviations

PCOS:

Polycystic ovary syndrome

m6A:

N6-methyladenosine

DEG:

Differentially expressed gene

HA:

Hyperandrogenemia

IR:

Insulin resistance

AGE:

Advanced glycation end product

GEO:

Expression Omnibus

RNA-seq:

RNA sequencing

Limma:

Linear Models for Microarray Data

PPI:

Protein–protein interaction

GSVA:

Gene set variation analysis

RF:

Random forest

SVM:

Support vector machine

ROC:

Receiver operating characteristic

DCA:

Decision curve analysis

PCA:

Principal component analysis

ssGSEA:

Single sample gene set enrichment analysis

GC:

Granulosa cell

References

  1. Joham AE, et al. Polycystic ovary syndrome. Lancet Diabetes Endocrinol. 2022;10(9):668–80.

    Article  CAS  PubMed  Google Scholar 

  2. Lizneva D, Suturina L, Walker W, Brakta S, Gavrilova-Jordan L, Azziz R. Criteria, prevalence, and phenotypes of polycystic ovary syndrome. Fertil Steril. 2016;106(1):6–15.

    Article  PubMed  Google Scholar 

  3. Walters KA, Moreno-Asso A, Stepto NK, Pankhurst MW, Rodriguez PV, Rodgers RJ. Key signalling pathways underlying the aetiology of polycystic ovary syndrome. J Endocrinol. 2022;255(1):R1–26.

    Article  CAS  PubMed  Google Scholar 

  4. Siddiqui S, Mateen S, Ahmad R, Moin S. A brief insight into the etiology, genetics, and immunology of polycystic ovarian syndrome (PCOS). J Assist Reprod Genet. 2022;39(11):2439–73.

    Article  PubMed  Google Scholar 

  5. Carmina E, Chu MC, Longo RA, Rini GB, Lobo RA. Phenotypic variation in hyperandrogenic women influences the findings of abnormal metabolic and cardiovascular risk parameters. J Clin Endocrinol Metab. 2005;90(5):2545–9.

    Article  CAS  PubMed  Google Scholar 

  6. Desrosiers R, Friderici K, Rottman F. Identification of methylated nucleosides in messenger RNA from Novikoff hepatoma cells. Proc Natl Acad Sci U S A. 1974;71(10):3971–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR. Comprehensive analysis of mRNA methylation reveals enrichment in 3’ UTRs and near stop codons. Cell. 2012;149(7):1635–46.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Dominissini D, et al. Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature. 2012;485(7397):201–6.

    Article  CAS  PubMed  Google Scholar 

  9. Zaccara S, Ries RJ, Jaffrey SR. Reading, writing and erasing mRNA methylation. Nat Rev Mol Cell Biol. 2019;20(10):608–24.

    Article  CAS  PubMed  Google Scholar 

  10. Hu Y, et al. New sights in cancer: Component and function of N6-methyladenosine modification. Biomed Pharmacother. 2020;122: 109694.

    Article  CAS  PubMed  Google Scholar 

  11. Arguello AE, DeLiberto AN, Kleiner RE. RNA Chemical Proteomics Reveals the N(6)-Methyladenosine (m(6)A)-Regulated Protein-RNA Interactome. J Am Chem Soc. 2017;139(48):17249–52.

    Article  CAS  PubMed  Google Scholar 

  12. Huang H, Weng H, Chen J. The Biogenesis and Precise Control of RNA m(6)A Methylation. Trends Genet. 2020;36(1):44–52.

    Article  CAS  PubMed  Google Scholar 

  13. Liu H, Zheng J, Liao A. The regulation and potential roles of m6A modifications in early embryonic development and immune tolerance at the maternal-fetal interface. Front Immunol. 2022;13: 988130.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Sui X, Klungland A, Gao L. RNA m6A modifications in mammalian gametogenesis and pregnancy. Reproduction. 2023;165(1):R1–8.

    Article  CAS  PubMed  Google Scholar 

  15. Yang B, Wang JQ, Tan Y, Yuan R, Chen ZS, Zou C. RNA methylation and cancer treatment. Pharmacol Res. 2021;174: 105937.

    Article  CAS  PubMed  Google Scholar 

  16. Wei B, et al. N(6)-Methyladenosine RNA modification: a potential regulator of stem cell proliferation and differentiation. Front Cell Dev Biol. 2022;10: 835205.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Wang T, Kong S, Tao M, Ju S. The potential role of RNA N6-methyladenosine in Cancer progression. Mol Cancer. 2020;19(1):88.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Sun X, Lu J, Li H, Huang B. The role of m(6)A on female reproduction and fertility: from gonad development to ovarian aging. Front Cell Dev Biol. 2022;10: 884295.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Barrett T, et al. NCBI GEO: archive for functional genomics data sets–update. Nucleic Acids Res. 2013;41(Database issue):D991–5.

    CAS  PubMed  Google Scholar 

  20. Popovic M, Sartorius G, Christ-Crain M. Chronic low-grade inflammation in polycystic ovary syndrome: is there a (patho)-physiological role for interleukin-1? Semin Immunopathol. 2019;41(4):447–59.

    Article  PubMed  Google Scholar 

  21. Yang Y, Qiao J, Tang RX, Li MZ. Genotype and haplotype determination of interleukin (IL) 1 beta (g. -511C>T and g. +3954C>T) and IL-1RN in polycystic ovary syndrome. Fertil Steril. 2010;94(1):384–6.

    Article  CAS  PubMed  Google Scholar 

  22. Wang C, et al. LncRNA GAS5 is upregulated in polycystic ovary syndrome and regulates cell apoptosis and the expression of IL-6. J Ovarian Res. 2020;13(1):145.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Li Y, et al. Multi-system reproductive metabolic disorder: significance for the pathogenesis and therapy of polycystic ovary syndrome (PCOS). Life Sci. 2019;228:167–75.

    Article  CAS  PubMed  Google Scholar 

  24. Zhang S, Deng W, Liu Q, Wang P, Yang W, Ni W. Altered m(6) A modification is involved in up-regulated expression of FOXO3 in luteinized granulosa cells of non-obese polycystic ovary syndrome patients. J Cell Mol Med. 2020;24(20):11874–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Zhou L, et al. N6-methyladenosine Demethylase FTO induces the dysfunctions of ovarian granulosa cells by upregulating Flotillin 2. Reprod Sci. 2022;29(4):1305–15.

    Article  CAS  PubMed  Google Scholar 

  26. Haouzi D, Assou S, Monzo C, Vincens C, Dechaud H, Hamamah S. Altered gene expression profile in cumulus cells of mature MII oocytes from patients with polycystic ovary syndrome. Hum Reprod. 2012;27(12):3523–30.

    Article  CAS  PubMed  Google Scholar 

  27. Magoffin DA, San RG, Muderspach LI. Insulin-like growth factor binding proteins in a natural pre-ovulatory follicle from a woman with polycystic ovary syndrome. Hum Reprod. 1995;10(9):2248–52.

    Article  CAS  PubMed  Google Scholar 

  28. Kwon H, Choi DH, Bae JH, Kim JH, Kim YS. mRNA expression pattern of insulin-like growth factor components of granulosa cells and cumulus cells in women with and without polycystic ovary syndrome according to oocyte maturity. Fertil Steril. 2010;94(6):2417–20.

    Article  CAS  PubMed  Google Scholar 

  29. Yu XX, et al. Ascorbic acid induces global epigenetic reprogramming to promote meiotic maturation and developmental competence of porcine oocytes. Sci Rep. 2018;8(1):6132.

    Article  PubMed  PubMed Central  Google Scholar 

  30. Zheng H, et al. Decreased expression of programmed death Ligand-L1 by Seven in Absentia homolog 2 in Cholangiocarcinoma enhances T-Cell-mediated antitumor activity. Front Immunol. 2022;13: 845193.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Chen J, Fang Y, Xu Y, Sun H. Role of m6A modification in female infertility and reproductive system diseases. Int J Biol Sci. 2022;18(9):3592–604.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Kontur C, Jeong M, Cifuentes D, Giraldez AJ. Ythdf m(6)A Readers Function Redundantly during Zebrafish Development. Cell Rep. 2020;33(13):108598.

    Article  CAS  PubMed  Google Scholar 

  33. Li X, Tian G, Wu J. Novel circGFRα1 promotes self-renewal of female germline stem cells mediated by m(6)A writer METTL14. Front Cell Dev Biol. 2021;9: 640402.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Raffel GD, et al. Ott1 (Rbm15) is essential for placental vascular branching morphogenesis and embryonic development of the heart and spleen. Mol Cell Biol. 2009;29(2):333–41.

    Article  CAS  PubMed  Google Scholar 

  35. Luan YY, Zhang L, Peng YQ, Li YY, Liu RX, Yin CH. Immune regulation in polycystic ovary syndrome. Clin Chim Acta. 2022;531:265–72.

    Article  CAS  PubMed  Google Scholar 

  36. Hu C, Pang B, Ma Z, Yi H. Immunophenotypic profiles in Polycystic Ovary Syndrome. Mediators Inflamm. 2020;2020:5894768.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Yao X, Wang X. Bioinformatics searching of diagnostic markers and immune infiltration in polycystic ovary syndrome. Front Genet. 2022;13: 937309.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Qu J, et al. Discovery of immune-related diagnostic biomarkers and construction of diagnostic model in varies polycystic ovary syndrome. Arch Gynecol Obstet. 2022;306(5):1607–15.

    Article  CAS  PubMed  Google Scholar 

  39. Gao Q, et al. Exploration of molecular features of PCOS with different androgen levels and immune-related prognostic biomarkers associated with implantation failure. Front Endocrinol (Lausanne). 2022;13: 946504.

    Article  PubMed  Google Scholar 

  40. Na Z, Guo W, Song J, Feng D, Fang Y, Li D. Identification of novel candidate biomarkers and immune infiltration in polycystic ovary syndrome. J Ovarian Res. 2022;15(1):80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Rudnicka E, et al. Chronic low grade inflammation in pathogenesis of PCOS. Int J Mol Sci. 2021;22(7):3789.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. de Alencar JB, Alves HV, Elpidio LN, Visentainer JE, Sell AM. Polymorphisms of cytokine genes and Polycystic Ovary Syndrome: a review. Metab Syndr Relat Disord. 2016;14(10):468–74.

    Article  PubMed  Google Scholar 

  43. Zhai Y, Pang Y. Systemic and ovarian inflammation in women with polycystic ovary syndrome. J Reprod Immunol. 2022;151: 103628.

    Article  CAS  PubMed  Google Scholar 

  44. Collins FS, Varmus H. A new initiative on precision medicine. N Engl J Med. 2015;372(9):793–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This work was supported by National Key Research and Development Program of China (2016YFC1000205), National Natural Science Foundation of China (81901566), and Natural Science Foundation of Guangdong Province (2022A1515011732).

Author information

Authors and Affiliations

Authors

Contributions

SZ designed the study, downloaded the raw data, performed the process of the statistical analyses, and prepared the manuscript, table and figures. RH checked the table and figures and revised the manuscript. SQ designed the study and edited the manuscript. All authors contributed to the article and approved the submitted version.

Corresponding authors

Correspondence to Rui Hua or Song Quan.

Ethics declarations

Ethics approval and consent to participate

The study is an integrated analysis based on published datasets from Gene Expression Omnibus (GEO) database. Therefore, ethics approval is not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhou, S., Hua, R. & Quan, S. N6-methyladenosine regulator-mediated methylation modification patterns and immune infiltration characterization in Polycystic Ovary Syndrome (PCOS). J Ovarian Res 16, 73 (2023). https://doi.org/10.1186/s13048-023-01147-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13048-023-01147-9

Keywords