Integration of epigenomic and transcriptome analyses of neural tube defects reveals methylation driver lncRNAs and mRNAs

doi:10.21203/rs.3.rs-1511943/v1

Background: Abnormal genome-wide methylation during embryogenesis is associated with neural tube defects (NTDs) at birth. Long noncoding RNAs (lncRNAs) may be promising biomarkers for nervous system-related diseases. Therefore, we aimed to investigate the role of lncRNAs with aberrant methylation in the pathogenesis of NTDs.

Methods: Pregnant mice were given retinoic acid (dissolved in corn oil, 50 mg/kg) to build the NTDs model by gavage. After collecting brain tissues, reduced representation bisulfite sequencing (RRBS) and lncRNAs sequencing were conducted. Differentially expressed lncRNAs (DElncRNAs) and mRNAs (DEGs) between NTDs and control group were screened, and then integrated with RRBS data to obtain genes with aberrant methylation, followed by functional enrichment analysis. Subsequently, protein-protein interaction (PPI) network and lncRNA-miRNA-mRNA network were constructed. Finally, qRT-PCR was applied to determine the expression levels of identified hub lncRNAs.

Results: A total of 8 DElncRNAs as well as 213 DEGs with aberrant methylation between NTD group and normal group were screened. By bioinformatics analysis, several hub lncRNAs including Gm15521, Gm4681, Gm13974 and Gm40638, were identified. Function analysis showed these genes were mainly enriched in axon guidance pathway. The qRT-PCR assay revealed that the expression level of Gm15521, Gm4681 and Gm13974 in the NTDs group was significantly lower than those in the control group.

Conclusion: The study screened DElncRNAs with aberrant methylation in the NTDs and the identified genes could be potential biomarkers for prenatal diagnosis of NTDs. These findings will provide a reference for further study on the regulatory mechanism of non-coding RNAs in the NTDs.

Neural tube defects

DNA methylation

transcriptome

lncRNAs

Neural tube defects (NTDs) are a group of severe congenital central nervous system malformations that result from a failure in the morphogenetic process of neural tube closure during early embryonic development [1]. They include spina bifida, anencephaly, and encephalocele. As the second-most prevalent malformation, NTDs have been reported to affect approximately 300,000–4,000,000 infants worldwide each year, resulting in approximately 88,000 deaths and 8.6 million disability-adjusted life years [2]. Fetuses affected by NTDs are usually demised intrauterine or die shortly after birth, while surviving infants are often accompanied by varying degrees of disability [3]. Generally, NTDs are considered to be caused by a combined effect of environmental, genetic and metabolic factors, especially epigenetic alterations [4]. However, the pathogenesis of NTDs is not clear.

DNA methylation, one of the most important epigenetic mechanisms, is involved in chromatin organization and gene expression [5]. Previous studies have reported that aberrant genome-wide methylation during embryogenesis is associated with NTDs at birth. For example, hypomethylation of genomic DNA [6] and ZIC4 [7] may increase the risk of NTDs. Price et al. [8] also demonstrated that abnormal methylation of specific DNA in spina bifida may be related to neural tube defects. In addition, RNA-directed DNA methylation requires stepwise binding of silencing factors to long non-coding RNA (lncRNAs), which have been reported to play important roles in the regulation of DNA methylation [9]. For example, H19 lncRNA alters genome-wide DNA methylation by regulating S-adenosylhomocysteine hydrolase [10]. HOTAIR lncRNA inhibits miR-122 expression in hepatocellular carcinoma via DNA methylation [11]. However, the abnormal DNA methylation-mediated lncRNAs involved in NTDs pathogenesis have not been identified.

The development of high-throughput sequencing techniques for DNA methylation and transcriptome analysis make genome-wide detection of methylation driving events feasible [12]. Meanwhile, omics technologies offer new possibilities for prenatal diagnosis of NTDs [13]. In this study, after construction of NTDs mice model, representation bisulfite sequencing (RRBS) and lncRNAs sequencing was performed. Based on the integration analysis of these data, the differentially expressed lncRNAs (DElncRNAs) and mRNA (DEGs) between NTDs and control groups were screened. DElncRNAs or DEGs related to abnormal methylation were also identified. Moreover, lncRNA-mRNA network, protein-protein interaction (PPI) network and lncRNA-miRNA-mRNA network were constructed. The study is aimed to screen the lncRNAs associated with epigenetic alternations in NTDs to look for new insights into the epigenetic regulation of lncRNA expression by DNA methylation in the pathogenesis of NTDs.

Animal models and samples

Eighty Kunming mice (9–10 weeks of gestation) were purchased from Caven Biogle (Suzhou) Model Animal Research Co. Ltd and then included in this study. They were placed in cages with equal proportions of male and female at 18:00, and then the animals were mated overnight. The vaginal plug was examined the next morning. Female mice with a vaginal plug at 8:00 were designed as embryonic days (E0d), while at 16:00 were regarded as E0.5d. Finally, 30 pregnant mice were obtained, which were divided into NTDs and control groups, with 15 mice in each group.

NTDs model was established as described in the previous study [14]. On E7.0d-7.25d, the mice in the NTDs group were administered with 50 mg/kg (body weight) of retinoic acid (RA, sigma, R2625) dissolved in corn oil by a one-time gavage and those in the control group was given an equal amount of corn oil. On E16.5d, pregnant mice were killed by cervical dislocation. Their embryos were dissected out from uteri and were examined with a dissecting microscope to confirm whether NTDs had occurred. Meanwhile, brain tissues (anterior end of the neural tube) of control embryos and of NTDs embryos (n = 5 per group) were collected and then stored in frozen for further sequencing. The animal experiments were approved by the Animal Ethics Committee of Kunming Medical University and followed the Institutional Animal Care and Use guidelines.

LncRNA sequencing

Total RNA was extracted from tissues by using the TRIzol reagent (Invitrogen, CA, USA). The concentration and purity of RNA were detected by using Nanodrop2000. After removing rRNA, the RNA library was constructed using the RNA Library Prep Kit for Illumina® (NEB, UK) and then sequenced on the Illumina Hi seq 2500 platform (Illumina, San Diego, CA, USA) with 2×150 bp paired-end read to obtain raw data.

To ensure the accuracy of further bioinformatics analysis, the raw reads were filtered to obtain clean reads via deleting the low-quality reads. Then, clean reads were mapped to the Mus musculus (GRCm39) reference genome by using HISAT2 software [15]. Read counts of lncRNAs and mRNA were generated by using RSEM, and the expression levels of them were normalized using the fragments per kilobase of transcript per million mapped reads (FPKM).

RRBS analysis

Genomic DNA was extracted from each sample and digested with the methylation-insensitive enzyme MspI. DNA fragments sized 150-300 bp were selected with bisulfite conversion and PCR amplification to obtain DNA library, followed by sequencing on the Illumina Hiseq 2500 Platform.

To obtain clean reads, we used the Trimming method to remove the sequencing adapters and low-quality reads of the raw data. Then, clean reads were mapped with the Mus musculus (mm10) reference genome by using Bsmap software which was also applied to detect the methylation sites.

Differential expression analysis

To further explore the genes or methylation involved in the pathogenesis of NTDs, a series of differential expression analyses between the NTDs and the control group were performed. Based on the RNA sequencing data, DESeq2 software was applied to screen the DEGs and DElncRNAs within NTDs vs. control, with the criterion of an adjusted p-value < 0.1 and |log2FoldChange (FC)| > 1. Of which, p value was adjusted with the negative binomial distribution test. In addition, differentially methylated regions (DMR) between NTDs and control samples were also identified based on the RRBS data using metilene software. DMRs with adjusted p-value<0.05 were considered statistically significant.

Integrative analysis of DMR and DElncRNAs/DEGs

The data from RNA sequencing and RRBS sequencing were integrated, and genes annotated in DMR anchor regions (such as Promoter, Exon, Intron, CGI, CGI Shore, Repeat, TSS, TES, etc.) were obtained. Next, these genes were compared with DEGs and DElncRNAs to identify methylation-driven DEGs and DElncRNAs.

Construction of DElncRNA-DEG co-expression network

Pearson correlation coefficient (PCC) between methylation-driven DEGs and lncRNA was calculated by using cor function. Pairs with PCC higher than 0.7 and p value less than 0.05 were selected. Then, co-expression network was constructed via Cytoscape (version 3.6.1) [16]. Meanwhile, the topological properties of nodes in this network were analyzed by using CytoNCA (version 2.1.6) [17].

Functional enrichment analysis

To further explore the biological function of nodes in the co-expression network, we used the clusterprofiler package (version 3.16.0) [18] to assess the GO-BP terms and KEGG pathways of DEGs as well as DElncRNAs. An adjusted p value < 0.05 was considered as significant enrichment.

Construction of protein-protein interaction (PPI) network

The relationship between target genes of lncRNAs was predicted by using STRING software(version 11.0) [19]. In brief, the species was set as mice and PPIs with a combined score > 0.15 were screened to generate a PPI network.

Prediction of miRNA

miRWalk 2.0 was applied to predict the miRNA in the DElncRNA-DEG co-expression network[20]. miRNA-DEG interactions were verified using seven databases: miRWalk, Microt4, miRanda, miRDB, miRNAMap, RNA22 and Targetscan. Pairs appearing in at least five databases were selected to construct miRNA-DEG network.

Construction of lncRNA-miRNA-mRNA network

The miRNAs-related lncRNAs were predicted using DIANA-LncBase v.2 database[21], and the lncRNAs appearing in the DElncRNA-DEG network were screened to obtain DElncRNA-miRNA pairs. Based on the identified DElncRNA-miRNA and miRNA-DEG pairs, the lncRNA-miRNA-mRNA network was constructed.

qRT-PCR validation

To validate the accuracy of bioinformatics analysis, the expression levels of several identified lncRNAs (Gm15521, Gm4681 and Gm13974) were evaluated using the qRT-PCR. Total RNA was isolated from tissue samples using TRIzol® reagent (Invitrogen, USA) and reverse transcribed to cDNA using PrimeScript™ 1st Strand cDNA Synthesis kit (Takara, China). qRT-PCR was carried out by using TaKaRa Ex Taq® (Takara, China) on the 7900HT Fast qPCR System. The relative mRNA expression of lncRNAs was normalized using GAPDH with the method of 2^−ΔΔCt. The specific primers used in qRT-PCR are listed in Table S1.

Statistical analysis

Statistics analysis was performed using GraphPad Prism 5.0 (Graphpad Software, San Diego, CA) and all data were presented as mean ± SD. The differences between two groups were compared by the Student’s paired t-test. The p value < 0.05 was considered statistically significant.

Establishment of NTD mouse models

The brain of mouse embryos in the control group was completely closed, with a full and smooth appearance, visible optic and otic vesicles, and intact spinal surface without laceration (Fig. 1A). In the NTDs group, the stillbirths of mouse embryos were increased; some typical morphological malformations were observed, including cracks in the top of the skull, abnormality of the hindbrain and face, and spina bifida (Fig. 1B).

Results of lncRNA sequencing and RRBS

A total of 945,603,356 raw reads were detected using lncRNA sequencing. After quality control, 939221946 clean reads were obtained, including 476,257,370 reads in the control group and 462,964,576 reads in the NTD group (Table S2). Furthermore, 95.23% of the reads matched the Mus musculus reference genome, indicating that our data were of high quality and suitable for subsequent analysis.

In the RRBS, 188,563,561 raw data were obtained and 182,999,476 clean reads remained after filtration (90,565,820 reads in the control group and 92,433,656 reads in the NTD group, Table S3). In the control group, the percentages of CG, CHG, and CHH were 47.57%, 0.13%, and 0.05%, respectively. In the NTDs group, the distribution of C sites was 47.15% CG, 0.13% CHG, and 0.05% CHH, respectively.

Identification of methylation-driven DEGs and DElncRNAs

After differential expression analysis, 639 (26 upregulated and 613 downregulated) DElncRNAs and 1132 (618 upregulated and 514 downregulated) DEGs between the NTD and control groups were screened (Fig. 2A). The heatmap for DElncRNAs and DEGs showed that these genes could significantly distinguish the NTD group from the control group (Fig. 2B and 2C). After integration with RRBS data, 8 DElncRNAs and 213 DEGs with abnormal methylation levels were identified. The top five DElncRNAs and DEGs are listed in Table 1.

Table 1

Top five DElncRNAs and DEGs with abnormal methylation level.
	Gene ID	Gene name	Log2FC	P value	Padjust	Regulate
lncRNA	ENSMUSG00000113105	Gm47371	-1.16756	8.29E-16	2.62E-12	down
	ENSMUSG00000115536	Gm32857	-1.23758	1.57E-15	3.73E-12	down
	ENSMUSG00000104512	Gm37165	-1.7099	1.78E-15	4.00E-12	down
	ENSMUSG00000098066	Gm26944	-23.4564	9.22E-15	1.37E-11	down
	ENSMUSG00000093675	Gm20618	-1.45402	2.15E-14	2.60E-11	down
mRNA	ENSMUSG00000053838	Nudcd3	-1.04612	1.43E-32	6.79E-28	down
	ENSMUSG00000035696	Rnf38	-2.22203	2.05E-21	3.24E-17	down
	ENSMUSG00000041258	Zfp236	-9.70547	3.20E-20	3.79E-16	down
	ENSMUSG00000017639	Rab11fip4	-27.541	8.97E-20	8.49E-16	down
	ENSMUSG00000021144	Mta1	-25.9946	8.74E-18	5.91E-14	down

LncRNA–mRNA co-expression network

According to the PCC score between DElncRNAs and DEGs, 99 co-expression pairs were obtained, including 6 DElncRNAs and 28 DEGs. The co-expression network is shown in Fig. 3. The degree of each node was calculated. Nodes with a degree greater than 5 are listed in Table 2. We observed that lncRNA Gm4681 had the highest degree, followed by lncRNAs 4933406K04Rik, Gm13974, Gm40638, and Gm15521.

Table 2

The nodes with degree larger than 5 in the co-expression network.
Name	Degree	Betweenness	Closeness	Log2FC	Padjust	Regulate	type
Gm4681	23	280.2948	0.297297	-1.4828	3.50E-09	down	lncRNA
4933406K04Rik	20	167.5543	0.282051	-1.1272	0.022408	down	lncRNA
Gm13974	20	131.6011	0.282051	-1.03811	9.26E-04	down	lncRNA
Gm40638	19	186.1886	0.277311	-1.28794	1.21E-05	down	lncRNA
Gm15521	16	56.36119	0.264	-1.24576	2.00E-04	down	lncRNA
Usp46	5	8.605687	0.264	-1.73308	0.008168	down	mRNA
Trerf1	5	8.605687	0.264	-1.17501	0.014953	down	mRNA
Tdrd9	5	8.605687	0.264	-1.17423	3.14E-04	down	mRNA
Syn3	5	8.605687	0.264	-1.37689	8.94E-04	down	mRNA
Styx	5	8.605687	0.264	-1.11316	0.002523	down	mRNA
Rpgr	5	8.605687	0.264	-1.94071	4.18E-06	down	mRNA
Ppp1r9a	5	8.605687	0.264	-1.54827	6.20E-06	down	mRNA
Lrp1	5	8.605687	0.264	-1.05218	0.026977	down	mRNA
Klhl29	5	8.605687	0.264	-1.53648	6.98E-04	down	mRNA
Dapk1	5	8.605687	0.264	-1.95977	0.030613	down	mRNA
Cdk8	5	8.605687	0.264	-1.18471	0.047516	down	mRNA
Camta1	5	8.605687	0.264	-1.27631	0.003395	down	mRNA
Auts2	5	8.605687	0.264	-1.58481	0.009469	down	mRNA

Functional enrichment analysis of nodes in the co-expression network

To further explore the biological functions of the DElncRNAs and DEGs, functional enrichment analysis was conducted. The results showed that these DElncRNAs were significantly enriched in 1409 GO-BP terms and 19 KEGG pathways. In particular, 4933406K04Rik, Gm15521, Gm4681, Gm13974, and Gm40638 were mainly enriched in several GO-BP terms, such as positive regulation of cell projection organization, righting reflex, and reflex (Fig. 4A), whereas 4933406K04Rik, Gm13974, Gm15521, and Gm40638 were primarily involved in the axon guidance pathway (Fig. 4B).

For the lncRNA target genes, 327 GO-BP terms and four KEGG pathways were identified. We observed that these genes were significantly enriched in GO-BP terms, such as positive regulation of cell projection organization (Fig. 4C), and KEGG pathways, such as axon guidance (Fig. 4D).

PPI network analysis

A PPI network was constructed based on the PPI pairs between DEGs. As illustrated in Fig. 5, the PPI network was composed of 22 nodes and 25 edges. All DEGs were downregulated.

Establishment of lncRNA-miRNA-mRNA regulatory network

A total of 584 miRNA-mRNA pairs (408 miRNAs and 22 mRNAs) were predicted. Four common DElncRNAs and 70 miRNAs were identified in this study. After integrating the lncRNA–miRNA and miRNA–mRNA pairs, a lncRNA-miRNA-mRNA network containing 91 nodes (70 miRNAs, 4 lncRNAs, and 17 mRNA) with 237 edges was constructed (Fig. 6). Nodes with a degree greater than 10 are summarized in Table 3. It was demonstrated that lncRNA 4933406K04Rik had the highest degree in the network, followed by lncRNAs Gm13974 and Gm15521. These nodes may be hub abnormally methylated genes involved in the pathogenesis of NTDs.

Table 3

The nodes with degree greater than 10 in ceRNA network
name	Degree	Betweenness	Closeness	Log2FC	Padjust	Regulate	type
4933406K04Rik	44	2791.1	0.647482	-1.1272	0.022408	down	lncRNA
Gm13974	43	2495.893	0.625	-1.03811	9.26E-04	down	lncRNA
Gm15521	23	845.0656	0.494505	-1.24576	2.00E-04	down	lncRNA
Ppp1r9a	22	673.2955	0.552147	-1.54827	6.20E-06	down	mRNA
Sema5a	19	515.323	0.471204	-1.11811	0.0155	down	mRNA
Itsn1	18	545.051	0.514286	-1.4752	7.32E-06	down	mRNA
Klhl29	17	779.9576	0.542169	-1.53648	6.98E-04	down	mRNA
Zfhx3	16	269.1233	0.454545	-1.09265	8.83E-09	down	mRNA
Magi1	11	239.4414	0.481283	-1.28005	0.001264	down	mRNA
Camta1	10	133.242	0.511364	-1.27631	0.003395	down	mRNA

Validation analysis by qRT-PCR

To validate the accuracy of the sequencing analysis results, the expression levels of the hub lncRNAs were evaluated using qRT-PCR. The results indicated that the expression levels of lncRNAs Gm15521, Gm4681, and Gm13974 were significantly lower in the NTD group than in the control group (p < 0.05; Fig. 7), which was consistent with the results of bioinformatics analysis.

Animal models and samples

Eighty Kunming mice (9–10 weeks old) were purchased from Caven Biogle (Suzhou) Model Animal Research Co. Ltd (Suzhou, Jiangsu, China). The mice were placed in cages with equal proportions of males and females at 18:00, and then mated overnight. The vaginal plug was examined the following morning. Female mice with vaginal plugs at 8:00 were designated as embryonic days (E0d), whereas those at 16:00 were regarded as E0.5d. Finally, 30 pregnant mice were obtained, which were divided into NTD and control groups, with 15 mice in each group.

The NTD model was established as described previously [14]. On E7.0d-7.25d, mice in the NTD group received a one-time administration of 50 mg/kg (body weight) retinoic acid (RA) (R2625; Sigma) dissolved in corn oil by gavage, whereas the control group was administered an equal amount of corn oil. On E16.5d, the pregnant mice were killed by cervical dislocation. Embryos were dissected from the uteri and examined using a dissecting microscope to confirm the occurrence of NTDs. The brain tissues (the anterior end of the neural tube) of control and NTD embryos (n = 5 per group) were collected and frozen at -80°C until further sequencing. The animal experiments were approved by the Animal Ethics Committee of Kunming Medical University, and were conducted in accordance with the guidelines of Institutional Animal Care and Use.

LncRNA sequencing

Total RNA was extracted from each tissue using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), and the concentration and purity of RNA were detected using Nanodrop2000 (Thermo Fisher Scientific, USA). After rRNA was removed, an RNA library was constructed using the RNA Library Prep Kit for Illumina® (NEB, UK) and then sequenced on the Illumina HiSeq 2500 platform with 2×150 bp paired-end reads to obtain raw data.

To ensure the accuracy of further bioinformatics analysis, the raw reads were filtered by deleting low-quality reads to obtain clean reads. Next, the clean reads were mapped to the Mus musculus (GRCm39) reference genome using the HISAT2 software [15]. Read counts of lncRNAs and mRNA were generated using RSEM, and their expression levels were normalized using fragments per kilobase of transcript per million mapped reads (FPKM).

RRBS analysis

Genomic DNA was extracted from each sample and digested with the methylation-insensitive enzyme MspI. DNA fragments sized 150–300 bp were selected by bisulfite conversion and PCR amplification to obtain a DNA library, followed by sequencing on the Illumina HiSeq 2500 Platform (Illumina, San Diego, CA, USA).

To obtain clean reads, we used the trimming method to remove sequencing adapters and low-quality reads from raw data. Next, the clean reads were mapped to the Mus musculus (mm10) reference genome using the Bsmap software. Finally, methylation sites were detected using the Bsmap software.

Differential expression analysis

To further explore the genes or methylation involved in the pathogenesis of NTDs, a series of differential expression analyses between NTD and control samples was performed. Based on RNA sequencing data, DEGs and DElncRNAs between NTD vs. control samples were screened using the DESeq2 software, with the criteria of an adjusted p-value < 0.1 and |log2FoldChange (FC)| > 1. The p-value was adjusted using the negative binomial distribution test. In addition, differentially methylated regions (DMR) between NTD and control samples were identified based on RRBS data using the Metilene software. DMRs with an adjusted p-value < 0.05 were considered significant.

Integrative analysis of DMR and DElncRNAs/DEGs

Data from RNA sequencing and RRBS sequencing analyses were integrated, and genes annotated in DMR anchor regions (such as Promoter, Exon, Intron, CGI, CGI Shore, Repeat, TSS, and TES) were obtained. Next, these genes were compared with DEGs and DElncRNAs to identify methylation-driven DEGs and DElncRNAs, respectively.

Construction of DElncRNA-DEG co-expression network

The Pearson correlation coefficient (PCC) between methylation-driven DEGs and lncRNAs was calculated using the cor function. Pairs with PCC > 0.7 and p-values < 0.05 were selected. A co-expression network was constructed via Cytoscape (version 3.6.1) [16]. The topological properties of nodes in this network were analyzed using CytoNCA (version 2.1.6) [17].

Functional enrichment analysis

To further explore the biological function of nodes in the co-expression network, we used the clusterprofiler package (version 3.16.0) [18] to assess the GO-BP terms and KEGG pathways of DEGs and DElncRNAs. DEGs and DElncRNAs were considered significantly enriched if the adjusted p-value < 0.05.

Construction of protein-protein interaction (PPI) network

The relationship between the target genes of the lncRNAs was predicted using the STRING (version 11.0) software [19]. In brief, the species was set as mice, and PPIs with a combined score > 0.15 were screened to generate a PPI network.

Prediction of miRNA

Furthermore, miRWalk 2.0 [20] was used to predict miRNAs in the DElncRNA–DEG co-expression network. miRNA-DEG interactions were verified using seven databases: miRWalk, Microt4, miRanda, miRDB, miRNAMap, RNA22, and Targetscan. Pairs appearing in at least five databases were selected to construct the miRNA-DEG network.

Construction of lncRNA–miRNA–mRNA network

In addition, miRNA-related lncRNAs were predicted using the DIANA-LncBase v.2 [21] database, and the lncRNAs appearing in the DElncRNA–DEG network were screened to obtain DElncRNA–miRNA pairs. Based on the identified DElncRNA–miRNA and miRNA–DEG pairs, a lncRNA–miRNA–mRNA network was constructed.

qRT-PCR validation

To validate the accuracy of the bioinformatics analysis results, the expression levels of several identified lncRNAs (Gm15521, Gm4681, and Gm13974) were evaluated using qRT-PCR. Total RNA was isolated from tissue samples using TRIzol® reagent (Invitrogen) and reverse-transcribed to cDNA using the PrimeScript™ 1st Strand cDNA Synthesis kit (Takara, China). qRT-PCR was carried out using TaKaRa Ex Taq® (Takara) on the 7900HT Fast qPCR System. The relative mRNA expression of lncRNAs was normalized to that of GAPDH using the 2^−ΔΔCt method. The primers used for qRT-PCR are listed in Table S1.

As the most common human birth defects, NTDs may lead to severe neurological handicaps and even death. Despite recent advances in elucidating the etiology of NTDs, lncRNAs involved in the pathogenesis of NTDs remain unidentified, especially genes associated with DNA methylation. In the present study, RRBS and lncRNAs sequencing was performed after the construction of the NTDs mouse model. Several hub lncRNAs (such as 4933406K04Rik, Gm15521, Gm4681, Gm13974, and Gm40638) with abnormal methylation were screened via a series of bioinformatics analyses. Enrichment analysis showed the lncRNAs were mainly enriched in positive regulation of cell projection organization (GO-BP term) and axon guidance (KEGG pathway). We also explored the relationship between lncRNA, miRNA, and mRNA via the construction of the lncRNA-miRNA-mRNA network. Furthermore, experimental verification indicated that the expression levels of lncRNAs Gm15521, Gm4681 and Gm13974 were significantly down-regulated in the NTDs group.

LncRNAs have the potential to identify characteristics and hallmarks of various diseases [22]. LncRNA-dependent mechanisms regulate transcription and CpG DNA methylation. Methylation modification modulates lncRNA function and stability by altering local RNA structure and affecting the nervous system [23, 24]. In this study, several NTDs-related lncRNAs and mRNAs with abnormal methylation were identified, including Gm47371, Gm32857, Nudcd3, and Rnf38. Nudcd3 was reported to be associated with diabetes mellitus-related atherogenesis [25]. Rnf38 was dysregulated in various cancers and might be a potential therapeutic target, such as osteosarcoma [26], colorectal cancer [27] and cervical carcinoma [28]. However, these hub genes have not been reported in central nervous system diseases, including NTDs. Hence, the role of these genes in the pathogenesis of NTDs needs to be further investigated.

To further elucidate the biological role and possible mechanisms of lncRNAs and mRNAs in NTDs pathogenesis, GO-BP and KEGG pathway analyses were conducted. We observed that both lncRNA and mRNA were enriched in GO-BP terms such as positive regulation of cell projection organization and KEGG pathways such as axon guidance. Axon guidance is a key process in forming the connection between neurons and targets [29]. In brief, several axon guidance-related molecules serve a crucial role in the regulation of neuronal connections by inducing cytoskeletal rearrangements during development [30]. Previous studies indicated that axon guidance mechanism provided potential contribution to neural dysfunction [31]. Ramya et al. [32] also revealed that axon guidance pathway was critical for brain development in neural stem cells exposed to different treatments. All these studies indicated that these key genes might involve in the pathogenesis of NTDs through regulating several pathways, such as axon guidance.

The regulatory axis and connectivity network of mRNA-miRNA-lncRNA play an active role in various pathological conditions, such as neurodegenerative diseases [33]. In this study, we predicted the miRNA-lncRNA and miRNA-mRNA pairs, followed by the construction of a regulatory network, consisting 70 miRNAs, 4 lncRNAs and 17 mRNAs. In this network, key lncRNAs and mRNAs connected with more nodes were screened. Ppp1r9a had the highest degree value, followed by Sema5a. Ppp1r9a is reported to play an important role in the formation of neurites and the movement of the spine [34]. The highly connected neighbor genes of Ppp1r9a had de novo mutations in individuals with neuropsychiatric disorders, indicating it may be regarded as a novel candidate gene for this disease [35]. Sema5a is overexpressed in the early stages of brain development, which is involved in axonal, neuronal and synaptic-related processes [36]. Previous studies showed that Sema5a was associated with autism spectrum disorder and mental health development disorders [37]. We also found Ppp1r9a targeted by lncRNAs Gm15521 and Gm13974; Sema5a targeted by lncRNA 4933406K04Rik. The qRT-PCR showed that the expression level of Gm15521 and Gm13974 was lower in the NTDs sample than those in the control sample. Nevertheless, the specific function of these lncRNAs and mRNAs has not been reported in the NTDs pathogenesis.

This study screened lncRNAs and mRNAs with methylation alteration in NTDs pathogenesis. There were some limitations. Our findings are obtained via bioinformatics analysis, and these results are preliminary and lack additional functional validation experiments. In addition, there are few studies on NTDs-related lncRNA sequencing, we still need to explore key molecules in the regulatory network to elucidate the pathogenesis of NTDs.

In summary, our study demonstrated that several hub lncRNAs, especially Gm15521, Gm4681, Gm13974 and Gm40638, were differentially expressed and with abnormal methylation in the NTDs mouse. These identified genes may have important roles in the pathogenesis of NTDs. However, further studies are needed to validate our findings.

NTDs: neural tube defects

lncRNAs: long noncoding RNAs

RRBS: reduced representation bisulfite sequencing

DElncRNAs: differentially expressed lncRNAs

DEGs: differentially expressed mRNAs

PPI: protein–protein interaction

FC: fold change

PCC: Pearson correlation coefficient

Ethics approval and consent to participate

All animal experiments were approved by the ethics committee of Kunming Medical University (approval number: KMMU2020446) and were conducted in accordance with the guidelines of Institutional Animal Care and Use.

Consent for publication

Not applicable.

Availability of data and materials

The datasets used and analyzed during the current study are available from the first author on reasonable request.

Competing interests

The authors declare that they have no competing interests.

Funding

The work was supported by the Yunnan Province Science and Technology Program (No. 202101AY070001-121) and Guangzhou Science and Technology Program (No. 202102010129).

Authors’ contributions

Jing Xu, Min Chen, and Yingting Li designed the study, analyzed the data, and performed the molecular and animal experiments. Jing Xu prepared the figures 1-5., Jing Xu, Min Chen, and Luting Zhang wrote the manuscript. All authors reviewed the manuscript.

Acknowledgements

[Please acknowledge anyone who contributed towards the article who does not meet the criteria for authorship including anyone who provided professional writing services or materials. You should obtain permission to acknowledge from all those mentioned in this section.]

Avagliano L, Massa V, George TM, Qureshy S, Bulfamante GP, Finnell RH. Overview on neural tube defects: From development to physical characteristics. Birth Defects Res. 2019;111(19):1455–1467. doi.org. 10.1002/bdr2.1380.
Zhou Y, Crider KS, Yeung LF, et al. Periconceptional folic acid use prevents both rare and common neural tube defects in China. Birth Defects Res. 2022. doi.org. 10.1002/bdr2.1983.
Zaganjor I, Sekkarie A, Tsang BL, et al. Describing the prevalence of neural tube defects worldwide: a systematic literature review. PloS one. 2016;11(4):e0151586. doi.org. 10.1371/journal.pone.0151586.
Finnell RH, Caiaffa CD, Kim S-E, et al. Gene Environment Interactions in the Etiology of Neural Tube Defects. Front Genet. 2021;12:659612. doi.org. 10.3389/fgene.2021.659612.
Gong F, Miller KM. Histone methylation and the DNA damage response. Mutat Res Rev Mutat Res. 2019;780:37–47. doi.org. 10.1016/j.mrrev.2017.09.003.
Wang L, Wang F, Guan J, et al. Relation between hypomethylation of long interspersed nucleotide elements and risk of neural tube defects. The American journal of clinical nutrition. 2010;91(5):1359–1367. doi.org. 10.3945/ajcn.2009.28858.
Huang Y, Lin S, Wang C, et al. Neural Tube Defects and Hypomethylation in Relation to Polycyclic Aromatic Hydrocarbon Exposure. Front Cell Dev Biol. 2020;8:582661. doi.org. 10.3389/fcell.2020.582661.
Price EM, Peñaherrera MS, Portales-Casamar E, et al. Profiling placental and fetal DNA methylation in human neural tube defects. Epigenetics Chromatin. 2016;9:6. doi.org. 10.1186/s13072-016-0054-8.
Lu J, Huang Y, Zhang X, Xu Y, Nie S. Noncoding RNAs involved in DNA methylation and histone methylation, and acetylation in diabetic vascular complications. Pharmacol Res. 2021;170:105520. doi.org. 10.1016/j.phrs.2021.105520.
Zhou J, Yang L, Zhong T, et al. H19 lncRNA alters DNA methylation genome wide by regulating S-adenosylhomocysteine hydrolase. Nature communications. 2015;6(1):10221. doi.org. 10.1038/ncomms10221.
Cheng D, Deng J, Zhang B, et al. LncRNA HOTAIR epigenetically suppresses miR-122 expression in hepatocellular carcinoma via DNA methylation. EBioMedicine. 2018;36:159–170. doi.org. 10.1016/j.ebiom.2018.08.055.
Selamat SA, Chung BS, Girard L, et al. Genome-scale analysis of DNA methylation in lung adenocarcinoma and integration with mRNA expression. Genome research. 2012;22(7):1197–1211. doi.org. 10.1101/gr.132662.111.
Huang W, Gu H, Yuan Z. Identifying biomarkers for prenatal diagnosis of neural tube defects based on "omics". Clin Genet. 2021. doi.org. 10.1111/cge.14087.
Yu J, Mu J, Guo Q, et al. Transcriptomic profile analysis of mouse neural tube development by RNA-Seq. IUBMB Life. 2017;69(9):706–719. doi.org. 10.1002/iub.1653.
Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–360. doi.org. 10.1038/nmeth.3317.
Shannon P, Markiel A, Ozier O, et al. Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Research. 2003;13(11):2498–2504. doi.org. 10.1101/gr.1239303.
Tang Y, Li M, Wang J, Pan Y, Wu FX. CytoNCA: a cytoscape plugin for centrality analysis and evaluation of protein interaction networks. Biosystems. 2015 Jan;127:67–72. doi: 10.1016/j.biosystems.
Yu G, Wang L, Han Y, He Q. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–287. doi.org. 10.1089/omi.2011.0118.
Szklarczyk D, Franceschini A, Wyder S, et al. STRING v10: protein–protein interaction networks, integrated over the tree of life. Nucleic Acids Research. 2014;43(Database issue):D447-452. doi.org. 10.1093/nar/gku1003.
Dweep H, Gretz N. miRWalk2.0: a comprehensive atlas of microRNA-target interactions. Nat Methods. 2015;12(8):697. doi.org. 10.1038/nmeth.3485.
Paraskevopoulou M, Vlachos I, Karagkouni D, et al. DIANA-LncBase v2: indexing microRNA targets on non-coding transcripts. Nucleic Acids Res. 2016;44(D1):D231-238. doi.org. 10.1093/nar/gkv1270.
Tang J, Fang X, Chen J, Zhang H, Tang Z. Long Non-Coding RNA (lncRNA) in Oral Squamous Cell Carcinoma: Biological Function and Clinical Application. Cancers. 2021;13(23). doi.org. 10.3390/cancers13235944.
Zhou X, Lin W, Ren Y, et al. Comparison of DNA Methylation in Schwann Cells before and after Peripheral Nerve Injury in Rats. BioMed research international. 2017;2017:5393268. doi.org. 10.1155/2017/5393268.
Khavari B, Cairns MJ. Epigenomic Dysregulation in Schizophrenia: In Search of Disease Etiology and Biomarkers. Cells. 2020;9(8). doi.org. 10.3390/cells9081837.
Huang Q, Deng G, Wei R, Wang Q, Zou D, Wei J. Comprehensive Identification of Key Genes Involved in Development of Diabetes Mellitus-Related Atherogenesis Using Weighted Gene Correlation Network Analysis. Front Cardiovasc Med. 2020;7:580573. doi.org. 10.3389/fcvm.2020.580573.
Zhou J, Tang Z-Y, Sun X-L. RNF38 inhibits osteosarcoma cell proliferation by binding to CRY1. Biochem Cell Biol. 2021;99(5):629–635. doi.org. 10.1139/bcb-2021-0093.
Long Y, Zhao Q, Huang Y. RNF38 enhances 5-Fluorouracil resistance in colorectal cancer by activating the Wnt pathway. J BUON. 2021;26(4):1246–1251.
Guo W-L, Li N, Ma J-L, Chen X-M, Shi F-Y. Inhibiting microRNA-301b suppresses cell growth and targets RNF38 in cervical carcinoma. Kaohsiung J Med Sci. 2020;36(11):878–884. doi.org. 10.1002/kjm2.12273.
Erdogan B, St Clair RM, Cammarata GM, Zaccaro T, Ballif BA, Lowery LA. Investigating the impact of the phosphorylation status of tyrosine residues within the TACC domain of TACC3 on microtubule behavior during axon growth and guidance. Cytoskeleton (Hoboken). 2020;77(7):277–291. doi.org. 10.1002/cm.21622.
Pasterkamp RJ, Burk K. Axon guidance receptors: Endocytosis, trafficking and downstream signaling from endosomes. Prog Neurobiol. 2021;198:101916. doi.org. 10.1016/j.pneurobio.2020.101916.
Kang D-S, Yang YR, Lee C, et al. Netrin-1/DCC-mediated PLCγ1 activation is required for axon guidance and brain structure development. EMBO Rep. 2018;19(11). doi.org. 10.15252/embr.201846250.
Ramya S, Shyamasundar S, Bay B, Dheen S. Maternal Diabetes Alters Expression of MicroRNAs that Regulate Genes Critical for Neural Tube Development. Frontiers in molecular neuroscience. 2017;10:237. doi.org. 10.3389/fnmol.2017.00237.
Moreno-García L, López-Royo T, Calvo A, et al. Competing Endogenous RNA Networks as Biomarkers in Neurodegenerative Diseases. International journal of molecular sciences. 2020;21(24). doi.org. 10.3390/ijms21249582.
Konopaske GT, Subburaju S, Coyle JT, Benes FM. Altered prefrontal cortical MARCKS and PPP1R9A mRNA expression in schizophrenia and bipolar disorder. Schizophr Res. 2015;164(1–3):100–108. doi.org. 10.1016/j.schres.2015.02.005.
Uddin M, Pellecchia G, Thiruvahindrapuram B, et al. Indexing Effects of Copy Number Variation on Genes Involved in Developmental Delay. Sci Rep. 2016;6:28663. doi.org. 10.1038/srep28663.
Ho WY, Chang J-C, Tyan S-H, et al. FUS-mediated dysregulation of Sema5a, an autism-related gene, in FUS mice with hippocampus-dependent cognitive deficits. Hum Mol Genet. 2019;28(22):3777–3791. doi.org. 10.1093/hmg/ddz217.
Wang Q, Liu Z, Lin Z, et al. Germline Mutations in Associated With Infantile Spasms. Front Genet. 2019;10:605. doi.org. 10.3389/fgene.2019.00605.

No competing interests reported.

TablesS1.docx

Integration of epigenomic and transcriptome analyses of neural tube defects reveals methylation driver lncRNAs and mRNAs

Status:

Version 1

Abstract

Figures

Background

Materials And Methods

Results

Discussion

Conclusion

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1