PCOS patients who underwent laparoscopic surgery for hydrosalpinx at the Department of Obstetrics and Gynecology, West China Second University Hospital, from January 2022 to February 2022, were enrolled in this study. Patients were diagnosed with PCOS according to the criteria of the European Society for Human Reproduction and Embryology/American Society for Reproductive Medicine (ESHRE/ASRM) (Rotterdam criteria) (2). Two or more of the following three criteria should be met including hyperandrogenism, oligo-/anovulation, and polycystic ovaries on transvaginal ultrasound (TVS). Patients with any known significant disorders that may overlap or confound the diagnosis of PCOS, such as congenital adrenal hyperplasia, functional hypothalamic amenorrhea (FHA), Cushing syndrome, thyroid diseases, hyperandrogenism due to non-PCOS causes, androgen-secreting ovarian tumors or adrenal tumors, and premature ovarian insufficiency (POI), were excluded. Patients who underwent laparoscopic surgery for hydrosalpinx without any clinical manifestation of PCOS were enrolled as the control group. Abdominal adipose tissue and subcutaneous adipose tissue were obtained, and samples were stored at −80°C until use. The Ethical Review Board of West China Second University Hospital, Sichuan University, approved this study [Approval number: No. 2021(033)].

The RNA‐seq dataset for abdominal adipose tissue and subcutaneous adipose tissue from the PCOS and control patients was uploaded to the Gene Expression Omnibus (GEO) database with accession number PRJNA1010231. Total RNA was harvested using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the instructions provided by the manufacturer. Total amounts and integrity of RNA were assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). RNA libraries were prepared for sequencing using the NEBNext^® Ultra RNA Library Prep Kit for Illumina^® (NEB, USA). Illumina NovaSeq 6000 was used, and the qualified libraries were pooled and sequenced on Illumina platforms with the PE150 strategy of Novogene Bioinformatics Technology Co., Ltd. (Beijing, China). Raw data of fastq format were first processed through in-house perl scripts. Then, we used the following parameters to trim the data: fastp -n 15 -q 5 -u 50 -l 150. Reads containing adapter, reads containing N base, and low-quality reads from raw data were removed to obtain clean data. Concurrently, Q20, Q30, and GC content clean data were calculated, and all of the downstream analyses were based on the high-quality clean data. Paired-end clean reads were mapped to the reference human genome (GRCh38) using Hisat2 (version 2.0.5) (http://github.com/infphilo/hisat2), and alignment files were generated in Binary Alignment Map (BAM) format. To quantify the gene expression level, featureCounts (version 1.5.0-p3) was adopted to count the read numbers mapped to each gene. Then, FPKM, expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced, of each gene was calculated based on the length of the gene, and read counts were mapped to this gene.

The WGS dataset for abdominal adipose tissue and subcutaneous adipose tissue from the PCOS and control patients was uploaded to the GEO database with accession number PRJNA1010231. Genomic DNA was extracted from the adipose samples using the standard phenol–chloroform protocol. DNA libraries were generated using the NEB Next^® Ultra™ DNA Library Prep Kit for Illumina (NEB, USA) according to the manufacturer’s instructions. DNA concentration was measured using Qubit^® DNA Assay Kit in Qubit^® 3.0 Flurometer (Invitrogen, USA). The sequencing libraries were analyzed for size distribution using NGS3K/Caliper and quantified using real-time PCR (3 nM). Clustering of the index-coded samples was performed on a cBot Cluster Generation System using Illumina PE Cluster Kit (Illumina, USA). Illumina NovaSeq 6000 was used, and the qualified libraries were pooled and sequenced on Illumina platforms with the PE150 strategy of Novogene Bioinformatics Technology Co., Ltd. (Beijing, China). We used Fastp (version 0.19.7) to perform quality control analysis of the generated raw data. A paired read was discarded if either one read contains adapter contamination, more than 10% of bases are uncertain in either one read, or if the proportion of low quality (Phredquality <5) bases is over 50% in either one read. Sequencing data is mapped to the reference human genome (GRCh38), alignment files in BAM format were obtained, and VCF files were generated by an intelligent web portal (GTX.Digest system) for genomics data interpretation (17). The filter parameters of single-nucleotide variations (SNVs) are shown as follows: genotype quality score ≥30; gene quality score ≥2.

We employed REDItools software for detecting candidate sites (18). Subsequently, mismatches between the initial reads and the reference genome were collected. Multiple test correction was conducted to adjust the p-value using the p.adjust function in R, and RNA-editing sites with a false discovery rate >0.05 were removed. RNA editing was annotated using the Annovar software (19). The Ensembl and National Center for Biotechnology Information (NCBI) databases (Homo_sapiens.GRCh38.100.gtf) provided the reference gene set for annotation. We also extracted the variants with A/G and C/T mismatch or T/C and G/A mismatch corresponding to plus and minus strand (20). The RNA-editing events with an alternative allele frequency (AAF) of 1% or higher and less than 100% were included. To eliminate SNPs, we download the SNP data from the dbSNP database (https://ftp.ncbi.nih.gov/snp/organisms/human_9606/BED/) and Ensembl database (https://ftp.ensembl.org/pub/release-111/variation/vcf/homo_sapiens/), and bedtools was used to remove all SNPs in our data from the databases. We also compared the detected RNA-editing sites with the known human RNA-editing sites from the REDIportal V2.0 database. Repeat sequences were obtained from the University of California, Santa Cruz (http://genome.ucsc.edu). Alu sequences were extracted using a Python program, and subsequently, the editing sites within the Alu sequences were analyzed (21). Visualization of RNA-editing sites of PCOS and control groups in different chromosomes was performed using the online web tool shinyCircos-V2.0 (https://venyao.xyz/shinycircos/) (22, 23).

We obtained essential (Ess) and nonessential (Noness) genes of humans from the OGEE (Online GEne Essentiality) database (https://v3.ogee.info/static/files/gene_essentiality.txt.gz) (24). We investigated the effect of gene essentiality of the detected RNA-editing genes because essential genes are those that are functionally critical for the survival of an organism or cell and capture the intricate nature of biological systems. Building upon this distinction, we discovered substantial differences in the numbers of edited genes between Ess and Noness categories.

In this study, we utilized the bamCoverage function of deeptools (version 3.5.5) to convert aligned bam files into bw files, therefore facilitating visualization of gene expression abundance in the integrative genomics viewer (IGV) (version 2.17.4) (http://www.broadinstitute.org/software/igv/download) (25). Additionally, we loaded bam files from WGS and RNA-seq to display the read mapping results across the entire genome and examined the RNA-editing sites of genes.

To illustrate and cluster the potential differences between PCOS and control groups, we adopted the partial least squares-discriminant analysis (PLS-DA) method (26). PLS-DA is a supervised dimensionality reduction method and commonly used in metabolomics since it can more effectively highlight metabolites that differentiate between conditions. R package ropls (version 1.34.0) was used to perform PLS-DA, and R package ggplot2 (version 3.3.6) was used to visualize the result.

The DESeq2 R package (version 1.20.0) (27) was utilized to conduct differential expression analysis among distinct groups. DESeq2 offers statistical routines to ascertain differential expression in digital gene expression data by employing a model grounded on the negative binomial distribution. To control the false discovery rate (FDR), the raw p-values underwent adjustment using the Benjamini–Hochberg (BH) method. The threshold for significantly differentially expressed genes was established as p.adjust ≤ 0.05 and |log2(foldchange)| ≥ 1.

The functional enrichment analysis was conducted using R package clusterProfiler (version 3.16.0) (28) to identify significantly enriched Gene Ontology (GO) biological processes terms. We not only performed GO analysis on the deferential RNA-editing sites but also on the differentially expressed genes between the abdominal or subcutaneous adipose tissue of PCOS and control patients. The BH method was employed for adjusting the p-value to identify significantly enriched GO biological processes terms, and terms with p.adjust <0.05 were considered significantly enriched.

Wilcoxon tests or Kruskal–Wallis tests were used in identifying differential RNA-editing sites and genes and RNA editing frequencies between two or more groups. The χ2 and F-tests were used as the statistical analysis for categorical variable, such as proportion of essential and nonessential genes. A value of p < 0.05 was used as the significance cutoff. For gene expression data, raw p-values were adjusted using the BH method, and adjusted p < 0.05 was considered significant in detecting functional enriched pathways with differentially edited or differentially expressed genes. The association between RNA editing level and gene expression level was calculated using Spearman’s rank correlation to obtain the correlation coefficient (r) and p-value (29). Statistical analyses were performed using R or GraphPad Prism 9 software.

This study enrolled five PCOS patients and three control patients ( Supplementary Table S1 ). Paired subcutaneous and abdominal adipose tissues were collected and used for WGS and RNA sequencing, and RNA-editing sites were identified ( Supplementary Figures S1A, B ). In subcutaneous adipose tissue of PCOS patients, an average of 8,213 RNA-editing events and 3,387 RNA-editing genes were identified, 7,898 sites and 3,329 genes were edited in the abdominal adipose tissue of PCOS, 4,575 sites and 2,366 genes in the subcutaneous adipose tissue of the control, and 5,267 sites and 2,911 genes in the abdominal adipose tissue of the control ( Supplementary Figures S2A, B ). Compared with the healthy controls, we found more RNA-editing sites and edited genes in the subcutaneous and abdominal adipose tissue of PCOS patients, but the differences in editing genes did not reach statistical significance. In addition, we found that approximately half of our detected RNA-editing sites match those already cataloged in the REDIportal V2.0 database ( Supplementary Figure S2C ). Circos plot was drawn to show the distributions of RNA-editing events in different chromosomes of different groups ( Figure 1A ). As shown in the figure, chromosome 1 contains the most RNA-editing events because it is the longest and, therefore, contains genes encoding the largest number of transcripts capable of undergoing RNA editing. The main types of nucleotide substitutions in humans are A‐to‐I and uridine‐to‐cytidine (U‐to‐C) followed by cytidine‐to‐uridine (C‐to‐U) and guanine-to-adenosine (G-to-A) ( Figure 1B ).

It is well known that A-to-I editing is catalyzed by ADAR enzymes (30). Hence, we investigated whether the differences in RNA editing are related to the differential expression of the ADAR genes. First, we calculated the expression levels of the ADAR gene using the RNA-seq data. The total expression of the ADARs is higher in PCOS than in the control, although it did not reach statistical significance ( Figure 2A ). Statistically significant correlations were also observed for ADAR gene expression and the number of editing sites (r = 0.531, p = 0.0347) ( Figure 2B ).

In this study, we found that RNA editing in adipose tissue mainly occurs in UTR3, exonic regions, introns, and noncoding regions ( Figure 3A ). Our findings demonstrated that nonsynonymous editing was more than synonymous editing in the adipose tissue of the PCOS and healthy control (all p < 0.05) ( Figure 3B ). These findings indicated that RNA editing may play an important role in the adipose tissue of PCOS. We also found that the number of RNA-editing events in the Alu regions in PCOS remains higher than in the control probably due to more total RNA-editing sites in PCOS than in the control ( Figure 3C ). In addition, the proportion of RNA-editing events in the Alu repeating region is also higher in the subcutaneous and abdominal adipose tissues of PCOS than in the control ( Figure 3D ).

Next, we identified the Ess and Noness genes of humans from the OGEE database. Ess genes are crucial for the survival of organisms, whereas Noness genes are not indispensable and rely on external conditions rather than intrinsic properties for their necessities. Then, we counted the Ess and Noness genes at the editing site of the CDS, and our results demonstrated that the number of Noness genes among all editing genes or nonsynonymous editing genes in humans is significantly higher than the Ess genes in all groups (all p < 0.05) ( Supplementary Figures S3A, B ). The higher number of editing sites in Noness genes is due to the fact that there are typically more Noness genes than Ess genes. To explore whether these nonsynonymous editing genes differ in the necessities of genes, we compared the editing frequency of Ess and Noness genes in the nonsynonymous editing genes. There was no significant difference in the editing frequency of Ess and Noness genes in all groups (all p > 0.05 in Supplementary Figure S3C ).

We also investigated the RNA editing level, also known as RNA editing frequency, in abdominal and subcutaneous adipose tissues of PCOS and control patients. There were no significant differences in RNA editing level between abdominal and subcutaneous adipose tissues of PCOS and control patients ( Supplementary Figure S4A ). Our findings revealed a slightly higher editing frequency of nonsynonymous genes than synonymous genes in subcutaneous adipose tissue of the PCOS group ( Supplementary Figure S4B ). In summary, we verified that RNA editing tends to present in Noness genes and that similar results exist in both abdominal and subcutaneous adipose tissues of PCOS and control patients.

It is also found that RNA editing differs in different stages, chromosomes, and editing types. In humans, RNA-editing sites are concentrated on chromosome 1 ( Figure 1C ). We explored the genes with nonsynonymous RNA editing and performed function enrichment analysis. First, Venn plots was drawn to identify tissue or PCOS-specific RNA-edited genes ( Supplementary Table S2 and Supplementary Figure S5 ). We found that only 2,034 genes in the abdominal adipose tissue of PCOS, but not in the control patients, were mainly enriched in immune responses, including B-cell tolerance induction, positive regulation of tumor necrosis factor (TNF) secretion, and tolerance induction dependent upon immune response pathways ( Figure 4A ). More specifically, major histocompatibility complex, class I, G (HLA-G), leukocyte immunoglobulin-like receptor A5 (LILRA5), leucine-rich repeat kinase 2 (LRRK2), Wnt family member 5A (WNT5A), frizzled class receptor 5 (FZD5), and LILRA2 genes were only edited in the abdominal adipose tissue of PCOS but not in the control. Similarly, only 2,077 genes edited in the subcutaneous adipose tissue of PCOS, but not in the control patients, were mainly enriched in interleukin-12 (IL-12) biosynthetic process and stem cell population maintenance ( Figure 4B ). Genes enriched in inflammation, such as lymphotoxin beta (LTB), TNF receptor-associated factor 6 (TRAF6), and nuclear factor kappa B subunit 1 (NFKB1), and genes enriched in stem cell population maintenance, such as delta-like canonical Notch ligand 1 (DLL1), notch receptor 1 (NOTCH1), NOTCH2, replication timing regulatory factor 1 (RIF1), FZD7, terminal uridylyl transferase 4 (TUT4), hook microtubule-tethering protein 3 (HOOK3), DNA ligase 4 (LIG4), forkhead box P1 (FOXP1), and forkhead box O1 (FOXO1), were merely edited in the subcutaneous adipose tissue of PCOS. In addition, genes differently edited between the abdominal and subcutaneous adipose tissues of PCOS patients were also enriched in immune processes ( Supplementary Figure S6 ).

RNA sequencing data were also analyzed to investigate the differentially expressed genes between the adipose tissue of PCOS and control patients. PLS-DA plot and heatmap of the cluster analysis of different samples revealed great homogeneity within groups ( Supplementary Figure S7 ). A total of 1,030 differentially expressed genes were identified between the subcutaneous adipose tissue of PCOS and control patients ( Supplementary Table S3 and Figures 5A, B ). GO enrichment analysis revealed that genes differentially expressed between the subcutaneous adipose tissue of PCOS and control patients were mainly involved in the immune responses such as Toll-like receptor 2 signaling pathway, regulation of B-cell activation, macrophage activation, lymphocyte activation involved in immune response, positive regulation of inflammatory response, and lymphocyte differentiation ( Figure 5C ). This is in accordance with our previous findings that genes only edited in the subcutaneous adipose tissue of PCOS but not in control patients were mainly enriched in immune responses ( Figure 4B ). In addition, we also drew a heatmap for the differentially expressed immune-related genes ( Figure 5D ).

We also identified 4,306 differentially expressed genes between abdominal adipose tissue of PCOS and control patients. However, instead of immune processes, genes differentially expressed between the abdominal adipose tissue of the PCOS and control patients were mostly enriched in signal transduction, cell structure organization, metabolic process, etc. ( Supplementary Table S4 and Supplementary Figure S8 ). We also explore the differences between the abdominal and subcutaneous adipose tissue of PCOS patients ( Supplementary Table S5 and Supplementary Figure S9A ). The upregulated genes in the subcutaneous adipose tissue of PCOS, compared with the abdominal adipose tissue of PCOS, were also mainly enriched in immune responses ( Supplementary Figure S9B ), and the downregulated genes were enriched in cellular process, cell death, etc. ( Supplementary Figure S9C ). Additionally, some known adipokine genes, such as chemerin-encoding gene retinoic acid receptor responder 2 (RARRES2), were also detected in our datasets ( Supplementary Tables S3 and S4 ).

Finally, we aimed to study the relationship between RNA editing level and RNA expression level of the differentially expressed genes and focused on genes of inflammation and immune response pathways. Our results found that there was a positive correlation between RNA editing level and RNA expression level in some genes, such as ataxia–telangiectasia mutated (ATM) (r = 0.53, p = 0.036), CD8 subunit alpha (CD8A) (r = 0.76, p = 0.00067), lymphoid enhancer-binding factor 1 (LEF1) (r = 0.55, p = 0.027), mucosa-associated lymphoid tissue lymphoma translocation protein 1 (MALT1) (r = 0.56, p = 0.025), and Toll-like receptor 7 (TLR7) (r = 0.51, p = 0.046) ( Supplementary Figure S10 ). Next, we used IGV software to visualize these RNA-editing sites and the same regions covered by DNA-seq data, such as ATM and MALT1 ( Supplementary Figure S11 ). Our data revealed RNA-editing events occurring at position chr18:58,753,294 in MALT1, and the RNA-seq results also showed significant changes in expression levels near this site. Moreover, the editing level at this site was notably higher in PCOS patients than those in normal controls. However, further experimental validation is required to determine whether this editing event can actually influence gene expression.

Moreover, we also looked into all the genes that are both differentially edited and differentially expressed between the PCOS and control groups. Our results indicated that 445 differentially expressed genes are exclusively edited in the abdominal adipose tissue of PCOS, and 94 differentially expressed genes are only edited in the subcutaneous adipose tissue of PCOS ( Figures 6A, B ). Heatmaps were drawn to directly visualize the RNA expression level of differentially expressed genes that were exclusively edited in PCOS ( Figures 6C, D ). Interestingly, using functional enrichment analysis ( Figures 6E, F ), we found three subcutaneous PCOS-specific edited differentially expressed genes, ATM, TUT4, and YTH N6-methyladenosine RNA-binding protein C2 (YTHDC2), enriched in the oocyte development and differentiation pathway, which may also indicate the important role of adipose tissue in the pathogenesis of PCOS.