The Alteration of m6A Modification at the Transcriptome-Wide Level in Human Villi During Spontaneous Abortion in the First Trimester.

A growing number of studies have demonstrated that N6 methyladenine (m6A) acts as an important role in the pathogenesis of reproductive diseases. Therefore, it is essential to profile the genome-wide m6A modifications such as in spontaneous abortion. In this study, due to the trace of human villi during early pregnancy, we performed high-throughput sequencing in villous tissues from spontaneous abortion (SA group) and controls with induced abortion (normal group) in the first trimester. Based on meRIP-seq data, 18,568 m6A peaks were identified. These m6A peaks were mainly located in the coding region near the stop codon and were mainly characterized by AUGGAC and UGGACG motif. Compared with normal group, the SA group had 2,159 significantly upregulated m6A peaks and 281 downregulated m6A peaks. Biological function analyses revealed that differential m6A-modified genes were mainly involved in the Hippo and Wnt signaling pathways. Based on the conjoint analysis of meRIP-seq and RNA-seq data, we identified thirty-five genes with differentially methylated m6A peaks and synchronously differential expression. And these genes were mainly involved in the Wnt signaling pathway, phosphatase activity regulation, protein phosphatase inhibitor activity, and transcription inhibitor activity. This study is the first to profile the transcriptome-wide m6A methylome in spontaneous abortion during early pregnancy, which provide novel insights into the pathogenesis and treatment of spontaneous abortion in the first trimester.

Introduction

Spontaneous abortion (SA) is considered to be one of the most common and severe complications during early pregnancy, which affects 10–15% of pregnant women (Rossen et al., 2018). The etiology of SA is multifactorial, which mainly includes endocrine abnormalities, immune abnormalities, abnormal uterine anatomy, prethrombotic state, chromosome abnormality and infection factors. There may still be other unknown factors contributing to SA, so further investigation is needed. To date, the role of some epigenetic modifications (DNA methylation, histone modification, and non-coding RNA) in SA has been well identified (Liu et al., 2018; Wang et al., 2019; Chen et al., 2021). As the most abundant epigenetic modification of mRNA in eukaryotic cells, m6A modification affects the stability (Wang et al., 2014; Ke et al., 2017; Huang et al., 2020), translation (Meyer et al., 2015; Wang et al., 2015; Shi et al., 2017), and splicing of mRNA (Xiao et al., 2016; Pendleton et al., 2017). Although previous studies have also illuminated that m6A modification plays important roles in the regulation of immune function and inflammatory response, the relationship between m6A methylation and SA remains to be elucidated.

As known, the m6A modification includes three main components: 1) “writers”, the methyltransferase complex, such as METTL3 (methyltransferase-like 3), METTL14 (methyltransferase-like 14) (Wang et al., 2016), and METTL16 (methyltransferase-like 16) (Mendel et al., 2018); 2) ‘readers’, RNA binding proteins, including YTHDF1/2/3 (YTH-family proteins 1/2/3), and YTHDC1/2 (YTH domain containing proteins 1/2) (Xu et al., 2015; Gao et al., 2019), and IGF2BP1/2/3 (insulin-like growth factor 2 mRNA binding proteins 1/2/3) (Huang et al., 2020); 3) “erasers”, demethylases, including ALKBH5 (alkB homolog 5) (Zheng et al., 2013) and FTO (fat mass and obesity-associated protein) (Jia et al., 2011). Recently, studies have shown that YTHDF2 can regulate oocyte maturation in female mice and negatively regulate the JAK-STAT pathway to affect the development of mouse nervous system (Ivanova et al., 2017). Wang et al. found that FTO can inhibit the expression level of myogenin, and thus inhibit the differentiation of muscle cells. Li et al. found that ALKBH5 might inhibit the invasion of trophoblast cells in patients with recurrent pregnancy loss, inhibit trophoblast invasion and thus affect mRNA stability. Furthermore, m6A modification has also been confirmed to be associated with embryonic stem cell differentiation (Aguilo et al., 2015; Geula et al., 2015), hematopoietic system development (Zhang et al., 2017), myogenesis (Wang et al., 2017), and early embryonic development (Sui et al., 2020). Xia’s lab mapped the transcriptome-wide m6A profiles of major fetal tissues including human placenta (Xiao et al., 2019), but the m6A methylome in human villi during spontaneous abortion in the first trimester has not been characterized.

In this study, due to the trace of early villous tissue, we performed high-throughput sequencing to determine the transcriptome-wide m6A methylome in human villi from patients with spontaneous abortion and controls with induced abortion in the first trimester. Using MeRIP-seq data, we further identified differential m6A peaks in villous tissue based on the comparison of spontaneous abortion and induced abortion samples. Then, we identified differentially expressed genes using RNA-seq data. Finally, the conjoint analysis of MeRIP-seq and RNA-seq revealed some genes with differentially methylated m6A peaks and synchronously differential expression, which might provide an alternative strategy for the therapy and prevention of spontaneous abortion.

Materials and Methods

Samples Collection

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.

MeRIP Sequencing

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-m6A 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 m6A antibody-enriched mRNAs and input mRNAs were prepared. Sequencing was carried out using an Illumina Hiseq 4000 platform according to the manufacturer’s instructions.

Data Analysis

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 m6A 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 m6A 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 m6A 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 Construction of Hub Gene Network

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.

Statistical Analyses

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

Data Availability

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.

Results

Profiles of m6A Modification in Human Villi With Spontaneous Abortion in the First Trimester

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 m6A 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 m6A peaks, which corresponded to transcripts of 2,087 genes, and 281 significantly downregulated m6A peaks, which corresponded to transcripts of 311 genes (| log2 (fold change) | > 0.585 and p < 0.05). The top 20 altered m6A peaks were listed in Table 1. Then, we investigated the distribution of m6A peaks in the SA and normal group and found that m6A 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, m6A peaks in the SA group showed a distinct pattern from m6A peaks in the normal group with a relative decrease in the number of m6A 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 m6A peaks per gene, we found that most genes only had one corresponding m6A peak (1,878/236 genes with upregulated and downregulated peaks, respectively) (Figure 1D). Furthermore, dysregulated m6A 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 m6A peaks were mainly characterized by AUGGAC and UGGACG motif (Figure 1F).

FIGURE 1

Characteristics of m6A methylation in human villi during spontaneous abortion in the first trimester. (A) Volcano plots showing the significantly altered m6A peaks. (B) Accumulation of the region of average m6A peaks along with all transcripts in SA group and normal group. (C) Pie charts showing the distribution of m6A peaks in SA group and normal group. (D) The distribution of altered m6A peaks per gene. (E) The distributions of altered m6A peaks in all chromosomes. (F) The top two m6A motifs enriched from the altered m6A peaks.

TABLE 1

Top 20 altered m6A peaks in human villi during spontaneous abortion in the first trimester.

Chr	Peak start	Peak end	Peak region	Gene name	Regulation	Log2(FC)	P value
chr4	158,171,401	158,172,272	intron	FAM198B	up	5.830	5.76E-03
chr4	119,033,176	119,033,301	3′ UTR	MYOZ2	up	5.829	4.81E-03
chr4	453,467	453,542	intron	ABCA11P	up	5.778	8.64E-03
chr1	169,376,663	169,377,063	exon	NME7	up	5.586	8.54E-03
chrX	101,488,487	101,488,537	non-coding	ARMCX4	up	5.580	7.82E-03
chr4	37,020,627	37,022,548	TTS	LOC100508631	up	5.563	8.18E-03
chr16	69,356,783	69,357,056	exon	TMED6	up	5.263	1.13E-02
chr8	18,221,796	18,221,846	5′ UTR	NAT1	up	5.198	1.48E-02
chr9	88,462,616	88,475,585	intron	NXNL2	up	5.165	1.59E-02
chr4	13,615,416	13,615,466	exon	BOD1L1	up	5.157	1.38E-02
chr4	372,742	375,592	3′ UTR	MIR571	down	−4.054	7.14E-03
chr11	78,658,228	78,658,303	exon	NARS2	down	−3.808	1.70E-02
chr16	12,572,889	12,573,014	3′ UTR	MIR4718	down	−3.588	3.27E-02
chr3	37,819,313	37,819,463	3′ UTR	ITGA9-AS1	down	−3.415	1.37E-03
chr7	83,135,139	83,135,189	exon	PCLO	down	−3.388	2.05E-02
chr1	241,595,581	241,595,642	TTS	CHML	down	−3.307	2.46E-02
chr1	147,757,292	147,758,409	3′ UTR	GJA5	down	−3.233	4.10E-08
chrX	72,204,881	72,204,956	3′ UTR	PIN4	down	−3.199	4.94E-02
chr3	37,819,488	37,819,538	3′ UTR	ITGA9-AS1	down	−3.167	2.36E-03
chr22	30,663,347	30,663,422	TTS	DUSP18	down	−3.159	2.39E-02

3′UTR, 3′untranslated region; 5′UTR, 5′untranslated region; TTS, transcription termination site.

GO and KEGG Pathway Enrichment Analysis of Differentially m6A-Modified mRNA

To investigate the biological significance of m6A 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 m6A 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 m6A 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 m6A 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.

FIGURE 2

GO and KEGG pathway enrichment analyses of differentially methylated mRNA. (A) The top 5 GO terms of genes with significantly upregulated m6A peaks. (B) The top 5 GO terms of genes with significantly downregulated m6A peaks. (C) The top 10 KEGG pathways of genes with significantly upregulated m6A peaks. (D) The top 10 KEGG pathways of genes with significantly downregulated m6A peaks.

FIGURE 3

Four major KEGG pathways. (A) Phosphatidylinositol signaling system. (B) Hippo signaling pathway–multiple species. (C) Hippo signaling pathway. (D) Wnt signaling pathway. The colors in the graph from red to green indicate the change from high to low in m6A dysregulation on genes. The colors in the graph from yellow to blue indicate the change from high to low in m6A dysregulation on compounds.

Overview of mRNA Expression Profiles and Conjoint Analysis of meRIP-Seq and RNA-Seq

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 m6As and host gene expression level. We found that the expression level of host gene of differential m6As was lower in both SA group and normal group compared with constitutive m6As (Supplementary Figure S4). Dividing all differentially methylated m6A peaks with all differentially expressed mRNAs into four groups (hyper-up, hyper-down, hypo-up, and hypo-down), we identified 34 hypermethylated m6A peaks in mRNAs that were significantly upregulated (28; hyper-up) or downregulated (6; hyper-down), while 2 hypomethylated m6A 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 m6A 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.

FIGURE 4

Conjoint analysis of MeRIP-seq and RNA-seq data. (A) Volcano plots showing the differentially expressed genes in villous tissues of SA group compared with normal group. (B) Heatmap plots showing the differentially expressed genes in villous tissues of SA group compared with normal group. (C) Four-quadrant plots showing the distribution of genes with significant changes in both the m6A modification and mRNA levels. (D) The top 5 GO terms of genes with significant changes in both the m6A modification and mRNA levels.

TABLE 2

The top 20 differentially expressed mRNAs in human villi during spontaneous abortion in the first trimester.

Gene name	Log2(FC)	Regulation	Location	Strand	P-value
LILRB4	6.994372489	up	chr19:54643889-54670359	+	2.56E-05
KRT6A	6.408677648	up	chr12:52487174-52493257	−	0.000353
TMEM176B	6.402794137	up	chr7:150791285-150801360	−	0.000225
TMEM176A	6.289546685	up	chr7:150800403-150805120	+	0.000221
WNT10A	6.231355116	up	chr2:218880363-218899581	+	0.000217
CTD-2020K17.3	6.217088894	up	chr17:45238028-45241734	−	0.000581
CXCL10	6.158743623	up	chr4:76021117-76023497	−	0.000178
TNC	5.941216251	up	chr9:115019578-115118257	−	1.42E-05
CD300E	5.778269717	up	chr17:74609887-74623738	−	2.10E-06
CXCL9	5.640406283	up	chr4:76001275-76007488	−	0.000375
HBE1	−13.27930895	down	chr11:5268345-5505617	−	1.55E-06
HBZ	−9.7258281	down	chr16:152687-154503	+	1.59E-05
SLC4A1	−8.109924704	down	chr17:44248385-44268141	−	3.30E-05
SPTB	−7.906942374	down	chr14:64746283-64879883	−	0.000189
COX4I2	−7.653108088	down	chr20:31637888-31645006	+	0.000574
LINC02484	−4.950756206	down	chr4:34120894-34269747	−	7.19E-10
ADAMTS18	−4.478881307	down	chr16:77247813-77435114	−	0.000323
RAPGEF4	−4.44120067	down	chr2:172735274-173052893	+	6.25E-10
GOLGA2P7	−4.061744557	down	chr15:84199311-84230136	−	0.000712
AGTR1	−3.91772292	down	chr3:148697784-148743008	+	4.73E-05

Discussion

In this study, we performed high-throughput sequencing to reveal the m6A transcriptome-wide map in human villi during spontaneous abortion in the first trimester. Using the MeRIP-seq data, we found 2,398 genes corresponding to 2,440 altered m6A peaks, which were highly enriched by the m6A consensus motif (RRACH). The enrichment was not perfectly overlapping the RRACH motif described in literature, which might be related to the parameters (-mask -rna -len 6) used for motif search. These genes were mainly involved in lysine degradation, Hippo signaling pathway, ubiquitin-mediated proteolysis, and glycosphingolipid and glycosaminoglycan biosynthesis. Through conjoint analysis of meRIP-seq and RNA-seq data, we identified 35 genes with differentially methylated m6A peaks and synchronously differential expression, which revealed the relationship between m6A methylation and gene expression. These genes were enriched in the Wnt signaling pathway, phosphatase activity regulation, protein phosphatase inhibitor activity, and transcription inhibitor activity. It was recently reported that m6A inhibition through targeted strategies was effective in counteracting different diseases, such as myeloid leukaemia. This could provide a background for the development of therapeutics and for further investigations in the future (Bedi et al., 2020; Garbo et al., 2021; Moroz-Omori et al., 2021; Yankova et al., 2021).

The cellular mechanisms underlying SA are the proliferation and apoptosis of cytotrophoblasts and human decidual cells (Cinar et al., 2012). Studies have shown that at 6-8 weeks in the first trimester of pregnancy, the expression level of TIMP-1 (tissue inhibitor of MMP-2) in villous tissue is significantly decreased, which leads to the abnormal invasion of trophoblast cells, and thus leads to spontaneous abortion (Kesanakurti et al., 2013). As reported, the high expression of e-cadherin (E-cad) in villous tissue affects the invasion of trophoblast, making it difficult for placenta implantation, resulting in spontaneous abortion (Li et al., 2017). MiR126, located in the region of epidermal growth factor-like domain 7 (EGFL7), negatively regulates vascular endothelial growth factor (VEGF), which reduces shallow implantation of trophoblasts, and finally leads to spontaneous abortion (Schmidt et al., 2007).

Basing on the MeRIP-seq data, we identified some differentially methylated mRNAs which were closely linked to many important pathways. KEGG pathway enrichment analysis results indicated that genes with upregulated m6A modification sites were involved in the regulation of lysine degradation. Some studies have confirmed the relationship between lysine metabolism and early embryo development. Studies illuminated that lysine deprivation during low-protein diets could adversely affect early embryo development (Van Winkle et al., 2020). Lysine was specific to LSD1, a demethylase, which regulated the expression and appropriate timing of key developmental regulators during early embryonic development (Foster et al., 2010). Our results indicated that m6A modification might affect early embryonic development by regulating lysine degradation. For these genes with upregulated m6A modification sites, another related pathway was Hippo signaling pathway. Hippo signaling plays a critical role in early embryonic development as low Hippo activity is required for trophoblast differentiation and high Hippo activity permits inner cell mass formation (Wu and Guan, 2021). During murine preimplantation embryogenesis, Hippo signaling pathway is known to play a significant role in lineage segregation and henceforth the formation of blastocysts (Sasaki, 2015). Our results suggested that modulating m6A modifications of the Hippo signaling pathway might be a possible therapy for in human villi during spontaneous abortion in the first trimester in the future.

In addition, another related pathway was ubiquitin-mediated proteolysis. Studies suggested that ubiquitin-mediated proteolysis could be used to regulate Hippo signaling and thus participate in early embryonic development (Ma et al., 2018). Therefore, m6A modification might affect ubiquitin-mediated proteolysis to regulate Hippo signaling and thus regulate early embryonic development. Genes with downregulated m6A modification sites were mainly enriched in glycosphingolipid and glycosaminoglycan biosynthesis. Glycosphingolipids (GSLs) were a class of ceramide-based glycolipids essential for embryo development in mammals (Yamashita et al., 1999; Russo et al., 2016), whether glycosaminoglycan (GAG) biosynthetic was important for mouse embryonic stem cells (mESCs) (Nairn et al., 2007). Based on our results, we hypothesized that m6A modification might influence the expression level of some genes related to glycosphingolipid and glycosaminoglycan biosynthesis and thus influence embryo development.

Through the conjoint analysis of MeRIP-seq and RNA-seq data, thirty-five genes with differentially methylated m6A peaks and synchronously differential expression in spontaneous abortion were discovered (Supplementary Table S1). Among these genes, IGFBP3, C/EBPβ may be regulated by m6A modification of mRNAs. IGFBP3 (Insulin-like growth factor binding protein 3) with high expression level suggested betted oocyte maturation and early embryo development (Wang et al., 2006). IGFBP3 is highly expressed in the endometrium and at the maternal-fetal interface, which promoted the matrix metalloproteinases 2 (MMP2) expression and cell migration in both human endometrial stromal cells (HESCs) and primary human decidual stromal cells (HDSCs) (Luo et al., 2020). Combined with our results, m6A modification may promote the expression of IGFBP3 and thus promote the high expression of MMP2 and the low expression of its tissue inhibitors TIMP-1, and finally lead to spontaneous abortion. CCAAT/enhancer binding protein β (C/EBPβ) is the earliest marker of enveloping layer (EVL) and is essential for EVL differentiation in zebrafish (Zhang et al., 2021). Studies indicated that C/EBPβ transcription factor could inhibit the mRNA decay of IL-8 and thus repress the inflammatory response (Zhang et al., 2010). Moreover, C/EBPβ is also a biomarker of endometrial receptivity and plays a conserved functional role during embryo implantation (Kannan et al., 2010). Combined with our results, altered m6A modification may influence the expression of C/EBPβ, regulate embryo implantation and thus influence early embryo development. However, detailed molecular mechanisms are still unknown and further exploration deserves careful consideration in the future.

Conclusion

Here, we systematically investigated the whole-transcriptome m6A profile of human villous tissues during spontaneous abortion in the first trimester, revealing a dynamic m6A methylation landscape in spontaneous abortion for the first time. Based on the conjoint analysis of MeRIP-seq and RNA-seq data, many genes with differentially methylated m6A peaks and synchronously differential expression were discovered. It indicated a potential link between m6A methylation and mRNA expression, and might provide an alternative therapeutic strategy for spontaneous abortion. In addition, the m6A modification profile might provide novel insights into the pathogenesis and treatment of spontaneous abortion during early pregnancy.