Affiliations:
Center for Human Genetic Research, Massachusetts General Hospital, Boston, Massachusetts, United States of America,
Harvard Medical School, Boston, Massachusetts, United States of America

Figures

Abstract

Background

Neural tube defects (NTDs) are one of the most common birth defects caused by a combination of genetic and environmental factors. Currently, little is known about the genetic basis of NTDs although up to 70% of human NTDs were reported to be attributed to genetic factors. Here we performed genome-wide copy number variants (CNVs) detection in a cohort of Chinese NTD patients in order to exam the potential role of CNVs in the pathogenesis of NTDs.

Methods

The genomic DNA from eighty-five NTD cases and seventy-five matched normal controls were subjected for whole genome CNVs analysis. Non-DGV (the Database of Genomic Variants) CNVs from each group were further analyzed for their associations with NTDs. Gene content in non-DGV CNVs as well as participating pathways were examined.

Results

Fifty-five and twenty-six non-DGV CNVs were detected in cases and controls respectively. Among them, forty and nineteen CNVs involve genes (genic CNV). Significantly more non-DGV CNVs and non-DGV genic CNVs were detected in NTD patients than in control (41.2% vs. 25.3%, p<0.05 and 37.6% vs. 20%, p<0.05). Non-DGV genic CNVs are associated with a 2.65-fold increased risk for NTDs (95% CI: 1.24–5.87). Interestingly, there are 41 cilia genes involved in non-DGV CNVs from NTD patients which is significantly enriched in cases compared with that in controls (24.7% vs. 9.3%, p<0.05), corresponding with a 3.19-fold increased risk for NTDs (95% CI: 1.27–8.01). Pathway analyses further suggested that two ciliogenesis pathways, tight junction and protein kinase A signaling, are top canonical pathways implicated in NTD-specific CNVs, and these two novel pathways interact with known NTD pathways.

Conclusions

Evidence from the genome-wide CNV study suggests that genic CNVs, particularly ciliogenic CNVs are associated with NTDs and two ciliogenesis pathways, tight junction and protein kinase A signaling, are potential pathways involved in NTD pathogenesis.

Funding: This work was jointly supported by the Ministry of Science and Technology of P. R. China, the National "973" Program on Population and Health (2013CB945404 to TZ; 2010CB529601 to BLW), the National Nature Science Fund (81100841 to XLC), the Beijing Excellent Scientist Fund (20061A0303200108 to XLC), Chinese Returned Oversea Scientist Fund to XLC by Beijing Science and Technology, and the Shanghai Science Fund (09JC1402400 and 09ZR1404500 to BLW). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Background

Neural tube defects (NTDs) are common and severe birth defects. They arise between the third and fourth week of embryogenesis because of partial or complete failure of neural tube closure [1]. Severe NTDs, such as craniorachischisis and anencephaly, are directly related to high morbidity and mortality. Less severe forms, such as open spina bifida, lead to life-long disabilities and impose a tremendous burden on affected families [1]. Approximately 1 in 1,000 pregnancies worldwide is affected by NTD [2]. The prevalence of NTDs is high in China, with approximately 27.4 per 10,000 pregnancies reported [3]. In particular, the Shanxi Province in China has the highest NTD occurrence at the rate of 138.7–199.4/10,000 [4], [5].

NTDs are considered to represent a complex disorder as both environmental and genetic factors contribute to NTDs [6]. It has been estimated that up to 70% of human NTDs can be attributed to genetic factors [7]. The recurrence risk of NTDs for siblings of indexed cases is 2–5%, which is approximately 50-fold higher than that of general population [8]. The robust NTD phenomena in gene knockout mouse models imply the existence of high penetrance NTD genes [9], [10]. A number of candidate genes from mouse NTDs were screened in humans, but so far only the MTHFR polymorphism and pathogenic mutations of the planar cell polarity (PCP) core gene were shown to be consistently associated with human NTDs [11], [12], [13], [14], [15], [16], [17], [18], [19]. These PCP gene mutations only account for 0.75–22.2% of different isolated NTDs, with low mutation rate in spinal and a higher rate in cranial NTDs [19], [20]. On the other hand, conventional karyotyping has revealed microscopic chromosomal abnormalities in 2.5–10.26% of fetal and newborn with NTDs [21], [22], [23], [24], suggesting that genomic imbalances are one of the genetic bases of NTDs. No specific NTD candidate gene has been implicated by cytogenetic studies due to low resolution of conventional karyotyping. With the advent of chromosomal microarray technology, we anticipate detection of smaller NTD-specific copy number variants (CNVs) that could reveal novel candidate genes and pathways associated with NTDs.

Methods

Study Design

A hospital-based case-control study was conducted from April 2006 to December 2008 in the Lvliang region of Shanxi Province. All pregnant women who enrolled this study provided their written informed consent to participate in this study. The study protocol was reviewed and approved by the Ethics Board of Capital Institute of Pediatrics.

Samples Recruitment

Collection of NTD-affected embryos took place during routine prenatal checkups in multiple local county hospitals. Detailed records of gestational ages, general development of the embryo, and embryo B-ultrasound data were kept by obstetricians. The epidemiological method was described in detail in our previous publication [25]. Medical abortions were carried out by obstetricians to terminate NTD-affected pregnancies. Pathologic diagnosis of NTDs was performed by an experienced pathologist according to the International Classification of Disease, Tenth Revision, codes Q00.0, Q05.9 and Q01.9 (http://apps.who.int/classifications). The embryos aborted for non-medical reasons from same region were also enrolled and used as matched controls. Routine prenatal checkups, questionnaire interviews and autopsies were completed for controls, and any embryos that harbored any pathological malformation or intrauterine growth retardation were excluded from control group.

Array-comparative Genomic Hybridization (Array-CGH)

Brain tissue (25 mg) from the cerebrum was used for DNA extraction (DNeasy Blood & Tissue kit, Qiagen, Valencia, CA). For the cases with anencephaly, residue brain tissue in the cavity of the skull or spinal tissue was collected for DNA extraction. Array-CGH was performed according to previously published methods [26] using the Agilent 244 K Oligonucleotide CGH microarray platform (Agilent Technologies Inc., Palo Alto, CA) in the Genetic Diagnostic Laboratory at Children’s Hospital Boston.

Workflow for the Identification of CNVs

Deletions and duplications were identified using at least five consecutive probes by the DNA Genomic Workbench Standard Edition 5.0.14 software (Agilent Technologies Inc., Palo Alto, CA). We applied a minimal size cutoff (>30 Kb) for CNV detection. A CNV was defined as non-DGV when it does not overlap, or overlaps less than 20% with reported CNVs in the Database of Genomic Variants (DGV, http://projects.tcag.ca/variation/). These CNVs were also compared with the CNVs detected in normal Chinese populations [27], [28]. CNVs that involve one or more Refseq genes or exonic region of a gene were defined as genic CNVs [29]. CNVs containing cilia genes were termed as ciliogenic CNVs.

Three primer pairs targeting the center and two sides of a CNV were designed for real-time quantitative PCR (Applied Biosystems, Carlsbad, CA). A normal male was used as a reference control and the ß-actin gene was used as an internal control. The assays were done in triplicate for each CNV.

Ingenuity Pathway Analysis

We uploaded the gene list generated from all non-DGV genic CNVs to the online Ingenuity Pathway Analysis platform (IPA, Ingenuity Systems Analysis, www.ingenuity.com) for identifying associated networks, related disorders and canonical pathways. Each identifier was mapped to its corresponding object in Ingenuity’s Knowledge Base. Fischer’s exact test was used to calculate a p-value that determined the probability of the genes involved in each pathway. A gene list containing 223 NTD candidate genes (Table S5) after literature review [9], [10], [30], [31], [32] was also analyzed using the same setting.

MKS1, MKS3(TMEM67) Mutation Screening

Meckel-Gruber Syndrome (MIM 249000) is a major cause of systemic NTD [24]. To rule out the involvement of Meckel-Gruber Syndrome, we did the mutation screening for MKS1 (MIM 609883) and MKS3/TMEM67 (MIM 609884) in 38 systemic NTD cases presenting encephalocele and/or urinary development malformation [33]. The coding exons and intronic flanking regions of MKS1/MKS3 were amplified and sequenced. Sequences were analyzed using Mutation Surveyor V3.30 (SoftGenetics, State College, PA). The sequences of primers designed for 42 amplicons were listed in the Table S6.

Statistical Analysis

Data were analyzed using SPSS version 10.0 statistics software (SPSS Inc., Chicago, IL). Normal distribution of CNV sizes was verified by the Kolmogorov-Smimov method. Fisher’s exact test or the χ2 test was used for categorical variables, Student’s t test was used for continuous variables, and Wilcoxon rank-sum test was used for rank categorical variables. Logistical regression was used to determine adjusted odds ratios and their associated 95% confidence intervals. A value of p<0.05 (two-tailed) was considered significant.

Results

Non-DGV CNVs were Enriched in NTD Samples

A total of 85 cases (NTD-affected embryos) and 75 controls were collected. No case had the positive familial NTD history, suggesting they all were sporadic NTDs. The detailed categories of NTD cases were listed on Table 1. No known pathogenic or novel, rare mutation of MKS1/MKS3 genes was detected.

A total of 1057 CNVs in 85 cases and a total of 821 CNVs in 75 controls were detected. The per sample CNV load either by count or by size was not different between cases and controls (CNV count per sample: 12.89 in cases vs. 12.07 in controls, p>0.05; CNV size per sample: 283 kb in cases vs. 330 kb in controls, p>0.05). Stratifying CNVs by size or grouping samples by CNV count did not show any significant difference among cases and controls either (Figure 1A and 1B).

Figure 1. Summary of whole genomic copy number variants (CNVs) in the two groups; Wilcoxon rank-sum test was used for statistical analysis.

Figure 1A: Distribution of all called genomic CNVs by size. Figure 1B: Distribution of all called genomic CNVs by frequency.

doi:10.1371/journal.pone.0054492.g001

After filtering out benign CNVs reported in DGV and Chinese-specific common CNVs (four in cases and five in controls) [27], [28], 81 non-DGV CNVs were identified, 55 in NTD cases and 26 in controls. We compared this non-DGV CNV list with the clinical diagnostic CNV dataset in Genetic Diagnostic laboratory, Boston Children’s Hospital and a Chinese CNV dataset consisting of about 300 mental retardation patients and their unaffected relatives. None of the non-DGV CNVs identified from NTD cases was present in these two datasets, suggesting that these non-DGV CNVs are very rare and/or NTD-specific. Out of 81 non-DGV CNVs, 59 CNVs (40 in NTD cases and 19 in controls) involve one or more genes (termed “genic CNVs”). Six different-sized genic CNVs validated by PCR are shown in Figure S1 (Table S1 for detailed primer information).

The average count of non-DGV CNVs per sample was 1.85-fold higher in NTDs than in controls (0.65 vs. 0.35) while the average count of non-DGV genic CNVs per sample was 1.88-fold higher in case than in controls (0.47 vs. 0.25). The proportions of samples with non-DGV CNVs and genic CNVs were significantly higher in cases than in controls (41.2% vs. 25.3%, p<0.05; 37.6% vs. 20%, p<0.05, see Table 2). Non-DGV genic CNVs had a lower p value than non-DGV CNVs (0.01 vs. 0.03, labeled *** in Table 2). Logistical regression analysis showed that non-DGV CNVs and genic CNVs lead to an increased risk for NTDs (non-DGV CNVs, OR = 2.26, CI = 1.10–4.74; non-DGV genic CNVs, OR = 2.65, CI = 1.24–5.87). Furthermore, for samples without non-DGV CNVs, we compared the whole genomic common CNV burden in cases and controls. Neither the count nor the size of genome-wide common CNVs was different in cases and controls (CNV count: 11.61 vs. 11.92, p>0.05; CNV size: 292 kb vs. 344 kb, p>0.05). The above analysis indicated that it is the non-DGV genic CNVs, not common CNVs that are associated with increased risk of NTDs. We also compared the non-DGV CNV size in NTDs and controls and no difference was identified (Table S2).

Enrichment of Ciliogenic CNVs in NTDs

The total number of genes involved in NTD-specific CNVs was about 3-fold that in controls (138 vs. 44, see Table S3 and S4). Importantly, nine genes (CALR3, CTNNA3, CYP27B1, LMO1, RAB19, RABL5, SLC gene family, ZIC, WNT) in NTD-specific CNVs are homologous to genes that have been reported as NTD candidate genes in mouse models (labeled yellow in Table S5), suggesting that a subset of genes in NTD-specific CNVs is directly contributing to the pathogenesis of NTDs in some samples. After referring to recently published literature reporting 828 ‘ciliome’ genes [34], forty one out of 138 genes from NTD-specific CNVs were identified as cilia genes. In controls, only twelve cilia genes were influenced by non-DGV CNVs.

IPA core analysis showed that bio-functional networks, related diseases/disorders, and canonical pathways were all different between gene lists from cases and controls (Table 3). In cases, the top-ranked network was related to “Cell Morphology, Developmental Disorder, Skeletal and Muscular Disorders”. Twenty-eight genes from cases (see Figure 2, green indicates deletion, red indicates duplication) involved in this network. The top related diseases and disorders involved were “Genetic and Neurological Disease”. The canonical pathway analysis revealed that “tight junction signaling” and “protein kinase A signaling” were two significant pathways involved in genes from NTD-specific CNVs (Figure 2, labeled as CP in the oval regions). The PAR3-PAR6-aPKC complex and CTNNB1 were shared components of the two canonical pathways.

We further identified that 12 cilia genes (GLIS3, CENPN, CDC123, WDFY2, C16orf61, SLC17A5, CTNNA3, CLDN15, CDC16, ATP5G1, PARD3 and PARD6G) from 10 cases participate this top network (see Figure 2 for their locations in the network). Both the core gene of the network-HNF4A and shared components (PARD3, PAR6D6G and CTNNB1, see Figure 2) belong to cilia genes. In order to further evaluate the involvement of cilia genes in NTDs, we compared the frequency of non-DGV CNVs containing cilia genes (termed as non-DGV ciliogenic CNVs) in NTD cases and controls. In NTDs, 21 cases carry non-DGV ciliogenic CNVs, which was significant higher than that in controls. (24.7% vs. 9.3%, p<0.05, see Table 2). Non-DGV ciliogenic CNVs led to a 3.19-fold (95% CI: 1.27–8.01) increased risk for NTDs. In order to further explore the association of non-DGV ciliogenic CNVs with subtypes of NTD (Table 4), NTD cases were stratified into cranial or spinal NTDs, isolated or systemic NTDs according to patient phenotypes. Table 4 shows that the frequency of ciliogenic CNVs was not significantly higher in cranial NTDs than in spinal NTDs (21.6% vs. 29.4%, p>0.05). Meanwhile, the frequency of ciliogenic CNVs was not significantly higher in systemic NTDs compared with isolated NTDs (28.6% vs. 13.6%, p>0.05). Interestingly, the systemic NTDs with urinary/adrenal developmental anomalies showed significantly more non-DGV ciliogenic CNVs compared with those without such anomalies (42.3%% vs 18.9%, p<0.05, one-tailed Fisher exact test). The common urinary/adrenal abnormalities included renal agenesis/hypoplasia, adrenal gland agenesis/hypoplasia and polycystic kidney. This analysis revealed a possible role of cilia genes in some specific systemic NTDs. Table 5 listed detailed phenotypes of 11 systemic NTD cases with non-DGV ciliogenic CNVs and abnormal urinary/adrenal development.

We performed an IPA for 223 known NTD candidate genes. The top network and related diseases/disorders were similar to what was identified with the NTD-specific genes. “Embryonic Development, Neural System Development and Function, Organ Development” and “Amino Acide Metabolism, Molecular Transport, Small Molecular Biochemistry” were the top networks, and “Developmental and Neurological Diseases” were the top related diseases/disorders. “One carbon pool by folate” and “Methionine Metabolism” were the canonical pathways of “Amino Acid Metabolism”. For the network “Embryonic Development, Neural System Development and Function, Organ Development”, the canonical pathways of known importance to neural tube closure, including sonic hedgehog signaling (SHH), bone morphogenic protein (BMP) signaling, Wnt/β-catenin signaling work were identified in this top network (Fig. 3, labeled as CP). Also we identified tight junction signaling and protein kinase A signaling interact with SHH and BMP signaling (Fig. 3, labeled as CP) by sheared molecules PRKAC.

IPA analysis was performed for 361 integrated genes consisted of known NTD candidate genes and 138 genes from NTD-specific CNVs. We observed that “Organ Development, Tissue Development, Embryonic Development” were in the top networks (Fig. 4), and “Developmental and Neurological Disorders” was the top related diseases/disorders. 8 genes from NTD-specific CNVs (APLF, GLIS3, PARP12, POLG, PYCR2, SUPT7L, TAF2, WNT4) involved in this top network. The SHH, BMP and Wnt/β-catenin signaling and protein kinase A signaling were identified in this network (Fig. 4, labeled as CP).

Figure 4. Top bio-functional networks in the integrated list of 361 genes.

“Organ Development, Tissue Development, Embryonic Development” was identified as top network in the integrated list of 361 genes. The interactions exist between the novel pathway and the known pathway of NTDs.

doi:10.1371/journal.pone.0054492.g004

Discussion

Association of Non-DGV Ciliogenic CNVs with NTDs

CNVs are copy number gains or losses of DNA segments ranging from one kilobase to multiple megabases in length that are important underlying factors for human genetic and phenotypic diversity [35]. Genomic CNVs contribute significantly to both rare and common disorders [36], [37], [38], [39], [40]. Previous studies showed that 6.5% of prenatally detected NTDs are associated with chromosome abnormalities using traditional cytogenetic methods [22]. These cytogenetic studies suggested that genomic imbalances could play a significant role in the pathogenesis of NTDs. The rates of chromosomal abnormalities detected by cytogenetic methods differ significantly in different type of NTDs, e.g. 2.3% in anencephaly, 7.1% in encephalocele, and 10.2% in meningomyelocele [22]. It is consistent across several studies that encephalocele have higher chromosomal abnormality rate than anencephaly, and spina bifida have higher chromosomal abnormality rate than craniorachischisis [41], [42], [43]. Among all the detected chromosomal abnormalities, trisomy was the most common chromosomal abnormality associated with NTDs [43]. Those reported chromosomal abnormalities detected by karyotyping are very large (>5 Mb) and helped little to narrow down causal genes for NTDs. Recent microarray-based genomic profiling technologies enabled us to detect smaller CNVs with great sensitivity and specificity. This study is based on the hypothesis that submicroscopic genomic CNVs are genetic risk factors for some NTDs and detection of small CNVs will help us to pinpoint novel NTD candidate genes and related developmental pathways.

To our knowledge, this is the first such study to evaluate the impact of genomic CNVs on the occurrence of NTD in Chinese patients. As a result, we detected many small CNVs in NTD samples. Not surprisingly most of these CNVs were reported in the DGV database and are presumably benign polymorphisms. A subset of non-DGV CNVs in NTD samples is novel and not shared with CNVs observed in unrelated clinical patient cohorts, suggesting that they are potentially NTD-specific, important for identifying causal genes for NTDs. Due to the lack of parental samples, we were not able to evaluate their de novo status. Instead, we compared the collective frequencies of non-DGV, particularly genic CNVs in cases and controls. Enrichment of ciliogenic CNVs was identified as a risk factor for human NTD. The significant enrichment of genic CNVs in NTD compared to controls is very similar to what was discovered in schizophrenic patients [29]. Xu et al. compared the distribution of genic CNVs in 200 schizophrenia individuals and 159 matched controls. Significant higher enrichment of rare genic CNVs was seen in familial schizophrenia than controls (relative enrichment = 2.7, p = 10−3). But unlike what was found in schizophrenia patients, we failed to identify any recurrent CNVs in NTDs that pose high penetrant effect on the phenotype.

Computerized functional analyses for genes encompassed in common/rare CNVs have been previously performed to reveal their potential roles in disease-susceptibility [29], [36], [44]. Common CNVs often encompass the “environmental sensor” genes. These genes are not necessarily critical for early embryonic development, but rather enable environmental adaptation. Such adaptive genes include olfactory receptors, immune and inflammatory response genes, cell signaling molecules and cell adhesion molecules [45]. In contrast, genes disrupted by rare CNVs in neurodevelopmental patients were significantly involved in psychological disorders and learning behaviors [36], [44]. Our study revealed that different functional pathways are involved with genic CNVs detected in case and control samples. In control samples, the top diseases were related to immunological and infectious diseases, which are similar to what were found in common CNVs. By contrast, genes related to neurological and genetic disorders were overrepresented in NTD groups, reflecting the neurodevelopmental nature of the neural tube closure process.

The signaling pathways known to be essential for neural tube closure include PCP signaling (non-canonical Wnt pathway), SHH, BMP, Wnt and Notch signaling and inositol phosphorylation metabolism [7]. Recently primary ciliogenesis processes have been proved to be important for primary neurulation via interaction with the SHH pathway [46], [47]. NTDs are frequent phenotypes in ciliary dysfunction disorders [48]. The overlap of developmental pathways and phenotypes implies a potential role for cilia-related genes in human NTDs. In this study, we identified the enrichment of ciliogenic CNVs in NTD cases. Moreover, two canonical pathways (tight junction signaling and protein kinase A signaling) identified from NTD-specific CNVs are primary cilio-related pathways [48], [49] and show interaction with known NTD pathways. Our results support the notion that some genes encompassed by NTD-specific CNVs are functionally related to ciliogenesis pathways and are potentially responsible for the NTDs.

Among cilia genes identified in NTD-specific CNVs, GLIS3 (MIM 610192) is a member of the GLI-similar zinc finger protein family and a nuclear transcription regulator. Mutations and a large deletion in the 5-UTR region of the GLIS3 gene have been reported in children with neonatal diabetes mellitus and congenial multiple malformation including polycystic kidney [50]. In addition, the association of GLIS3 CNVs with diabetes was reported [51]. Interestingly, diabetes and gestational diabetes are recognized risk factors for NTDs [52]. Case 1 in this study was an embryo with open spina bifida and bilateral adrenal gland hypoplasia, atrial septal defect. It carried a single rare duplication covering the GLIS3 5-UTR and exon 1 (see Table 5 for phenotype, see Figure S1 for microarray picture and validation). Functional network analysis demonstrated that GLIS3 participates in the network involved in NTD pathogenesis (see Figure 2 and 4 for its role in top networks, blue arrowheads). PARD3-PARD6 (MIM 606745) has been proved to be essential for epithelial polarity during embryogenesis and neurulation [53]. This polarity complex also localizes to cilia and regulates ciliogenesis, tight junction during kidney morphogenesis [54]. Zebrafish model revealed that PARD3 is critical in orchestrating the midline crossing during neurulation [55]. PARD3 depletion in zebrafish results in intracellular cilia and hydrocephalus [56]. human PARD3/PARD6 knock-out in MDCK cell resulted in severe morphogenetic defect on three-dimensional culture condition, forming multiple lumen which is similar to polycystic kidney phenotype in human [57]. The case carrying PARD3/PARD6G CNVs (case 8 and case 2 in Table 5, see Figure S1 for microarray picture and validation) presented with adrenal gland hypoplasia in addition to craniorachischisis. In case 8, PARD3 was the only gene in the sole non-DGV CNV detected. With the available maternal samples, we excluded the maternal inheritance of PARD3 deletion in case 8. CTNNA3 (MIM 607667) recruits E-cadherin and beta-catenin to cell-cell contacts, and plays role in cell-cell adhesion and tight junction. Its homolog ctnnbip1 was reported as candidate gene of mouse NTDs [9]. The mice with mutant ctnnbip1 presented with exencephaly (equivalent to anencephaly in human). Case 5 carried CTNNA3 and ACTR3B duplication, and exhibited open spina bifida, kidney agenesis and polycystic kidney phenotype. ACTR3B duplication is possible genetic factor to explain why the NTD phenotype in case 5 was different from ctnnbip1-mutated mouse. The maternal inheritance of CTNNA3 and ACTR3B duplication was excluded by its mother’s blood.

Complex Genetic Heterogeneity of NTDs in Shanxi Province

We enrolled NTDs cases from Shanxi Province, which has the highest NTD prevalence in the world. This is a region that was not supplemented with folic acid at the time of sample collection. Previous studies supported that low folate levels, and polymorphism of folate-metabolic genes were risk factors for NTDs in this region [58], [59], [60]. Genes regulating early embryonic development were also identified as risk factors of NTDs in this region [61], [62], [63]. The significant involvement of copy number variants in NTD revealed by this study further suggests the complex genetic heterogeneity underlying NTD pathogenesis in the Shanxi region. We noticed two unique features pertaining to our NTD cohort. 1) None of the cases harbored CNVs larger than 3 Mb, the majority of CNVs were between 100–500 kb in size and no trisomy was detected in this not previously screened NTD cohort. Previous studies identified large chromosomal abnormalities including trisomy in a significant percentage of NTD cases among patients from Western countries [21], [22], [23], [41], [42]; this phenomenon may indicate underlying genetic differences between our cohort collected in Shanxi region and cases collected for previous studies. 2) The NTD cases collected for this study consisted of many cases with additional malformations outside the nervous system. For example, 63 cases were categorized as systemic NTDs. Our finding of the enrichment of cilia genes in NTD related CNVs is consistent with the pleiotropic roles of cilia genes in the development of different systems. Among systemic NTD subtype, the borderline significance between systemic NTDs accompanying abnormal urinary/adrenal gland development and ciliogenic CNVs may reflect their unique genetic basis. Adrenal hypoplasia/agenesis have been reported as associated morphologic anomalies in fetus with NTDs, and 29–57% NTD cases were found to have low adrenal weights [64], [65], [66]. Although the role of cilia genes on adrenal gland development is unknown, knock out of cilia genes (PARD3/PARD6) in mammary stem/progenitor cells can result in severe morphogenetic defects, forming multiple lumens in the mammary gland [67]. Future studies with larger sample size will help to verify the association between systemic NTDs, urinal/adrenal malformation and cilia genes.

In conclusion, we detected the enrichment of genic CNVs, particular ciliogenic CNVs in NTD cases relative to controls, suggesting the association of ciliogenic CNVs with NTD. Pathway analysis further revealed that two ciliogenesis pathways including tight junction and protein kinase A signaling are involved with the pathogenesis of human NTDs.

Supporting Information

Examples of six rare genic CNVs and their validations by PCR. A–F show three heterozygous deletions and three heterozygous duplications. A–C represent the 139 kb deletion in the PARD3 region (chr10:34875595-35015198), the 30 kb deletion in the NRG3 region (chr10:84631526-84661716), and the 42 kb deletion in the SLIT2 region (chr4:20030784-20086137). D–F show three heterozygous duplications with different sizes. They are a 636 kb duplication in the PRKX region (chrx:3528099-4164677), a 107 kb duplication in the GLIS3 region (chr9:4142060-4249876), and a 58 kb duplication in the PARD6G region (chr18:76009131-76067279). G–I show the validations of the CNVs by long-range PCR and real-time quantitative PCR. G and H present bilateral breakpoints of A and B validated by long-range PCR following sequencing. Sequences were blasted (USCS hg18) online to get an exact joint site. I shows the real-time PCR validation for C–F. A normal male without these CNVs was used as control sample. CNVs were calculated using log 10 of CNV ratio (the height of case/mean height of control sample). The minus index represents deletions and the plus index represents duplications. For each CNV, three pairs of primers pairing to the head, middle and tail regions were used to amplify three segments.

Acknowledgments

We are grateful to all the families for their NTD-affected pregnancies enrolled in this study. We also extend our thanks to the obstetricians for their cooperation and support of the collection of samples at the local hospitals in Shanxi Province, and the pathologists in the Department of Pathology at the Capital Institute of Pediatrics for their diagnostic work. We appreciate Allan Zheng for his help in editing the manuscript.