Transcriptome-wide N6-methyladenosine modification profiling of long non-coding RNAs in patients with recurrent implantation failure.

We performed genome-wide profiling of m 6 A-modified lncRNAs in three control participants and three patients with RIF group. We have submitted the data to the Gene Expression Omnibus (GEO accession number GSE205398 ). We identified 10,504 and 10,512 m 6 A peaks in the RIF and control groups, respectively. Of these peaks, 8761 were common to both the RIF and control groups, 1743 were unique to the RIF group and 1751 were unique to the control group (Fig. 1 A). As shown in Fig. 1 B, the 10,504 m 6 A peaks in the RIF group were mapped to 7766 lncRNAs, and the 10,512 m 6 A peaks in the control group were mapped to 7720 lncRNAs in the control group. Of these lncRNAs, 6782 were common to the RIF and control groups, 984 were unique to the RIF group and 938 were unique to the control group. Among the 8761 common peaks, there were 1868 significantly differentially methylated peaks (fold change ≥ 2 and P  ≤ 0.0001). The 1868 differentially methylated m 6 A peaks were located across 623 lncRNAs. Among these lncRNAs, 158 were hypomethylated and 465 were hypermethylated (Fig. 1 C). All methylated lncRNAs were clustered (fold-change ≥ 2.0 and P  ≤ 0.00001), indicating differential expression between the RIF and control groups (Fig. 1 D). We searched for motifs that were enriched in the 50-bp area around m 6 A peak by DREME software [ 33 ]. As shown in Fig. 1 E, the identified m 6 A peaks contain a conserved RRACH motif, where R represents purine, A is m 6 A and H is a non-guanine base). This finding suggests the existence of m 6 A methylation mechanism. Fig. 1 The characteristics of m 6 A methylation in RIF. A A Venn diagram of m 6 A modification sites identified in lncRNAs from the RIF group (green) and the control group (blue). B A Venn diagram of m 6 A-modified lncRNAs from the RIF group (green) and the control group (blue). C The general number of differentially methylated peaks and associated lncRNAs; orange and blue represent hypermethylated and hypomethylated lncRNAs, respectively. D Cluster analysis of the m 6 A-modified lncRNA genes in the RIF and control groups. Red and green represent up- and downregulated peaks, respectively (fold-change ≥ 2.0, P  < 0.00001). E The top five m 6 A motifs enriched in the RIF and control groups The characteristics of m 6 A methylation in RIF. A A Venn diagram of m 6 A modification sites identified in lncRNAs from the RIF group (green) and the control group (blue). B A Venn diagram of m 6 A-modified lncRNAs from the RIF group (green) and the control group (blue). C The general number of differentially methylated peaks and associated lncRNAs; orange and blue represent hypermethylated and hypomethylated lncRNAs, respectively. D Cluster analysis of the m 6 A-modified lncRNA genes in the RIF and control groups. Red and green represent up- and downregulated peaks, respectively (fold-change ≥ 2.0, P  < 0.00001). E The top five m 6 A motifs enriched in the RIF and control groups As shown in Fig. 2 A, compared with the control group, there were 1443 significantly upregulated m 6 A peaks and 425 significantly downregulated m 6 A peaks in the RIF group (fold change ≥ 2 and P  ≤ 0.0001). They are listed in Supplementary Table 1. The top 10 hypermethylated lncRNAs and hypomethylated lncRNAs are listed in Tables 2 and 3 , respectively. We divided the m 6 A-modified lncRNAs into six categories based on the positional relationship between the m 6 A methylated lncRNA and mRNA: bidirectional, exon sense overlap, intron sense overlap, intron antisense, natural antisense and intergenic [ 34 ]. We found that 51.02% of the differentially m 6 A-methylated lncRNAs were in the exon sense overlap group (Fig. 2 B). We found that the majority of lncRNAs had one m 6 A peak. We explored the global role of m 6 A modification on the regulation of gene expression, especially the distribution pattern on different chromosomes (Fig. 2 D). The differentially m 6 A-methylated peaks could be mapped to all chromosomes, except the Y chromosome; chromosomes 1, 16 and 19 were particularly well represented. These results suggest that chromosomal specificity regarding m 6 A modification could regulate the complexity of gene expression and contribute to RIF. Table 2 Top ten hypermethylated lncRNAs Chromosome txStart txEnd transcript_id GeneName Foldchange chr2 3688901 3689022 NR_045659 COLEC11 96.7 chr2 219256098 219256280 ENST00000490872 SLC11A1 93.8 chr17 18394001 18394380 ENST00000584127 LGALS9C 75.5 chr17 72599881 72600260 ENST00000569279 CTD-2006K23.1 75.4 chr12 10365056 10365220 ENST00000539408 GABARAPL1 73.3 chr19 11908261 11908500 ENST00000589227 CTC-499B15.7 73.1 chr3 183240361 183240924 ENST00000487643 KLHL6 72.3 chr2 233215941 233216019 ENST00000427961 ECEL1P3 71.2 chr19 47235821 47235957 ENST00000600716 CTB-174O21.2 67.9 chr2 166933597 166933662 ENST00000599041 AC010127.3 67.0 Table 3 Top ten hypomethylated lncRNAs Chromosome txStart txEnd transcript_id GeneName Foldchange chr22 17155996 17156936 ENST00000338526 ANKRD62P1-PARP4P3 169.3 chr14 90314776 90314840 ENST00000550103 EFCAB11 140.1 chr9 80864161 80864460 ENST00000536374 CEP78 112.9 chr7 101032148 101032361 ENST00000413033 RP5-1106H14.1 108 chr9 21395981 21396345 ENST00000569618 RP11-354P17.15 102.7 chr3 77679221 77679620 ENST00000470802 ROBO2 100.1 chr8 104029261 104029817 ENST00000518264 NPM1P52 97 chr1 71927243 71927500 ENST00000587066 ZRANB2-AS2 89.8 chr13 78233547 78233560 ENST00000450718 SPTLC1P5 87.6 chr6 27235891 27236080 ENST00000604325 XXbac-BPGBPG24O18.1 85 Top ten hypermethylated lncRNAs Top ten hypomethylated lncRNAs Fig. 2 Distribution of differentially methylated m 6 A sites. A The volcano plot indicates the distribution of differential m 6 A peaks the in RIF. Red and blue represent up- and downregulated peaks, respectively (fold-change ≥ 2.0, P  < 0.00001 with Fisher’ s exact test). B The positional relationship between differentially m 6 A methylated lncRNAs and mRNAs. Orange, DodgerBlue, olive, light grey, blue and green represent exon sense-overlap, intergenic, intronic antisense, natural antisense, intron sense-overlap and bidirectional, respectively. C The distribution of altered m 6 A peaks per lncRNA. Blue and orange represent down- and upregulated peaks, respectively. D The chromosomal distribution of all differentially methylated sites within lncRNAs. Blue and orange represent down- and upregulated peaks, respectively Distribution of differentially methylated m 6 A sites. A The volcano plot indicates the distribution of differential m 6 A peaks the in RIF. Red and blue represent up- and downregulated peaks, respectively (fold-change ≥ 2.0, P  < 0.00001 with Fisher’ s exact test). B The positional relationship between differentially m 6 A methylated lncRNAs and mRNAs. Orange, DodgerBlue, olive, light grey, blue and green represent exon sense-overlap, intergenic, intronic antisense, natural antisense, intron sense-overlap and bidirectional, respectively. C The distribution of altered m 6 A peaks per lncRNA. Blue and orange represent down- and upregulated peaks, respectively. D The chromosomal distribution of all differentially methylated sites within lncRNAs. Blue and orange represent down- and upregulated peaks, respectively We explored the physiological and pathological significance of m 6 A modification in patients with RIF by performing GO and KEGG pathway analyses for the genes associated with differentially methylated lncRNA. For GO analysis, we determined the top 10 significantly enriched genes associated with hyper- and hypomethylated lncRNAs related to biological processes, cellular components and molecular functions. The genes associated with hypermethylated lncRNAs are significantly enriched in neuron projection, regulation of GTPase activity and nucleoside-triphosphatase regulator activity (Fig. 3 A). The genes associated with hypomethylated lncRNAs are significantly enriched in regulation of chromosome segregation, the Golgi stack, and peptidase activity (Fig. 3 B). For KEGG pathway analysis, we found that the genes associated with hypermethylated lncRNAs in patients with RIF are significantly associated with choline metabolism; the PLD signalling pathway; and valine, leucine and isoleucine degradation (Fig. 3 C). The genes associated with hypomethylated lncRNAs are involved in the p53 signalling pathway, transcriptional misregulation in cancer and natural killer (NK) cell–mediated cytotoxicity (Fig. 3 D). We found that the PLD (Fig. 4 A) and p53 (Fig. 4 B) signalling pathways [ 35 ] are the most significant enriched pathways. Fig. 3 A and B The top 10 GO terms for genes associated with hyper- and hypomethylated lncRNAs, respectively. Orange indicates biological process, green indicates molecular function and blue indicates cellular component. C and D The top 5 KEGG pathways of genes associated with hyper- and hypomethylated lncRNAs, respectively Fig. 4 Diagrams of the ( A ) phospholipase D and ( B ) p53 signalling pathways. These images have been reproduced from a previous publication [ 35 ] according to an open access licence A and B The top 10 GO terms for genes associated with hyper- and hypomethylated lncRNAs, respectively. Orange indicates biological process, green indicates molecular function and blue indicates cellular component. C and D The top 5 KEGG pathways of genes associated with hyper- and hypomethylated lncRNAs, respectively Diagrams of the ( A ) phospholipase D and ( B ) p53 signalling pathways. These images have been reproduced from a previous publication [ 35 ] according to an open access licence To explore the mRNAs regulated by lncRNAs, we screened 5 of 454 hypermethylated lncRNAs with a fold-change > 70 and 5 of 151 hypomethylated lncRNAs with a fold-change > 60 (Table 4 ). We constructed a ceRNA network based on lncRNA–miRNA–mRNA association analysis (Fig. 5 ). The network consists of the top five miRNAs (defined by exact probability from lncRNASNP2) combined with a screened lncRNA and the top five mRNAs (defined by the cumulative weighted context +  + score) bound to the miRNAs, including 10 lncRNAs, 49 miRNAs and 232 mRNAs (Supplementary Table 2). From this ceRNA network, it is clear that lncRNAs regulate miRNAs and mRNAs. In a previous study, the authors showed that lncRNAs could act as an miRNA sponge to interact with miR-548 s and increase ALDH1A3 mRNA expression. This action promotes glucose metabolism and cell proliferation in hepatocellular carcinoma [ 36 ]. We speculate that LINC01152 could act as a sponge for hsa-miR-6801-3p to affect LIF expression. Table 4 LncRNAs screened for lncRNA-miRNA-mRNA analysis transcript_id GeneName Foldchange Regulation NR_045659 COLEC11 96.7 hypermethylate ENST00000490872 SLC11A1 93.8 hypermethylate ENST00000539408 GABARAPL1 73.3 hypermethylate ENST00000487643 KLHL6 72.3 hypermethylate ENST00000599041 AC010127.3 67.0 hypermethylate ENST00000550103 EFCAB11 140.1 hypomethylate ENST00000536374 CEP78 112.9 hypomethylate ENST00000413033 RP5-1106H14.1 108 hypomethylate ENST00000569811 GOLGA8B 78.1 hypomethylate ENST00000538693 LINC01152 73.5 hypomethylate Fig. 5 The lncRNA–miRNA–mRNA regulatory network in RIF. The pink circular nodes represent upregulated lncRNAs, the purple circular nodes represent downregulated lncRNAs, the yellow triangular nodes represent miRNAs and the blue circular nodes represent mRNAs LncRNAs screened for lncRNA-miRNA-mRNA analysis The lncRNA–miRNA–mRNA regulatory network in RIF. The pink circular nodes represent upregulated lncRNAs, the purple circular nodes represent downregulated lncRNAs, the yellow triangular nodes represent miRNAs and the blue circular nodes represent mRNAs We evaluated the association between the expression of lncRNA and m 6 A methylation by conducting a conjoint analysis of the m 6 A-seq and RNA-Seq data. We identified 72 lncRNAs (Supplementary Table 3) with differential m 6 A methylation and expression levels (Fig.  6 C). We divided these genes into four categories: 27 hypermethylated and self-regulated lncRNAs, 19 hypermethylated and self-downregulated lncRNAs, 33 hypomethylated and self-downregulated lncRNAs and 1 hypomethylated and self-upregulated lncRNA. Among the 72 lncRNA in Supplementary Table 3, we randomly selected two (ENST00000416361 and ENST00000471664). We determined their m 6 A methylation abundance using the Integrative Genomics Viewer (IGV) (Fig.  6 A, B) and validated the findings with MeRIP-qPCR (Fig.  7 A, B). Finally, RT-qPCR revealed that the relative mRNA expression of ENST00000416361 and ENST00000471664 in endometrial tissue from the RIF group and Control group. was consistent with the sequencing results (Fig.  7 C, D). Fig. 6 Validation of the m 6 A methylation abundances of lncRNAs. IGV plots of ( A ) ENST00000416361 and ( B ) ENST00000471664 m 6 A methylation and expression abundances in the RIF and control groups. Each plot shows the MeRIP-qPCR primer binding regions. At the bottom of the IGV plot, blue and red represent primer positions and the middle region is the pcr product region. IP refers to the methylation level, and input refers to the expression level. C The four-quadrant plot shows the inclusion of differentially methylated DEGs(Differentially Expressed Genes) for m 6 A peaks Fig. 7 The m 6 A methylation abundances of ( A ) ENST00000416361 and ( B ) ENST00000471664 was validated in endometrial tissue with MeRIP-qPCR. The relative mRNA expression of ( C ) ENST00000416361 and ( D ) ENST00000471664 in endometrial tissue from the RIF group and Control group was determined with RT-qPCR. Statistical significance was determined using Student's t -test (* P  < 0.05, ** P  < 0.01, *** P  < 0.001, **** p  < 0.0001) Validation of the m 6 A methylation abundances of lncRNAs. IGV plots of ( A ) ENST00000416361 and ( B ) ENST00000471664 m 6 A methylation and expression abundances in the RIF and control groups. Each plot shows the MeRIP-qPCR primer binding regions. At the bottom of the IGV plot, blue and red represent primer positions and the middle region is the pcr product region. IP refers to the methylation level, and input refers to the expression level. C The four-quadrant plot shows the inclusion of differentially methylated DEGs(Differentially Expressed Genes) for m 6 A peaks The m 6 A methylation abundances of ( A ) ENST00000416361 and ( B ) ENST00000471664 was validated in endometrial tissue with MeRIP-qPCR. The relative mRNA expression of ( C ) ENST00000416361 and ( D ) ENST00000471664 in endometrial tissue from the RIF group and Control group was determined with RT-qPCR. Statistical significance was determined using Student's t -test (* P  < 0.05, ** P  < 0.01, *** P  < 0.001, **** p  < 0.0001)

This study was approved by the Ethics Committee of First Hospital of Lanzhou University (No: LDYYLL2019-45). The participants were recruited from the First Hospital of Lanzhou University. The inclusion criteria were (1) < 40 years old; (2) a regular menstrual cycle (25–35 days); and (3) normal basal serum sex hormone levels, namely < 10 mIU/mL for follicle-stimulating hormone (FSH), < 10 mIU/mL for luteinising hormone (LH) and < 50 pg/mL for oestradiol (ES), measured on days 2–3 of the menstrual cycle. The exclusion criteria were: (1) chromosomal anomaly based on a peripheral blood test; (2) tested positive for anticardiolipin antibody or lupus anticoagulant; (3) polycystic ovary syndrome (PCOS); (4) uterine anomalies (congenital uterine anomaly, fibroid, polyps and intrauterine adhesions); (5) abnormal blood glucose level or thyroid dysfunction function test; and/or (6) use steroid hormone in the preceding 2 months. According to the RIF definition, women with RIF had failed to achieve clinical pregnancy after at least four high-quality embryos had been transferred for a minimum of three cycles [ 1 ]. The control participants, who became pregnant during the cycle after sampling, underwent IVF/ICSI-ET due to fallopian tube or male factors. Six women were recruited for this study: three for the RIF group and three for the control group. Endometrial samples were obtained on day LH + 7 via pipe suction curettage, placed in cryopreservation tubes and stored at − 80ºC. These samples were from the middle secretory phase, as confirmed by histology. Total RNA was isolated from the endometrial samples with the TRIzol Reagent (Invitrogen, Gibco-BRL, Bethesda, MD, USA) according to the manufacturer’s instructions. A NanoDropND-1000 was utilised to determine the RNA concentration; the 260/280 ratio was 1.8–2.1. The quality of the RNA was assessed based on the 18S/28S ribosomal RNA (rRNA) band intensity in 1% agarose gels containing ethidium bromide. In addition, following the immunoprecipitation, a 2100 Bioanalyzer (Agilent) was used to assess the RNA quality. The m 6 A-Seq service was provided by CloudSeq Inc. (Shanghai, China). In accordance with manufacturer's instructions, total RNA was immunoprecipitated using the m 6 A-IP Kit (GenSeq Inc.). In brief, RNA was randomly fragmented to about 200 nt using the RNA Fragmentation Reagent, and Protein A/G beads were rotated for 1 h to couple protein to antibody. Bead-linked antibodies were incubated with RNA fragments for 4 h at 4 °C. After incubation, the RNA/antibody complexes were washed several times. Then, captured RNA was eluted from the complexes and purified. The RNA libraries for the IP and input samples were constructed using the Low Input Whole RNA Library Prep Kit (GenSeq), following the manufacturer’s instructions. After evaluating the quality of the libraries with a 2100 bioanalyzer (Agilent), sequencing was carried out on a NovaSeq instrument (Illumina). High-throughput RNA-seq was performed by Cloud-Seq Biotech. Briefly, rRNA was removed from total RNA with the rRNA Removal Kit (GenSeq, Inc.). Then, the rRNA-depleted samples were used to construct libraries with the Low Input RNA Library Prep Kit (GenSeq, Inc.) according to the manufacturer’s instructions. The library quality and quantity were analysed using a 2100 BioAnalyzer (Agilent). The libraries were sequenced with a Novaseq 6000 sequencer (Illumina) with 150 base pair (bp) paired-end reads. Three differentially expressed genes with differentially methylated m 6 A peaks were used for MeRIP-qPCR validation. RNA was extracted and quantified as described above. A small proportion of fragmented RNA was used as the input control. The remaining RNA was incubated with magnetic beads coupled with an m 6 A antibody. Then, m 6 A RNA was eluted from the magnetic beads, according to the manufacturer’s instructions. The m 6 A IP samples and input samples were collected for reverse transcription (RT)-qPCR. The primers for MeRIP-qPCR are shown in Table  1 . Table 1 Primer sequences used for RT-qPCR transcript_id Primer sequences (5’-3’) Forward Reverse ENST00000416361 GGGGAAATGTGGGGAAAAGA GATCAACAGCATCCCAAGGC ENST00000471664 AGTGACTCTTCTTGGACCAGG ATTGTGCCTCTCCAATCTGC Primer sequences used for RT-qPCR Briefly, the Novaseq 6000 reads were subjected to quality control with Q30 [ 23 ]. After 3′ adaptor-trimming and removing low-quality reads with the cutadapt software (v1.9.3), the high-quality trimmed reads were used to analyse lncRNAs and mRNAs. They were aligned to the reference genome (UCSC hg19) with the hisat2 software (v2.0.4) [ 24 ]. Then, the HTSeq software (v0.9.1) [ 25 ] was used to obtain the gene level (mRNA) and transcript level (lncRNA) raw counts, and edgeR (v3.16.5) [ 26 ] was used for normalisation. Differentially expressed mRNAs were identified based on the p-value and fold change. Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed based on the differentially expressed mRNAs. Methylated sites on RNAs were identified with the MACS software [ 27 ]. Differentially methylated sites were identified using diffReps [ 28 ]. The lncRNA target genes were predicted, and GO and KEGG analyses were performed based on these target genes [ 29 ].The miRNA-target gene prediction software programs lncRNASNP2 [ 30 ] and TargetScan [ 31 ] were used to predict miRNAs and mRNAs, which were combined for the screened lncRNAs. The ceRNA network was plotted using Cytoscape (v3.9.1) [ 32 ].

The success rate of assisted reproductive technology exceeds 50%, but some patients fail after recurrent implantation, which leads to large financial and psychological burdens. Recurrent implantation failure (RIF) is usually defined as a failure to achieve clinical pregnancy after at least four high-quality embryos are transferred in women under the age of 40 years within at least three fresh or frozen cycles [ 1 ]. Low endometrial receptivity is considered to be the main cause of RIF [ 2 ]. During a normal menstrual cycle, 6–8 days after ovulation is considered to be the window of implantation when endometrial receptivity is highest [ 3 ]. At this time, pinopodes, which are membranous protuberances, can be observed under a scanning electron microscope; they are considered to be a morphological marker of endometrial receptivity [ 4 ]. There are several molecular change that are associated with endometrial receptivity, have been detected in women with RIF, including mucin 1, integrin β3, homeobox A10 and leukaemia inhibitor factor ( LIF ) [ 5 – 7 ]. However, the precise aetiology and pathogenesis of RIF remains unclear. N 6 ‐methyladenosine (m 6 A) was first discovered in 1974. It is the most common RNA modification in eukaryotes [ 8 ]. Researchers have confirmed that the formation, removal and other functions of m 6 A are accomplished by three types of regulatory proteins [ 9 ]. m 6 A modification is catalysed by the methyltransferase ‘writer’ complex that comprises methyltransferase-like 3 ( METTL3 ); methyltransferase-like 14 ( METTL14 ); Wilms’ tumour 1-associating protein ( WTAP ); and other proteins with methyltransferase capability, including RNA-binding motif protein 15 ( RBM15 ) and zinc finger CCCH-type containing 13 ( ZC3H13 ) [ 10 ]. The demethylase complex includes ALKB homolog 5 ( ALKBH5 ) and fat mass and obesity associated gene ( FTO ); it removes m 6 A, thus acting as an ‘eraser’. Finally, m 6 A-modified RNA can be recognised and regulated by m 6 A-binding protein complexes, including YTH domain family proteins 1–3 ( YTHDF1–3 ), and other proteins, which acting as ‘readers’ [ 11 ]. A 2021 study confirmed downregulated METTL16 and WTAP messenger RNA (mRNA) levels and upregulated ALKBH5 and IGF2BP2 mRNA levels in the endometrium of patients with infertility [ 12 ]. Of note, it has been proved that global m 6 A methylation and METTL3 expression are significantly increased in the endometrial tissues from women with RIF [ 13 ]. This evidence indicates that m 6 A modification is critical in the endometrium of patients with RIF. Long noncoding RNAs (lncRNAs) are a class of noncoding RNAs that are longer than 200 nucleotides (nt). For a long time, they were regarded as transcriptional noise because of their low expression and non-protein coding features. Recently, lncRNAs have been indicated to play an important role in embryo implantation. For example, TUNAR might participate in embryo implantation by regulating the adhesion of blastocysts to endometrial epithelium and the proliferation and decidualisation of endometrial stromal cells (ESCs) [ 14 ]. CECR3 and five other lncRNAs together with their competing endogenous RNA (ceRNA) sub-networks, might be involved in immunological activity, growth factor binding, vascular proliferation, apoptosis, steroid biosynthesis in the uterus and preparing endometrium for embryo implantation [ 15 ]. Therefore, the regulation mediated by lncRNAs may affect embryo implantation. In addition, researchers have demonstrated that lncRNAs present abundant m 6 A methylation [ 16 , 17 ]. The m 6 A-modified lncRNA profile has been established in many diseases, including coronary artery disease [ 18 ], gastric cancer [ 19 ], ameloblastoma [ 20 ] and systemic lupus erythematosus [ 21 ]. Recently, m 6 A-modified lncRNAs have been shown to inhibit human trophoblast cell proliferation and to induce miscarriage [ 22 ]. However, the m 6 A-modified lncRNA methylome in patients with RIF has not yet been determined. In this study, we used m 6 A-modified RNA immunoprecipitation sequencing (m 6 A-seq) to establish the m 6 A-methylation transcription profiles in RIF for the first time. We identified differentially methylated peaks and then constructed a ceRNA was constructed by using bioinformatics to reveal a regulatory relationship between lncRNAs, microRNAs (miRNAs) and mRNAs. Our findings provide new insights into the pathogenic mechanism of RIF.

RIF is a challenging clinical dilemma. m 6 A modifications play a vital role in reproductive diseases, but prior to our study, the m 6 A-modified lncRNA profile of RIF had not been characterised. We established the transcriptome-wide lncRNA m 6 A modification profile in RIF by using m 6 A-Seq. We identified 1868 differentially methylated lncRNA m 6 A peaks in endometrial tissue from patients with RIF, including 1443 significantly upregulated m 6 A peaks and 425 significantly downregulated m 6 A peaks. Furthermore, we explored differential expression patterns of differentially methylated lncRNAs as well as GO and KEGG pathway analyses to reveal the potential functions of differentially methylated transcripts. We also constructed a ceRNA network to reveal the regulatory relationships between lncRNAs, miRNAs and mRNAs. Our findings provides insights regarding the directions for the diagnosis and treatment of RIF in the future. m 6 A plays a crucial role in human fertility. As early as 2013, researchers identified that ALKBH5 , a genes that encodes an m 6 A ‘eraser’, impairs fertility by affecting apoptosis of spermatocytes in meiotic metaphase [ 37 ]. Subsequently, several m 6 A regulators have been discovered to play an important role in reproductive diseases, such as endometriosis [ 38 ], polycystic ovary syndrome [ 39 ] and primary ovarian insufficiency [ 40 ]. In one study, the authors used a colourimetric m 6 A quantification strategy to examine the m 6 A level in control subjects and patients with RIF [ 13 ]. They found a significantly increased global m 6 A level in endometrial tissues from patients with RIF compared with the controls. This finding consistent with our study (Fig. 1 D) and indicates that m 6 A is important for the pathogenesis of RIF. Non-coding RNAs are important in RIF. The HOX family members are critical regulators in endometrial decidualisation [ 41 ]. In patients with RIF, HOXA11-AS , a lncRNA in the HOX gene family, is elevated and there is impaired PKM2 splicing [ 42 ]. LINC02190 is upregulated in the endometrium of patients with RIF and decreases the adhesion rate of Ishikawa and JAR cells [ 43 ]. Given the abundant published evidence, the role of lncRNAs in RIF cannot be ignored. Consistently, our study has provided novel findings regarding m 6 A modified of lncRNA in RIF. As shown in Fig. 2 B, most differentially m 6 A-modified lncRNAs are exon sense overlap. In a previous study, the authors found that the majority of differentially m 6 A-modified lncRNAs in colorectal cancer are also exon sense overlap [ 44 ]. However, there have been different findings from other studies. Most differentially methylated m 6 A sites within lncRNAs are located in intergenic region in coronary artery disease [ 18 ] and gastric cancer [ 19 ]. Thus, the dominant category of m 6 A-modified lncRNAs can differ depending on the disease. Of note, exon sense overlap lncRNAs are relevant to cell proliferation, invasion and metastasis, and thus the differentially m 6 A-modified exon sense overlap mRNAs we identified may play in RIF. We found that the differentially methylated lncRNAs are related to several important biological pathways, especially choline metabolism, the PLD signalling pathway and amino acid metabolism. Choline metabolism is related to the PLD signalling pathway. PLD hydrolyses phosphatidylcholine into phosphatidic acid, a lipid signalling molecule, and choline [ 45 ], a micronutrient and a methyl donor [ 46 ]. Arzu Yurci et al. [ 47 ] reported lower choline levels in patients with RIF compared with individuals with normal fertility based on magnetic resonance spectroscopy. These results suggest that choline may participate in the pathogenesis of RIF through m 6 A methylation. Moreover, genes associated with hypermethylated lncRNAs are related to the metabolism of amino acids such as valine and histidine. Increasing evidence has demonstrated that amino acid metabolism is a key regulator in RIF. In previous studies, researchers found that valine is significantly downregulated in women with repeated implantation success compared with women with RIF [ 48 ]. In a recent systematic review, the authors found that alanine, aspartate and glutamate metabolism has the greatest impact on RIF [ 49 ]. Based on our results, we suggest that m 6 A methylation influences the outcome of embryo implantation via altered valine and histidine metabolism. A previous study confirmed that lncRNA hypomethylation can promote cancer cell cycle progression [ 50 ]. Consistently, a recent review suggested that m 6 A methylation can promote cancer progression [ 51 ]. According to our GO analysis, the hypomethylated lncRNAs are mostly related to the cell cycle, implying that changes in m 6 A methylation affect cell proliferation in patients with RIF. The genes associated with hypomethylated lncRNAs are mainly involved in immunological processes, including NK cell–mediated cytotoxicity and the chemokine signalling pathway (Fig. 3 D). The ability to kill other cells is an important effector mechanism of the immune system; NK cells are the major mediators of this activity [ 52 ]. NK cells are recruited and activated by ovarian hormones and play pivotal roles throughout pregnancy. Decidual natural killer (dNK) cells release chemokines that induce trophoblast invasion, tissue remodelling, embryonic development and placentation. NK cells can also mediate cytotoxicity and carry out immune defence in case of an in utero infection [ 53 – 55 ]. Researchers have suggested an association between RIF and abnormally elevated uterine NK cell levels and cytotoxicity [ 54 ]. It has been proposed that intrauterine administration of peripheral blood mononuclear cells (PBMCs) modulates the maternal immune response through a cascade of chemokines to favour implantation. This therapy could significantly improve the clinical pregnancy rate, a finding confirmed with a meta-analysis [ 56 , 57 ]. Based on our results, we hypothesise that m 6 A methylation associated with immunological processes influences the pregnancy outcome in patients with RIF. Genes associated with hypomethylated lncRNAs are also involved in the p53 signalling pathway. p53 is a well-known tumour suppressor, but little is known about its other functions. With the aid of the p53MH algorithm, researchers identified LIF as a p53-regulated gene [ 58 ]. p53 is crucial for embryonic implantation through upregulation of uterine LIF transcription [ 59 ]. Recently, it has been proposed that the absence of PARP-1 and PARP-2 increases p53 signalling and the population of senescent decidual cells [ 60 ]. In a case–control survey with 100 patients with RIF and 100 normal pregnancies, the researchers suggested that p53 gene polymorphisms could be a genetic predisposing factor for RIF [ 61 ]. Based on the results of our study, m 6 A might participate in the p53 signalling pathway and contribute to the pathogenic mechanism of RIF. We constructed a ceRNA network to analyse the regulatory relationship between lncRNAs, miRNAs and mRNAs (Fig. 4 ). The potential target genes of several lncRNAs in this ceRNA are related to reproductive disease. For example, LINC01152 could act as a sponge for hsa-miR-6801-3p to affect LIF expression (Fig. 4 ). LIF activates signal transducer and activator of transcription 3 ( STAT3 ) precursors via uterine LIF receptors, allowing successful blastocyst implantation [ 62 ]. In addition to LIF , PRX , CD69 and HLA-F are potential target genes in the ceRNA network. The expression of the PRX is significantly decreased in ectopic endometrium compared with eutopic or control endometrium [ 63 ]. In an endometrial immune assessment, the sub-fertile population exhibited increased CD69 activation [ 64 ]. HLA-F fluctuates during the menstrual cycle, with high levels during the implantation window. HLA-F protein expression correlates with the number of CD56-positive NK cells in the mid-secretory endometrium [ 65 ]. The functions of abovementioned genes in decidualisation and RIF require further investigation. However, it is clear that the network provides insights to further understand the role of lncRNAs in the occurrence and development of RIF. We verified the m 6 A methylation abundance of three lncRNAs (ENST00000416361 and ENST00000471664) by using MeRIP-qPCR. ENST00000416361 is upregulated in coronary artery disease and is related to inflammation and lipid metabolism [ 66 ]. ENST00000471664 plays a significant role in GLIS3 knockout mouse pancreatic islets [ 67 ]. However, these studies are limited to bioinformatics. There have been many studies on the involvement of m 6 A methylation in the expression of lncRNAs [ 68 , 69 ]. We suspect that this process might influence lncRNA expression in RIF, although this speculation requires further investigation. There are several limitations to our study. First, we only used a small number of samples. To make up for this deficiency, we used MeRIP-qPCR to ensure the accuracy of our m 6 A-seq results. Second, we performed a descriptive bioinformatics analysis. It is necessary to explore the relationship of m 6 A and lncRNAs in RIF in subsequent experiments.

In summary, we revealed the m 6 A lncRNA methylation landscape by using high-throughput sequencing;. This endeavour has provided new perspectives on the pivotal roles of epigenetic changes in RIF. We generated a ceRNA network that demonstrated a regulatory relationship between lncRNAs, miRNAs and mRNAs. Additional experiments should be carried out to verify the differential methylated lncRNAs in patients with RIF.

Supplementary Material 1. Supplementary Material 2. Supplementary Material 3. Supplementary Material 1. Supplementary Material 2. Supplementary Material 3.