The Alteration of m6A Modification at the Transcriptome-Wide Level in Human Villi During 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 m 6 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 6 A peaks, which corresponded to transcripts of 2,087 genes, and 281 significantly downregulated m 6 A peaks, which corresponded to transcripts of 311 genes (| log2 (fold change) | > 0.585 and p < 0.05). The top 20 altered m 6 A peaks were listed in Table 1 . Then, we investigated the distribution of m 6 A peaks in the SA and normal group and found that m 6 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 6 A peaks in the SA group showed a distinct pattern from m 6 A peaks in the normal group with a relative decrease in the number of m 6 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 6 A peaks per gene, we found that most genes only had one corresponding m 6 A peak (1,878/236 genes with upregulated and downregulated peaks, respectively) ( Figure 1D ). Furthermore, dysregulated m 6 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 6 A peaks were mainly characterized by AUGGAC and UGGACG motif ( Figure 1F ). Characteristics of m 6 A methylation in human villi during spontaneous abortion in the first trimester. (A) Volcano plots showing the significantly altered m 6 A peaks. (B) Accumulation of the region of average m 6 A peaks along with all transcripts in SA group and normal group. (C) Pie charts showing the distribution of m 6 A peaks in SA group and normal group. (D) The distribution of altered m 6 A peaks per gene. (E) The distributions of altered m 6 A peaks in all chromosomes. (F) The top two m 6 A motifs enriched from the altered m 6 A peaks. Top 20 altered m 6 A peaks in human villi during spontaneous abortion in the first trimester. 3′UTR, 3′untranslated region; 5′UTR, 5′untranslated region; TTS, transcription termination site. To investigate the biological significance of m 6 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 6 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 6 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 6 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 . GO and KEGG pathway enrichment analyses of differentially methylated mRNA. (A) The top 5 GO terms of genes with significantly upregulated m 6 A peaks. (B) The top 5 GO terms of genes with significantly downregulated m 6 A peaks. (C) The top 10 KEGG pathways of genes with significantly upregulated m 6 A peaks. (D) The top 10 KEGG pathways of genes with significantly downregulated m 6 A peaks. 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 m 6 A dysregulation on genes. The colors in the graph from yellow to blue indicate the change from high to low in m 6 A dysregulation on compounds. 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 6 As and host gene expression level. We found that the expression level of host gene of differential m 6 As was lower in both SA group and normal group compared with constitutive m 6 As ( Supplementary Figure S4 ). Dividing all differentially methylated m 6 A peaks with all differentially expressed mRNAs into four groups (hyper-up, hyper-down, hypo-up, and hypo-down), we identified 34 hypermethylated m 6 A peaks in mRNAs that were significantly upregulated (28; hyper-up) or downregulated (6; hyper-down), while 2 hypomethylated m 6 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 6 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. 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 m 6 A modification and mRNA levels. (D) The top 5 GO terms of genes with significant changes in both the m 6 A modification and mRNA levels. The top 20 differentially expressed mRNAs in human villi during spontaneous abortion in the first trimester.

Here, we systematically investigated the whole-transcriptome m 6 A profile of human villous tissues during spontaneous abortion in the first trimester, revealing a dynamic m 6 A 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 m 6 A peaks and synchronously differential expression were discovered. It indicated a potential link between m 6 A methylation and mRNA expression, and might provide an alternative therapeutic strategy for spontaneous abortion. In addition, the m 6 A modification profile might provide novel insights into the pathogenesis and treatment of spontaneous abortion during early pregnancy.

In this study, we performed high-throughput sequencing to reveal the m 6 A 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 m 6 A peaks, which were highly enriched by the m 6 A 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 m 6 A peaks and synchronously differential expression, which revealed the relationship between m 6 A 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 m 6 A 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 m 6 A 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 m 6 A modification might affect early embryonic development by regulating lysine degradation. For these genes with upregulated m 6 A 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 m 6 A 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, m 6 A modification might affect ubiquitin-mediated proteolysis to regulate Hippo signaling and thus regulate early embryonic development. Genes with downregulated m 6 A 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 m 6 A 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 m 6 A peaks and synchronously differential expression in spontaneous abortion were discovered ( Supplementary Table S1 ). Among these genes, IGFBP3, C/EBPβ may be regulated by m 6 A 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, m 6 A 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 m 6 A 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.

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, m 6 A 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 m 6 A modification plays important roles in the regulation of immune function and inflammatory response, the relationship between m 6 A methylation and SA remains to be elucidated. As known, the m 6 A 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, m 6 A 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 m 6 A profiles of major fetal tissues including human placenta ( Xiao et al., 2019 ), but the m 6 A 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 m 6 A 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 m 6 A 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 m 6 A peaks and synchronously differential expression, which might provide an alternative strategy for the therapy and prevention of spontaneous abortion.

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 6 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 6 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 6 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 6 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 6 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.