Identifying key m6A-methylated lncRNAs and genes associated with neural tube defects via integrative MeRIP and RNA sequencing analyses

Objective: N6-methyladenosine (m6A) is a common post-transcriptional modification of messenger RNAs (mRNAs) and long non-coding RNAs (lncRNAs). However, m6A-modified lncRNAs are still largely unexplored. This study aimed to investigate differentially m6A-modified lncRNAs and genes involved in neural tube defect (NTD) development. Methods: Pregnant Kunming mice (9–10 weeks of age) were treated with retinoic acid to construct NTD models. m6A levels and methyltransferase-like 3 (METTL3) expression were evaluated in brain tissues of the NTD models. Methylated RNA immunoprecipitation sequencing (MeRIP-seq) and RNA sequencing (RNA-seq) were performed on the NovaSeq platform and Illumina HiSeq 2,500 platform, respectively. Differentially m6A-methylated differentially expressed lncRNAs (DElncRNAs) and differentially expressed genes (DEGs) were identified, followed by GO biological process and KEGG pathway functional enrichment analyses. Expression levels of several DElncRNAs and DEGs were evaluated by quantitative reverse transcription-polymerase chain reaction (qRT-PCR) for validation. Results: m6A levels and METTL3 expression levels were significantly lower in the brain tissues of the NTD mouse model than in controls. By integrating MeRIP-seq and RNA-seq data, 13 differentially m6A-methylated DElncRNAs and 170 differentially m6A-methylated DEGs were identified. They were significantly enriched in the Hippo signaling pathway and mannose-type O-glycan biosynthesis. The qRT-PCR results confirmed the decreased expression levels of lncRNAs, such as Mir100hg, Gm19265, Gm10544, and Malat1, and genes, such as Zfp236, Erc2, and Hmg20a, in the NTD group. Conclusion: METTL3-mediated m6A modifications may be involved in NTD development. In particular, decreased expression levels of Mir100hg, Gm19265, Gm10544, Malat1, Zfp236, Erc2, and Hmg20a may contribute to the development of NTD.


Introduction
Neural tube defects (NTDs) are common congenital abnormalities caused by the failure of the neural tube to close during embryogenesis (Yadav et al., 2021). The prevalence of NTDs is estimated to be 18.6 per 10,000 live births (Finnell et al., 2021). Babies with NTDs are more likely to be stillborn, die shortly after birth, or develop different degrees of disability (Huang et al., 2021). The etiology of NTDs is complex and is associated with interactions between genetic factors and diverse environmental factors (Avagliano et al., 2019;Kakebeen and Niswander, 2021). However, the molecular mechanisms underlying NTDs have not yet been fully elucidated.
An N 6 -methyladenosine (m 6 A) modification is a dynamic and reversible process modulated by methyltransferase "writers" (such as methyltransferase-like 3 (METTL3)) and demethylase "erasers" (such as alkB homolog 5 (ALKBH5)) (Zhang W. et al., 2021). m 6 A modifications are crucial for the regulation of RNA metabolism, including RNA stability, translation, alternative splicing, and translocation (Jiang et al., 2021). Moreover, m 6 A modifications have functions in embryonic development and neurodevelopmental diseases (Yen and Chen, 2021). m 6 A is the most prevalent messenger RNA (mRNA) and long non-coding RNA (lncRNA) modification (Tang et al., 2021). lncRNAs are a group of RNA transcripts longer than 200 nucleotides without open reading frames (Iyer et al., 2015). They have been implicated in the development of neurodevelopmental and neuropsychiatric disorders (Aliperti et al., 2021). Furthermore, m 6 A-modified lncRNAs are involved in various diseases. For instance, the lncRNA DNA methylation-deregulated and RNA m6A readercooperating lncRNA (DMDRMR) interacts with the m 6 A reader insulin-like growth factor 2 mRNA-binding protein 3 (IGF2BP3) to stabilize target genes, like cyclin-dependent kinase 4 (CDK4), in an m 6 A-dependent manner, thus exerting an oncogenic effect in clear cell renal cell carcinoma (Gu et al., 2021). m 6 A modifications of the lncRNA ZNFX1 Antisense RNA 1 (ZFAS1) and RAB22A, member RAS oncogene family (RAB22A) via METTL14 contributes to the development of atherosclerosis (Gong et al., 2021). Additionally, m 6 A-modified lncRNAs play pivotal roles in obstructive nephropathy (Liu P. et al., 2020), intervertebral disc degeneration (Li et al., 2022), and muscle development (Xie et al., 2021). However, the key m 6 A-modified lncRNAs and their regulatory mechanisms in NTD development have not yet been thoroughly investigated.
In the present study, we constructed a mouse model of NTD and investigated the overall m 6 A levels and METTL3 expression. We then performed m 6 A-modified RNA immunoprecipitation sequencing (MeRIP-seq) and RNA sequencing (RNA-seq) to compare brain tissues of NTD embryos and control embryos and conducted comprehensive bioinformatics analyses to identify key m 6 A-modified lncRNAs and genes associated with NTDs. Moreover, the expression levels of m 6 A-modified lncRNAs and genes were experimentally validated. These results are expected to improve our understanding of the molecular mechanisms underlying NTDs.

Animal models and samples
The Ethics Committee of the Third Affiliated Hospital of Guangzhou Medical University approved this study (2022-041, date of approval: 1 June 2022). Equal numbers of male and female Kunming mice (9-10 weeks) (Caven Biogle (Suzhou) Model Animal Research Co. Ltd., Suzhou, China) were mated overnight. The vaginal plug was examined the following day, and 25 pregnant mice were used for subsequent experiments. The day (08:00) a vaginal plug was observed was regarded as the embryonic day 0 (E0d), and 16:00 was considered E0.5d. NTD models were established as described previously (Yu et al., 2017). On E7.0d-7.25d, the mice in the NTD group (N = 18) were administered corn oil-dissolved retinoic acid (50 mg/kg of body weight) (R2625; Sigma, St. Louis, MO, USA) by one-time gavage. Mice in the control group (N = 7) were administered an equal amount of corn oil. On E16.5d, the mice were sacrificed by cervical dislocation, and embryos were taken from the uteri. NTDs were confirmed using a dissecting microscope. Brain tissues (anterior end of the neural tube) of NTD and control embryos were collected and frozen for storage.

m 6 A quantification
Total RNA was isolated from the brain tissues of the NTD and control groups using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). Using an m 6 A RNA Methylation Quantification Kit (Abcam, Cambridge, MA, USA), m 6 A levels were colorimetrically quantified by determining absorbance at 450 nm.
Quantitative reverse transcriptionpolymerase chain reaction (qRT-PCR) m 6 A methyltransferase METTL3 expression in brain tissues of the NTD and control groups was detected by real-time qRT-PCR. Total RNA was extracted from brain tissues using TRIzol reagent (Invitrogen). Reverse transcription for cDNA synthesis was performed using the PrimeScript ™ first strand cDNA Synthesis Frontiers in Genetics frontiersin.org Kit (Takara, Beijing, China). Real-time qRT-PCR was conducted using Power SYBR Green PCR Master Mix (Thermo Fisher Scientific, Waltham, MA, USA) and the 7900HT Fast qPCR System (Applied Biosystems, Foster City, CA, USA). Cycling conditions were as follows: initial denaturation at 95°C for 10 min and 40 cycles of 95°C for 15 s and 60°C for 60 s, followed by a melt curve analysis from 60°C to 95.0°C in increments of 0.5°C per 10 s. The internal control was glyceraldehyde-3-phosphate dehydrogenase (GAPDH), and the relative expression levels of lncRNAs and genes were calculated using the 2 −ΔΔCT method.
Methylated RNA immunoprecipitation (IP) sequencing (MeRIP-seq) and differential methylation analysis Total RNA was extracted from the brain tissues of three NTD embryos and three control embryos using TRIzol reagent (Invitrogen) and was treated with DNase I (Roche Diagnostics, Mannheim, Germany) to remove residual DNA. RNA was fragmented and immunoprecipitated with a mixture of beads and an anti-m 6 A antibody (Abcam) for 6 h at 4°C. The mixture was then immunoprecipitated by incubation with beads resuspended in IP reaction buffer (fragmented RNA, 5′ IP buffer, and RNasin Plus RNase Inhibitor) at 4°C for another 2 h. Then, the immunoprecipitated RNA was eluted and used for m 6 A MeRIP library construction using the SMARTer Stranded Total RNA-Seq Kit v2 (Pico Input Mammalian, Takara/ Clontech), following the manufacturer's protocols. Sequencing was performed on the NovaSeq platform. The raw data have been deposited in the NCBI Sequence Read Archive (SRA) database under accession number PRJNA879256.
Raw reads were filtered using Trimmomatic (v0.36) (Bolger et al., 2014), followed by a quality assessment to ensure the reliability of subsequent analyses. The data were aligned to the reference genome (mm10_gencode) using STAR (v2.5.2a) (Dobin et al., 2013). Uniquely mapped reads were used for subsequent analyses. Peak calling was used to detect regions significantly enriched in RNA (peaks), which were candidate m 6 A-methylated sites. Peak calling was performed using the MetPeak package (Cui et al., 2016) Differentially methylated sites on RNAs (m 6 A peaks) between NTD and control samples were analyzed using MeTDiff (Cui et al., 2018) in R based on MeRIP-seq data. The cutoff values were |fold change| > 1 and p < 0.05. RNAs with differentially methylated m 6 A sites were classified as mRNAs, lncRNAs, miRNAs, pseudogenes, or others to better understand the function of m 6 A methylation in various pathological processes. Moreover, the R package ChIPseeker (v1.24.0) (Yu et al., 2015) was used to annotate differentially methylated m 6 A sites and to analyze their distribution in functional regions according to their locations in RNA transcripts (i.e., 5′ untranslated region (UTR), 3′UTR, first exon, other exon, first intron, other intron, and distal intergenic regions). mRNAs with differentially methylated m 6 A sites in the 3′UTR region were analyzed.

RNA sequencing (RNA-seq)
Total RNA was isolated from the brain tissues of five NTD embryos and five control embryos using TRIzol reagent (Invitrogen), following RNA concentration detection using the NanoDrop 2000. Ribosomal RNA (rRNA) was removed, and the RNA was fragmented. The RNA library was established using the Illumina TruSeq RNA Sample Prep Kit and then loaded on an Illumina HiSeq 2,500 platform for 150 bp paired-end sequencing. The raw data have been deposited in the NCBI Sequence Read Archive (SRA) database under accession number PRJNA879256.
Raw reads were subjected to data filtration to remove lowquality reads. Read alignment to the reference genome (Mus_musculus.GRCm39) was conducted using HISAT2 (v2.1.0) (Kim et al., 2015). lncRNA and mRNA read counts were generated using RSEM (Li and Dewey, 2011), and mRNA or lncRNA expression levels were quantified.
To explore the key molecules in the pathogenesis of NTD, differentially expressed genes (DEGs) and differentially expressed lncRNAs (DElncRNAs) between NTD and control samples were identified using the DESeq2 package (v1.22.2) in R based on RNAseq data. The p-value was adjusted using the Benjamini-Hochberg (BH) method (Benjamini and Hochberg, 1995). The cutoff values were adjusted p < 0.05 and |fold change| > 1.
Integrative lncRNA/mRNA analysis with differentially methylated m 6 A sites and DElncRNAs/DEGs Data from MeRIP-seq and RNA-seq analyses were integrated. Then, m 6 A-methylated DElncRNAs and DEGs in which m 6 A levels were negatively correlated with expression levels were identified by analyzing the m 6 A levels of lncRNAs/ mRNAs with differentially methylated m 6 A sites and DElncRNA/DEG expression levels.

Construction of a lncRNA-mRNA coexpression network
The "cor" function in R was used to calculate the Pearson correlation coefficients (PCCs) between m 6 A-methylated DElncRNAs and DEGs. The false discovery rate (FDR)adjusted p-value was used for multiple testing correction. Pairs with a PCC >0.95 and FDR <0.05 were selected. The lncRNA-mRNA co-expression network was established using Cytoscape (version 3.6.1) (Shannon et al., 2003). The topological properties of this co-expression network, including node degree, Frontiers in Genetics frontiersin.org betweenness, and closeness, were analyzed using the CytoNCA plugin (version 2.1.6) (Tang et al., 2015).

Functional enrichment analysis
Gene Ontology (GO) biological process (BP) (The Gene Ontology Consortium, 2019) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000) pathway enrichment analyses were performed using the clusterProfiler package (version 3.16.0) to elucidate the functions of m 6 A-methylated DElncRNAs and DEGs in the co-expression network (Yu et al., 2012). A value of p < 0.05 (adjusted by the BH method) was considered significant.

Construction of a protein-protein interaction (PPI) network
Based on the STRING (version 11.0) database (Szklarczyk et al., 2015), the interactions between DEGs co-expressed with DElncRNAs were predicted, and a PPI network was constructed. The species was set as mice, and the PPI score was 0.15 (Yang et al., 2020;Zhao et al., 2022).

Validation of key DElncRNAs and DEGs using qRT-PCR
qRT-PCR was performed to detect the expression of several identified DElncRNAs and DEGs in the NTD and control groups to validate the reliability of the bioinformatics analysis. The DElncRNAs were Gm5165, Mir100hg, Gm19265, Gm10544, A730017L22Rik, and Malat1. The DEGs included Zfp236, Erc2, Nudcd3, and Hmg20a. qRT-PCR was conducted as described in 2.3.

Statistical analysis
All data are presented as the mean ± standard deviation (SD). The differences between the NTD and control groups were analyzed using Student's t-tests in GraphPad Prism 5.0 (GraphPad Software, San Diego, CA, USA). p < 0.05 was considered statistically significant.

Results
Overall m 6 A levels and METTL3 expression were decreased in NTDs We first established NTD mouse models by the intragastric administration of corn oil-dissolved retinoic acid. In the control group, the brains of mouse embryos were completely closed, the appearance was full and smooth, the optic vesicle and otic vesicle were visible, and the spinal surface was intact without laceration. In the NTD group, the rate of stillbirth increased. Some typical morphological malformations were observed, including cracks in the top of the skull, abnormalities of the hindbrain and face, and spina bifida ( Figure 1A). These data suggested that the retinoic acid-induced NTD mouse model was successfully established. We then analyzed overall m 6 A levels in NTD and control brain samples. Compared with levels in the control group, the NTD group had significantly lower m 6 A levels, indicating that methylation reactions could be compromised in NTD embryos ( Figure 1B). Furthermore, METTL3 expression in the NTD group was remarkably lower than that in the control group ( Figure 1C), indicating that METTL3 might be essential for m 6 A modification in NTDs.

Identification of differentially methylated m 6 A sites based on MeRIP-seq data
The NTD and control groups (N = 3 per group) were subjected to MeRIP-seq. Supplementary Table S1 provides a statistical summary of the raw and clean reads obtained by MeRIP-seq. After sequence alignment to the reference genome, the rates of uniquely mapped reads were higher than 70% (Supplementary Table S2), indicating that the data quality was sufficiently high for subsequent analyses. Based on MeRIPseq data, 468 hypomethylated m 6 A sites and 7,487 hypermethylated m 6 A sites were found in the NTD group compared to the control group (Figure 2A). According to RNA categories, we found that 50, 402, 3, and 13 hypomethylated m 6 A sites were located in lncRNAs, mRNAs, pseudogenes, and other regions, respectively ( Figure 2B), and 310, 7,128, 50, and 32 hypermethylated m 6 A sites were located in lncRNAs, mRNAs, pseudogenes, and other regions, respectively ( Figure 2C). Figure 2D shows the distribution of differentially methylated m 6 A sites in functional regions after peak annotation. When the mRNA 3′UTR had an m 6 A modification, 81 and 3,502 mRNAs had hypomethylated and hypermethylated m 6 A sites, respectively.

Screening DElncRNAs and DEGs based on RNA-seq data
The NTD and control groups (N = 5 per group) were subjected to RNA-seq. In general, 945, 603, 356 raw reads were detected using RNA-seq. Supplementary Table S3 shows a statistical summary of the raw and clean reads obtained from RNA-seq after quality control. After sequence alignment to the reference genome, the rates of uniquely mapped reads were higher than 84.97% (Supplementary Table S4), indicating that the data quality was Frontiers in Genetics frontiersin.org sufficiently high for subsequent analyses. After a differential expression analysis, 639 DElncRNAs (26 upregulated and 613 downregulated) and 1,132 DEGs (618 upregulated and 514 downregulated) were identified between the NTD and control groups ( Figure 3A). As shown in a heatmap, the identified DElncRNAs ( Figure 3B) and DEGs ( Figure 3C) could distinguish NTD samples from control samples.

Identification of m 6 A-related DElncRNAs and DEGs
Further integrative analyses of lncRNAs/mRNAs with differentially methylated m 6 A sites and DElncRNAs/DEGs yielded 13 differentially m 6 A-methylated DElncRNAs (corresponding to 20 transcripts) and 170 differentially m 6 A-methylated DEGs (corresponding to 201 transcripts).

Analysis of lncRNA-mRNA co-expression and PPI networks
Based on PCC scores, 171 co-expression pairs were obtained, including nine differentially m 6 A-methylated DElncRNAs and 58 differentially m 6 A-methylated DEGs. Figure 4A shows the co-expression network. By analyzing the topological properties of nodes in the co-expression network, it was observed that lncRNAs, such as Mir100hg, A730017L22Rik, Gm10544, Gm5165, and Gm19265, had higher degrees than those of DEGs. In addition, a PPI network was established based on interactions between differentially m 6 A-methylated DEGs and included 49 nodes and 131 interaction pairs ( Figure 4B). All DEGs in the network were downregulated in NTD.

Functional enrichment analyses of nodes in the lncRNA-mRNA co-expression network
Functional enrichment analyses of the differentially m 6 A-methylated DElncRNAs and DEGs were performed. In total, 1250 GO-BP terms and 23 KEGG pathways were significantly enriched for differentially m 6 A-methylated DElncRNAs. For instance, A730017L22Rik, Gm19265, Gm29260, Gm5165, and Mir100 hg were significantly enriched in developmental cell growth ( Figure 4C). Gm10545, Gm16096, Gm29260, and Mir100 hg were involved in mannosetype O-glycan biosynthesis. A730017L22Rik, Gm19265, and Gm5165 were involved in the Hippo signaling pathway ( Figure 4D). Moreover, the differentially m 6 A-methylated FIGURE 1 m 6 A levels and METTL3 expression were decreased in NTDs. (A) The dissecting microscope images for embryo in the control group and NTD groups on E16.5d. Arrows indicated the site of defects. (B) A quantitative m 6 A analysis was conducted to explore m 6 A enrichment in the brain tissues of NTD and control embryos (N = 5 per group). (C) qRT-PCR was carried out to investigate METTL3 expression in the brain tissues of NTD and control embryos (N = 3 per group). Compared to control group, *p < 0.05 and ***p < 0.001.

Frontiers in Genetics frontiersin.org
DEGs were associated with 125 GO-BP terms, including neuron migration ( Figure 4E), and one KEGG pathway, mannose-type O-glycan biosynthesis.

Validation by qRT-PCR
The expression levels of several DElncRNAs and DEGs were evaluated using qRT-PCR. The NTD group had dramatically lower expression levels of lncRNA, including Mir100hg, Gm19265, Gm10544, and Malat1, than those in the control group (p < 0.05, Figure 5A). Similarly, expression levels of key genes, such as Zfp236, Erc2, and Hmg20a, were lower in the NTD group than in the control group (p < 0.05, Figure 5B). These data were in line with the results of the bioinformatics analysis. There were no significant differences in Gm5165, A730017L22Rik, or Nudcd3 expression between the NTD and control groups. Discussion m 6 A modifications contribute to the regulation of gene expression during embryonic development (Li C. et al., 2021).
METTL3 has been identified as a critical methyltransferase in m 6 A modification and has functions in various biological processes (Liu S. et al., 2020). The METTL3-mediated m 6 A modification is essential for mammalian embryonic development (Sui et al., 2020). A previous study has revealed that levels of m 6 A modification and METTL3 expression are lower in ethionine-induced NTD than in controls (Zhang L. et al., 2021). Consistent with these findings, we found that the retinoic acid-induced NTD mouse model had dramatically decreased overall m 6 A levels and METTL3 expression levels than those in control mice. These findings suggested that METTL3mediated m 6 A modifications may be involved in NTD development. m 6 A modifications can modulate the expression of RNAs, including lncRNAs (Meyer and Jaffrey, 2014). There is increasing evidence that m 6 A-mediated epitranscriptomic changes can modulate lncRNAs in the developing cortex of mouse brains . We identified differentially m 6 A-methylated DElncRNAs associated with NTDs via an integrative analysis of MeRIP-seq and RNA-seq data, including Mir100hg, Gm19265, Gm10544, and Malat1. The neurogenic lncRNA Mir100 hg is a miR-125b and let-7 precursor. Both are implicated in neural development Frontiers in Genetics frontiersin.org (Rybak et al., 2008;Le et al., 2009). Malat1 is extremely abundant in brain tissues and is associated with neurological disorders, such as stroke, Alzheimer's disease, and retinal neurodegeneration Meng et al., 2019). In addition to lncRNAs, we identified m 6 A-methylated DEGs, including Zfp236, Erc2, and Hmg20a. Hmg20a is mainly expressed in hypothalamic astrocytes and is crucial for neuronal integrity preservation (Lorenzo et al., 2021). Erc2 participates in synaptic and neuronal functions (Lenihan et al., 2017). Nevertheless, the specific functions of Gm19265, Gm10544, and Zfp236 in neuronal development and nervous system-related diseases have not been reported.
Our results showed that these lncRNAs and genes are downregulated in NTD, implying that their dysregulation may be associated with NTD development. Further studies are needed to investigate the roles of these lncRNAs and genes in NTDs.
Notably, the differentially m 6 A-methylated DElncRNAs and DEGs were significantly enriched in multiple GO terms and KEGG pathways. O-Mannose-linked glycans are highly enriched in the brain and are vital for nervous system functions (Morise et al., 2014). They are also involved in brain development and remyelination (Gao et al., 2019). Our analysis revealed that Gm10545, Gm16096, Gm29260, and Mir100 hg were involved in mannose-type O-glycan biosynthesis, providing a mechanism by which Frontiers in Genetics frontiersin.org m 6 A-methylated DElncRNAs affected the development of NTDs. Moreover, A730017L22Rik, Gm19265, and Gm5165 were implicated in the Hippo signaling pathway, which has established roles in neuronal cell differentiation, neuroinflammation, and neuronal death (Cheng et al., 2020;Li X. et al., 2021). The Hippo/YAP signaling pathway has been implicated in the development of neurological diseases (Bao et al., 2017). Our findings revealed the potential role of the Hippo signaling pathway in NTD development. Moreover, the identified differentially m 6 A-methylated DEGs were significantly enriched in GO-BP terms, including neuron migration, and one KEGG pathway, mannose-type O-glycan biosynthesis. Neuronal migration is a fundamental step in nervous system formation (Stockinger et al., 2011). These pathways and GO functions provide insight into the roles of key m 6 A-methylated DElncRNAs and DEGs in the development of NTD.
In conclusion, our study is the first to identify differentially m 6 A-methylated DElncRNAs and DEGs associated with NTDs. Our findings demonstrated that METTL3-mediated m 6 A modifications may be involved in the development of NTD. In particular, decreased expression levels of lncRNAs,  qRT-PCR assay verified the lncRNA and gene expression levels in the brain tissues of NTD and control embryos. (A) Expression levels of lncRNAs, such as Mir100hg, Gm19265, Gm10544, and Malat1 (N = 3 per group). (B) Expression levels of genes, such as Zfp236, Erc2, and Hmg20a (N = 3 per group). Compared to control group, *p < 0.05, **p < 0.01, and ***p < 0.001.
Frontiers in Genetics frontiersin.org 09 such as Mir100hg, Gm19265, Gm10544, and Malat1, as well as key genes, such as Zfp236, Erc2, and Hmg20a, may be crucial in NTD development. However, the sample size was small, and the expression levels of the identified m 6 A-methylated lncRNAs and genes were not validated in neural tube cells isolated from fetal brains. Additional functional experiments are needed to reveal the regulatory mechanisms of key DElncRNAs and DEGs in NTDs.

Data availability statement
The raw data of MeRIP-seq and RNA-seq have been deposited in the NCBI Sequence Read Archive (SRA) database under accession number PRJNA879256, respectively.

Ethics statement
The animal study was reviewed and approved by The Animal Ethics Committee of the Third Affiliated Hospital of Guangzhou Medical University.