Fetal DNA hypermethylation in tight junction pathway is associated with neural tube defects: A genome-wide DNA methylation analysis

ABSTRACT

Neural tube defects (NTDs) are a spectrum of severe congenital malformations of fusion failure of the neural tube during early embryogenesis. Evidence on aberrant DNA methylation in NTD development remains scarce, especially when exposure to environmental pollutant is taken into consideration. DNA methylation profiling was quantified using the Infinium HumanMethylation450 array in neural tissues from 10 NTD cases and 8 non-malformed controls (stage 1). Subsequent validation was performed using a Sequenom MassARRAY system in neural tissues from 20 NTD cases and 20 non-malformed controls (stage 2). Correlation analysis of differentially methylated CpG sites in fetal neural tissues and polycyclic aromatic hydrocarbons concentrations in fetal neural tissues and maternal serum was conducted. Differentially methylated CpG sites of neural tissues were further validated in fetal mice with NTDs induced by benzo(a)pyrene given to pregnant mice. Differentially hypermethylated CpG sites in neural tissues from 17 genes and 6 pathways were identified in stage 1. Subsequently, differentially hypermethylated CpG sites in neural tissues from 6 genes (BDKRB2, CTNNA1, CYFIP2, MMP7, MYH2, and TIAM2) were confirmed in stage 2. Correlation analysis showed that methylated CpG sites in CTNNA1 and MYH2 from NTD cases were positively correlated to polycyclic aromatic hydrocarbon level in fetal neural tissues and maternal serum. The correlation was confirmed in NTD-affected fetal mice that were exposed to benzo(a)pyrene in utero. In conclusion, hypermethylation of the CTNNA1 and MYH2 genes in tight junction pathway is associated with the risk for NTDs, and the DNA methylation aberration may be caused by exposure to benzo(a)pyrene.

KEYWORDS:

• Benzo(a)pyrene
• DNA methylation
• neural tube defects
• tight junction

Introduction

Neural tube defects (NTDs) are a spectrum of severe congenital malformations of fusion failure of the neural tube during early embryogenesis.^1,2 As a multifactorial disorder, NTDs arise from a complex interaction of genetic and environmental factors. The genetic component is estimated to be higher than 50% and over 200 genes associated with NTDs have been identified in mouse models. However, few causative genes have been confirmed in human to date. Environmental risk factors include deficient folate intake, exposure to polycyclic aromatic hydrocarbons (PAHs), anticonvulsant therapy, and diabetes mellitus, among others.^1,3 It is well established that periconceptional folic acid supplementation decreases the risk of NTDs by 50–75%.⁴ However, the prevalence of NTDs remains relatively stable after a dramatic decline caused by folic acid fortification,⁵ and the proportion of NTD cases that can be attributed to known risk factors is lower than one-third.⁶ The etiologies and mechanisms underlying NTD development remain to be elucidated.

Epigenetic alterations are thought to be involved in the development of complex diseases, including NTDs. Aberration in DNA methylation has been proposed as the mechanism underlying folic acid prevention of NTDs. Folate is a substrate for the activation of methyl groups in one-carbon metabolism and plays a critical role in DNA synthesis and the generation of S-adenosylmethionine, the major methyl donor for many important methylation reactions, including DNA methylation. Studies have shown a decreased methylation capacity in NTDs.^{4, 7} The only well-established genetic risk factor for human NTDs is the 677C>T variant in the 5,10-methylene tetrahydrofolate reductase gene (MTHFR). The MTHFR 677C>T variant leads to global DNA hypomethylation, which is more pronounced under low folic acid conditions.⁸ However, the exact mechanism by which alterations in DNA methylation would cause NTDs is still unknown.

Besides vague effect of global DNA methylation, disruption of methylation at specific genes has been linked to NTDs, such as imprinted genes,^9,10 transposon genes,⁷ DNA repair genes,¹¹ planar-cell polarity genes,¹² folate receptor genes,¹³ and HOX genes.¹⁴ Hypermethylation or hypomethylation of these specific genes was identified in NTD cases by targeted gene methylation profiling. However, apart from these findings, DNA methylation landscapes of NTDs have remained a mystery until now. A genome-wide DNA methylation-profiling method is required to allow a discovery-based, less biased, and whole-genome study.

We hypothesized that aberrant DNA methylation of key pathway genes might result in increased risk for abnormal embryonic development and NTDs. We tested this hypothesis by profiling genome-wide DNA methylation in fetal neural tissues of NTD cases and non-malformed controls in two stages and evaluating the association between methylation changes and NTD risk. We further examined the correlation of CpG site methylation in fetal neural tissues and PAH levels in NTD fetus and mother. Finally, we assessed the association by examining the methylation profile in fetal mice with NTDs induced by benzo(a)pyrene (BaP) that was given to the pregnant dam.

Materials and methods

Study population

The study subjects were recruited from five rural counties of Shanxi Province (Taigu, Pingding, Xiyang, Shouyang, and Zezhou) in northern China from 2011 to 2014. Cases included terminated fetus following prenatal diagnosis of an NTD, which were ascertained through a population-based birth defect surveillance program.¹⁵ Controls were induced fetus with no congenital malformations and were matched to cases by sex, mother's county of residence, and date of mother's last menstrual period. Information on obstetric characteristics was collected through face-to-face interviews and by viewing medical records. Fasting blood samples from NTD case women and control women were collected at delivery or termination of NTD-affected pregnancies. Spinal cord, brain, and liver tissues of fetuses were collected at pregnancy termination by experienced pathologists. Blood and tissue samples were stored at −80°C until analysis. The study protocol was approved by the institutional review board of Peking University, and written informed consent was obtained from the mothers before investigation.

Animal experiment

Virgin female and male ICR mice of 8–9 weeks old and weighing 25–28 g were used for this experiment. These mice were maintained on a 12-h light/dark cycle and were allowed free access to food and water. Presence of the female vaginal plug in the next morning after pairing was designated as embryonic day (ED) 0.5. Pregnant mice were randomly divided into four groups and were treated with BaP dissolved in corn oil intraperitoneallly (i.p.) at doses ranging from 250 to 350 mg kg⁻¹ daily beginning on ED 7.0 for four consecutive days. The control group was treated with 15 ml kg⁻¹ of vehicle on the same days. Embryos were harvested on ED 10.5. Alive fetuses were carefully inspected for visible external malformations under a dissecting microscope. NTD-affected embryos were classified as showing distinct evidence of failed closure of the neural tube. From these embryos, the anterior part of the neural tube was dissected and stored at −80°C until methylation analysis. All animal procedures were approved by the Institutional Animal Care and Use Committee of Peking University.

Methylation assay

Human methylation study was performed in two stages. In stage 1, tissue samples from 10 NTD cases (5 spina bifida and 5 anencephaly) and 8 non-malformed control fetuses recruited between 2011 and 2012, were used to perform methylation analysis with Infinium HumanMethylation450 BeadChip assay (450K; Illumina, San Diego, CA, USA). Stage 2 was conducted to confirm the results of stage 1. In this stage, tissue samples from 20 NTD cases (10 spina bifida and 10 anencephaly) and 20 non-malformed control fetuses, recruited from 2013 to 2014, were analyzed using a Sequenom MassARRAY system.

DNA was extracted from frozen neural tissues using the QIAamp DNeasy Blood and Tissue Kit (QIAGEN, Hilden, Germany). Bisulfite modification of 1μg DNA was conducted using an EZ DNA Methylation Kit (Zymo Research, CA, USA) according to the manufacturer's instruction. The bisulfite conversion reaction was performed in duplicate for each sample to minimize potential bias caused by variable conversion efficiency, and pooled bisulfite-treated DNA was used for subsequent array analysis. Genome-wide DNA methylation was performed using the Illumina 450K array to interrogate approximately 480,000 CpG sites across the genome, according to Illumina's standard protocol.

Illumina GenomeStudio software (Illumina) was used to extract signal intensities for each probe, perform initial control quality checks, and estimate β scores and detection P values. β score was defined as the proportion of total signal from the methylation-specific probe or color channel [β = intensity of the methylated allele /(intensity of the unmethylated allele + intensity of the methylated allele + 100)]. Detection P value was defined as the 1-P value computed from the background model characterizing the chance that the target sequence signal was distinguishable from the negative controls. Pooled t-tests with multiple comparison tests were used to identify differentially methylated CpG sites. Differentially methylated CpGs were defined if they had a false discovery rate <0.05 and a β difference in methylation >0.2 (20% methylation) between cases and controls. Hierarchical clustering of the methylation data was performed with the top 500 significant sites to cluster the phenotypes (case vs. control). To determine the overall similarity between samples, principal component analysis that project β values into an explanatory principal components was performed using R software (http://www.r-project.org/). Clustering of DNA methylation data based on phenotype-associated probes was identified by ANOVA analysis.

Pathway analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID) (http://david.abcc.ncifcrf.gov) was used to analyze biological features associated with genes hypo- or hyper-methylation in neural tissues from NTD cases compared with controls. The gene lists were obtained with the criteria of a P value < 0.05 and Δβ-value > |0.2| (between cases and controls) due to the input limit of genes to the DAVID tool. The Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used to identify functional gene groups and ontology terms that are significantly overrepresented.

Validation of array targets

Target validation of differentially methylated CpGs in human was performed using the Sequenom EpiTYPER (Sequenom, San Diego, USA) in stage 1 and stage 2. PCR primers were designed using the Sequenom EpiDesigner software (http://www.epidesigner.com/) and the sequences were provided in Table 1. About 20 ng of bisulphite DNA was amplified by PCR, followed by reverse transcription, fragmentation, and analysis on a mass spectrometer (Sequenom, Inc., San Diego, USA). EpiTYPER Analyzer software version 1.0 was used to translate mass signal patterns into quantitative DNA methylation levels of different CpG sites. Samples of cases and controls were determined simultaneously under uniform conditions to avoid bias during bisulfite treatment, PCR amplification, and microarray detection.

Table 1. Characteristics of neural tube defect cases and controls in stage 1 (n = 18)^a.

Characteristic	Cases (n = 10)	Controls (n = 8)	P^b
Gestational week, mean (SD)	25.1 (7.6)	25.6 (8.9)	0.907
Gender, n (%)^c
Male	4 (44.4)	3 (37.5)	0.772
Female	5 (55.6)	5 (62.5)

^a The total number may not be equal to the number of cases or controls due to missing or unknown data.

^b Independent samples test for gestational week; Pearson's χ² test for gender.

^c One case with missing information on gender.

Methylation analysis of mouse tissues was also performed using the Sequenom EpiTYPER (Sequenom, San Diego, USA). The three groups with BaP exposure of three doses were combined and then divided into two groups according to the phenotype of neural tube. Thus, the neural samples for methylation analysis were divided into three groups: controls, normal-neural-tube fetus with BaP exposure, and NTD fetus with BaP exposure. DNA was extracted from frozen neural tube using TIANamp Genomic DNA Kit (TIANGEN Biotech, Beijing, China). A total of 500 ng of genomic DNA from each sample was bisulfite-treated with the EZ DNA Methylation Kit (Zymo Research, CA, USA). PCR primers were designed using the Sequenom EpiDesigner software (http://www.epidesigner.com/) and the sequences were provided in Table S2.

Table 2. Distribution of genomic regions for significant CpG sites passing the stringent filtering criteria^a.

Genomic region of CpG sites	Total CpG sites (%)	Number of hypermethylated CpG sites (%)	Number of hypomethylated CpG sites (%)
TSS1500	68917 (14.2)	915 (7.4)	970 (8.9)
TSS200	52238 (10.8)	410 (3.3)	615 (5.6)
5′UTR	42669 (8.8)	1016 (8.2)	914 (8.4)
1stExon	22727 (4.7)	187 (1.5)	262 (2.4)
Gene body	161582 (33.3)	5142 (41.5)	3798 (34.8)
3′UTR	17487 (3.6)	413 (3.3)	355 (3.3)
Others/open sea	119579 (24.6)	4300 (34.7)	3997 (36.6)
CpG Island (CGI)^b	150180 (31.0)	1433 (11.6)	717 (6.6)
N Shelf (0–2 kb from CGI)	24821 (5.1)	596 (4.8)	538 (4.9)
N Shore (0–2 kb from CGI)	62807 (12.9)	1599 (12.9)	1258 (11.5)
S Shelf (0–2 kb from CGI)	22275 (4.6)	489 (3.9)	415 (3.8)
S Shore (0–2 kb from CGI)	49171 (10.1)	1156 (9.3)	965 (8.8)
Intergenic	175945 (36.3)	7110 (57.4)	7018 (64.3)
Total	485199	12383	10911

^a CpG loci were categorized in the context of: (1) CpG islands, CpG shores, CpG shelves or others; and (2) genomic location of TSS1500 (1.5 kb within transcriptional start site), TSS200 (200 bp within transcriptional start site), 5′UTR (5′untranslated region) and 1st exon, gene body and 3′UTR (3′untranslated region) or others/open sea. For significant hypermethylation/hypomethylated sites, we selected CpG sites with the following criteria: (1) the mean difference in methylation level between case and control tissues is >20%; (2) pooled t-test P-value < 0.05.

^b CGI defined as DNA regions with at least 500 bp, a GC percentage ≥55% and with a CpG ratio (observed CpG/expected CpG) >0.65.

PAH analysis

The methods used for the analysis of thirteen parent PAHs have been described in detail elsewhere.^{16, 17} After concentration, each sample was spiked with internal standards and analyzed with an Agilent 7890A-5975C gas chromatograph and mass spectrometer equipped with a HP-5MS capillary column (30 m × 0.25 mm × 0.25 μm). Triglycerides and cholesterol contents were measured using the oxidative method (INTEC Technology Co., China) with an automatic biochemical analyzer (Olympus AU400, Japan). During chemical analysis, the case-control status of the samples was masked from the analysts. PAH concentrations were expressed on a lipid weight basis and reported as ng/g lipid.

Statistical analysis

Differences in proportions of population characteristics between groups were examined with Pearson's χ² test. Independent samples t-test and paired t-test were used to identify CpG sites that were differentially methylated between cases and controls with Sequenom MassARRAY in sample sets for validation in stage 1 and 2. The comparison of PAH concentrations in fetal liver between cases and controls was performed by independent samples t-test. In the mouse study, because the sample size of the NTD subgroup of different BaP treatments was small, we merged data of different concentrations of BaP treatments to increase statistical power. Differentially methylated CpG sites in neural tissues among fetal mice with and without BaP exposure during ED 7.0–10.5 were detected by one-way ANOVA. Correlation analysis of differentially methylated CpG sites and PAH levels in fetal neural tissues and maternal blood was performed by Spearman's correlation coefficient. Statistical analyses were conducted using SPSS 18.0. A two-tailed P value < 0.05 was considered statistically significant.

Results

DNA methylation pattern in NTD cases compared with controls

The characteristics of the NTD cases and controls in stage 1 are presented in Table 1. No significant differences were observed between case and control groups regarding gestational week and gender.

Using stringent criteria of CpG sites with a false discovery rate <0.05 and a β difference in methylation level >20%, a total of 23,294 (4.8%) differentially methylated CpG sites, out of 485,199 CpG sites across the entire genome, were identified between cases and controls. Among all significant CpG sites, 12,383 (53%) were significantly hypermethylated, and 10,911 (47%) were significantly hypomethylated in neural tissues of NTD cases. Analysis of genomic location of CpG sites showed that differentially methylated probes decreased in transcription start site (TSS) 1500 (7.4% for hypermethylated CpG sites and 8.9% for hypomethylated CpG sites versus 14.2% of total CpG Sites) and TSS 200 (3.3% and 5.6% vs. 10.8%), and increased in gene body (41.5% and 34.8%, respectively, vs. 33.3%) and open sea (34.7% and 36.6%, respectively, vs. 24.6%), compared with their representation on the microarray (Table 2). In addition, analysis of neighborhood location of CpG sites revealed a significant decrease in CpG island (11.6% for hypermethylated and 6.6% for hypomethylated CpG sites vs. 31.0%), and a significant enrichment of differentially methylated probes in intergenic regions (57.4% and 64.3%, respectively, vs. 36.3%), compared with their representation on the microarray (Table 2).

Fig. 1 shows the results of hierarchical cluster analysis of the top 500 significant CpG sites (based on |β difference|, cut off value about 0.45) that distinguished lesion neural tissues in NTD cases from normal neural tissues in controls. Excellent separation of cases and controls was observed. To gain more statistical insight into the data, a principal component analysis was performed to visualize the difference in DNA methylation between case and control samples. As expected, principal component analysis almost entirely separated NTD cases from controls, while the two NTD subgroups, namely spina bifida and anencephaly clusters, did not separate (Fig. S1).

Figure 1. Hierarchical cluster analysis of the top 500 significantly differentially methylated CpG sites in neural tissues of NTD cases and controls of Population 1. Rows represent probes and columns represent samples. Cells are colored according to level of methylation (blue = low methylation, red = high methylation). The top 500 significant CpG sites can clearly distinguish lesion tissues from normal neural tissues.

Identification of functionally related gene groups

A total number of 97 gene ontology term clusters were generated when analyzing the hypermethylated CpG sites in NTD cases compared with controls by KEGG pathway analysis (Table 3). These gene ontology terms include several KEGG pathways involving Wnt signaling pathway, tight junction, regulation of actin cytoskeleton, axon guidance, aldosterone-regulated sodium reabsorption, metabolism of xenobiotics by cytochrome P450, and heparan sulfate biosynthesis with the enrichment score of 1.80–3.64. Wnt signaling pathway and tight junction were the only two pathways identified of significantly differential methylation between NTD cases and controls after Bonferroni and Benjamini correction for multiple comparisons.

Table 3. KEGG pathway analysis of the differentially methylated genes in neural tissues from NTD cases as compared to controls^a.

Pathway_ID	Term	P-value	Gene symbol	Fold enrichment	Bonferroni	Benjamini
hsa04810	Regulation of actin cytoskeleton	4.29E-03	FGF19, FGFR2, FGFR4, APC2, INS-IGF2, BCAR1, RDX, BDKRB2, PFN3, KRAS, PAK2, TIAM2, TIAM1, PAK1, FGF1, PIK3R1, ARHGEF7, CHRM3, CHRM2, RRAS2, CHRM1, ITGA7, CYFIP2, MYLK, MYH10	1.80	0.53	0.12
hsa04310	Wnt signaling pathway	1.43E-04	PRKCA, TBL1XR1, CTBP1, CTBP2, APC2, PPP2R5C, VANGL2, MMP7, TP53, SMAD3, FZD3, MAPK10, DAAM2, CCND1, SFRP2, PPP2CA, FRAT1, SIAH1, CAMK2A, NFATC1, WNT8B	2.37	0.02	0.02
hsa04530	Tight junction	1.86E-04	PRKCA, PRKCZ, CLDN19, VAPA, CLDN6, MYH2, MPP5, MYH7, CLDN20, CDK4, CTNNA1, CLDN14, CTNNA2, EPB41L2, KRAS, PPP2CA, RRAS2, PPP2R2B, PPP2R2C, TJP2, PTENP1, MYH10	2.45	0.03	0.02
hsa04360	Axon guidance	2.43E-02	PLXNB1, PLXNA2, PLXNB2, DPYSL5, KRAS, PAK2, FYN, RGS3, SEMA3B, PAK1, SEMA4D, NFATC1	1.85	0.99	0.38
hsa04960	Aldosterone-regulated sodium reabsorption	1.24E-03	PRKCA, SGK1, KRAS, INS-IGF2, HSD11B1, NR3C2, NEDD4L, SCNN1A, PIK3R1, KCNJ1	3.64	0.19	0.07
hsa00980	Metabolism of xenobiotics by cytochrome P450	4.47E-02	CYP1B1, CYP2C19, ADH4, ALDH1A3, EPHX1, MGST1	2.24	1.00	0.51
hsa00534	Heparan sulfate biosynthesis	2.66E-02	EXTL3, EXT1, HS3ST3B1, GLCE	3.44	0.99	0.37

^a KEGG, Kyoto Encyclopedia of Genes and Genomes.

P-value was obtained with modified Fishers exact statistic and adjusted for multiple testing with Bonferroni and Benjamini methods

Candidate gene selection and validation of findings in stage 1

The selection of candidate genes was based on both differentially methylated CpG sites generated with microarray analysis and differentially methylated pathway analysis. From the list of genes generated with KEGG pathway analysis (Table 3) and the list of 12,383 differentially hypermethylated CpG sites (Table 2), the following criteria for candidate gene selection were applied: (1) without a single significantly hypomethylated CpG site in whole gene; (2) with ≥2 hypermethylated CpG sites in one gene; and (3) within the top 20 rates of hypermethylated CpG sites in total CpG sites determined in the microarray analysis. The 20 candidate marker genes were FYN, ALDH1A3, CYP1B1, MGST1, EXTL3, BDKRB2, CHRM1, CYFIP2, FGF19, FGFR2, PIK3R1, TIAM2, CTNNA1, MYH2, PPP2R2B, CAMK2A, MMP7, TBL1XR1, VANGL2, and WNT8B (Table 4).

Table 4. Selection of differentially hypermethylated genes in neural tissues from NTD cases as compared to controls for validation using Sequenom MassARRAY in stage 1.

Pathway	Gene symbol	Gene full name	Location
Axon guidance	FYN	FYN proto-oncogene, Src family tyrosine kinase	6q21
Metabolism of xenobiotics by cytochrome P450	ALDH1A3	aldehyde dehydrogenase 1 family member A3	15q26.3
	CYP1B1	cytochrome P450 family 1 subfamily B member 1	2p22.2
	MGST1	microsomal glutathione S-transferase 1	12p12.3-p12.1
Heparan sulfate biosynthesis	EXTL3	exostosin like glycosyltransferase 3	8p21
Regulation of actin cytoskeleton	BDKRB2	bradykinin receptor B2	14q32.2
	CHRM1	cholinergic receptor muscarinic 1	11q13
	CYFIP2	cytoplasmic FMR1 interacting protein 2	5q33.3
	FGF19	fibroblast growth factor 19	11q13.1
	FGFR2	fibroblast growth factor receptor 2	10q26
	PIK3R1	phosphoinositide-3-kinase regulatory subunit 1	5q13.1
	TIAM2	T-cell lymphoma invasion and metastasis 2	6q25.2
Tight junction	CTNNA1	catenin alpha 1	5q31.2
	MYH2	myosin, heavy chain 2, skeletal muscle, adult	17p13.1
	PPP2R2B	protein phosphatase 2 regulatory subunit B, beta	5q32
Wnt signaling pathway	CAMK2A	calcium/calmodulin dependent protein kinase II alpha	5q32
	MMP7	matrix metallopeptidase 7	11q22.2
	TBL1XR1	transducin (beta)-like 1 X-linked receptor 1	3q26.32
	VANGL2	VANGL planar cell polarity protein 2	1q23.2
	WNT8B	Wnt family member 8B	10q24

Sequenom MassARRAY platform was used for validation of the differential methylation of 20 selected genes found in stage 1. Table S3 shows the percentage methylation of the analyzed CpG sites in neural tissues from NTD cases and controls. The results are in very good agreement with the array data, except for 3 genes that were not examined. The mean methylation levels in most CpG sites in the 17 genes from NTD tissues were higher than those from control tissues, although a few CpG sites only showed hypermethylated trend in NTD tissues.

Validation in an independent sample set in stage 2

Further, the differential methylation of 17 genes was independently validated in a separate sample set in stage 2 using Sequenom MassARRAY technology. The characteristics of the NTD cases and controls in stage 2 are presented in Table S4. Gestational week and gender had no significant differences between cases and controls. Table S5 shows the percentage methylation of the analyzed CpG sites from Table S3 in neural tissues from NTD cases and controls. The methylation levels in CpG sites in the 12 genes were higher in NTD tissues than in control tissues (Table S5). Of the 12 genes, 6 genes (BDKRB2, CTNNA1, CYFIP2, MMP7, MYH2, and TIAM2) showed multiple (3–5) CpG sites with hypermethylation, and the other genes only had 1 or 2 hypermethylated CpG sites.

Table 5. Differentially methylated CpG sites in neural tissues of fetal mice with and without BaP exposure during ED 7.0–10.5^a.

	Controls (1)		Normal neural tube fetus with BaP exposure (2)		NTD fetus with BaP exposure (3)		Difference between (3) and (1)	Difference between (3) and (2)	Difference between (1) and (2)
CpG sites	Mean	SD	Mean	SD	Mean	SD	P-value	P-value	P-value	P for trend
Ctnna1_CpG_1	0.83	0.09	0.89	0.10	0.92	0.08	0.034	0.460	0.172	0.033
Ctnna1_CpG_2	0.40	0.10	0.45	0.06	0.48	0.06	0.027	0.533	0.132	0.026
Ctnna1_CpG_3	0.80	0.11	0.87	0.05	0.87	0.05	0.061	0.980	0.065	0.057
Ctnna1_CpG_4	0.40	0.10	0.45	0.06	0.48	0.06	0.027	0.533	0.132	0.026
Ctnna1_CpG_5	0.63	0.09	0.69	0.11	0.75	0.08	0.008	0.212	0.152	0.008
Ctnna1_CpG_6.7	0.57	0.10	0.62	0.10	0.69	0.07	0.013	0.137	0.313	0.013
Ctnna1_CpG_8	0.40	0.11	0.44	0.08	0.52	0.06	0.007	0.076	0.363	0.008
Ctnna1_CpG_9	0.51	0.10	0.58	0.11	0.68	0.11	0.002	0.081	0.171	0.003
Ctnna1_CpG_10	0.46	0.16	0.57	0.12	0.65	0.09	0.003	0.186	0.079	0.003
Ctnna1_CpG_11.12	0.60	0.09	0.62	0.14	0.67	0.09	0.151	0.353	0.627	0.151
Ctnna1_CpG_13	0.63	0.13	0.76	0.14	0.80	0.11	0.007	0.555	0.032	0.006
Myh2_3_CpG_1	0.77	0.10	0.71	0.17	0.81	0.08	0.429	0.094	0.346	0.429
Myh2__3_CpG_2	0.50	0.03	0.53	0.09	0.57	0.09	0.072	0.335	0.406	0.072
Myh2__3_CpG_4	0.34	0.03	0.38	0.02	0.42	0.07	0.001	0.072	0.061	0.001
Myh2__3_CpG_5	0.54	0.05	0.61	0.04	0.64	0.08	0.001	0.344	0.013	0.001
Myh2__3_CpG_6	0.33	0.04	0.36	0.07	0.38	0.06	0.076	0.521	0.261	0.076

^a Independent Samples Test.

Correlation of CpG site methylation in fetal neural tissues with PAH concentration in fetal liver tissues and in maternal blood

Our previous study showed that the risk of NTDs increased with PAH concentration in maternal serum and placenta.^16,18 To gain more knowledge about the relationships among PAH exposure, gene methylation and NTD risk, correlation analysis of differentially methylated CpG sites in fetal neural tissues and PAH concentrations in fetal neural tissues and maternal serum was performed in the combined cohort of cases in stage 1 and 2 (Table S6). The results showed that differentially methylated CpG sites in CTNNA1 and MYH2 were positively correlated with PAH concentrations in NTD fetal neural tissues, and also in NTD maternal serum. Levels of PAHs in fetal liver from NTD cases as compared to controls were presented in Table S7. Although the differences between cases and controls did not reach statistical significance, it is apparent that the case group had higher levels of low and high molecular weight PAHs and of total PAHs. In addition, our previous study showed higher levels of PAHs in serum of mothers of NTD cases than in serum of control mothers.¹⁶

Differentially methylated CpG sites of neural tissues in NTD-fetal mice

To validate the findings of the correlation analysis in the human study, the methylation analysis was conducted in neural tissues of fetal mice with and without BaP exposure from ED 7.0 through ED 10.5. NTDs in ED 10.5 mice embryos induced by BaP are shown in Fig. S2. As shown in Table 5, differentially methylated CpG sites in Ctnna1 and Myh2 were significantly higher in NTD mice with BaP exposure than those in non-exposed controls. The methylated levels of CpG sites in Ctnna1 and Myh2 in normal-neural-tube fetal mice with BaP exposure tend to be higher than those in controls and be lower than those in NTD mice, although statistical significance was not reached. Methylation levels of most CpG sites in Ctnna1 and Myh2 showed increasing trends for the three groups, except for the CpG site of Ctnna1 _CpG_3, Ctnna1 _CpG_11.12, Myh2_3_CpG_1, Myh2_3_CpG_2, and Myh2_3_CpG_6.

Discussion

Abnormal DNA methylation has been implied in the etiology of NTDs and proposed as the mechanism underlying folic acid prevention of NTDs. In this study, to obtain a full picture of the methylation pattern, genome-wide DNA methylation was profiled in fetal neural tissues. Differentially hypermethylated CpG sites in neural tissues from 17 genes and 6 pathways were identified initially, of which the CpG sites from 6 genes (BDKRB2, CTNNA1, CYFIP2, MMP7, MYH2, and TIAM2) were validated in an independent cohort of subjects. Correlation analysis showed that hypermethylation of the CTNNA1 and MYH2 genes from the tight junction pathway is associated with the risk for NTDs and with levels of PAHs in fetal neural tissues of NTD cases. The hypermethylated CpG sites in Ctnna1 and Myh2 were further found in neural tissue of mice with NTDs induced by BaP.

Although epigenetic mechanisms are thought to be involved in the complex etiology of NTDs,⁴ little progress has been made in this area. Rochtus et al. reported methylation data for the HOXB7 gene from a genome-wide DNA methylation analysis using the Illumina 450K array using blood samples as surrogate for brain tissue of NTD cases. Different methylation patterns between blood and brain tissue hampered the interpretation of their results. Global DNA hypomethylation in fetal brain tissue was found to be associated with NTD-affected pregnancy.¹⁹ Global DNA methylation status was not a specific biomarker and could not indicate the site-specific alteration of methylation in genes that are critical to neural tube development. Our study for the first time profiled genome-wide DNA methylation of NTDs in lesion neural tissues of fetal NTD cases. Out of the differentially hypermethylated CpG sites identified in 17 genes and 6 pathways using the microarray, the hypermethylated CpG sites of genes regulating actin cytoskeleton (BDKRB2, CYFIP2, and TIAM2), tight junction (CTNNA1, MYH2), and Wnt signaling pathway (MMP7) were validated in an independent case-control sample matched by gestational age and fetal gender. The findings of our study were inconsistent with a previous study,²⁰ which is the only one that reported genome-wide DNA methylation in human NTDs. The authors showed that neither aberrant CpG site in specific genes nor genome-wide methylation were characteristic of placental chorionic villi, spinal cord, brain, and muscle tissue in NTD cases, except for the kidneys of spina bifida cases. Given the great effect of gestational age and fetal gender on DNA methylation, failure to match the two variables between cases and controls may confound Price et al.'s results, which may partly contribute to the difference between their results and ours.

To date, no hypermethylated CpG sites in genes that regulate actin cytoskeleton and tight junction pathways have been reported that could be related to the risk of human NTDs. The dynamic nature of the actin cytoskeleton is critical for regulating cellular processes, such as division, polarity, adhesion, migration, and morphology.^21,22 The neural fold seems highly dependent on the actin cytoskeleton and sensitive to interference in cytoskeletal dynamics.²³ Animal studies have shown that the architecture and dynamics of the actin cytoskeleton must be precisely regulated during neural tube closure. Several lines of knockout mice lacking cytoskeleton-associated proteins, such as paladin,²⁴ Nap1,²⁵ p190RhoGap,²⁶ Mena/profilin,²⁷ and Vinculin²⁸ have shown NTDs. Tight junction is defined as a circumferential ring at the apex of epithelial cells that seals adjacent cells to one another and plays a crucial role in the formation and physiological regulation of epithelial cell sheets.²⁹ Downregulation of tight junction components occurs in the neuroepithelium during neural tube closure and in pre-migratory neural crest cells during the initiation of epithelial-to-mesenchymal transition.^30,31 Whole-genome choroidal neovascularization analysis revealed that tight junction is involved with the pathogenesis of human NTDs.³² It is well known that the Wnt signaling pathway, especially the planar cell polarity (PCP) signaling pathway, plays an essential role in neural tube development.^33,34 Genetic defects in the PCP pathway have been strongly associated with NTDs in animal models and human studies.³⁵ Aberrant methylation of the FZD3 gene, a PCP component, was observed in brain tissue in human spina bifida and showed an association with human NTD development.¹² Our findings expand the knowledge of the relation between these pathways and NTD development.

Prenatal exposure to PAHs has been associated with some adverse birth outcomes, such as low birth weight and preterm birth.^36,37 Our group previously reported that an elevated risk of NTDs was associated with maternal exposure to PAHs in the Shanxi province of northern China,^16,18 from which the subjects of this study were recruited. Aberrant DNA methylation occurring in a population with a high-level of PAH exposure may imply a role of DNA methylation changes, caused by exposure to PAH, in the development of NTDs. However, no human study examined the effect of prenatal PAH exposure on DNA methylation in relation to NTDs. A cohort study of the Columbia Center for Children's Environmental Health reported that prenatal PAH exposure in African-American and Dominican women was associated with lower global methylation and hypermethylation of interferon γ in umbilical cord white blood cells.^38,39 In the present study, differentially methylated CpG sites in CTNNA1 and MYH2 were positively correlated with PAH concentrations in NTD fetal neural tissues, and also in maternal serum of NTD cases. Further, the methylation analysis was performed in neural tissues of fetal mice with and without BaP exposure during ED 7.0–10.5. Consistent with the findings of the human study, the methylated CpG sites in CTNNA1 and MYH2 were significant higher in NTD mice with BaP exposure than those in controls without BaP exposure.

α-, β-, and γ-catenins (CTNNA1 encodes α-catenin1) are proteins associated with the cytoplasmic domain of cadherins, a family of transmembrane cell adhesion molecules.⁴⁰ α-catenin is highly expressed in neural tissues⁴¹ and links cadherin complexes to the actin cytoskeleton, and is necessary for functional cell-cell adhesion by these complexes, which is required for normal embryonic development.⁴⁰ Besides binding to the cadherin /catenin complex, β-catenin is involved in the Wnt signaling pathway and plays a key role in the development of the neural tube.⁴² In the absence of PCP function, a non-canonical Wnt-frizzled-disheveled signaling cascade, cells fail to interdigitate mediolaterally resulting in an unfused neural tube.¹ Fluctuations in α-catenin can impact upon both the nuclear accumulation of β-catenin and Wnt signaling and result in dysregulation of Wnt/β-catenin signaling.^43,44 Therefore, we speculate that downregulation of CTNNA1 expression might affect the development of the neural tube through regulation of the Wnt/β-catenin signaling, although cardiac-specific α-E-catenin conditional knockout mice exhibit abnormalities in cardiac muscle,⁴⁵ but not NTDs.

Myosin II is a multifunctional motor protein that controls cell behavior, such as cell adhesion, cell migration, and cell division.⁴⁶ Vertebrate neural tube formation involves two dramatic morphogenetic events, convergent extension, and bending and rolling of the neural plate. Myosin II is essential for both of these processes. It has been demonstrated that Myosin II and its upstream regulator ROCK are critical for both normal neural cell cortical tension and shape and for autonomous convergent extension of the neural tissue.^47-49 Moreover, depletion of Myosin II results in disruption of the polarized cytoskeleton, loss of the polarized protrusive activity characteristic of intercalating cells, eventual loss of cell-cell and cell-matrix adhesion, and failure of blastopore closure.⁵⁰ No phenotype information for Myh2 knockout mice is available in the literature. Our findings suggest that hypermethylation of MYH2 might be related to NTDs by alteration of cell polarity and cell behavior during neurulation.

The main limitation of this study is the lack of information regarding the expression levels (mRNA) of the identified hypermethylated CpG site-genes. Fresh neural tissue samples of NTD cases are extremely difficult to collect, as most fetuses with an NTD have been deceased for several hours before leaving the maternal body. We have not yet been able to secure human neural tissues fresh enough to obtain mRNA.

In conclusion, through genome-wide DNA methylation profiling of fetal neural tissues from NTD cases and unaffected controls, we demonstrated that hypermethylation of the CTNNA1 and MYH2 genes in the tight junction pathway is associated with the risk for NTDs; the changes in DNA methylation may be related to exposure to environmental pollutants, such as BaP in the present study. These findings provide novel mechanistic insights into the etiology of NTDs involving environmental pollutants and epigenetics. Future studies are needed to determine whether maternal exposure to PAHs affect the expression of genes in the tight junction pathway in fetal neural tissue, and elucidate the impact of abnormal expression of CTNNA1 and MYH2 genes on the development of NTDs.

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Funding

This work was supported by grants from the National Natural Science Foundation of China (Grant No. 31371523 and 81472987); the Ministry of Education of China (Grant No. 20130001110064); Beijing Natural Science Foundation (Grant No. 7162094); and the National Key Research and Development Program, Ministry of Science and Technology, P.R. China (Grant No. 2016YFC1000501).