Skip to main content

Association between rare variants in specific functional pathways and human neural tube defects multiple subphenotypes



Neural tube defects (NTDs) are failure of neural tube closure, which includes multiple central nervous system phenotypes. More than 300 mouse mutant strains exhibits NTDs phenotypes and give us some clues to establish association between biological functions and subphenotypes. However, the knowledge about association in human remains still very poor.


High throughput targeted genome DNA sequencing were performed on 280 neural tube closure-related genes in 355 NTDs cases and 225 ethnicity matched controls,


We explored that potential damaging rare variants in genes functioning in chromatin modification, apoptosis, retinoid metabolism and lipid metabolism are associated with human NTDs. Importantly, our data indicate that except for planar cell polarity pathway, craniorachischisis is also genetically related with chromatin modification and retinoid metabolism. Furthermore, single phenotype in cranial or spinal regions displays significant association with specific biological function, such as anencephaly is associated with potentially damaging rare variants in genes functioning in chromatin modification, encephalocele is associated with apoptosis, retinoid metabolism and one carbon metabolism, spina bifida aperta and spina bifida cystica are associated with apoptosis; lumbar sacral spina bifida aperta and spina bifida occulta are associated with lipid metabolism. By contrast, complex phenotypes in both cranial and spinal regions display association with various biological functions given the different phenotypes.


Our study links genetic variant to subphenotypes of human NTDs and provides a preliminary but direct clue to investigate pathogenic mechanism for human NTDs.


Neural tube closure (NTC) is a highly orchestrated process wherein the neural plate bends to form the neural folds, which then meet and fuse to close the neural tube. In mammals, NTC initiates sequentially at different levels of the body axis [1]. Primary neurulation in mice initiates at the hindbrain/cervical boundary, called closure I, and then NTC proceeds by zippering together of the neural folds in a bidirectional manner into the hindbrain and along the spinal region [2]. Additional closure points are initiated at the rostral extremity of the forebrain (closure 3), and can also occur at the midbrain-forebrain boundary (closure 2). Fusion progresses along the spine, culminating in final closure at the posterior neuropore, at the level of the second sacral segment. Formation of the spinal cord at lower sacral and caudal levels is accomplished by a different process called “secondary” neurulation in which condensed mesenchyme hollows out to form a tube. Neuromesodermal progenitors biomechanically couple zippering point and drive caudal NTC [3]. In human embryos, neurulation events have been described that correspond to closure 1 and 3 in the mouse, whereas closure 2 may not occur in human embryos [4].

Perturbations of NTC can lead to neural tube defects (NTDs), a common and severe birth defect. Failure to close the neural tube in the cranial region (exencephaly or called anencephaly after degradation of the exposed neural tissue) leads to death before or at birth. Infants born with caudal NTDs (spina bifida) have increased risk of mortality, and those that survive often face life-long disabilities and neurological, cognitive, urologic, and gastrointestinal complications. A fundamental principle of neurulation is that the levels of the body axis that undergo primary neurulation are susceptible to “open” NTDs (for example, anencephaly, spina bifida aperta and craniorachischisis). By contrast, defective secondary neurulation leads to “closed” forms of spina bifida (also called “dysraphism” condition/spina bifida occulta), which represent the failure of the emerging spinal cord to become separated from other tissue derivatives of the tail bud. Failure of closure 1 leads to the most severe NTD, craniorachischisis, in which the neural tube is open throughout the midbrain, hindbrain and entire spinal region. If closure 1 is completed but closure of the cranial neural tube is incomplete, this leads to anencephaly. Failure of closure 3 is uncommon but, when present, yields split face with anencephaly. In the spinal region, failure of final closure at the posterior neuropore yields open spina bifida (Spina bifida aperta/Myelomeningocele), in which the upper limit can vary as to the axial level, depending on precisely when the progression of zippering becomes arrested.

In addition to the complex morphogenetic movements described above, numerous cellular and molecular functions must be tightly regulated to coordinate proliferation, differentiation, apoptosis, patterning and cell shape changes. Indeed, animal models have revealed over 300 genes that regulate the process of NTC, as disruption in the function of each gene can cause NTD. Animal models have also revealed how some of these genes contribute to NTC [1, 5, 6]. Narrowing and lengthening of the neural plate requires a process called convergent extension, wherein cell intercation mediolaterally in order to elongate the tissue along the rostrocaudal axis. Planar cell polarity (PCP) pathway is well documented to function in convergent extension [7,8,9] and mutations in PCP genes lead to failure to undergo closure 1, resulting in craniorachischisis [10]. Signaling through the BMP, Sonic hedgehog (Shh), FGF, Wnt pathways and cilia-related genes coordinate patterning of the neural tissue [11, 12]. Moreover, Shh signaling emanating from the notochord regulates dorsal lateral hinge point formation [13] and disruption of Shh signaling can lead to both spina bifida and exencephaly [1]. Apical constriction converts columnar cells into wedge-shaped cells, which involves cytoskeleton proteins such as Shroom3, Abl, and Mena [14]. Normal cytoskeletal function is required for cranial and spinal NTC and deficiency leads to exencephaly and spina bifida [1, 15]. Cell cycle is crucial for NTC and the balance between continued proliferation and neuronal differentiation may be critical for successful closure [1]. Live imaging has shown programmed cell death in cranial regions during neurulation [16] as well as dynamic cell behaviors and cell extensions during fusion [2, 15]. Prevention of neural fold fusion leads to hindbrain to forebrain exencephaly and thoracolumbosacral spina bifida [17,18,19]. Epigenetic modifiers, such as DNA methylation, chromatin modification and nucleosome assembly, contribute to proper closure of the mammalian neural tube [5]. Additionally, genes functioning in folate one carbon metabolism and glucose metabolism, as well as multiple vitamins and minerals are essential for NTC [20,21,22]. Despite the wealth of knowledge from animal models, there remains a large gap in knowledge relative to the association of these pathways with NTD phenotypes in humans.

Direct mutation screening of candidate genes has been carried out in cohorts of patients, largely using case-control association studies. NTDs are considered to be a complex disease with polygenic and multi-factorial etiology, and recent studies indicate that rare allele variation can have a large effect size of causation with respect to complex diseases [23]. Mutations in PCP pathway core genes VANGL1 and VANGL2 were explored first in craniorachischisis cases and then other NTD patients [24, 25], and the contribution of PCP gene mutations to human NTDs is well-established in several cohorts [26,27,28]. Another major emphasis has been on the evaluation of folate-related genes. A mutation screening study on rare variations within 31 folate-related genes in 480 NTD case-control population uncovered ethnic-specific risk signatures for spina bifida [29]. However, this has given inconsistent results between cohorts and populations, indicating that very few of the tested genes have a major causative effect [30]. We and others have screened a subset of retinoid-related genes in a range of NTD cases and found loss of function rare variants in retinoic acid degradative enzyme encoded gene CYP26B1, CYP26A1 and retinoic acid receptor encoded genes [31,32,33].

To gain a broader understanding of the genetics of NTDs, in the present study we sequenced 280 NTC-related genes in 355 NTDs cases with multiple phenotypes and 225 ethnicity-matched controls in a high throughput manner. Our purpose was to identify rare mutations in a broad set of genes found to be critical for NTC in animal models and to then determine whether there is a relationship with human NTD phenotypes. This has allowed us to establish associations between molecular functions and some NTD sub-phenotypes in human.

Subjects and methods


As described in our previous paper [33, 34], high-throughput DNA sequencing were performed on genomic DNA samples collected from 355 subjects with NTD in a Han Chinese population ranging in age from gestational week (GW) 12 to 10-years old and from multiple local hospitals in six provinces in China. Two hundred twenty-five ethnicity-matched controls were collected from non-medically related terminations and were free of any NTDs. The enrolled pregnant women were diagnosed by local clinicians using ultrasonography. Individuals with NTDs that had been assessed by clinical geneticists were enrolled and were placed into at least one of the following diagnostic groups: craniorachischisis, anencephaly, encephalocele, spina bifida (aperta, cystica, or occulta).

The study was approved by the Committee of Medical Ethics of the Capital Institute of Pediatrics (Beijing, China). Written informed consent was obtained from the parents. We carried out the study in accordance with The Code of Ethics of the World Medical Association (Declaration of Helsinki) for experiments involving humans and in accordance with the approved guidelines.

Genomic DNA sequencing

Genomic DNA was extracted using a Truseq DNA Sample preparation kit (Illumina Inc., San Diego, CA), libraries constructed with Agilent Custom SureSelect Enrichment Kit, and run on an Agilent Custom enrichment array (Probe Code: BI426526171). Sequence reads were aligned to the UCSC human genome GRCh37/hg19 using BWA (v0.5.9) [35]. The average depth of coverage in the present sequencing was 22.5X. Based on cDNA sequence, nucleotide numbering uses + 1 as the A of the ATG translation initiation codon in the reference sequence, with the initiation codon as codon 1. Single-nucleotide variants (SNVs) were called using GATK (Samtools pileup (MAPQ30))(version 0.1.17) [36] and Varscan [37] (minimum coverage = 1, minimum alterative allele reads = 1, minimum variation frequency > 0.03) and short indels (insertions and deletions) were called using Varscan with a standard (minimum coverage = 2, minimum alternative allele reads = 2, minimum variation frequency > 0.1). Genotypes were called using Bayescall [38]. Variants were annotated with ANNOVAR [39]. Rare variants were filtered out using the dbSNP in NCBI, the 1000 Genomes Project, and the shared variants in NTDs cases and controls were also excluded.


A wide range of NTD phenotypes in the current cohort

This study focuses on 355 NTD cases and 225 ethnicity-matched controls enrolled in China during 2005–2011. The NTD phenotypes were described as craniorachischisis (anencephaly continuous with exposed spinal cord), anencephaly (lack of brain and cranial vault subsequent to failure to close the cranial neural tube), encephalocele (meningeal sac, often containing brain tissue, protruding from the skull), spina bifida aperta (exposed spinal cord, usually called meningomyelocele), spina bifida cystica (spinal cord defect covered by meningeal sac), spina bifida occulta (skin-covered lesion involving two or more vertebrae, also called spina dysraphism) [40]. As shown in Fig. 1, 74 out of the 355 NTD cases showed cranial deficits only. This included 17 affected with anencephaly and 57 cases with encephalocele. Amongst encephalocele cases, 2 were located in the frontal and occipital regions, 35 were located in the occipital region, and 9 were in the frontal region. Of the remaining 11 encephalocele cases, detailed clinical information was not available (labeled encephalocele no details).

Fig. 1
figure 1

Category for neural tube defects cases in the present study. Phenotypes are arranged as order of rostral to caudal from left to right. The numbers in lower frames indicate summarized numbers of phenotypes

The number of cases with NTD limited to the spinal region was 178. The majority of these were spina bifida aperta (89 cases) which occurred in lumbar-sacral segments (50/89), thoracic-lumbar-sacral segments (20/89) and thoracic-lumbar segments (11/89) (Fig. 1). Of the 26 spina bifida occulta cases, the majority involved lumbar-sacral segments (13/26) and thoracic-lumbar segments (5/26). Many of spina bifida cystica cases did not have detailed information on the segments involved (26/44), hence we recruited 10 lumbar-sacral spina bifida cystica. Spina bifida cystica was not found in cervical segments (Fig. 1).

More complex phenotypes affecting both the cranial and spinal regions were also observed. This included 19 cases of craniorachischisis and 57 cases of anencephaly combined with spina bifida, the majority being anencephaly with upper segment spina bifida aperta (from thoracic to occipital). Four cases exhibited anencephaly and lower level spina bifida aperta encompassing the lumbar or sacral segments (Fig. 1). The remaining 27 cases were encephalocele combined with spina bifida aperta (Fig. 1; in 2 cases there was no detailed clinical information).

A summary of the phenotypes in the current cohort is 74 cases with only cranial malformations, 178 cases with only spinal phenotypes, and 103 cases with complex phenotypes in both cranial and spinal regions. Viewing all cases together, there are 177 with cranial NTDs and 281 with spinal NTDs (Fig. 1).

High throughout targeted sequencing of neural tube closure (NTC)-related genes

To gain an understanding of the genetic contribution to a broad set of NTD phenotypes, we undertook a population-based case-control mutation screening study to reveal potential disease-causative rare variants by sequencing a selected set of 280 genes implicated in NTC from human and animal studies [1, 5, 6, 30, 40]. These 280 genes encompassed more than twenty-one signaling pathways and key developmental pathways (upplementary Table 1). For these 280 genes we sequenced the coding regions (CDS), splice junctions, and 2000 bases upstream of the genes for a total of 4.38 Mb of the human genome. The CDS length sequenced for each pathway is shown in Supplementary Table 1. We applied criteria in calling cohort-specific variants by eliminating from further analysis shared rare variants found in both NTD and control samples. SIFT and Polyphen2 were used to identify putatively damaging SNVs [41, 42]. All missense mutations, together with splicing mutations, frameshift mutations or nonsense mutations were defined as putatively damaging rare variants (PDRVs) in the present study. By these criteria we identified NTDs case-specific 791 PDRVs in 213 genes and 387 common PDRVs in 150 genes only in controls, all PDRVs were heterozygous.

PDRVs in the pathways of chromatin modification, apoptosis, retinoid metabolism and lipid metabolism are enriched in NTDs

To explore the possible genetic contribution of each signaling pathway to human NTDs, we first separately analyzed the number of PDRVs in NTDs and controls within each pathway. The data indicate that the genes which function in chromatin modification, apoptosis, retinoid-related and lipid metabolism harbor significantly more NTD-specific PDRVs (Fig. 2). Regarding genes that function in chromatin modification, we found 75 PDRVs in NTD cases which was significantly more than in controls (75/355 in NTDs vs. 27/225 in controls; P = 0.019; two-sided Fisher’s exact test) (Fig. 2a). For genes that function more specifically in chromatin remodeling (ACTL6A, CECR2, SMARCA4, SMARCC1 and ATRX) there were 22 PDRVs in 355 NTD cases, which was slightly but significantly more than in controls (5/225; P = 0.041, two-sided Fisher’s exact test).

Fig. 2
figure 2

Four pathways significantly associated with human neural tube defects. In the present study, PDRVs in genes functioning in four pathways have been explored being associated with the current NTD cohort, they are (a) chromatin modification; (b) apoptosis; (c) retinoid metabolism related genes; and (d) lipid metabolism. The occurrences of missense variants are shown as green font; of splicing variants are shown as blue; of frameshift variants are shown as purple. The background color of different phenotypes are same as in Fig. 1. (e) Polygenic or digenic genetic PDRVs were observed in the present NTD cohort. The columns represent the number of NTD cases who carries certain numbers of PDRVs in different genes

Defective apoptosis can lead to cranial NTD in mouse [16]. Here we analyzed 14 genes implicated in the apoptosis pathway and found 28 PDRVs in NTD cases, significantly enriched relative to controls (28/355 vs. 4/225, P = 0.0023, two-sided Fisher’s exact test) (Fig. 2b). Interestingly, the intrinsic pro-apoptotic genes implicated in exencephaly in mice TP53 [43], APAF1 [44], CASP9 [45] and CASP3 [46] harbored 16 PDRVs in human NTD cases, considerably increased relative to controls (16/355 vs. 1/225, P = 0.0044, two-sided Fisher’s exact test) (Fig. 2b). By contrast, we did not find a PDRV in the anti-apoptotic gene BCL10 in our cohort (Fig. 2b), although its mutant show exencephaly in mice [47].

Retinoic acid is a small lipophilic molecule derived from vitamin A via the retinoid metabolism pathway and gene mutations in this pathway cause NTD phenotypes in mice [48]. As we previously reported [33], in six retinoid metabolism genes we found 12 occurrences of PDRVs in NTD cases versus 0/225 in controls (P = 0.0046, two-sided Fisher’s exact test) (Fig. 2c) and PDRVs were predominantly in the CYP26B1 gene (9/12), which acts to degrade retinoic acid. In the three genes sequenced in the lipid metabolism pathway, there were significantly more PDRVs in NTD cases (36/355 in NTDs vs. 10/225 in controls, P = 0.025) (Fig. 2d). Of note, a frameshift variant in the APOB gene was found in control (NM_000384: exon26:c.10373delT; p.M3458fs) but overall the APOB gene harbored more PDRVs in NTD cases (28/355) versus controls (10/225). PDRVs in the LEPR gene were only found in NTD cases (8/355 vs. 0/225, P = 0.026, Fig. 2d). In total, in the current NTD cohort, polygenic PDRVs were carried by 137 cases, while 97 cases were digenic, by contrast, no PDRV in the targeted genes was found in 48 cases (Fig. 2e). These confirm polygenic or digenic genetic features for human NTDs, however, it is also necessary to extend the target list in the future.

In summary, evaluation of our current set of 280 genes and 21 pathways showed that four signaling pathways were significantly enriched for NTD-specific PDRVs. Furthermore, PDRVs in genes of these four pathways were associated with 34% (120/355) of NTD cases. These data suggest the importance of PDRVs in these pathways in the pathogenesis of NTDs and highlight these biological functions as critical for human NTC.

Craniorachischisis is associated with PDRVs in PCP signaling and new pathways of chromatin modification and retinoid metabolism

We next asked whether we could discern possible relationships between more specific NTD phenotypes and variants in molecular pathways. This was possible as our present cohort encompasses a range of NTD sub-phenotypes including cranial, caudal, and complex NTD as shown in Fig. 1. We started with an analysis of craniorachischisis, the most severe NTD phenotype and which has been associated with defects in genes within the PCP pathway in both mice [5] and human [49]. Our sequencing study encompassed 18 PCP genes, including the core PCP genes VANGL1/2, DVL1/2/3, CELSR1 and PRICKLE1/2 [50]. The results showed that there were 69 PDRVs in PCP genes in the total cohort of 355 NTD cases versus 34 PDRVs in 225 controls, thus not reaching statistical significance when viewed across the range of NTDs cases (Fig. 3a). However, as expected, PDRVs in PCP genes were significantly associated with craniorachischisis (9/19 in craniorachischisis cases vs. 34 in all 225 controls; P = 0.0210; two-sided Fisher’s exact test, Fig. 3b). Within this smaller cohort of 9 craniorachischisis cases, there were 9 different occurrences of PDRVs in PCP pathway genes (yellow background in Fig. 3c highlights these 9 cases). Thus, our data support previous findings that the PCP pathway genetically contributes to human craniorachischisis. It is of note that some craniorachischisis do not carry PDRVs in the set of PCP pathway genes sequenced, either suggesting the presence of variants in untranslated regions of the sequenced PCP genes or in other PCP genes that were not sequenced or alternative mechanism(s) underlying this pathology. In support of the latter possibility, we found a relationship between PDRVs in genes involved in chromatin modification and retinoid metabolism that were associated with the craniorachischisis phenotype (Fig. 3b, pink and blue highlight in Fig. 3c). Interestingly, the degradation enzyme for retinoic acid encoded by the CYP26B1 gene harbors significant number of PDRVs in craniorachischisis cases (2/19 in craniorachischisis vs. 0/225 in all controls; P = 0.007; Fig. 3b and c). Therefore, our sequencing of even a limited number of NTD genes reveals new candidate genes related to chromatin modification and retinoid metabolism associated with craniorachischisis. Furthermore, our sequencing uncovered additional genes and pathways (Fig. 3c) that may play causal roles and are worthy of consideration in other cases of this most severe NTD, craniorachischisis.

Fig. 3
figure 3

Planar cell polarity pathways and craniorachischisis. In the present cohort PDRVs occur in multiple NTD phenotypes but are significantly associated with craniorachischisis (a). The background color of different phenotypes are same as in Fig. 1. Interestingly, except for planar cell polarity pathway, craniorachischisis phenotype is statistically associated with PDRVs in chromatin modification and retinoid metabolism genes (b), fisher exact test was used to calculated P value, and numbers in brackets means occurrences of PDRVs in each gene. In (c) we list all 19 craniorachischisis in the present study and the PDRVs we found in these cases. Yellow background means genes in planar cell polarity pathway; purple background represents genes functioning in chromatin modification; sky blue background represents retinoid metabolism related genes. In all panels, numbers or gene symbols in green represent missense PDRVs occurred; in blue represent splicing PDRVs; in purple represents frameshift PDRVs and in red represents nonsense PDRVs

Cranial NTD phenotypes are associated with PDRVs in chromatin modification, retinoid and glucose metabolism, apoptosis and one-carbon metabolism

We next addressed whether particular pathways may be more closely associated with cranial NTD sub-phenotypes. Statistical analysis of the 177 NTD cases that include a cranial phenotype (Fig. 1) indicated that PDRVs in genes in chromatin modification, apoptosis and retinoid metabolism were associated with prevalence of cranial NTD phenotypes with or without spinal NTD phenotypes (P = 0.024; P = 0.011 and P = 0.0016, respectively; Fig. 4a). Regarding cranial-only NTD phenotypes, PDRVs were enriched in genes involved in chromatin modification and retinoid metabolism (P = 0.0054 and P = 0.0011; Fig. 4b). In 17 cases affected with only anencephaly (Fig. 4c), PDRVs were enriched in chromatin modification genes as well as in glucose metabolism genes, particularly the case D102 carries two different PDRVs in the INSR gene (INSR: NM_000208: c.1882C > G:p.P628A and c.4028G > A:p.R1343Q).

Fig. 4
figure 4

Summary on PDRVs burden genes related functional pathways in cranial NTDs. (a) Summary of PDRVs in all 177 NTDs cases who display NTDs phenotype in cranial regions with or without spinal regions. (b) Summary in cases who have only cranial NTD phenotypes. (c-f) show summaries of PDRVs burden genes and their functional pathways in cases affected with only anencephaly (c); only encephalocele (d); only frontal encephalocele (e); only occipital encephalocele (f). Numbers in brackets means occurrences of PDRVs in each gene, numbers or gene symbols in green represent missense PDRVs; in blue represent splicing PDRVs; in purple represents frameshift PDRVs. Gene symbols linking with colon represents two concurrent PDRVs in case

In all 57 encephalocele cases, PDRVs in genes in chromatin modification, apoptosis, retinoid metabolism and one carbon metabolism were found significantly associated with the phenotype prevalence (Fig. 4d). Interestingly, three pathways were related to encephalocele in frontal regions (Fig. 4e), whereas for occipital encephalocele the PDRVs were statistically accumulated only in genes functioning in retinoid metabolism (Fig. 4f), correlating with a role for retinoid signaling in hindbrain patterning [51]. In general, it appears that PDRVs in genes in chromatin modification can contribute to diverse cranial phenotypes with or without spinal phenotypes, including anencephaly and encephalocele. In contrast, PDRVs in genes in retinoid metabolism are associated with all types of encephalocele, but not anencephaly, at least in our cohort of patients. Finally, PDRVs in genes in apoptosis and one carbon metabolism are significantly associated with frontal encephalocele or encephalocele only (Fig. 4d and e).

Spina bifida phenotypes are associated with PDRVs in apoptosis, retinoid and lipid metabolism, and cytoskeleton genes

Spina bifida remains a common NTD phenotype worldwide and we enrolled a total of 281 cases with spina bifida, of which 178 had only spinal NTDs (Fig. 1). When viewed as the entire cohort of 281 spina bifida with or without cranial phenotype, there was significant enrichment for PDRVs in apoptotic pathway genes (23/281 vs. 4/225 controls; P = 0.0023, two-sided Fisher’s exact test) and lipid metabolism genes (32/281 vs. 10/225 controls; P = 0.0093, two-sided Fisher’s exact test) (Fig. 5a). However, when the 178 cases with only spinal NTDs were considered, the association with lipid metabolism genes was not significant. Instead, these 178 cases were significantly associated with PDRVs in apoptosis and retinoid metabolism genes (15/178 and 4/178, respectively in spina bifida cases vs. 4/225 and 0/225 in controls, Fig. 5b). In addition, PDRVs in apoptosis genes are related to spina bifida aperta and spina bifida cystica (P = 0.0079 and P = 0.030, Fig. 5c and d). Taken together, PDRVs in apoptosis genes seem to be highly associated with cases of spina bifida aperta and spina bifida cystica.

Fig. 5
figure 5

Summary on PDRVs burden genes related functional pathways in spinal NTDs. Summary in cases affected with spinal NTD phenotype whatever who are with or without cranial NTD (a), or affected with only spinal NTD phenotypes (b), or with only spina bifida aperta (c), or with spina bifida cystica (d). Numbers in brackets means occurrences of PDRVs in each gene, numbers or gene symbols in green represent missense PDRVs; in blue represent splicing PDRVs; in purple represents frameshift PDRVs. Gene symbols linking with colon represents two concurrent PDRVs in case

Because we had a large cohort of spina bifida cases that varied in position of the affected segments as well as severity, we next evaluated our data within specific caudal NTD sub-phenotypes. Lumbar sacral spina bifida is thought to result from failure of primary neurulation at the caudal end and failure of secondary neurulation [1]. Our present cohort included 73 cases with only lumbar sacral spina bifida and these showed an association only with lipid metabolism genes (Fig. 6a). PDRVs in the APOB gene occurred in a significant number of lumbar sacral spina bifida cases (10/73 vs. 10/225 controls, P = 0.018, two-sided Fisher’s exact test), although Apob knockout mice exhibit only exencephaly [52]. We then separated the lumbar sacral spina bifida cases into spina bifida aperta, spina bifida cystica, and spina bifida occulta. Lumbar sacral spina bifida aperta is the most widespread NTD phenotype in human. In 50 cases, we found PDRVs in genes functioning in cytoskeleton and lipid metabolism (8 PDRVs in the APOB gene) were significantly enriched (P = 0.044 and 0.012, respectively; Fig. 6b). The regulation of the actin cytoskeleton is critical in murine NTC [5] and we found PDRVs in many genes that control the actin cytoskeleton including the SHROOM family of proteins. Our data identified a frameshift variant in the SHROOM3 gene in case D183 (NM_020859:c.1782delC:p.N594fs; indicated in orange font in Fig. 6b) and another case (TJ-QS39) which carries a SHROOM4 NM_020717:c.2011G > C: p.A671P variant and a MARCKSL1:NM_023009:c.226G > A: p.A76T variant (Fig. 6b).

Fig. 6
figure 6

Summary on PDRVs burden genes related functional pathways in lumbosacral spina bifida. Summary in cases affected with only lumbar sacral spina bifida (a), or affected with only spinal NTD phenotypes (b), or with only spina bifida aperta (c), or with spina bifida cystica (d). Numbers in brackets means occurrences of PDRVs in each gene, numbers or gene symbols in green represent missense PDRVs; in blue represent splicing PDRVs; in purple represents frameshift PDRVs. Gene symbols linking with colon represents two concurrent PDRVs in case

In the 10 lumbar sacral spina bifida cystica cases, we found PDRVs in 3 apoptotic genes (P = 0.0038; Fig. 6c). For the 13 lumbar sacral spina bifida occulta cases, PDRVs in genes in lipid metabolism were weakly correlated (P = 0.041, Fig. 6d). Although the number of NTD sub-phenotype cases is too few to draw strong inferences, it is interesting that lumbar sacral spina bifida occulta was related to PDVRs in LEPR gene (2/13 vs. 0/225, P = 0.0037), whereas PDRVs in LEPR were not associated with lumbar sacral spina bifida aperta.

Within our enrolled cases there were a few other spinal sub-phenotypes encompassing other segments. In eleven cases affected by thoracic lumbar spina bifida aperta there were 3 PDRVs in the lipid metabolism genes LEPR and APOB (P = 0.029; Fig. S1A). Interestingly, the case D126 carried two APOB PDRVs (APOB:NM_000384:c.11230C > G: p.L3744V and c.7777A > G: p.I2593V). Within our two cases of thoracic spina bifida aperta-only, one case carried two PDRVs in genes involved in nuclear pore complex function (NUP98:NM_139132:c.4837C > T: p.R1613C and c.5026C > A: p.H1676N) (Fig. S1B). Our five cases of thoracic lumbar spina bifida occulta harbored PDRVs in genes in chromatin modification (4/5 vs. 27/225 controls) and neural development (3/5 vs. 12/225 controls; P = 0.013 and P = 0.008, respectively; Fig. S1C). Three cases of lumbar spina bifida cystica revealed one PDRV in the retinoid metabolism gene CYP26B1 (1/3 vs. 0/225, P = 0.017; Fig. S1D). In the one case of thoracic spina bifida cystica there was a PDRV in apoptosis gene CASP9 (P = 0.042; Fig. S1E).

Complex NTDs affecting both cranial and spinal regions are associated with PDRVs in a range of functional pathways

Our present cohort also encompasses complex phenotypes which affect both cranial and spinal regions. In the 40 cases with anencephaly and occipital cervical spina bifida aperta, the sequencing data revealed 7 PDRVs in lipid metabolism genes (P = 0.011, Fig. 7a). This included the case D98 which carries two heterozygous APOB missense variants (APOB: NM_000384 c.11230C > G: p.L3744V and c.7777A > G: p.I2593V)(Fig. 7a). Interestingly, we had one case with anencephaly and thoracic lumbar spina bifida aperta (the neural tube was closed in cervical region), who carried two missense mutations in the MAP3K4 gene involved in MAPK signaling (MAP3K4: NM_005922: c.877G > A: p.D293N and c.2711C > G: p.T904S)(Fig. 7b). In mouse Map3k4 expresses along the edges of neural folds and the mutant mice display exencephaly from forebrain to hindbrain, curly tail and lower spina bifida phenotypes or combination of exencephaly and spina bifida/curly tail [53, 54]. In two cases with anencephaly and spina bifida cystica, missense mutations were found In two Notch signaling pathway genes; one in the E3 ubiquitin-protein ligase MIB2, which positively regulates the Delta-mediated Notch signaling, and the other in JAG1, a ligand for multiple Notch receptors (2/2 vs. 16/225 controls, P = 0.028; Fig. 7c). In 27 cases of encephalocele and spina bifida there were 3 PDRVs in apoptotic genes (P = 0.035; Fig. 7d). In an unique case with frontal occipital encephalocele and occipital cervical spina bifida aperta, there was a PDRV in the pro-apoptotic gene APAF1 (APAF1:NM_181868:c.1004C > G: p.P335R; P = 0.042; Fig. 7e). Notably, in 2 cases of frontal encephalocele and occipital cervical thoracic spina bifida apertas, there were 4 PDRVs in four genes functioning in ECM and adhesion (4/2 vs. 46/225 controls; P = 0.010; Fig. 7f). This included concurrent PDRVs in integrin alpha6 encoded gene ITGA6 and a basement membrane protein heparin sulfate proteoglycan 2 encoded gene HSPG2 in one individual; and concurrent PDRVs in two cell adhesion genes NCAM1 and CDON in another individual. This suggests a possible synergistic effect of ECM and adhesion on the closure events disrupted in cases of frontal encephalocele and occipital cervical thoracic spina bifida aperta. In 5 cases of frontal encephalocele and occipital cervical spina bifida aperta, we observed 3 PDRVs in the cytoskeleton genes SHROOM2, SHROOM3, and VCL encoding VINCULLIN (P = 0.040; Fig. 7g). There was a relationship to PDRVs in neural development genes with the complex phenotypes of occipital encephalocele and spina bifida represented by occipital encephalocele and holorachischisis (spina bifida aperta in all spinal regions, 3 PDRVs in 7 cases vs. 12/225) and occipital encephalocele and occipital, cervical and thoracic spina bifida aperta (2/3 cases vs. 12/225) (P = 0.016 and P = 0.028, respectively; Fig. 7h and i). TRPM6, encodes a channel protein crucial for magnesium homeostasis and mutations in Trpm6 in mice result in exencephaly and spina bifida occulta [55]. We found 2 PDRVs in TRPM6 in two cases of occipital encephalocele and occipital, cervical and thoracic spina bifida aperta and 1 PDRV in one case of occipital encephalocele and holorachischisis (P = 0.0004 and P = 0.0343; respectively, Fig. 7h and i).

Fig. 7
figure 7

Summary on PDRVs burden genes related functional pathways in cranial and spinal complex phenotypes. Summary of PDRVs burden genes in cases affected with anencephaly and spinal NTD phenotypes (a-c), or encephalocele and spinal NTD phenotypes (d-i). Numbers in brackets means occurrences of PDRVs in each gene, numbers or gene symbols in green represent missense PDRVs; in blue represent splicing PDRVs; in purple represents frameshift PDRVs. Gene symbols linking with colon represents two concurrent PDRVs in case


The association of craniorachischisis and PCP gene variants should be highlighted in human NTD genetics [34, 49, 50]. In the present study, our data indicated that in addition to PCP genes, variants in retinoic acid related gene CYP26B1 and chromatin modification genes are associated with craniorachischisis prevalence. This suggests that genetic variants in PCP genes are not unique causative factors for craniorachischisis in human and reveals new genes for consideration for this severe NTD. Indeed, mouse embryos deficient for glycine decarboxylase (Gldc) exhibit craniorachischisis as well [22]. The CYP26 proteins along with retinaldehyde dehydrogenase control the regional levels of retinoic acid. Cyp26b1, together with Cyp26a1 and Cyp26c1, corrects Nodal expression during gastrulation, when retinaldehyde dehydrogenase is not expressed [56]; mice lacking all three CYP26 genes manifest duplication of neural tube which is more pronounced than Cyp26a−/c- double mutant mice [56], emphasizing the crucial role of Cyp26b1. Chromatin modification genes are necessary for NTC, although mouse embryos carrying mutations in these genes mainly exhibit exencephaly [1]. Nonetheless, multiple chromatin modification genes such as Dnmt3b and Cecr2 show expression patterns distributed at the leading edge of the neural folds along the rostrocaudal axis during neurulation [57, 58]. Kat2a is also expressed widely along the cranial and spinal neural tube [59]. Jarid2, a regulator of histone methyltransferases complex, is expressed at the midbrain-hindbrain boundary and distributed in spinal regions as well [60]. These suggest that these molecules potentially play roles in both spinal and cranial regions. Moreover, in our present results, all variants in chromatin modification genes and CYP26B1 gene are concurrent with other variants in neural tube-related genes (Fig. 3c), so we hypothesize that the craniorachischisis phenotype is also resulted from the synergistic effect of NTC-related polygenic functional variants.

Rare variants in chromatin modification genes, apoptosis genes and retinoid metabolism genes significantly contribute to cranial NTD phenotypes, including anencephaly and encephalocele (Fig. 4), further emphasizing the importance of these signaling pathways in cranial neurulation. The morphology of the elevating neural folds differs between cranial and spinal regions in mammals. In the mouse embryonic midbrain, the neural folds are initially bi-convex with their tips orientated away from the midline. Then, dorsolateral bending occurs, generating bi-concave neural folds and orienting the tips towards the midline for fusion [50, 61]. In addition, the emigration of cranial neural crest may enable dorsolateral bending [4]. In contrast, spinal region does not exhibit a bi-convex elevation phase and the neural folds remain straight except for focal bending sites at the midline and dorsolaterally; and neural crest emigration in the spinal region does not begin until several hours after NTC is complete [62]. In addition, during cranial neurulation, a functional actin cytoskeleton, emigration of the cranial neural crest, spatio-temporally regulated apoptosis, and a balance between cell proliferation and the onset of neuronal differentiation are all required for normal dorsolateral bending [50]. Moreover, closure of the cranial neural tube is essential not only for maintenance of brain development but also for initial formation of much of the skull, with contributions from both cranial mesenchyme and cranial neural crest [4, 63]. Coordinate regulation of transcriptional networks by epigenetic regulators is absolutely crucial for proper cranial NTC [5].

In our present NTD cohort, variants in one carbon metabolism are significantly associated with encephalocele, in particular frontal encephalocele (Fig. 4d and e). A previous study of 239 newborns with spina bifida and 241 non-malformed controls reported a Hispanic genetic risk profile pointing to alterations in genes functioning in purine biosynthesis, whereas that in non-Hispanic whites implicated in homocysteine metabolism [29]. In glycine cleavage system genes, rare variants are associated with multiple phenotypes including spina bifida and anencephaly [22, 64]. Therefore, except for a possible causation of different ethnicities, association of one carbon metabolism to NTD subphenotypes is still to be determined in a larger number of populations.

Variants in the lipid metabolism genes APOB and LEPR appear to contribute to lumbar sacral spina bifida in the present study (Fig. 6). NTC in lumbar sacral regions is attributed to secondary neurulation. Therefore, our data provides clues that lipid metabolism is essential for this process in human. APOB encoded apolipoproteins B is a major structural component of very low density lipoproteins, intermediate density lipoprotein, low density lipoprotein, chylomicrons and lipoproteins. In Apob knockout mice or mice carrying a dysfunctional truncated Apob protein, homozygous embryos at gestational day 9.5 were either runted or appeared to be nonviable, and the few viable homozygous embryos at gestational day 10.5 appeared to be exencephalic. Moreover, even heterozygous mice had significant reductions in total plasma cholesterol and HDL cholesterol compared with wildtype mice [52, 65, 66]. Thus, APOB can have widespread effects and heterozygous allelic variants can reduce biological function. Similar to our results, LEPR rs1137100 G allele showed a significant increase in the risk of an NTD-affected offspring when inherited from the mother (2.43 fold), and an increased risk for lower spina bifida aperta (Lumbar 1 and lower, 3.20 fold) [67]. A study in Ireland of 520 spina bifida cases and their families and 994 controls found over-transmission of the LEPR rs1805134 minor C allele associated with spina bifida (relative risk of 1.5 fold) [68]. Taken together with our data, there is a significant association of genetic variants in lipid metabolism genes and human lumbar sacral spina bifida aperta.

Variants in apoptosis genes are associated with cranial and spinal NTDs, particularly in the frontal encephalocele, spina bifida aperta and spina bifida cystica (Figs. 4 and 5). Variants in both pro-apoptotic (APAF1, TP53, CASP3 and CASP9 etc.) and anti-apoptotic genes (NFKB1, IKBKB etc.) were found in the present study, consistent with the idea that too little or too much cell death can disrupt the fine regulation of NTC. Indeed, mice carrying mutations in apoptosis genes exhibit exencephaly and spina bifida [1, 69] associated with excess neural cells [69]. On the other hand, excessive apoptosis in the lumbar sacral neuroepithelium has been observed with retinoic acid-induced spina bifida [70]. In human NTDs cases, more apoptotic cells were observed in the central nervous tissue [71]. Higher levels of p53 were observed in anencephaly cases [71]. Cleaved caspase 3 levels were elevated in encephalocele cases and cleaved caspase 8 levels were higher in spina bifida cases relative to controls [71]. Even under harsh environment, failure of NTC is probably due to p53 stabilization and excessive apoptosis [72]. Altogether, abnormal apoptosis might be a common causal factor for cranial and spinal NTDs and abnormal apoptosis might contribute to a wide variety of human NTDs.

Adhesion is key process for fusion of the neural folds and genetic disturbance that alters the adhesion of neural folds should prevent NTC [10]. In two frontal encephalocele and occipital cervical thoracic spina bifida aperta cases, two pairs of rare variants in ECM and adhesion genes were found (Fig. 7f). Integrin Itga6/Itga3 double mutant mice or Hspg2 mutant embryos exhibit exencephaly [73, 74]. We observed PDRVs in both ITGA6 and HSPG2 genes in one case. In a study in American Caucasian simplex lumbar sacral spina bifida aperta families, a nucleotide polymorphism in NCAM1 gene may influence NTD risk [75]. Similarly, we found variants in NCAM1 and CDON adhesion genes in some cases.


Through sequencing NTC-related genes we explored associations of PDRVs in signaling pathways to subphenotypes of human NTDs. Our studies add a significant body of data relative to the genetic causes of human NTDs. Our identification of rare variants in a wide range of genes in a relatively large cohort of NTD cases and controls provides robust insight into plausible critical genetic factors that regulate human neurulation. Moreover, our studies contribute new insight into NTD phenotype-genotype relationships in human. These results lay a foundation for a more detailed understanding of human NTC and how variations in genetic factors may lead to NTDs.

Availability of data and materials

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


  1. Copp AJ, Greene ND. Genetics and development of neural tube defects. J Pathol. 2010;220(2):217–30.

    CAS  PubMed  PubMed Central  Google Scholar 

  2. Massarwa R, Niswander L. In toto live imaging of mouse morphogenesis and new insights into neural tube closure. Development. 2013;140(1):226–36.

    CAS  PubMed  PubMed Central  Google Scholar 

  3. Galea GL, Cho YJ, Galea G, Mole MA, Rolo A, Savery D, Moulding D, Culshaw LH, Nikolopoulou E, Greene NDE, et al. Biomechanical coupling facilitates spinal neural tube closure in mouse embryos. Proc Natl Acad Sci U S A. 2017;114(26):E5177–86.

  4. Copp AJ. Neurulation in the cranial region--normal and abnormal. J Anat. 2005;207(5):623–35.

    PubMed  PubMed Central  Google Scholar 

  5. Wilde JJ, Petersen JR, Niswander L. Genetic, epigenetic, and environmental contributions to neural tube closure. Annu Rev Genet. 2014;48:583–611.

    CAS  PubMed  PubMed Central  Google Scholar 

  6. Harris MJ, Juriloff DM. An update to the list of mouse mutants with neural tube closure defects and advances toward a complete genetic perspective of neural tube closure. Birth Defects Res A Clin Mol Teratol. 2010;88(8):653–69.

    CAS  PubMed  Google Scholar 

  7. Goto T, Keller R. The planar cell polarity gene strabismus regulates convergence and extension and neural fold closure in Xenopus. Dev Biol. 2002;247(1):165–81.

    CAS  PubMed  Google Scholar 

  8. Wallingford JB, Rowning BA, Vogeli KM, Rothbacher U, Fraser SE, Harland RM. Dishevelled controls cell polarity during Xenopus gastrulation. Nature. 2000;405(6782):81–5.

    CAS  PubMed  Google Scholar 

  9. Davidson LA, Keller RE. Neural tube closure in Xenopus laevis involves medial migration, directed protrusive activity, cell intercalation and convergent extension. Development. 1999;126(20):4547–56.

    CAS  PubMed  Google Scholar 

  10. Copp AJ, Greene ND, Murdoch JN. The genetic basis of mammalian neurulation. Nat Rev Genet. 2003;4(10):784–93.

    PubMed  Google Scholar 

  11. Henrique D, Abranches E, Verrier L, Storey KG. Neuromesodermal progenitors and the making of the spinal cord. Development. 2015;142(17):2864–75.

    CAS  PubMed  Google Scholar 

  12. Stamataki D, Ulloa F, Tsoni SV, Mynett A, Briscoe J. A gradient of Gli activity mediates graded sonic hedgehog signaling in the neural tube. Genes Dev. 2005;19(5):626–41.

    CAS  PubMed  PubMed Central  Google Scholar 

  13. Ybot-Gonzalez P, Cogram P, Gerrelli D, Copp AJ. Sonic hedgehog and the molecular regulation of mouse neural tube closure. Development. 2002;129(10):2507–17.

    CAS  PubMed  Google Scholar 

  14. Wallingford JB, Niswander LA, Shaw GM, Finnell RH. The continuing challenge of understanding, preventing, and treating neural tube defects. Science. 2013;339(6123):1222002.

    PubMed  PubMed Central  Google Scholar 

  15. Rolo A, Savery D, Escuin S, de Castro SC, Armer HE, Munro PM, Mole MA, Greene ND, Copp AJ. Regulation of cell protrusions by small GTPases during fusion of the neural folds. Elife. 2016;5:e13273.

    PubMed  PubMed Central  Google Scholar 

  16. Yamaguchi Y, Shinotsuka N, Nonomura K, Takemoto K, Kuida K, Yosida H, Miura M. Live imaging of apoptosis in a novel transgenic mouse highlights its role in neural tube closure. J Cell Biol. 2011;195(6):1047–60.

    CAS  PubMed  PubMed Central  Google Scholar 

  17. Holmberg J, Clarke DL, Frisen J. Regulation of repulsion versus adhesion by different splice forms of an Eph receptor. Nature. 2000;408(6809):203–6.

    CAS  PubMed  Google Scholar 

  18. Pyrgaki C, Liu A, Niswander L. Grainyhead-like 2 regulates neural tube closure and adhesion molecule expression during neural fold fusion. Dev Biol. 2011;353(1):38–49.

    CAS  PubMed  PubMed Central  Google Scholar 

  19. Ting SB, Wilanowski T, Auden A, Hall M, Voss AK, Thomas T, Parekh V, Cunningham JM, Jane SM. Inositol- and folate-resistant neural tube defects in mice lacking the epithelial-specific factor Grhl-3. Nat Med. 2003;9(12):1513–9.

    CAS  PubMed  Google Scholar 

  20. Blom HJ, Shaw GM, den Heijer M, Finnell RH. Neural tube defects and folate: case far from closed. Nat Rev Neurosci. 2006;7(9):724–31.

    CAS  PubMed  PubMed Central  Google Scholar 

  21. Davidson CM, Northrup H, King TM, Fletcher JM, Townsend I, Tyerman GH, Au KS. Genes in glucose metabolism and association with spina bifida. Reprod Sci. 2008;15(1):51–8.

    CAS  PubMed  Google Scholar 

  22. Leung KY, Pai YJ, Chen Q, Santos C, Calvani E, Sudiwala S, Savery D, Ralser M, Gross SS, Copp AJ, et al. Partitioning of one-carbon units in Folate and methionine metabolism is essential for neural tube closure. Cell Rep. 2017;21(7):1795–808.

    CAS  PubMed  PubMed Central  Google Scholar 

  23. Manolio TA, Collins FS, Cox NJ, Goldstein DB, Hindorff LA, Hunter DJ, McCarthy MI, Ramos EM, Cardon LR, Chakravarti A, et al. Finding the missing heritability of complex diseases. Nature. 2009;461(7265):747–53.

    CAS  PubMed  PubMed Central  Google Scholar 

  24. Kibar Z, Torban E, McDearmid JR, Reynolds A, Berghout J, Mathieu M, Kirillova I, De Marco P, Merello E, Hayes JM, et al. Mutations in VANGL1 associated with neural-tube defects. N Engl J Med. 2007;356(14):1432–7.

    CAS  PubMed  Google Scholar 

  25. Lei YP, Zhang T, Li H, Wu BL, Jin L, Wang HY. VANGL2 mutations in human cranial neural-tube defects. N Engl J Med. 2010;362(23):2232–5.

    CAS  PubMed  Google Scholar 

  26. De Marco P, Merello E, Piatelli G, Cama A, Kibar Z, Capra V. Planar cell polarity gene mutations contribute to the etiology of human neural tube defects in our population. Birth Defects Res A Clin Mol Teratol. 2014;100(8):633–41.

    PubMed  Google Scholar 

  27. Lei Y, Zhu H, Yang W, Ross ME, Shaw GM, Finnell RH. Identification of novel CELSR1 mutations in spina bifida. PLoS One. 2014;9(3):e92207.

    PubMed  PubMed Central  Google Scholar 

  28. Yang XY, Zhou XY, Wang QQ, Li H, Chen Y, Lei YP, Ma XH, Kong P, Shi Y, Jin L, et al. Mutations in the COPII vesicle component gene SEC24B are associated with human neural tube defects. Hum Mutat. 2013;34(8):1094–101.

    CAS  PubMed  Google Scholar 

  29. Marini NJ, Hoffmann TJ, Lammer EJ, Hardin J, Lazaruk K, Stein JB, Gilbert DA, Wright C, Lipzen A, Pennacchio LA, et al. A genetic signature of spina bifida risk from pathway-informed comprehensive gene-variant analysis. PLoS One. 2011;6(11):e28408.

    CAS  PubMed  PubMed Central  Google Scholar 

  30. Greene ND, Stanier P, Copp AJ. Genetics of human neural tube defects. Hum Mol Genet. 2009;18(R2):R113–29.

    CAS  PubMed  PubMed Central  Google Scholar 

  31. Tran PX, Au KS, Morrison AC, Fletcher JM, Ostermaier KK, Tyerman GH, Northrup H. Association of retinoic acid receptor genes with meningomyelocele. Birth Defects Res A Clin Mol Teratol. 2011;91(1):39–43.

    CAS  PubMed  Google Scholar 

  32. Rat E, Billaut-Laden I, Allorge D, Lo-Guidice JM, Tellier M, Cauffiez C, Jonckheere N, van Seuningen I, Lhermitte M, Romano A, et al. Evidence for a functional genetic polymorphism of the human retinoic acid-metabolizing enzyme CYP26A1, an enzyme that may be involved in spina bifida. Birth Defects Res A Clin Mol Teratol. 2006;76(6):491–8.

    CAS  PubMed  Google Scholar 

  33. Li H, Zhang J, Chen S, Wang F, Zhang T, Niswander L. Genetic contribution of retinoid related genes to neural tube defects. Hum Mutat. 2018;39(4):550–562.

  34. Chen S, Zhang Q, Bai B, Ouyang S, Bao Y, Hi H, Zhang T. MARK2/Par1b insufficiency attenuates DVL gene transcription via histone Deacetylation in lumbosacral spina bifida. Mol Neurobiol. 2017;54:6304–6316.

  35. Li H, Durbin R. Fast and accurate short read alignment with burrows-wheeler transform. Bioinformatics. 2009;25(14):1754–60.

    CAS  PubMed  PubMed Central  Google Scholar 

  36. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R. Genome project data processing S: the sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.

    PubMed  PubMed Central  Google Scholar 

  37. Koboldt DC, Chen K, Wylie T, Larson DE, McLellan MD, Mardis ER, Weinstock GM, Wilson RK, Ding L. VarScan: variant detection in massively parallel sequencing of individual and pooled samples. Bioinformatics. 2009;25(17):2283–5.

    CAS  PubMed  PubMed Central  Google Scholar 

  38. Kao WC, Stevens K, Song YS. BayesCall: a model-based base-calling algorithm for high-throughput short-read sequencing. Genome Res. 2009;19(10):1884–95.

    CAS  PubMed  PubMed Central  Google Scholar 

  39. Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res. 2010;38(16):e164.

    PubMed  PubMed Central  Google Scholar 

  40. Copp AJ, Stanier P, Greene ND. Neural tube defects: recent advances, unsolved questions, and controversies. Lancet Neurol. 2013;12(8):799–810.

    PubMed  PubMed Central  Google Scholar 

  41. Kumar P, Henikoff S, Ng PC. Predicting the effects of coding non-synonymous variants on protein function using the SIFT algorithm. Nat Protoc. 2009;4(7):1073–81.

    CAS  PubMed  Google Scholar 

  42. Adzhubei IA, Schmidt S, Peshkin L, Ramensky VE, Gerasimova A, Bork P, Kondrashov AS, Sunyaev SR. A method and server for predicting damaging missense mutations. Nat Methods. 2010;7(4):248–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  43. Sah VP, Attardi LD, Mulligan GJ, Williams BO, Bronson RT, Jacks T. A subset of p53-deficient embryos exhibit exencephaly. Nat Genet. 1995;10(2):175–80.

    CAS  PubMed  Google Scholar 

  44. Yoshida H, Kong YY, Yoshida R, Elia AJ, Hakem A, Hakem R, Penninger JM, Mak TW. Apaf1 is required for mitochondrial pathways of apoptosis and brain development. Cell. 1998;94(6):739–50.

    CAS  PubMed  Google Scholar 

  45. Kuida K, Haydar TF, Kuan CY, Gu Y, Taya C, Karasuyama H, Su MS, Rakic P, Flavell RA. Reduced apoptosis and cytochrome c-mediated caspase activation in mice lacking caspase 9. Cell. 1998;94(3):325–37.

    CAS  PubMed  Google Scholar 

  46. Urase K, Kouroku Y, Fujita E, Momoi T. Region of caspase-3 activation and programmed cell death in the early development of the mouse forebrain. Brain Res Dev Brain Res. 2003;145(2):241–8.

    CAS  PubMed  Google Scholar 

  47. Ruland J, Duncan GS, Elia A, del Barco BI, Nguyen L, Plyte S, Millar DG, Bouchard D, Wakeham A, Ohashi PS, et al. Bcl10 is a positive regulator of antigen receptor-induced activation of NF-kappaB and neural tube closure. Cell. 2001;104(1):33–42.

    CAS  PubMed  Google Scholar 

  48. Duester G. Retinoic acid synthesis and signaling during early organogenesis. Cell. 2008;134(6):921–31.

    CAS  PubMed  PubMed Central  Google Scholar 

  49. Juriloff DM, Harris MJ. A consideration of the evidence that genetic defects in planar cell polarity contribute to the etiology of human neural tube defects. Birth Defects Res A Clin Mol Teratol. 2012;94(10):824–40.

    CAS  PubMed  Google Scholar 

  50. Nikolopoulou E, Galea GL, Rolo A, Greene ND, Copp AJ. Neural tube closure: cellular, molecular and biomechanical mechanisms. Development. 2017;144(4):552–66.

    CAS  PubMed  Google Scholar 

  51. Rhinn M, Dolle P. Retinoic acid signalling during development. Development. 2012;139(5):843–58.

    CAS  PubMed  Google Scholar 

  52. Huang LS, Voyiaziakis E, Markenson DF, Sokol KA, Hayek T, Breslow JL. apo B gene knockout in mice results in embryonic lethality in homozygotes and neural tube defects, male infertility, and reduced HDL cholesterol ester and apo A-I transport rates in heterozygotes. J Clin Invest. 1995;96(5):2152–61.

    CAS  PubMed  PubMed Central  Google Scholar 

  53. Abell AN, Rivera-Perez JA, Cuevas BD, Uhlik MT, Sather S, Johnson NL, Minton SK, Lauder JM, Winter-Vann AM, Nakamura K, et al. Ablation of MEKK4 kinase activity causes neurulation and skeletal patterning defects in the mouse embryo. Mol Cell Biol. 2005;25(20):8948–59.

    CAS  PubMed  PubMed Central  Google Scholar 

  54. Chi H, Sarkisian MR, Rakic P, Flavell RA. Loss of mitogen-activated protein kinase kinase kinase 4 (MEKK4) results in enhanced apoptosis and defective neural tube development. Proc Natl Acad Sci U S A. 2005;102(10):3846–51.

    CAS  PubMed  PubMed Central  Google Scholar 

  55. Walder RY, Yang B, Stokes JB, Kirby PA, Cao X, Shi P, Searby CC, Husted RF, Sheffield VC. Mice defective in Trpm6 show embryonic mortality and neural tube defects. Hum Mol Genet. 2009;18(22):4367–75.

    CAS  PubMed  PubMed Central  Google Scholar 

  56. Uehara M, Yashiro K, Takaoka K, Yamamoto M, Hamada H. Removal of maternal retinoic acid by embryonic CYP26 is required for correct nodal expression during early embryonic patterning. Genes Dev. 2009;23(14):1689–98.

    CAS  PubMed  PubMed Central  Google Scholar 

  57. Okano M, Bell DW, Haber DA, Li E. DNA methyltransferases Dnmt3a and Dnmt3b are essential for de novo methylation and mammalian development. Cell. 1999;99(3):247–57.

    CAS  PubMed  Google Scholar 

  58. Banting GS, Barak O, Ames TM, Burnham AC, Kardel MD, Cooch NS, Davidson CE, Godbout R, McDermid HE, Shiekhattar R. CECR2, a protein involved in neurulation, forms a novel chromatin remodeling complex with SNF2L. Hum Mol Genet. 2005;14(4):513–24.

    CAS  PubMed  Google Scholar 

  59. Xu W, Edmondson DG, Evrard YA, Wakamiya M, Behringer RR, Roth SY. Loss of Gcn5l2 leads to increased apoptosis and mesodermal defects during mouse development. Nat Genet. 2000;26(2):229–32.

    CAS  PubMed  Google Scholar 

  60. Takeuchi T, Yamazaki Y, Katoh-Fukui Y, Tsuchiya R, Kondo S, Motoyama J, Higashinakagawa T. Gene trap capture of a novel mouse gene, jumonji, required for neural tube formation. Genes Dev. 1995;9(10):1211–22.

    CAS  PubMed  Google Scholar 

  61. Morriss-Kay GM. Growth and development of pattern in the cranial neural epithelium of rat embryos during neurulation. J Embryol Exp Morphol. 1981;65(Suppl):225–41.

    PubMed  Google Scholar 

  62. Erickson CA, Weston JA. An SEM analysis of neural crest migration in the mouse. J Embryol Exp Morphol. 1983;74:97–118.

    CAS  PubMed  Google Scholar 

  63. Jiang X, Iseki S, Maxson RE, Sucov HM, Morriss-Kay GM. Tissue origins and interactions in the mammalian skull vault. Dev Biol. 2002;241(1):106–16.

    CAS  PubMed  Google Scholar 

  64. Narisawa A, Komatsuzaki S, Kikuchi A, Niihori T, Aoki Y, Fujiwara K, Tanemura M, Hata A, Suzuki Y, Relton CL, et al. Mutations in genes encoding the glycine cleavage system predispose to neural tube defects in mice and humans. Hum Mol Genet. 2012;21(7):1496–503.

    CAS  PubMed  Google Scholar 

  65. Farese RV Jr, Ruland SL, Flynn LM, Stokowski RP, Young SG. Knockout of the mouse apolipoprotein B gene results in embryonic lethality in homozygotes and protection against diet-induced hypercholesterolemia in heterozygotes. Proc Natl Acad Sci U S A. 1995;92(5):1774–8.

    CAS  PubMed  PubMed Central  Google Scholar 

  66. Homanics GE, Smith TJ, Zhang SH, Lee D, Young SG, Maeda N. Targeted modification of the apolipoprotein B gene results in hypobetalipoproteinemia and developmental abnormalities in mice. Proc Natl Acad Sci U S A. 1993;90(6):2389–93.

    CAS  PubMed  PubMed Central  Google Scholar 

  67. Suazo J, Pardo R, Castillo S, Martin LM, Rojas F, Santos JL, Rotter K, Solar M, Tapia E. Family-based association study between SLC2A1, HK1, and LEPR polymorphisms with myelomeningocele in Chile. Reprod Sci. 2013;20(10):1207–14.

    PubMed  PubMed Central  Google Scholar 

  68. Carter TC, Pangilinan F, Troendle JF, Molloy AM, VanderMeer J, Mitchell A, Kirke PN, Conley MR, Shane B, Scott JM, et al. Evaluation of 64 candidate single nucleotide polymorphisms as risk factors for neural tube defects in a large Irish study population. Am J Med Genet A. 2011;155A(1):14–21.

    PubMed  Google Scholar 

  69. Cecconi F, Piacentini M, Fimia GM. The involvement of cell death and survival in neural tube defects: a distinct role for apoptosis and autophagy? Cell Death Differ. 2008;15(7):1170–7.

    CAS  PubMed  Google Scholar 

  70. Wei X, Li H, Miao J, Zhou F, Liu B, Wu D, Li S, Wang L, Fan Y, Wang W, et al. Disturbed apoptosis and cell proliferation in developing neuroepithelium of lumbo-sacral neural tubes in retinoic acid-induced spina bifida aperta in rat. Int J Dev Neurosci. 2012;30(5):375–81.

    CAS  PubMed  Google Scholar 

  71. Wang L, Lin S, Yi D, Huang Y, Wang C, Jin L, Liu J, Zhang Y, Ren A. Apoptosis, expression of PAX3 and P53, and Caspase signal in fetuses with neural tube defects. Birth Defects Res. 2017;109(19):1596–604.

    CAS  PubMed  Google Scholar 

  72. Li H, Zhang J, Niswander L. Zinc deficiency causes neural tube defects through attenuation of p53 ubiquitylation. Development. 2018;145(24):dev169797.

  73. De Arcangelis A, Mark M, Kreidberg J, Sorokin L, Georges-Labouesse E. Synergistic activities of alpha3 and alpha6 integrins are required during apical ectodermal ridge formation and organogenesis in the mouse. Development. 1999;126(17):3957–68.

    PubMed  Google Scholar 

  74. Arikawa-Hirasawa E, Watanabe H, Takami H, Hassell JR, Yamada Y. Perlecan is essential for cartilage and cephalic development. Nat Genet. 1999;23(3):354–8.

    CAS  PubMed  Google Scholar 

  75. Deak KL, Boyles AL, Etchevers HC, Melvin EC, Siegel DG, Graham FL, Slifer SH, Enterline DS, George TM, Vekemans M, et al. SNPs in the neural cell adhesion molecule 1 gene (NCAM1) may be associated with human neural tube defects. Hum Genet. 2005;117(2–3):133–42.

Download references


We sincerely thank the participants for their cooperation and unwavering dedication to this study.


This work was supported by the following grants: the National Natural Science Foundation of China (Grant No. 81701441, No.81771584);National Key Basic Research Program of China (2016YFC1000502).

Author information

Authors and Affiliations



TZ, LN and HL wrote the paper. JZ and XY conceived the idea and evaluated the patients. FW and HL analyzed the genetic results. HW supervised the diagnosis and treatment. All authors approved the fnal version of the manuscript.

Corresponding authors

Correspondence to Jizhen Zou or Huili Li.

Ethics declarations

Ethics approval and consent to participate

This study was approved by the ethics committee of the Capital Institute of Pediatrics (SHERLLM 2009002) and was performed in accordance with the ethical standards laid down in the 1964 Declaration of Helsinki and its later amendments. Informed consent Written informed consent was obtained from the patients and/or their guardians prior to sample collection and open sural nerve biopsy.

Consent for publication


Competing interests

All authors declare that they have no confict of interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Table S1

. The summary for genomic DNA sequencing targeted genes.

Additional file 2: Figure S1.

Summary on PDRVs burden genes related functional pathways in other spinal phenotypes. Summary of PDRVs burden genes in cases affected with spinal NTD phenotypes in other segments. Numbers in brackets means occurrences of PDRVs in each gene, numbers or gene symbols in green represent missense PDRVs. Gene symbols linking with colon represents two concurrent PDRVs in case.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zou, J., Wang, F., Yang, X. et al. Association between rare variants in specific functional pathways and human neural tube defects multiple subphenotypes. Neural Dev 15, 8 (2020).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: