Villous tissues from patients with spontaneous abortion and controls with induced abortion used in this research were obtained with written informed consent from all participants. All tissues were approved by the Medicine Ethics Committee of Shenzhen Second People’s Hospital (Approval number, 20210517001-FS01). The inclusion criteria were: (1) Patients who were clinically diagnosed with spontaneous abortion for the first time or healthy women who underwent voluntary induced abortion; (2) Age between 18 and 35; (3) Gestational ages between 6-8 weeks. The exclusion criteria were: (1) Patients with fetal chromosomal or congenital abnormalities; (2) Patients with abnormal uterine structures; (3) Patients with polycystic ovary syndrome, endometriosis, and thyroid disease; (4) Patients with vaginitis. Three spontaneous abortion samples and three induced abortion samples were obtained from elective terminations of apparently normal pregnancies. These samples were further used for MeRIP-seq. The villous tissue was rinsed in precooled normal saline three times, while the tissue with a diameter of 1 mm was cut with ophthalmic scissors. After the tissue was mixed in Trizol solution in a volume ratio of 1:10, it was ground to homogenate with a freeze grinder, and the left tissue was frozen in liquid nitrogen.

Total RNA from each sample was isolated using TRIzol reagent (Invitrogen) and fragmented into ∼100-nucleotide-long fragments by zinc acetate. Next, Affinity-purified anti-m⁶A polyclonal antibodies (Abcam) were used for immunoprecipitation to analyze approximately 300 μg of fragmented RNA. After stringent washing with a high-salt buffer (400 mM NaCl, 0.05% NP-40, 10 mM Tris–HCl), competitive buffer (150 mM NaCl, 0.05% NP-40, 10 mM Tris–HCl, 0.25 mg ml−1 mix of adenosine, uridine, guanosine and cytidine), high-detergent buffer (150 mM NaCl, 0.5% NP-40,10 mM Tris–HCl), and immunoprecipitation buffer (150 mM NaCl, 0.05% NP-40, 10 mM Tris–HCl), bound RNA was eluted by competition with 1 mg ml−1 N6-methyladenosine (Selleckchem) and used for library construction using the NEBNext Ultra RNA Library Prep Kit v2 for Illumina. After removal of ribosomal RNA using the Epicentre Ribo-zero rRNA Removal Kit (Epicentre), total RNA from each tissue was fragmented and a library was constructed using the NEBNext Ultra RNA Library Prep Kit v2 for Illumina (New England Biolabs) as input RNA. RNA-seq libraries of m⁶A antibody-enriched mRNAs and input mRNAs were prepared. Sequencing was carried out using an Illumina Hiseq 4000 platform according to the manufacturer’s instructions.

Trimmomatic (v.0.27) (Bolger et al., 2014) was used for quality control of paired-end sequencing data. Reads that mapped to rRNA and tRNA sequences (obtained from the UCSC gene annotation (hg38)) using bowtie2 (v.2.3.4) (Langmead and Salzberg, 2013) were discarded, and the remaining reads were aligned to the GRCh38 using hisat2-align (v.2.1.0) (Kim et al., 2015). Unique reads with high mapping quality were retained using Picard (v.2.16.0) and SAMtools (v.1.7.0). MACS2 (v.2.1.1) (Gaspar, 2018) was used to identify m⁶A peaks with the parameter ‘--nomodel’ and ‘-q 0.05’. ExomePeak2 package (v.1.2.0) (Meng et al., 2014) was used for the identification of differentially methylated peaks. The GFF annotation file was referred to determine the strand information of m⁶A peaks. The findMotifsGenome.pl Perl script from the Homer software suite was used for motif search with the “-mask -rna -len 6” parameters. Genomic locations were split into CDS, 5′UTR, 3′UTR, promoter (2 kb upstream and 100 bp downstream of the TSS), and intergenic regions. The distribution of m⁶A peaks on mRNA was analyzed using the R package Guitar (v.1.7.0) (Cui et al., 2016). Differentially expressed genes were identified using the RNA-seq data (the corresponding MeRIP-seq input library data) by the R package DEseq2 (v.1.32.0) (Love et al., 2014). The R package clusterProfiler (v.4.1.3) (Wu et al., 2021) was used to calculate the biological significance of differentially methylated genes and differentially expressed genes through the Gene Ontology (GO) database and the latest Kyoto Encyclopedia of Genes and Genomes (KEGG) database.

The STRING (v11.5) (https://string-db.org/cgi/input.pl) (Szklarczyk et al., 2021) has been widely applied to construct a protein-protein interaction (PPI) network. Based on those DEGs, the “Multiple proteins” option was selected. The minimum required interaction score was set as “high confidence (0.700)” and a PPI network was constructed. And then, the cytoHubba (Chin et al., 2014) was employed to identify hub genes. The eccentricity algorithm was selected and twenty top-ranked genes were chosen as hub genes. Finally, Cytoscape (v3.9.0) (Demchak et al., 2014) was used to visualize the hub gene network.

The t-test was used for comparing the statistical significance between two groups. For each analysis, p < 0.05 was considered as statistically significant.

The raw sequencing and processed data reported in this study have been deposited in the Sequence Read Archive (SRA) and Gene Expression Omnibus (GEO) database and are accessible at https://dataview.ncbi.nlm.nih.gov/object/PRJNA786693 and https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE193052, respectively.

Due to the trace of early villous tissue, we performed a micro MeRIP-seq analysis of villous tissues from three patients with spontaneous abortion (SA group) and three controls with induced abortion (normal group) in the first trimester. We detected a total of 18,568 m⁶A peaks in these two groups by R package exomePeak2. As shown in Figure 1A, compared with normal group, SA group had 2,159 significantly upregulated m⁶A peaks, which corresponded to transcripts of 2,087 genes, and 281 significantly downregulated m⁶A peaks, which corresponded to transcripts of 311 genes (| log2 (fold change) | > 0.585 and p < 0.05). The top 20 altered m⁶A peaks were listed in Table 1. Then, we investigated the distribution of m⁶A peaks in the SA and normal group and found that m⁶A peaks in the SA and normal group were primarily enriched in the coding sequence (CDS) near the stop codon and the whole CDS region, respectively (Figure 1B). However, m⁶A peaks in the SA group showed a distinct pattern from m⁶A peaks in the normal group with a relative decrease in the number of m⁶A peaks in the coding sequence (CDS) (42.08 vs. 45%) and 5′ untranslated region (5′UTR) (1.64 vs. 2%) and a relative increase in the 3′ untranslated region (3′UTR) (19.93 vs. 19%) (Figure 1C). By analyzing the distribution of m⁶A peaks per gene, we found that most genes only had one corresponding m⁶A peak (1,878/236 genes with upregulated and downregulated peaks, respectively) (Figure 1D). Furthermore, dysregulated m⁶A peaks were found in all chromosomes, except chrY, and were mainly found in chr1, chr2, chr3, chr4, chr5, chr6, chr7, chr17 (Figure 1E). Moreover, the m⁶A peaks were mainly characterized by AUGGAC and UGGACG motif (Figure 1F).

To investigate the biological significance of m⁶A modification in villous tissues of patients with spontaneous abortion and controls with induced abortion in the first trimester, we performed GO and KEGG pathway enrichment analyses of differentially methylated mRNAs. GO ontology was classified into three categories: biological process (BP), cellular component (CC), and molecular function (MF). The top five significantly enriched BPs, CCs, and MFs of genes with upregulated and downregulated m⁶A peaks were shown in Figures 2A,B, respectively. The results in Figure 2A indicated that GO terms such as the regulation of Wnt signaling pathway, Rho protein signal transduction, transcription coregulator activity, and transcription corepressor activity were significantly enriched, and GO terms such as regulation of RNA biosynthetic process, regulation of nucleic acid-templated transcription and DNA-binding transcription factor activity were significantly enriched in Figure 2B. For KEGG pathway enrichment analysis, we found that genes with upregulated m⁶A peaks in villous tissues of patients with spontaneous abortion in the first trimester were significantly associated with the lysine degradation, ubiquitin-mediated proteolysis, herpes simplex virus 1 infection, adherens junction, phosphatidylinositol signaling system, inositol phosphate metabolism, and Hippo signaling pathway - multiple species (Figure 2C). Genes with downregulated m⁶A peaks were significantly associated with herpes simplex virus 1 infection, glycosphingolipid biosynthesis-ganglio series, Hippo signaling pathway, tryptophan metabolism, and Wnt signaling pathway (Figure 2D). The enrichment of genes in the four major pathways is shown in Figure 3.

Through RNA-seq (meRIP-seq input library), we detected the transcriptome profiles of villous tissues from patients with spontaneous abortion and controls with induced abortion in the first trimester. R package DESeq2 was used to detect differentially expressed genes (DEGs). Compared with normal group, SA group had 254 significantly upregulated genes and 133 significantly downregulated genes (| log2 (fold change) | > 0.585 and p < 0.05; Figures 4A,B). The MA plot was visualized for these DEGs (Supplementary Figure S1). The top 20 DEGs are listed in Table 2. The top 5 significantly enriched BPs, CCs, and MFs of genes with upregulated and downregulated expressed genes and top 10 KEGG pathways were displayed in Supplementary Figure S2. Based on these DEGs, we construct a PPI network through the STRING database. The hub genes selected from the PPI network are visualized in Supplementary Figure S3. According to the eccentricity scores, we identified twenty hub genes with highest confidence scores from the network, and found that most of hub genes were related to immune response and embryonic development. Then, we conducted conjoint analysis of the MeRIP-seq and RNA-seq data and explored the relationship between differential m⁶As and host gene expression level. We found that the expression level of host gene of differential m⁶As was lower in both SA group and normal group compared with constitutive m⁶As (Supplementary Figure S4). Dividing all differentially methylated m⁶A peaks with all differentially expressed mRNAs into four groups (hyper-up, hyper-down, hypo-up, and hypo-down), we identified 34 hypermethylated m⁶A peaks in mRNAs that were significantly upregulated (28; hyper-up) or downregulated (6; hyper-down), while 2 hypomethylated m⁶A peaks in mRNAs that were significantly upregulated (1; hypo-up) or downregulated (1; hypo-down) (Figure 4C, Supplementary Table S1). Finally, we performed GO and KEGG pathway enrichment analysis to explore the biological function of those genes (35) with differentially methylated m⁶A peaks and differential expression. The top 5 significantly enriched BPs, CCs, and MFs indicated that these genes were mainly enriched in the ossification (BP category), nuclear outer membrane (CC category), and core promoter sequence-specific DNA binding (MF) (Figure 4D). However, no KEGG pathways were significantly enriched.