Genome-wide profiling of m6A modification and its role in the transcriptome of sequelae of pelvic inflammatory disease.

Several limitations in our experimental and methodological approach should be considered. First, the MeRIP-seq technology itself, while enabling genome-wide profiling, has inherent constraints: its resolution is limited to 100 to 200 nucleotides, precluding single-nucleotide mapping, and it provides relative rather than absolute quantification of m6A (e.g., vs LC-MS/MS). Furthermore, a technical bias exists where highly expressed transcripts are more readily detected, which we mitigated through input normalization but cannot be entirely eliminated. Specifically regarding immunoprecipitation, we acknowledge that the absence of an IgG control is a limitation, as it could further help account for nonspecific antibody binding. Future studies should employ techniques with higher resolution for validation, include IgG controls to ensure IP specificity, and utilize LC-MS/MS to obtain an absolute quantification of global m6A levels.

The sequelae of pelvic inflammatory disease (SPID) are legacy lesions of pelvic inflammatory disease (PID). PID refers to inflammatory lesions in the female upper reproductive tract. Untimely treatment or poorly controlled inflammation can prolong the course of the disease, thus causing SPID. SPID pathologically manifests as tissue adhesion, hyperplasia, scar tissue formation, and chronic aseptic inflammation in the later stage. [ 1 ] Moreover, SPID is clinically characterized by symptoms such as persistent lower abdominal pain, lumbosacral pain, abnormal leucorrhea, and infertility. [ 2 ] According to epidemiological surveys, PID is mostly prevalent among young sexually mature women, with a global incidence of ~2% to 12%. [ 3 , 4 ] It has been shown that women with a history of PID have a 20-fold greater risk of SPID than those without a history of PID; SPID patients have a greatly increased risk of secondary infertility and chronic pelvic pain compared with women of normal childbearing age. [ 5 ] Therefore, SPID causes great harm to the reproductive health and quality of life of women of childbearing age. N6-methyladenosine (m6A) methylation is the methylation of the sixth N atom of adenine base A of mRNA, which is currently considered to be the most common methylation in prokaryotes. The m6A methylation process involves 3 main classes of proteins: methyltransferases, demethylases, and recognition proteins. [ 6 , 7 ] Owing to the coordinated action of these 3 types of proteins, the m6A methylation process exhibits dynamic and reversible characteristics. [ 8 ] m6A methylation is a key process that regulates mRNA splicing, translocation, translation, stability, and advanced structure. [ 7 ] Moreover, m6A methylation-mediated regulation of gene expression plays an important role in the physiological and pathological processes of organisms. [ 9 ] In recent years, as m6A detection technology has progressed, the relationship between m6A modification and the pathogenesis of various diseases has been gradually explored. An increasing number of studies have confirmed that m6A methylation is closely related to various systemic diseases of the body. [ 10 – 13 ] Moreover, m6A methylation has also been reported to be involved in mediating the occurrence and development of diverse inflammatory diseases, such as acute pancreatitis, [ 14 ] myocarditis, [ 15 ] and osteoarthritis. [ 16 ] Wang Xiaotong reported that METTL3 -mediated m6A methylation of SIRT1 mRNA affects the ectopic implantation of embryonic stem cells by regulating inflammation and autophagy, thereby hindering the progression of endometriosis. [ 17 ] Wang Huamin et al also suggested that m6A methylation could activate the immune response by promoting the function of dendritic cells in T-cell activation and the translation of CD40, CD80, and other proteins. [ 18 ] m6A methylation is important in the regulation of inflammation and the immune response. [ 19 ] Hence, there may be a potential link between m6A methylation and the development of SPID. This study aimed to establish an expression profile of m6A methylation in SPID patients and to explore the pathogenesis of SPID at the mRNA level on the basis of m6A methylation levels. The m6A methylation level and mRNA expression level in SPID patients were analyzed via m6A-seq and RNA-seq at the whole-transcript level. The sequencing results were verified via MeRIP-qPCR and RT-qPCR. Regulating m6A methylation sites is expected to be a therapeutic direction for SPID.

Conceptualization: Min Zhang, Yunxia Xu. Funding acquisition: Yunxia Xu. Investigation: Min Zhang, Yanling Jiang. Methodology: Min Zhang. Resources: Yunxia Xu. Supervision: Yunxia Xu. Writing – original draft: Min Zhang. Writing – review & editing: Min Zhang, Yunxia Xu.

A total of 10 subjects who visited the Department of Gynecology of the First Affiliated Hospital of Anhui University of Traditional Chinese Medicine were included in this study, comprising 5 SPID patients (SPID group) and 5 healthy volunteers (control group) as biological replicates. Total RNA was extracted from each subject’s blood sample, and subsequent MeRIP-seq and RNA-seq library construction and sequencing were performed independently for each sample. The recruitment period was from August 1, 2022, to August 1, 2023. Total RNA was extracted from the patient’s peripheral blood via TRIzol reagent. RNA was immunoprecipitated via an M6A-specific antibody (Sigma-Aldrich, cat. no ABE572). m6A sequencing was performed by Shanghai Bohao Biotechnology Co., Ltd. (Shanghai, China). m6A RNA immunoprecipitation was performed via the GenSeqTM M6A-MERIP Kit (GenSeq Inc., Cyberjaya, Malaysia) following the manufacturer’s instructions. The input samples without immunoprecipitation and m6A IP samples were used for RNA-Seq library generation. Library quality was evaluated via an analyzer. The detailed quality metrics for each individual library are provided in Table S1, Supplemental Digital Content, https://links.lww.com/MD/Q951 . Library sequencing was performed on an Illumina NovaSeq instrument with a 150-bp paired-end reading. The raw data were unspliced and subjected to quality control via Cutadapt (v2.5), and clean reads were obtained. Sequencing data quality analysis was conducted via FastQC ( https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ ), and the sequencing weight distribution, base content distribution, and information such as the repeat sequence fragment ratio were obtained. The remaining reads were then aligned to the human Ensemble genome GRCh38 via the HISAT2 aligner (v2.1.0) under the parameter “--rna-strandness RF.” The enriched region was identified via the exomePeak R software package (v2.13.2) with the parameters “PEAK_CUTOFF_PVALUE = 0.05, PEAK_CUTOFF_FDR = 0.05, FRAGMENT_LENGTH = 200,” that is, peak calling. The peak was then annotated and analyzed. Differential m6A peaks were identified via exome peak R encapsulation with the following parameters: “PEAK_CUTOFF_PVALUE = 0.05, PEAK_CUTOFF_FDR = 0.05, FRAGMENT_LENGTH = 200.” m6A-RNA-associated genomic features were visualized via the Guitar R package (v1.16.0). The m6A peak with an False Discovery Rate (FDR)  < 0.05 was selected, and motif analysis of the enriched region was conducted via HOMER (v4.10.4) with the parameters “-len 6 -rna.” Gene expression was normalized via fragments per kilobase million reads, followed by gene expression statistics across genomes via feature counts. The obtained counts were subjected to standardized processing through DESeq2, and differential expression analysis was conducted with FDR  1 as selection criteria. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed via the clusterprofile R package (v3.6.0). GO and KEGG analyses were used for functional enrichment analysis of differentially expressed m6A-methylated mRNAs and differentially expressed mRNAs. [ 20 ] The significance level of each entry was calculated via the Fisher test. Significantly enriched terms were defined as those with an FDR < 0.05. MeRIP-qPCR and RT-qPCR were used to detect the m6A methylation level and expression level of the candidate genes, respectively. Total RNA was extracted from the peripheral blood of 5 SPID patients via TRIzol (Thermo, 15596018). An immunocapture solution was prepared, and RNA was released/recovered. cDNA was obtained via enrichment and RT (42°C for 20 minutes and 75°C for 5 minutes) on a PCR apparatus. The fluorescence quantitative PCR was performed with cDNA as a template for fluorescence quantification. The primer sequences of each detection indicator are shown in Table 1 . The reaction was carried out under the following conditions: 95°C for 60 seconds, followed by 40 cycles of 95°C for 20 seconds and 60°C for 45 seconds. Primers for each detection index. bp = base pair. SPSS Statistics 24.0 and GraphPad Prism 9 software were used for statistical processing and analysis. The significance level was calculated through the Fisher test. The paired t test was used for the analysis of comparisons between groups. A P -value of <.05 was considered to indicate a statistically significant difference.

To validate the quality and saturation of our sequencing data, we performed sequencing saturation analysis, as shown in Figures S1 and S2, Supplemental Digital Content, https://links.lww.com/MD/Q951 , which depict the saturation curves for all Input and IP samples, respectively. The peak m6A methylation in the peripheral venous blood of SPID patients and controls was compared. The difference in m6A methylation between the SPID group and the control group is shown in Figure 1 A. There were 14,065 and 12,204 m6A methylation-modified genes in the SPID group and the control group, respectively, of which 5698 were unique modifier genes in the SPID group. Moreover, compared with the control group, the SPID group had 1861 nonmodifier genes with m6A methylationmodifications. Additionally, there were a total of 24,987 peaks in the control group and 28,578 peaks in the SPID group, of which 8364 peaks were unique to the SPID group (Fig. 1 B). Distribution of m6A-modified epitranscriptomes. (A) Histogram of peaks common to and unique to the SPID group and the control group; (B) histogram of m6A-modified genes common to and endemic to the SPID group and the control group; (C) pie chart of the m6A-modified distribution in the SPID group and the control group; (D) overall distribution of m6A modifications on mRNAs in the SPID group and the control group; (E) histogram of the number of peaks on different chromosomes for the SPID group and the control group; (F) MeRIP-seq motif analysis in the SPID group. In panels A, B, E, the control group is represented in orange and the SPID group in green. MeRIP = methylated RNA immunoprecipitation, SPID = sequelae of pelvic inflammatory disease. As shown in Figures 1 C, 55.9% of the m6A methylation peaks of mRNAs were distributed in the 3’ untranslated region, and 24.6% were in the coding sequence (CDS) in the control group. In contrast, the m6A methylation of mRNAs in the SPID group was increased by 3% in the 3’ UTR and decreased by 1.6% in the CDS region. This difference is presented more intuitively in Figure 1 D. Additionally, an analysis of the distribution of m6A peaks in each gene revealed that ~54% of the genes had separate m6A methylation sites. The distribution of m6A methylation at different chromosomal loci was analyzed. The results suggested that chromosome 1 had the highest m6A methylation frequency, with 3073 m6A methylation peaks, followed by chromosome 2, with 2464 m6A methylation peaks, and then chromosome 12, with 1909 m6A methylation peaks (Fig. 1 E). Through further comparison, the SPID group presented more m6A peaks on each chromosome than did the control group did, and there was no significant difference in the distribution of m6A peaks on chromosomes between the SPID group and the control group. Through analysis of the mRNA sequences corresponding to the m6A methylation peaks, the m6A motifs in SPID were measured and are shown in Figure 1 F, where the most significant mRNA methylation occurred at the GGAC motif, which was consistent with the m6A modification classical site. Variance analysis of 2 groups of gene sequences was performed with |log2 FC| > 1 and FDR < 0.05 as screening standards. The differentially expressed genes (DEGs) were subjected to visualization analysis via a volcano map. Among the 1776 genes with m6A methylation, 205 genes (65 hypermethylated genes and 140 demethylated genes) with significant differences in m6A methylation were screened. Table 2 shows the top 10 hypermethylated genes and the top 10 hypomethylated genes. A volcano map of differentially expressed m6A methylation genes was generated (Fig. 2 A). Top 10 hypermethylated genes and top 10 hypomethylated genes. Chr = chromosome, FC = fold change, Hyper = hypermethylation, Hypo = hypomethylation. Transcriptome information analysis of differentially expressed m6A-methylated genes in the SPID. (A) Volcano map of differentially expressed m6A methylation genes between the SPID group and the control group; (B) KEGG pathway analysis of differentially expressed m6A-modified mRNAs between the SPID group and the control group; (C) GO analysis of differentially expressed m6A-modified mRNAs between the SPID group and the control group. In the volcano plot (A), hypermethylated genes in the SPID group are shown in red, hypomethylated genes in blue, and nonsignificant genes in gray. FC = fold change, GO = Gene Ontology, KEGG = Kyoto Encyclopedia of Genes and Genomes, SPID = sequelae of pelvic inflammatory disease. As shown by the enrichment analysis results of the identified DEGs, the DEGs were involved mainly in biological processes such as neutrophil activation, neutrophil-mediated immunity, proteasomal protein decomposition, and neutrophil activation involved in the immune response. Moreover, the KEGG enrichment analysis results revealed that 180 pathways were significantly enriched, of which 84 signaling pathways were upregulated and 89 signaling pathways were downregulated. Most of the enriched signaling pathways were related to RNA splicing and translation. The results of the GO and KEGG enrichment analyses of the differentially expressed m6A-methylated genes are depicted in Figure 2 B and C. After normalization of gene expression via fragments per kilobase million reads, the distribution intensity and abundance of normalized expression after the alignment of mRNA sequencing sequences to the genome were depicted in the violin diagram (Fig. 3 A). Differentially expressed genes were screened for FDR  1 and then subjected to visualization analysis with a volcano map and heatmap (Figure 3 B and C). Among the 310 DEGs, 190 significantly DEGs were obtained, including 74 upregulated genes and 116 downregulated genes. Table 3 presents the top 10 upregulated genes and the top 10 downregulated genes. Top 10 upregulated genes and top 10 downregulated genes. Transcriptome information analysis of differentially expressed mRNAs in the SPID. (A) Fragments per kilobase million reads distribution violin diagram of the SPID group and the control group; (B) volcano map of differential mRNA expression genes between the SPID group and the control group; (C) cluster plot of differential mRNA expression genes in the SPID group and the control group; (D) KEGG pathway analysis of differential m6A-modified mRNAs in the SPID group and the control group; (E) GO analysis of differential m6A-modified mRNAs in the SPID group and the control group. In the volcano plot (B), upregulated genes are shown in red, downregulated genes in blue, and nonsignificant genes in gray. GO = Gene Ontology, KEGG = Kyoto Encyclopedia of Genes and Genomes, SPID = sequelae of pelvic inflammatory disease. According to the GO analysis results, 163, 163, and 161 GO terms were significantly enriched in biological processes, cellular components, and molecular functions, respectively (Fig. 3 E), especially in the negative regulation of immune system processes, angiogenesis, regulation of intracellular protein transport and other biological processes. Similarly, the KEGG analysis results revealed 50 significantly enriched pathways (Fig. 3 D). Among them, systemic lupus erythematosus, antigen processing and presentation, graft-versus-host disease, and Staphylococcus aureus infection were the most significant pathways. Through combined analysis of m6A-seq and RNA-seq data, 9281 genes not only presented differential m6A methylation but also presented altered mRNA expression levels. However, not all of these genes were significant. After screening with FDR  1 (for the same genes, the RNA-modified differentiation ratio of the largest peak), 6 genes with significant differences in both the m6A methylation level and the mRNA expression level were finally obtained, including NOV , ADAM23 , TET1 , NUF2 , ZNF92 , and ARHGAP11A . The detailed information is shown in Table 4 . Among these 6 genes, 3 had m6A hypermethylation and downregulated mRNA expression, 2 had m6A hypomethylation and upregulated mRNA expression, and 1 had m6A hypomethylation and downregulated mRNA expression. The intersection of m6A methylation and mRNA expression for each gene is shown via a Venn diagram (Fig. 4 C). Genes with differences in m6A methylation and mRNA expression. Diff = Different, FC = Fold change. Correlational analyses between RNA sequencing and MeRIP sequencing. (A) Four-quadrant plot of m6A modification-related gene expression in the SPID group and the control group; (B) cluster plot of differentially expressed genes in the SPID group and the control group; (C) Venn diagram of m6A modification-related gene expression in the SPID group and the control group; (D) cumulative differential mRNA abundance. In the Venn diagram (C), hypermethylated genes are shown in blue, hypomethylated genes are shown in green, upregulated genes are shown in red, and downregulated genes are shown in yellow. MeRIP = methylated RNA immunoprecipitation, SPID = sequelae of pelvic inflammatory disease. As confirmed by the cumulative frequency map, there was a correlation between the differential m6A methylation-modified genes and the mRNA expression levels (Fig. 4 D), indicating that the differential m6A methylation-modified transcripts presented different mRNA expression levels. The top 3 mRNAs ( NUF2 , TET1 , and ADAM23 ) were selected on the basis of the fold difference in m6A methylation through RNA sequencing (RNA-seq) and MeRIP-seq association analysis. To validate the sequencing results via MeRIP-qPCR, 5 SPID patient samples and 5 healthy control samples were collected to verify gene expression. The relative expression levels of m6A modifications in each gene were detected (Fig. 5 A–C). There was a significant difference in the relative m6A modification level of each gene between the SPID group and the control group ( P  < .01). Compared with the control group, the SPID group presented increased m6A methylation levels of the TET1 and ADAM23 genes, with TET1 showing greater upregulation. In contrast, the m6A methylation level of NUF2 was downregulated, which was consistent with previous sequencing results. PCR was performed to further confirm the differences in the expression of these 3 genes, and the results are shown in Figure 5 D–F. The expression of the 3 genes was opposite to their respective methylation levels, which was consistent with previous experiments. Selection and clinical validation of m6A-modified target genes. (A–C) Relative m6A levels of verified genes. (D–F) Relative quantity of verified genes; ** P  < .01, compared with the control group. In all panels, the control group is represented by black bars, and the SPID group is represented by gray bars. SPID = sequelae of pelvic inflammatory disease.

m6A-seq and RNA-seq libraries were constructed, and the changes in m6A methylation levels and mRNA expression levels in the peripheral blood of SPID patients were investigated through methylation-RNA immunoprecipitation combined with high-throughput sequencing; the results were analyzed via bioinformatics. The results revealed 12,204 and 14,065 m6A methylation-modifying genes in the control group and the SPID group, respectively. A total of 205 differentially expressed m6A modifier genes were further analyzed, including 65 hypermethylated genes and 140 demethylated genes. Compared with the control group, the SPID group presented significantly different levels of m6A methylation in the peripheral blood. However, there was no significant difference in the distribution of m6A methylations between the SPID group and the control group, and the m6A modifications of more than half of the genes were distributed in the 3’ untranslated region or the adjacent 3’ CDS. Consistent results have been reported in the research of Dominissini Dan et al. [ 21 ] They reported that m6A methylation sites preferentially occur in 2 distinct regions: around the stop codon and within the long internal exon; this distribution pattern is highly conserved between humans and mice. The results of the motif analysis revealed that the m6A methylation sites were mainly present in the RRACH motif, which was consistent with the classical m6A methylation sites. To further clarify the molecular mechanism of these differential genes in SPID, GO and KEGG enrichment analyses were performed on differential m6A methylation-modified mRNAs and differentially expressed mRNAs, respectively. First, different m6A methylation genes are involved primarily in neutrophil activation, neutrophil-mediated immunity, and neutrophil activation, which are involved in the immune response. Moreover, the KEGG analysis results revealed 180 significantly enriched pathways, especially those related to spliceosome and RNA transport. According to the results of the GO and KEGG enrichment analyses, m6A methylation may mediate the occurrence and development of SPID by affecting the inflammatory immune response and mRNA translation. GO and KEGG enrichment analyses of the differentially expressed mRNAs revealed that the differentially expressed mRNAs associated with SPID were involved mainly in the negative regulation of immune system processes, angiogenesis, the regulation of intracellular protein transport, and other pathways (such as systemic lupus erythematosus, antigen processing and presentation, graft-versus-host disease, and S aureus infection). Through comprehensive analysis of m6A-seq and RNA-seq data, 6 genes ( NUF2 , ARHGAP11A , TET1 , ZNF92 , ADAM23 , NOV ) whose mRNA expression levels were synchronously different when m6A methylation levels were different were ultimately selected. Among these 6 genes, 3 genes ( NUF2 , TET1 , ADAM23 ) were selected for validation in another 5 healthy individuals and 5 SPID patients on the basis of their molecular biological functions and fold changes in m6A methylation. According to the MeRIP-qPCR results, the m6A-methylated mRNA levels of TET1 and ADAM23 were upregulated in the serum of SPID patients, with TET1 showing more significant upregulation. The m6A-methylated mRNA level of NUF2 was downregulated, which was consistent with our sequencing results. Therefore, NUF2 , TET1 , and ADAM23 may play crucial roles in the occurrence and development of SPID, and they may be important targets for SPID treatment. However, the specific molecular mechanism by which m6A methylation of these genes regulates SPID is still unclear and needs further exploration. NUF2 is a part of the linker between kinetochores and tubulin subunits of the spindle, which is closely related to cell apoptosis and proliferation. [ 22 ] Disruption of NUF2 can induce defects in kinetochore attachment and the spindle checkpoint and may eventually induce cell death in mitotic cells. [ 23 ] NUF2 has been confirmed by several studies [ 24 – 26 ] to be significantly expressed in multiple tumor tissues, such as hepatocellular carcinoma, lung adenocarcinoma, and cervical cancer. Notably, our KEGG analysis of differentially expressed m6A methylation genes revealed that mTOR signaling pathways closely related to cellular metabolic processes (such as immunity and autophagy) were enriched. Studies [ 27 ] have shown through GEPIA, GSEA, and GSVA that NUF2 knockdown can reduce the expression of PIK3CA, PIK3CB, p-mTOR, and p-AKT, indicating that NUF2 is involved in the occurrence of ovarian cancer through the PI3K/AKT/mTOR signaling pathway. Moreover, Tricarico R et al reported that the NUF2 gene expression level is significantly negatively correlated with CD4 + T cells and B cells in intestinal tumors. [ 28 ] CD4 + T cells are the main “regulators” of the immune system, and CD4 + T cells can promote B cells to play an immune role. Our sequencing and PCR results demonstrated that the methylation level of NUF2 was downregulated and that the mRNA expression level was upregulated in SPID patients, which may lead to the activation of the PI3K/AKT/mTOR signaling pathway, thereby causing immune disorders in CD4 + T and B cells. As a DNA demethylase, TET1 plays an important role in the process of active DNA demethylation. Tricarico et al [ 29 ] reported that TET1 -mediated DNA demethylation inhibits the immune response in the intestine. However, the m6A methylation of TET1 mRNA has not been adequately studied. This study revealed that the m6A methylation level of TET1 in the serum of SPID patients was significantly different from that in the serum of healthy people, and the subsequent validation results increased the confidence of the sequencing data. Therefore, TET1 may not only play an important role in DNA demethylation but also play an irreplaceable role in m6A methylation of mRNAs. ADAM23 is a member of the disintegrin and metalloproteinase (ADAM) family. Hypermethylation of the promoter region of ADAM23 can silence ADAM23 expression. It has been reported [ 30 – 32 ] that hypermethylation in the promoter region of ADAM23 induces cell–cell and cell–matrix adhesion in various cancer tissues (such as gastric cancer, non-small cell lung cancer, and colorectal cancer). Verbisck N V reported [ 33 ] that this adhesive effect may arise from ADAM23 negatively regulating αvβ3 integrin activation. Our sequencing results revealed that the ADAM23 gene was hypermethylated and expressed at low levels in SPID patients. Therefore, we hypothesized that the pelvic tissue adhesion of SPID may be attributed to the hypermethylation of the promoter region of ADAM23 , thus inducing the activation of adhesive interactions. SPID is a common inflammatory disease among women of childbearing age. Our study revealed that m6A methylation may affect disease occurrence and development by regulating cell metabolism, immune activation, and adhesion activation. Overall, the discovery of m6A methylation-modifying genes may provide a new perspective for developing therapeutic drugs for SPID. However, it is worth noting that there are still challenges in the clinical application of m6A methylation modifier genes in SPID. The regulation of m6A methylation may provide a target for SPID treatment. However, this requires further investigation. Moreover, clarifying the complex regulatory mechanism of SPID related to drug m6A methylation may offer new strategies for SPID therapy.