Epitranscriptomic profiling in human placenta: N6-methyladenosine modification at the 5'-untranslated region is related to fetal growth and preeclampsia.

Intracellular mRNA levels are not always proportional to their respective protein levels, especially in the placenta. This discrepancy may be attributed to various factors including post-transcriptional regulation, such as mRNA methylation (N6-methyladenosine: m6A). Here, we conducted a comprehensive m6A analysis of human placental tissue from neonates with various birth weights to clarify the involvement of m6A in placental biology. The augmented m6A levels at the 5'-untranslated region (UTR) in mRNAs of small-for-date placenta samples were dominant compared to reduction of m6A levels, whereas a decrease in m6A in the vicinity of stop codons was common in heavy-for-date placenta samples. Notably, most of these genes showed similar expression levels between the different birth weight categories. In particular, preeclampsia placenta samples showed consistently upregulated SMPD1 protein levels and increased m6A at 5'-UTR but did not show increased mRNA levels. Mutagenesis of adenosines at 5'-UTR of SMPD1 mRNAs actually decreased protein levels in luciferase assay. Collectively, our findings suggest that m6A both at the 5'-UTR and in the vicinity of stop codon in placental mRNA may play important roles in fetal growth and disease.

appropriate‐for‐date

coding DNA sequence

continuous POI

Database for Annotation, Visualization and Integrated Discovery

differentially expressed genes

differentially methylated genes

embryonic stem

false discovery rate

fragments per kilobase of exon per million mapped fragments

Gene Ontology

transformed human embryonic kidney cell

heavy‐for‐date

horseradish peroxidase

immunoprecipitated

Ingenuity Pathway Analysis

intrauterine growth restriction

N6‐methyladenosine

methylated RNA immunoprecipitation followed by sequencing

next‐generation sequencing

phosphate‐buffered salineFBS, fetal bovine serum

Principal component analysis

preeclampsia

peak‐over‐input

peak‐over‐median

quantitative real‐time PCR

small‐for‐date

transforming growth factor beta

untranslated region

western blot

INTRODUCTION

Per the central dogma of molecular biology, information encoded in DNA is transcribed to mRNA, which, in turn, is translated into proteins.1 Regulation occurs at various stages of gene expression; one such stage is post‐transcriptional regulation, which is classified into two types: cis and trans. Regulation via RNA‐binding proteins and microRNA constitutes trans‐type post‐transcriptional regulation, whereas cis‐type post‐transcriptional regulation is characterized by modifications at the untranslated regions (UTRs) and modification of RNA bases after transcription.2 Over 100 different types of RNA modifications have been reported, including N6‐methyladenosine (m6A), 5‐methylcytosine (m5C), N1‐methyladenosine (m1A), and pseudouridine (Ψ).3 Of these, m6A is the most frequently observed chemical post‐transcriptional mRNA modification.4, 5 Methylated RNA immunoprecipitation followed by sequencing (MeRIP‐Seq) is an epitranscriptome‐wide assay to determine the presence of m6A.6, 7 Since the development of this assay in 2012, the functions of m6A have gradually become clearer.6, 7 Sequences containing m6A are usually located in the vicinity of a stop codon, especially within the 3′‐UTR, and they have a consensus sequence of RRACHR, where R is a purine and H is any base except for G.6, 7 m6A modifications are involved in post‐transcriptional regulation, especially in determining the stability and lifespan of mRNA.2 Furthermore, several m6A regulators have been reported, such as the methylating enzyme m6A writer protein (METTL3, METTL14, and WTAP),8 the demethylating enzyme m6A eraser protein (FTO, ALKBH5),9 and m6A reader proteins (YTH family), that recognize m6A.10, 11, 12, 13, 14, 15 The site of m6A in mRNA differs with cell and tissue types; moreover, m6A levels change in response to external stimuli,16, 17 thereby functioning as a dynamic type of modification that fine‐tunes gene expression.18 Notably, although the levels of mRNA (arising from gene expression) and protein are positively correlated, the correlation is weak, that is, it is not a perfect correlation.19, 20 Therefore, elucidation of the specifics of post‐transcriptional mRNA regulation can help clarify the role of genes involved in various cellular events.

The placenta is an essential organ in female members of the clade Eutheria (including humans), which forms a barrier between the body of a mother and a fetus and functions to provide nutrition, gas exchange, and elimination of waste products.21 Placental DNA generally tends to be less methylated than that in other organs,22 with retrotransposon‐derived genes playing a role in placental development and maintenance.23, 24 Moreover, a recent placental single‐cell analysis revealed that proliferation, differentiation, and regeneration are maintained in normal‐term placenta.25 Therefore, placental gene expression regulation is unique, and the gene expression profile of a postpartum placenta may reveal essential information regarding the long‐term functions of the placenta throughout pregnancy.

Recently, the number of small‐for‐date (SFD) children being born in perinatal clinics has increased, and medical advancements have greatly improved the prognosis of SFD infants26; however, such infants continue to be at an increased risk of long‐term illnesses such as type II diabetes mellitus, high blood pressure, and hyperlipidemia.27 Moreover, infants that are SFD owing to preeclampsia (PE) undergo oxidative stress and inflammation as a consequence of their maternal abnormal spiral arteries and abnormal villar neovascularization during the placental development.28 However, a transcriptome analysis revealed only two genes with different expression levels between normal placentas and the placentas from SFD infants, with no associated complications.29 Conversely, 98 genes with differential expression levels in placentas from SFD infants with PE compared to normal control placentas have been reported,29 along with 41 genes with differential expression levels in placentas from macrosomia infants compared to those from healthy infants.29, 30 A proteome‐wide analysis of the placentas of mice with intrauterine growth restriction (IUGR) after artificial fertilization identified 178 types of proteins with significant changes in expression levels compared to control placentas. Conversely, there were no significant differences in the mRNA expression levels of these 178 genes, thereby rendering them to be of considerable interest in analyzing post‐transcriptional regulation. Among these 178 proteins, those involved in post‐transcriptional and ‐translational regulation were significantly more frequent and significantly upregulated.20 In other words, these results suggest that the unique characteristics of placentas in SFD infants cannot be determined through transcriptome analysis; rather, these characteristics result from post‐transcriptional and translational regulation. In addition, expression levels of FTO in placentas, which are known to demethylate m6A,9 are positively correlated with birth weight31 and prenatal fetal head circumference at 34 weeks of gestation.32 Considering the aforementioned findings, we hypothesized that m6A in placental mRNA constitutes an important mechanism in post‐transcriptional regulation and is related to fetal development. We therefore conducted MeRIP‐Seq on human placental tissue samples obtained from mothers of infants of various birth weights, profiled the m6A epitranscriptome of these specimens, and investigated the relationship between fetal development and placental m6A.

MATERIALS AND METHODS

Aim, design, and setting of the study

In human placental tissue, m6A in mRNA obtained from chorionic villi may play an important role in fetal development. The present epitranscriptome‐wide study aimed to clarify the relationship between m6A modification and fetal development or perinatal disease in placenta samples collected from 17 patients (Set1) and 8 patients (Set2) who underwent a cesarean section. RNA‐Seq and comprehensive m6A analysis (MeRIP‐Seq) were performed with Set1. All subjects were recruited from the National Research Institute for Child Health and Development, Tokyo, Japan. All participants provided informed consent in accordance with the tenets of the Declaration of Helsinki, and ethical approval was obtained from the ethics committee of the National Research Institute for Child Health and Development Ethics Committee (reference number 2014/630, Japan).

Definition of fetal growth categories and preeclampsia

The appropriate‐for‐date (AFD) group was defined as having a birth weight greater than or equal to the 10th percentile and less than the 90th percentile of the Japanese population. The SFD group met the criteria that both birth weight and height were less than the 10th percentile, and the heavy‐for‐date (HFD) group was in the 90th percentile or higher. Gravity was also considered. The diagnostic criteria for PE was as follows: blood pressure greater than or equal to 140 mm Hg systolic or greater than or equal to 90 mm Hg diastolic on two occasions at least 4 hours apart after 20 weeks of gestation in a woman with a previously normal blood pressure. Proteinuria was considered as greater than or equal to 300 mg per 24 hours urine collection (or this amount extrapolated from a timed collection) or a protein/creatinine ratio greater than or equal to 0.3.33

Cell culture and RNA extraction from placental tissue and cell lines

Chorionic villi samples were washed in phosphate‐buffered saline (PBS) three times to eliminate maternal blood. Then, the samples (approximately 1 mg) were treated with TRIzol (Thermo Fisher Scientific, Waltham, MA, USA) and homogenized using scissors within 15 min after placental delivery. Total RNA was extracted using TRIzol reagent in accordance with the manufacturer's instructions.

The transformed human embryonic kidney cell line (HEK293T) was a gift from Dr Senji Shirasawa.34 HEK293T was incubated in Dulbecco's modified Eagle's medium (Thermo Fisher Scientific) containing 4.5 g/l d‐glucose and GlutaMAX supplemented with 10% heat‐inactivated fetal bovine serum (FBS) and penicillin/streptomycin at 37°C and 5% CO2. Total RNA was extracted using the RNeasy Mini Kit (Qiagen, Hilden, Germany) in accordance with the manufacturer's instructions. DNA was simultaneously eliminated using an RNase‐Free DNase Set (Qiagen).

Poly(A) enrichment of RNA from total RNA of a placenta and HEK293T cells was performed in two rounds using Dynabeads Oligo(dT)25 (Thermo Fisher Scientific). The quality of each total RNA and poly(A) RNA was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

JEG3 (human choriocarcinoma cell line) purchased from ECACC was incubated in Eagle's minimal essential medium (Nacalai Tesque, Kyoto, Japan) containing 10% inactivated FBS and penicillin/streptomycin at 37°C and 5% CO2. Approximately 10% of cells from sub‐confluent cultures were split off into a new dish, and these cells were cultured to sub‐confluent state. For stimulation, JEG3 cells were treated with or without either transforming growth factor beta (TGFβ)1 or TGFβ3 (5 ng/mL, recombinant human TGFβ1 (240‐B) (R&D Systems, Minneapolis, MN, USA) and TGFβ3 (243‐B3) (R&D Systems) for 24 hours. Subsequently, RNA and protein were extracted. Total RNA was also extracted using the RNeasy Mini Kit (Qiagen) in accordance with the manufacturer's instructions. DNA was eliminated using an RNase‐Free DNase Set (Qiagen).

Protein extraction, antibodies, and western blot analysis

After washing three times with cold PBS, the cells were incubated with 500 mM RIPA buffer (NIPPON GENE, Tokyo, Japan) by adding 1% NP‐40 and 1% protease inhibitor cocktail (Merck KGaA, Darmstadt, Germany) for 10 min on ice. Thereafter, they were collected with a cell scraper, homogenized with a 27G needle, centrifuged at 15 000 g at 4°C for 15 min, and the supernatant was collected.

Western blot (WB) was performed using human villi tissues. Primary antibodies were used in the dilutions as mentioned below. Following incubation with horseradish peroxidase (HRP)‐conjugated secondary antibodies, membrane blots were visualized using ECL start Western Blotting Detection Reagent (RPN3243) (GE Healthcare, Chicago, IL, USA) and imaged on ImageQuant LAS4000 (GE Healthcare). The band intensity of SMPD1 and ACTB was digitized by ImageQuant LAS4000 using the ImageQuant TL Analysis Toolbox (GE Healthcare). Then, the amounts of SMPD1 normalized by ACTB were compared to controls for each analysis.

Antibodies against β‐Actin (A2228, mouse monoclonal [WB 1:5000]) and SMPD1 (mAbcam74281, mouse monoclonal [WB 1:500]) were purchased from Merck and Abcam (Cambridge, UK), respectively. The secondary antibodies, horse anti‐mouse IgG (HRP‐linked Antibody #7076) and sheep anti‐mouse IgG (HRP‐Linked Whole Ab NA931), were purchased from Cell Signaling Technology (Danvers, MS, USA) and GE Healthcare, respectively.

RNA fragmentation

RNA samples were chemically fragmented into 200‐nucleotide‐long fragments through a 25 s incubation at 94°C in Nature fragmentation buffer (10 mM ZnCl2 and 10 mM Tris‐HCl; pH 7).6 Fragmentation was terminated with 0.5 M ethylenediaminetetraacetic acid (pH 8.0) and the reaction mixture was incubated on ice. Thereafter, fragmented RNA for immunoprecipitation was purified using an RNeasy MinElute Cleanup Kit (Qiagen) in accordance with the manufacturer's instructions. Finally, the size of each sample was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies).

Immunoprecipitation and reverse transcription

After two‐round calibration of 1 M IP buffer (1 M NaCl, 10 mM sodium phosphate, 0.05% Triton‐X), 100 µl of Dynabeads M‐280 Sheep anti‐Rabbit IgG (11203D) (Thermo Fisher Scientific) was linked to 2.5 µg of anti‐m6A polyclonal antibody (Cat. No. 202 003) (Synaptic Systems, Goettingen, Germany) by rotating slowly at 4°C in 1 M IP buffer (500 µl) overnight. The fragmented poly(A) RNA (approximately 1 µg) was incubated with 2.5 µg of the affinity‐purified anti‐m6A polyclonal antibody in 140 mM IP buffer (140 mM NaCl, 10 mM sodium phosphate, and 0.05% Triton‐X) with gradual rotation at 4°C for 1 hours. Thereafter, the complex was washed three times with 140 mM IP buffer and twice with both low‐ and high‐salt buffer (low: 50 mM NaCl, 0.1% Nonidet [R] P‐40, 10 mM Tris‐HCl [pH 7.5] 500 mM NaCl; high: 0.1% Nonidet [R] P‐40, 10 mM Tris‐HCl [pH 7.5]) to reduce non‐specific binding of RNA to the anti‐m6A antibody. m6A‐tagged fragmented mRNAs were competitively eluted from the beads with free N6‐methyladenosine (M2780) (Sigma‐Aldrich, St. Louis, MO, USA) in 140 mM IP buffer with an adequate amount of RNase Inhibitor added (Thermo Fisher Scientific). Furthermore, immunoprecipitated RNA was purified using an RNeasy MinElute Cleanup Kit (Qiagen) in accordance with the manufacturer's instructions. Reverse transcription of input and immunoprecipitated (IP) RNA samples was carried out using SuperScript III (Thermo Fisher Scientific) in accordance with the manufacturer's instructions.

RNA‐Seq and MeRIP‐Seq

Strand‐specific and 101‐bp paired‐end RNA‐Seq libraries were generated using SureSelect Strand Specific RNA (Agilent Technologies) in accordance with the manufacturer's protocol and sequenced on an Illumina HiSeq1500 or 2500 (Illumina, San Diego, CA. USA) platform to a depth of ~60 million reads each. Low‐quality reads were discarded using the Illumina CASAVA or bcl2fastq tool. After exclusion of the adaptor sequence and extremely short reads, the reads were mapped on the 1,000 genomes reference genome sequence (hs37d5)(ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz), using tophat‐2.0.7 35 and bowtie2‐2.0.6.36 Final bam files were generated after removal of completely matched pair‐end reads, using Picard (https://broadinstitute.github.io/picard/).

Here, we verified our m6A detection system. We performed next‐generation sequencing (NGS) following MeRIP of cultured HEK293T cells, whose epitranscriptome has been reported previously,7 to evaluate the reproducibility of the MeRIP‐Seq. To identify the m6A locations, a peak‐over‐input (POI) score was calculated for every 50‐nt window to determine the m6A levels. The m6A regions were defined as a window with a POI score > 3. These regions were annotated using HOMER, with hg19 as a reference (see method). The proportions of m6A regions were as follows: 3′‐UTR (29%), the CDS excluding the last exon (31%), and the last exon in some cases including the 3′‐UTR (25%) (Supplemental Figure S1A, left), equivalent to those previously reported.7 The m6A‐modified windows were concentrated in the vicinity of the stop codon in HEK293T cells (Supplemental Figure S1B, left). The regions of m6A had the same consensus sequence of RRACH as reported previously (Supplemental Figure S1C, left). We also evaluated the m6A detection procedure based on POI score using published MeRIP‐Seq data of a human embryonic stem (ES) cell line (H1), which is available from Gene Expression Omnibus, NCBI (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52600).37 The results were similar to our HEK293T cell and reported results 37; the distribution and consensus sequence of m6A‐modified regions (Supplemental Figure S1A, right and Supplemental Figure S1C, right). These results indicated that our MeRIP‐Seq experimental system has high reproducibility, concurrent with previous in vitro and in silico results. Hence, this system was used for profiling of the placental m6A epitranscriptome.

RNA‐Seq analysis

Transcript assembly was conducted using Cufflinks 2.2.1.35 The assemblies were merged using Cuffmerge. Differentially expressed genes (DEGs) related to the fetal growth phenotype groups were assessed using Cuffdiff. On DEG analysis, a q value (false discovery rate [FDR]‐corrected P value) <.05 was considered statistically significant. The Fragments Per Kilobase of exon per Million mapped fragments (FPKM) value of placenta obtained using Cufflinks and DEGs between each fetal growth category data are shown in E‐MTAB‐6507 in ArrayExpress.

DETECTION OF m6A REGIONS: PEAK‐OVER‐INPUT SCORING

Identification of m6A sites

A previously reported strategy was adopted to identify m6A regions in the IP sample relative to the input sample.12 Input and IP bam files were divided to 50‐nt windows using igvtools.38 In the input sample, selecting windows, with more than 15 counts of reads and corresponding windows were selected for the IP files. On normalizing read counts by each median read count of input and IP windows, a peak‐over‐median (POM) score was assigned for each window. Finally, a POI score was assigned to each m6A region by calculating the ratio of POM in the IP sample to that in the input sample. The following formula was used to determine the POI score:

Putative m6A sites were defined when the POI score was higher than 3. The windows were annotated using the HOMER annotation tool (http://homer.ucsd.edu/homer/motif/index.html), and the position of each m6A site was classified into four mutually exclusive mRNA structural regions including the 5′‐UTR, coding DNA sequence (CDS), last exon, 3′‐UTR, and the added intron.

Continuous POI score

Continuous POI (cPOI) score was defined as the summation of POI scores annotating the same mRNA structural regions for one gene, as m6A regions often overlapped and clustered.7 We calculated the cPOI score of each transcript at each region and compared the score among different tissues and placenta samples (E‐MTAB‐6507 in ArrayExpress). The same analysis was carried out for each isoform (E‐MTAB‐6507 in ArrayExpress).

m6A motif analysis

The m6A regions with a POI score greater than 10 were selected for consensus motif finding. Sequences of these windows were analyzed using HOMER 39 motif tools. The windows were aligned to the reference genome (hg19) and “–rna” option was used to fit RNA sequences using findMotifsGenome.pl.

Gene ontology (GO) analysis

GO analysis was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) 40 to characterize our list of genes. As a functional annotation, we examined the GO term of Biological processes. FDR was adopted as adjusted P value as a criterion of significant difference. Those with FDR less than 0.05 were regarded as significant GO terms.

Metagene analysis

Metagene analysis was performed using the ver 2.0.041 for R version 3.5.0. to include RNA‐related genomic features. Metagene plot including 5′‐UTR, CDS, and 3′‐UTR was generated using POI score and the input and IP BAM files.

m6A peak of representative transcripts

The mean POI score of each window for each fetal growth category was plotted using the ggplot2 package42 for R version 3.5.0.

Calculation of correlation coefficients

For gene expression data, FPKM with more than 0.3 were transformed to their Log2 values. For m6A modification data, the cPOI scores were added to 1.0 and were transformed to their Log2 values.

Correlation coefficients of FPKM and the cPOI score of each transcript at the 5′‐UTR and in the vicinity of the stop codon in 17 placenta samples were calculated using R version 3.5.0 (https://cran.r-project.org/). Correlation coefficients were calculated as well as the average score among the same groups.

Ingenuity Pathway Analysis (IPA)

Molecular and functional annotations were performed using IPA software (TOMY Digital Biology Co., Tokyo, Japan). Regarding epitranscriptome features of human normal placenta, the average cPOI value of three AFD samples was calculated at the 5′‐UTR and in the vicinity of the stop codon. Functional annotation was performed for each of the top 10% genes using IPA. To be considered significant, resulting P values of the molecular and cellular functions had to fall below 0.05.

To investigate mRNA expression and m6A modification in normal placenta, molecular annotation was performed based on the gene list with the top 10% average values of FPKM and cPOI from three AFD samples. Proportions of each molecular type are listed. Relatedly, differentially methylated genes (DMGs) were compared in AFD samples and the proportions listed for each of two genic regions.

Principal component analysis (PCA)

Using prcomp of R version 3.5.0, PCA was conducted at the RNA expression and m6A modification levels. FPKM ≤ 0.3 was regarded as 0, and the common logarithmic transformation was performed to FPKM. Logarithmic transformation to base 2 was performed to cPOI. Finally, PCA was conducted using data of 17 placenta samples.

RT‐qPCR

Total RNA was extracted using an RNeasy Mini Kit (Qiagen, Hilden, Germany) in accordance with the manufacturer's instructions. First‐strand cDNA was synthesized using SuperScript III Reverse Transcriptase (Thermo Fisher Scientific, Waltham, MA, USA). Quantitative real‐time PCR (RT‐qPCR) was performed with TB Green Premix Ex Taq II (Tli RNaseH Plus) (TaKaRa Bio Inc, Shiga, Japan) and an Applied Biosystems 3500 (Thermo Fisher Scientific). The relative gene expression was normalized against expression of a housekeeping gene, ACTB or 18SrRNA. PCR primers were designed as shown below and synthesized at FASMAC (Kanagawa, Japan). All the primer sequences were validated using UCSC In‐Silico PCR to ensure their specificity, and melting curve analysis was performed to confirm the amplification of single targets.

Primer information: SMPD1 primer set: Fwd; 5ʹ‐TCC GCC TCA TCT CTC TCA AT‐3ʹ, Rev; 5ʹ‐ATT CCA GCT CCA GCT CTT CA‐3ʹ; Luc2 primer set: Fwd; 5ʹ‐TGC AAA AGA TCC TCA ACG TG‐3ʹ, Rev; 5ʹ‐AAT GGG AAG TCA CGA AGG TG‐3ʹ; ACTB primer set: Fwd; 5ʹ‐CGT CTT CCC CTC CAT CGT‐3ʹ, Rev; 5ʹ‐AGG GTG AGG ATG CCT CTC TT‐3ʹ, 18SrRNA primer set: Fwd; 5ʹ‐AAA CGG CTA CCA CAT CCA AG‐3ʹ, Rev; 5ʹ‐CCT CCA ATG GAT CCT CGT TA‐3ʹ.

MeRIP‐qPCR

In the experiment using villi tissue from patients and cell line, the primer covering the region subjected to the m6A modification at the 5′‐UTR of SMPD1 was designed as the m6A‐positive region primer, and the exon‐exon junction not subjected to the m6A modification was designed as the m6A‐negative region primer. Based on the difference between the Ct values of both, the extent to which the m6A modification was enriched in the 5′‐UTR of SMPD1 within the same individual was calculated and further compared between samples using the ΔΔCt method.43 RT‐qPCR adopted the same method as described above.

Primer information: SMPD1 m6A‐positive region primer set: Fwd; 5ʹ‐AGC AGT CAG CCG ACT ACA GAG‐3ʹ, Rev; 5ʹ‐CAT TGT CGC GCT TCC TAC AC‐3ʹ; SMPD1 m6A‐negative region primer set: Fwd; 5ʹ‐TG AAG AGC TGG AGC TGG AAT‐3ʹ, Rev; 5ʹ‐CCA GGA TTA AGG CCG ATG TA‐3ʹ.

Plasmid construction and mutagenesis assays

The full‐length 5′‐UTR of SMPD1 (National Center for Biotechnology Information reference sequence NM_000543.4) was PCR amplified from human placental gDNA using primers containing HindIII and NcoI restriction site (Takara Bio Inc, Shiga, Japan) at the 5ʹ end; Fwd: 5ʹ‐AAG CTT AGC TGT CAG AGA TCA GAG G‐3ʹ and Rev: 5ʹ‐ CCA TGG TGT CGC GCT TCC TAC ACG GG‐3ʹ and then inserted upstream of the firefly luciferase of pGL4.53[luc2/PGK] (Promega, WI, USA). For mutagenesis assay, PCR was performed using the mutation insert primer: Rev; 5ʹ‐ CCA TGG TGT CGC GCT TCC TAC ACG GGG CTG GTT CGT CTG CC CG CCC CC‐3ʹ and the above‐mentioned Fwd primer; then, the PCR product was inserted into pGL4.53. All constructs were confirmed by DNA sequencing.

Dual‐luciferase reporter assay

For dual‐luciferase reporter assay, 200 ng of empty pGL4.53 (vector), pGL4.53 with wild‐type or mutant SMPD1‐5′‐UTR, and 20 ng of pNL1.1 [Nluc] Vector (NanoLuc luciferase control reporter vector) (Promega, WI, USA) were co‐transfected into JEG3 cells in a 96‐well plate. To detect translational efficacy, 2 μg/mL actinomycin D was added at 0, 2, or 8 hours before collection. Cells were collected to measure the relative luciferase activity using Nano‐Glo Dual‐Luciferase Reporter Assay System (Promega). The ratio of Firefly to NanoLuc luciferase activity was used as the relative luciferase activity, and all experiments were performed in triplicate.

Statistical analysis

The cPOI score for each transcript at each region in AFD, SFD PE+, PE−, and HFD were compared using a Mann‐Whitney U test with R version 3.3.0, with all P values less than 0.05 considered statistically significant. Hierarchical clustering analysis was performed with Euclidean distance and the complete linkage method, using the ggplot2 package for R version 3.3.0. A dendrogram and heatmap were generated. mRNA and protein levels, quantified using RT‐qPCR and WB respectively, were analyzed using the Student's t test with un‐paired and two‐sided methods, which were conducted using R version 3.3.0. P values less than .05 were considered statistically significant.

RESULTS

m6A in placental mRNA

MeRIP‐Seq was performed to identify m6A locations in an epitranscriptome‐wide analysis of placental mRNA for postpartum placentas of mothers of three infants with normal birth weights (AFD). The characteristics of AFD‐1, AFD‐2, and AFD‐3 are shown in Table 1. The transcriptome‐wide numbers of m6A mRNAs encoding proteins were 6727, 6509, and 6832 for AFD‐1, AFD‐2, and AFD‐3, respectively. For one placenta, the primary m6A regions were the 3′‐UTR (30%), the CDS of the last exon including the 3′‐UTR (26%), and the CDS excluding the last exon (28%), thereby constituting approximately 85% of the modified regions (Figure 1A); the percentage of 5′‐UTR was only 3%. The m6A‐modified windows were also concentrated in the vicinity of the stop codon in the placental tissue as well as in HEK293T and human ES cells (see method section and Supplemental Figure S1A). We also identified the same consensus sequence with the highest frequency for placental m6A: RRACH (Figure 1B and Supplemental Figure S1C).

Table 1

Characteristic of the placental samples (Set 1)

Sample name	Gestational age	Birth body weight (g)	Height (cm)	Placental weight (g)	Body/placental weight ratio	Maternal age (y.o.)	Parity	RIN	Fetal sex	Preeclampsia
AFD1	38w0d	2739	46.5	664	4.13	41	1	7.4	F	−
AFD2 (AFD‐2)	38w4d	3191	50	530	6.02	36	1	7.5	F	−
AFD3	37w3d	2842	47	560	5.08	39	1	7.8	M	−
AFD4 (AFD‐3)	37w5d	3146	50	560	5.62	35	1	7.3	M	−
AFD5	37w5d	2512	46	452	5.56	40	2	8.4	F	−
AFD6 (AFD‐1)	38w0d	3010	48	516	5.83	33	0	7.6	F	−
SFD1	37w6d	2206	44	385	5.73	32	1	6.8	F	−	SFD1
SFD2	36w0d	1490	39	370	4.03	32	0	6.6	F	+	PE1
SFD3	37w6d	1993	44	315	6.33	30	0	7.3	M	−	SFD2
SFD4	35w2d	1332	39.5	295	4.52	33	0	7.8	M	+	PE2
SFD5	31w5d	1078	37.8	275	3.92	43	2	7.8	M	+	PE3
SFD6	37w0d	1926	42	310	6.21	42	0	8	F	−	SFD3
HFD1	38w6d	3750	50	660	5.68	41	1	7.5	F	−
HFD2	38w0d	4141	52.5	845	4.90	36	1	7.9	F	−
HFD3	38w3d	3544	50	595	5.96	34	1	7.6	F	−
HFD4	38w5d	3520	52	630	5.59	42	0	8.2	F	−
HFD5	38w5d	3400	52	520	6.54	34	0	7.8	M	−

Abbreviations: AFD, appropriate‐for‐date; SFD, small‐for‐date; HFD, heavy‐for‐date; RIN, RNA Integrity Number.

Figure 1

Characteristics of placental mRNA m6A modification. A, Window distribution and proportions of m6A‐modified regions in mRNA obtained from one appropriate‐for‐date (AFD) placenta sample. The regions around the coding DNA sequence (CDS) and the stop codon (the last exon and the 3′‐UTR) constitute approximately 85% of the modified locations. B, Motif analysis of window sequences with a large amount of m6A‐modified poly(A) RNAs (windows that satisfied peak‐over‐input [POI] >10) using HOMER shows a significant concentration of the GGAC motif (P = 10−575). C, Illustration of the concept of continuous POI (cPOI) scores. We calculated and analyzed cPOI scores for m6A levels in two genic regions (5′‐untranslated region [UTR] and vicinity of the stop codon). D, Metagene analysis of m6A modification site for one normal placenta representative, HEK293T cells, and human embryonic stem (ES) cells. E, The number of modified genes that were common among the three appropriate‐for‐date (AFD) placenta samples at the 5′‐UTR (left) and in the vicinity of the stop codon (right). F, Left: Overlap analysis of the genes with the top 10% mRNA expression levels in AFD placentas and the top 10% genes with the highest m6A levels at the 5′‐UTR in the three AFD placenta samples were conducted. Right: A similar overlap analysis was conducted using the top 10% genes in the vicinity of the stop codon. There were virtually small overlaps in either of the regions. G, Gene Ontology (GO) analysis was conducted on the top 10% genes with m6A in the vicinity of stop codon that were shared by the three AFD placenta samples. GO terms with significant concentrations are shown, where an FDR adjusted P < .05 was considered significant

We quantified the enrichment of m6A as the POI score in each 50‐nt window (Figure 1C). To visualize the gene‐wide distribution of m6A, we plotted POI on mRNA using metagene profile (Figure 1D). Comparison of m6A distribution between placenta, HEK293T, and human ES cells revealed that the frequency of m6A at the 5′‐UTR in placenta was higher than that in HEK293T cells but not that in human ES cells. Therefore, we proceeded to analyze the detailed features of m6A modification in human placenta.

Differences in the degree of overlap and gene function characteristics among individual mRNAs containing m6A regions

The 5′‐UTR6, 16 and the domains in the vicinity of the stop codon11, 13, 15 are being widely studied for determining the functions of post‐transcriptional m6A regulations. In this study, we analyzed the characteristics of placental mRNAs with m6A. Fifty‐six percent of m6A‐modified regions existed in either the 3′‐UTR or the last exon, when the data were divided into 50‐nt windows (Figure 1A). This observation is consistent with previous reports 6, 7, 44; hence, we analyzed the m6A sites in the vicinity of a stop codon integrating both the 3′‐UTR and the last exon. m6A peaks are frequently detected as clustered regions.7 Thus, we digitized the data and reported the cPOI scores for m6A clusters at the 5′‐UTR and in the vicinity of the stop codon; mRNAs with high m6A levels were identified in this manner (Figure 1C). The numbers of genes with m6A‐modified mRNAs the three AFD placentas were 1034, 976, and 896 at the 5′‐UTR and 4836, 4784, and 5119 in the vicinity of the stop codon. In total, 645 genes with m6A‐modified mRNAs at the 5′‐UTR were found to be shared by all the three AFD samples, corresponding to 62%, 66%, and 72% of the m6A‐modified genes in each sample. Furthermore, 4124 genes with m6A‐modified mRNAs in the vicinity of the stop codon were shared by all three AFD samples, corresponding to 85%, 86%, and 81% of the m6A‐modified genes in each sample (Figure 1E). Of all the modified placental mRNAs, the overlap between the top 10% of m6A levels and of expression levels was 0% for mRNAs at the 5′‐UTR and 2.2% for those in the vicinity of the stop codon (Figure 1F). These results indicate that m6A modification is independent of the mRNA expression level. The targets of m6A modification among individuals were more common in the vicinity of stop codons than at the 5′‐UTR.

The most significantly enriched Database for Annotation, Visualization, and Integrated Discovery (DAVID) GO term of top 10% genes with highest m6A modification levels in the vicinity of the stop codon in common with three AFD placental mRNA was “transcription, DNA‐templated” (Figure 1G). Functional annotation using IPA also resulted in “Gene Expression: Transcription of RNA” (P = 3.0 × 10−11) as the top function, which was much stronger than the next highest “Cellular Growth and Proliferation, Connective Tissue Development and Function, Tissue Development: Proliferation of fibroblast cell lines” (P < 3.0 × 10−5; Table 2). Conversely, at the 5′‐UTR, GO analysis did not detect any significant GO terms in genes with highest m6A levels (top 10% cPOI scored genes), whereas functional annotation using IPA suggested “1. Cell Death and Survival, Organismal Injury and Abnormalities: Cell death of endothelial cells” (P < 1.08 × 10−5; Table 2) as the top out of four significant functions. An IPA molecular annotation analysis revealed that 15.6% and 22.7% of the top 10% genes with highest m6A levels at the 5′‐UTR and in the vicinity of the stop codon were “transcription regulator,” respectively (Table 3). It was obviously enriched compared with the fact that 6.0% of the top 10% genes with highest expression levels were transcriptional regulators. Transcriptional regulation was the top‐ranked function for the genes with high mRNA m6A levels in HEK293T cells in a previous study.7 In this study, we revealed that the mRNA m6A levels for these genes were concentrated in the vicinity of stop codons, and transcription regulators were markedly enriched in genes highly modified at the 5′‐UTR. These results were also confirmed in human ES cells and HEK293T cells (Supplemental Figure S1A).

Table 2

Functional annotation of genes with high m6A levels in placentas from the healthy birth weight group

Top Molecular and Cellular Functions
Methylated mRNA regions	Category	Function Annotation	P Value
5′‐UTR	1. Cell Death and Survival, Organismal Injury and Abnormalities	Cell death of endothelial cells	1.08E‐05
	2. Post‐Translational Modification	Ubiquitination	1.58E‐04
	3. Cellular Movement	Invasion of cells	2.45E‐04
	4. Lipid Metabolism, Small Molecule Biochemistry	Metabolism of membrane lipid derivative	4.97E‐04
Vicinity of stop codon	1. Gene Expression	Transcription of RNA	3.00E‐11
	2. Cellular Growth and Proliferation, Connective Tissue Development and Function, Tissue Development	Proliferation of connective tissue cells	3.03E‐05
	3. Cellular Development, Tissue Development	Differentiation of epithelial cells	4.05E‐05
	4. Carbohydrate Metabolism	Synthesis of polysaccharide	4.95E‐05
	5. Cancer, Cell Morphology, Organismal Injury and Abnormalities	Cell rounding of astrocytoma cells	1.98E‐03

This table shows the functional annotations of mRNAs for which m6A continuous peak‐over‐input (cPOI) scores ranked among the top 10% at the 5′‐untranslated region (UTR) and in the vicinity of the stop codon using Ingenuity Pathway Analysis for the three placentas from the appropriate‐for‐date (AFD) group.

Table 3

Molecular annotation of gene expression and m6A modification in healthy placentas using IPA

Molecular type(s)	Gene expression	m6A modification
		5′‐UTR	Vicinity of stop codon
	AFD_3 common (top10% 598genes)	AFD_3 common (top10% 64genes)	AFD_3 common (top10% 412genes)
Enzyme (%)	20.1	21.9	17.3
Transcription regulator (%)	6.0	15.6	22.7
Transporter (%)	9.5	10.9	4.6
Kinase (%)	2.7	4.7	5.1
Peptidase (%)	4.9	6.3	2.0
Transmembrane receptor (%)	3.0	1.6	2.9
Phosphatase (%)	1.7	3.1	2.0
G‐protein coupled receptor (%)	0.5	3.1	2.0
Growth factor (%)	0.8	1.6	0.2
Ion channel (%)	0.5	0.0	0.7
Ligand‐dependent nuclear receptor (%)	0.0	0.0	0.2
Translation regulator (%)	1.5	0.0	0.0
Cytokine (%)	0.5	0.0	0.0
MicroRNA (%)	1.2	0.0	0.0
Other (%)	47.1	31.3	40.2
Total	100	100	100

Molecular annotation with Ingenuity Pathway Analysis (IPA) was performed on the expressed gene and the gene subjected to m6A modification. Gene expression was considered as input in genes commonly expressed in three cases of appropriate‐for‐date (AFD) placentas as shown in Figure 1. Regarding the m6A‐modified gene, we annotated the genes (5ʹ‐untranslated region [UTR] and the vicinity of the stop codon for the top 10%, respectively) that were commonly modified in the previous three AFD cases. Moreover, numbers in parentheses () indicate candidate counts.

These results indicate that fine‐tuning of mRNA and protein levels of transcription regulator genes is achieved through mRNA m6A both at the 5′‐UTR and in the vicinity of the stop codon. This ability could eventually regulate various cell functions by controlling target gene activation in placental cells and in other tissues.

Potential m6A of placental mRNAs at the 5′‐UTR associated with fetal growth

Next, we aimed to detect fetal growth‐related m6A modifications in placental mRNAs. Three placentas of AFD fetuses, as well as placentas of six SFD and five HFD fetuses were collected, in addition to the previous three AFD placentas, to obtain a total of 17 placenta samples representing different birth weight groups (Table 1). The SFD group included three patients with PE. To investigate the relationship between mRNA m6A and fetal development, we calculated the cPOI scores for placental mRNAs, at the 5′‐UTR and in the vicinity of the stop codon, which can regulate the amount of RNA post‐transcriptionally. Metagene analysis was also conducted to check for variations among all the placenta samples. Indeed, placental m6A at the 5′‐UTR was more variable between individuals than did that at CDS and the 3′‐UTR (Figure 2A). The collective MeRIP RNA‐Seq data and metagene analysis data (data before normalization with input data) indicate that read enrichment around the initiation codon was more pronounced than that around the stop codon only in the placenta (Supplemental Figure S1B).

Figure 2

Characteristics of the 5′‐UTR and the vicinity of the stop codon of m6A‐modified placental mRNA. A, Metagene analysis of m6A modification site for all placenta samples. B, Correlation coefficients between gene expression (top) and m6A modification level using FPKM and the continuous peak‐over‐input (cPOI) scores at the 5′‐untranslated region (middle) or in the vicinity of the stop codon (bottom), respectively, among all placentas. FPKM values were processed with log10 and cPOI values were with log2. C, Heatmap and hierarchical clustering of m6A levels with all methylated transcripts among fetal birth weight categories. Hierarchical clustering was conducted for each birth weight category to generate heat maps based on the cPOI score at the 5′‐UTR and in the vicinity of the stop codon. The continuous POI (cPOI) scores were transformed to Log2 values after adding + 1 to original scores and are displayed as colors ranging from white to red as shown in the key. The range of cPOI score after logarithmic transformation in the color key was 0 to 8.5 at the 5′‐UTR and 0 to 11 in the vicinity of stop codon. The number of all the m6A‐modified mRNAs in each category was indicated in the top left of each heatmap

Next, RNA expression levels and m6A levels were compared among individual samples. Correlation tests showed that m6A levels showed much lower correlation between individuals than RNA expression did. Notably, m6A levels at the 5′‐UTR among individuals and birth weight groups had the lowest correlation (Figure 2B). Therefore, m6A modification levels at the 5′‐UTR may be used to more easily determine the characteristics of individual samples and fetal growth and obtain new information that could not be obtained using RNA expression analysis alone.

PCA analysis of gene expression patterns revealed unique patterns in placentas from SFD with PE cases, compared to those in other placentas (Supplemental Figure S2A, red circles); hence, samples from SFD fetuses were categorized as PE + and PE−, yielding four birth weight categories for the placenta samples: AFD, SFD PE+, SFD PE−, and HFD. The cPOI score of the 5′‐UTR and the domains in the vicinity of stop codon for each transcript were shown by heat map and hierarchically clustered for each of the three placenta categories, as compared with AFD group. Placentas from the same category localized in adjacent clusters, suggesting that each fetal growth category contains unique m6A at the 5′‐UTR (Figure 2C). Collectively, these results indicate that m6A at the 5′‐UTR might have some important features to understand the mechanism underlying PE and abnormal fetal growth.

Next, to identify the genes with significant differences in m6A levels in different birth weight groups, the Mann‐Whitney U test, a nonparametric test, was conducted using the mRNA cPOI scores for the two regions at the 5′‐UTR and in the vicinity of the stop codon for each category. Compared to the AFD group, all birth weight groups showed higher levels of m6A modifications at the 5′‐UTR (Figure 3A) but considerably lower levels of m6A modifications in the vicinity of stop codon (Figure 3B). The number of genes with significantly different m6A levels at the 5′‐UTR (DMGs) among birth weight category pairs AFD and SFD PE+, AFD and SFD PE−, and AFD and HFD were 86, 133, and 91, respectively and those in the vicinity of the stop codon were 384, 230, and 369 DMGs, respectively (Figure 3B). Notably, compared to the AFD group, the HFD group showed decreased mRNA methylation in 84% (309/369) of the genes; this proportion was markedly higher than that in the other two SFD groups (Figure 3B, right). The expression levels of genes involved in m6A (writer: METTL3/METTL14/WTAP; eraser: FTO/ALKBH5; reader: YTHDF1/YTHDF2) were determined for the 17 placenta samples used in this study; no significant differences were observed in their expression levels among all categories (Supplemental Table S1). However, the m6A eraser, FTO, is reported to be upregulated in the placentas of newborn infants with a heavy birth weight 31, 32; this is consistent with our results showing that the HFD samples had conspicuously decreased methylation. Further Cuffdiff‐based analysis to identify DEGs among the birth weight categories revealed virtually few overlaps with the aforementioned DMGs (Supplemental Figure S2B). Excluding the SFD PE + cases, comparisons of the AFD and SFD PE − and the AFD and HFD groups revealed that the mRNA levels of 95%‐97% of the DMGs were not significantly different. In other words, although gene expression levels were relatively similar, many candidate genes that were differentially methylated among all birth weight categories were identified. However, SFD PE + placentas had unique gene expression patterns compared to those in other placentas. The proportion of genes with differences in m6A levels owing to differences in gene expression levels was comparatively high in this category. Nevertheless, approximately 80% of the genes yielded results, not from differences in gene expression levels, but exclusively from differences in m6A modification levels (Supplemental Figure S2B). For instance, in the SFD PE + and AFD groups, the SAV1 and PTMS mRNAs were differentially methylated at the 5′‐UTR, but the gene expression levels of these two transcripts were not statistically different between the SFD PE + and AFD groups (Figure 3C and Supplemental Figure S3A). These two were the most differentially methylated genes (ranking top, Supplemental Data) between SFD PE + and AFD and differential methylated regions at the 5′‐UTR of both genes had “GGAC” motifs (Figure 3C left).

Figure 3

The difference of m6A modification has little relation to the difference of gene expression among fetus growth categories. A, Heatmap and hierarchical clustering of m6A levels with significant difference among fetal birth weight categories. mRNAs with significantly different m6A levels among fetal birth weight categories were analyzed using the Mann‐Whitney U test based on the continuous peak‐over‐input (cPOI) scores around both the 5′‐untranslated region (UTR) and in the vicinity of the stop codon (with a significance level of P < .05). The continuous POI (cPOI) scores were transformed to Log2 values after adding + 1 to original scores and are displayed as colors ranging from white to red as shown in the key. The range of cPOI score after logarithmic transformation in the color key was 0 to 8.5 at the 5′‐UTR and 0 to 11 in the vicinity of stop codon. B, Proportions of mRNAs with increased or decreased m6A modification levels for the regions at the 5′‐untranslated region (UTR) and in the vicinity of the stop codon in small‐for‐date preeclampsia (SFD PE+), SFD PE−, heavy‐for‐date (HFD), and appropriate‐for‐date (AFD) placentas are shown as bar charts. The Mann‐Whitney U test, a nonparametric test, was conducted using the cPOI score for each mRNA in each two region; at the 5′‐UTR and in the vicinity of the stop codon. The significance level was set to nominal P < .05. C, Two representatives of the most modified genes at the 5′‐UTR: SAV1 and PTMS. SFD PE + placenta showed considerably higher m6A modification at the 5′‐UTR but similar mRNA expression levels for SAV1 and PTMS. Left: Conceptual diagram of m6A modification at the 5′‐UTR of these genes. The numbers indicated the location of consensus motif “GGAC” from the translational start site; center: m6A modification level at the 5′‐UTR of these genes (Mann‐Whitney U test, P value = 1.1 × 10−2 [SAV1], 1.1 × 10−2 [PTMS]); right: mRNA expression of these genes (Student's t test, P value = .866 [SAV1], 0.476 [PTMS])

All whole comparisons between DMGs showed no significant GO terms. On the other hand, molecular annotation analysis using IPA indicated that the proportions of “transcription regulator” were higher in upregulated DMGs both at the 5′‐UTR (12.9%‐14.0%) and in the vicinity of stop codons (10.2%‐13.4%) than in upregulated DEGs (6.3%‐8.8%) (Supplemental Table S2). These data indicate that mRNA m6A regulation was not prominently involved in any specific placental biological function related to fetal growth. However, several placental functions could be controlled through post‐transcriptional gene regulation of “transcription regulator” in each fetal size category. Furthermore, methylation at 5′‐UTR of SFD PE − and demethylation in the vicinity of stop codons of HFD were representative.

SMPD1 might be involved in the pathophysiology of PE as its expression is regulated by m6A modification at the 5ʹ‐UTR

The etiology of PE is better understood than that of SFD and HFD. In the top Molecular and Cellular Functions annotated by IPA, the “homeostasis of sphingolipids” category related to the group of genes that exhibited changes in m6A levels in the SFD PE + group (Supplemental Table S3) included an increase in m6A levels of the 5′‐UTR of SMPD1 mRNA. SMPD1 is an enzyme involved in the synthesis of the sphingolipid ceramide; notably, its mRNA levels in PE placentas do not change, although the protein expression levels are known to increase.45 In the present study, we detected a significant increase in m6A levels in SMPD1 at the 5′‐UTR of the SFD PE + compared to AFD groups, although there was no significant difference in mRNA levels between them. SMPD1 ranked 76th out of the 86 transcripts that were differentially methylated at the 5′‐UTR between SFD PE + and AFD (Supplemental Data). A comparison of m6A density at the 5′‐UTR between SMPD1 and the average level of other genes in human placenta is shown in supplemental Table S4. As more than 90% of the transcripts were not methylated at the 5′‐UTR, m6A densities at the 5′‐UTR of SMPD1 (2.930‐3.395) were higher than the average m6A amount for all expressed transcripts (0.147‐0.209) in all the 17 placenta samples; however, when focusing only on the 5′‐UTR‐modified transcripts, m6A density of SMPD1 was almost equal to the average levels in all the 17 placenta samples (3.170‐3.289). Although the SMPD1 mRNA expression levels were not significantly different (P value = .76), the cPOI scores of SMPD1 at positions chr11: 6411775 and chr11: 6411825 were significantly increased in PE cases (Mann‐Whitney U test, P value = 4.7 × 10−2; Figure 4A‐C). Notably, the consensus motif “GGAC” repeated twice at the latter location, within nt −41 to −33 from the translational start site of SMPD1 mRNA (Figure 4A). Next, we aimed to verify the protein levels of SMPD1 in SFD PE + placentas. Unfortunately, the protein samples from the placentas used for MeRIP‐Seq could not be analyzed because of poor quality. Therefore, we newly collected five AFD and three SFD PE + placenta samples (Supplemental Table S5). SMPD1 RNA expression levels and levels of m6A at 5′‐UTR were simultaneously examined with protein levels in the new placenta samples. We confirmed that SMPD1 protein was significantly increased in the SFD PE + placentas (P value = 1.7 × 10−4) and m6A of SMPD1 5′‐UTR was also augmented in SFD PE+ (P value = 7.2 × 10−3) without any significant change in RNA expression (P value = .184; Figure 4D‐F).

Figure 4

Levels of SMPD1 protein, RNA expression, and m6A modification at the 5′‐UTR in human placenta. A, Conceptual diagram of m6A modification level around the 5′‐UTR of SMPD1 was plotted based on POI score in appropriate‐for‐date (AFD) and small‐for‐date preeclampsia (SFD PE+) specimens (Set1). B, The cPOI score of SMPD1 at the 5′‐UTR between two groups. (AFD 6 vs SFD PE + 3. Mann‐Whitney U test, P value = 4.7 × 10−2) (Set1). C, RNA expression levels of SMPD1 between two groups by FPKM value using Cufflinks (AFD 6 vs SFD PE + 3. Student's t test, P value = .76). The comparison of m6A modification between new AFD and SFD PE + placenta samples (Set 1). D, SMPD1 protein expression in AFD and SFD PE + placenta using western blotting. SMPD1 levels in SFD PE + were significantly higher than those in AFD (AFD 5 vs SFD PE + 3. Student's t test, P value = 1.7 × 10−4) (Set 2). E, m6A levels of SMPD1 at the 5′‐UTR in SFD PE + were significantly higher than those in AFD (AFD 5 vs SFD PE + 3. Student's t test, P value = 7.2 × 10−3) (Set 2). F, Gene expression level of SMPD1 in SFD PE + was similar to that in AFD (AFD 5 vs SFD PE + 3. Student's t test, P value = .184) (Set2). G, Relative luciferase activity and RNA expression after treatment with the transcription inhibitor actinomycin D (2 μg/mL) and at 2 and 8 h in JEG3 cells transfected with pGL4.53 with wild‐type or the m6A site mutation of SMPD1‐5′‐UTR (SMPD1_Wt and SMPD1_Mut). Data represent mean ± SE. Left: Luciferase reporter assay (n = 3, SMPD1_Wt vs SMPD1_Mut, Student's t test, P value = 3.6 × 10−3 [0 h], 1.9 × 10−3 [2 h], 2.1 × 10−3 [8 h]). Right: Luc2 RNA expression. (n = 3, SMPD1_Wt vs SMPD1_Mut, Student's t test, P value = .244 [0 h], 0.058 [2 h], 0.601 [8 h])

To confirm whether m6A at 5′‐UTR of SMPD1 can influence translation efficacy, we conducted luciferase reporter and mutagenesis assays. The replacement of “A” to “U” in the consensus motif “GGAC” sequence, which is repeated twice at the 3′ side of the 5′‐UTR of SMPD1, reduced the luciferase activity just after treatment with the transcription inhibitor actinomycin D (2 μg/mL) when inserted upstream of the firefly luciferase gene in the plasmid (Figure 4G, left). The mutation effect continuously changed luciferase activity at 2 and 8 hours of the actinomycin D treatment (Figure 4G, left). Compared to that in the JEG3 cells with wild‐type 5′‐UTR of SMPD1, the sequence‐mutated JEG3 cells showed a 78% decrease in luciferase activity (Figure 4G, left). On the other hand, RNA expression did not change significantly at all time comparisons (Figure 4G, right). These results suggest that the adenosines in the 5′‐UTR of SMPD1 are involved in translational efficiency. Collectively, our findings suggest that m6A modification of SMPD1 at the 5′‐UTR may influence protein translation post‐transcriptionally in both cultured cells and placenta.

We also identified genes other than SMPD1 that may be involved in certain illnesses of the placenta by regulating m6A modification of mRNAs. Notably, the involvement of these genes could not be identified by analyzing the RNA expression levels alone.

DISCUSSION

In this study, we profiled the placental epitranscriptome and observed the features of m6A modification by immunoprecipitation‐based MeRIP‐seq, which has not been previously accomplished. The m6A modification mechanism at the 5′‐UTR is different from that at the vicinity of the stop codon.2, 6, 7, 16 METTL3, METTL14, WTAP, and KIAA1429 are known members of the m6A writer complex.18 However, WTAP‐independent methylation sites only exist at the 5′‐UTR, and such methylations are reported to have a positive correlation with translation efficiency.17 m6A modification at the 5′‐UTR allows eIF4F to directly bind to m6A without the action of YTHDFs, which are m6A reader proteins, leading to cap‐independent ribosome recruitment that promotes translation initiation.46 In response to cellular stress, numerous additional domains show changes in m6A levels at the 5′‐UTR, compared to those in the vicinity of the stop codon.12, 16 Our results indicate that m6A at the 5′‐UTR might have some important features that will enable the elucidation of mechanisms underlying PE and abnormal fetal growth. Many reports focused on only methylated transcripts itself but did not focus on the methylated genic region. However, we precisely analyzed m6A alterations by genic regions. We showed that the increase in m6A levels at the 5′‐UTR of SMPD1 mRNA and of its protein levels were concomitantly observed in SFD PE + and found that disruption of m6A motif sequence at 5′‐UTR of SMPD1 significantly decreased the translation efficacy (Figure 4G). These results suggest that m6A levels at 5′‐UTR may be involved in the renewal of translation.45, 47 Especially in PE placentas, TGFβs are elevated.48, 49 It has been reported that the TGFβ family recruits the m6A writer protein complex to the mRNA via SMAD 2/3.50 Given this, to further verify the mechanism of regulation of m6A modification of SMPD1, we investigated the influence of TGFβ family on SMPD1 expression using a placenta‐derived cultured cell line. Interestingly, the results showed that TGFβ3 stimulation increased m6A levels at the 5′‐UTR of SMPD1 without any changes in SMPD1 mRNA levels (Supplemental Figure S4); however, TGFβ3 stimulation was not enough to markedly upregulate SMPD1 protein levels (data not shown). Thus, it may indicate that TGFβ3 could be partially involved in SMPD1 protein levels in PE placentas through m6A modification at the 5′‐UTR. We further revealed that the genes related to triacylglycerol levels had the strongest correlation with the 10 genes that exhibited increased m6A at the 5′‐UTR in HFD placentas, suggesting that aberrant lipid metabolism may be a characteristic of mothers of HFD infants (Supplemental Figure S3B,C). Notably, the mRNA expression levels of these genes of interest did not differ significantly (Supplemental Figure S3D). The observed increases in m6A at the 5′‐UTR of other genes may suggest that these modifications influence the translational efficiency of genes involved in illnesses in response to environmental conditions differing from the AFD conditions.

Various stresses such as oxidative stress and inflammation exist in the placentas with IUGR, such as SFD and PE.51, 52, 53, 54 It has been recently reported that one cellular response to hypoxia‐ischemia is mRNA stabilization and promotion of translation through m6A modification of the transcription factors JUN and MYC, demonstrating that mRNA degradation is not the only effect of m6A.55 We identified an increase in the proportion of transcriptional regulator genes accompanying changes in m6A levels, especially at the 5′‐UTR of upregulated DMGs among the different birth weight categories (Supplemental Table S2). Additionally, JUN is among the list of 5′‐UTR DMGs among these categories (Supplemental Data). Translational regulation via changes in m6A levels at the 5′‐UTR of transcription factor genes may be important in the environmental response of cells and in the regulation of placental function. In the future, using techniques such as first‐stage cell culture of patient‐derived villi,56 we intend to investigate details of the relationships between various stresses and m6A levels especially at the 5′‐UTR to further clarify the significance of m6A in placental mRNAs.

Several mRNAs showed different m6A levels but similar expression levels (Figure 3B). There was no difference in gene expression levels between AFD and SFD (both PE + and PE−) groups; however, YY1 and Notch2, which are known to have important functions in the placenta, were among the genes with significantly higher 5′‐UTR m6A levels in both the SFD groups, compared to those in the AFD group (Supplemental Data).57, 58 Conversely, FLT1, which is considered to have a high expression in PE placentas,59 also exhibited significantly higher expression in SFD PE + placentas in the present study (data not shown), than that in the AFD group; however, there was no difference in m6A modification at the 5′‐UTR or in the vicinity of the stop codon (Supplemental Data). Genes, including FLT1, that show extremely abnormal gene expression appear to be beyond the control of m6A; rather, the difference in their RNA expression likely contributes to the phenotype. Mice genetically altered to have a loss‐of‐function of either Fto or Alkbh5, which are m6A erasers, have no developmental defects except decreased birth weight and some male fertility impairments.60, 61 Moreover, new evidence suggested that m6A in Alkbh5 may partially influence postnatal mouse cerebellum development.62 However, demethylation of m6A does not appear necessary for proper fetus/placenta development. Conversely, a gene manipulation study in zebrafish, using mettl3, an m6A methylation enzyme, reported that m6A is necessary for development; in the study, homozygous mettl3 mutants were lethal to embryos.63 However, no studies have investigated whether m6A is necessary for human placental function; hence, this warrants further investigation in future studies.

In our experimental system, mRNA isoforms were analyzed without distinction. Even in an analysis considering the isoforms, our results were similar to the results depicted in Figure 2C (Supplemental Figure S5). However, the activity of each isoform is unknown. Furthermore, the effect of the primary isoform is a large one. Hence, for the purposes of this study, we did not consider isoforms in our analysis. We also found that approximately 1000 genes among long non‐coding RNAs underwent m6A modifications and numerous modifications in the CDS of many mRNAs (these data are available under accession number E‐MTAB‐6507 in ArrayExpress), which possibly affected splicing. These changes could also be related to placental functions and fetal development. Moreover, MeRIP‐seq was not cost‐effective, thereby forcing us to keep our specimen numbers small; therefore, there may be a lack of statistical power. In the future, we plan to increase the number of specimens or proceed with experiments using pure cells.

There is increasing evidence of differential gene expression between healthy and diseased individuals or disease models; however, in the present study, we identified differences in post‐transcriptional modification for multiple genes with unchanged gene expression between different fetal growth categories. Organs responsible for the pathology of diseases usually undergo many drastic changes. It is possible that many of these changes may be associated with genes showing normalized gene expression but different post‐transcriptional regulation. In such candidates, some genes may be found to play important roles in furthering our understanding of the pathology of diseases. Moreover, the concept of the epitranscriptome, including m6A modification, especially with respect to dividing the methylated genic region, is expected to be studied in various diseases. Our results strongly indicate that m6A modification of mRNA may be involved in placental function as well as in phenotypes and illnesses in various other tissues via translational efficiency.

CONFLICT OF INTEREST

The authors declare that they have no conflict of interest.

AUTHOR CONTRIBUTIONS

T. Kawai and K. Taniguchi designed research; K. Taniguchi and T. Kawai analyzed data; K. Taniguchi, T. Kawai, J. Kitawaki, J. Tomikawa, K. Nakabayashi, H. Sago, and K. Hata performed research; K. Taniguchi, T. Kawai, and K. Hata wrote the paper; K. Okamura and K. Nakabayashi contributed new reagents or analytic tools; K. Okamura developed software necessary to perform and record experiments.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.