Burden of rare deleterious variants in WNT signaling genes among 511 myelomeningocele patients

Abstract

Genes in the noncanonical WNT signaling pathway controlling planar cell polarity have been linked to the neural tube defect myelomeningocele. We hypothesized that some genes in the WNT signaling network have a higher mutational burden in myelomeningocele subjects than in reference subjects in gnomAD. Exome sequencing data from 511 myelomeningocele subjects was obtained in-house and data from 29,940 ethnically matched subjects was provided by version 2 of the publicly available Genome Aggregation Database. To compare mutational burden, we collapsed rare deleterious variants across each of 523 human WNT signaling genes in case and reference populations. Ten WNT signaling genes were disrupted with a higher mutational burden among Mexican American myelomeningocele subjects compared to reference subjects (Fishers exact test, P ≤ 0.05) and seven different genes were disrupted among individuals of European ancestry compared to reference subjects. Gene ontology enrichment analyses indicate that genes disrupted only in the Mexican American population play a role in planar cell polarity whereas genes identified in both populations are important for the regulation of canonical WNT signaling. In summary, evidence for WNT signaling genes that may contribute to myelomeningocele in humans is presented and discussed.

About the Authors:

Luke Hebert

Roles Conceptualization, Formal analysis, Investigation, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing

Affiliation: Department of Pediatrics, Division of Medical Genetics, McGovern Medical School at The University of Texas Health Science Center at Houston (UTHealth), Houston, TX, United States of America

Paul Hillman

Roles Conceptualization, Formal analysis, Investigation, Methodology, Software, Visualization, Writing – review & editing

Affiliation: Department of Pediatrics, Division of Medical Genetics, McGovern Medical School at The University of Texas Health Science Center at Houston (UTHealth), Houston, TX, United States of America

Craig Baker

Roles Conceptualization, Formal analysis, Investigation, Methodology, Software, Visualization, Writing – review & editing

Current address: Munroe-Meyer Institute, University of Nebraska Medical Center, Omaha, NE, United States of America

Affiliation: Department of Pediatrics, Division of Medical Genetics, McGovern Medical School at The University of Texas Health Science Center at Houston (UTHealth), Houston, TX, United States of America

Michael Brown

Roles Data curation, Formal analysis, Methodology, Software, Writing – review & editing

Affiliation: Human Genetics Center, Department of Epidemiology, Human Genetics and Environmental Sciences, School of Public Health at The University of Texas Health Science Center at Houston (UTHealth), Houston, TX, United States of America

Allison Ashley-Koch

Roles Data curation, Writing – review & editing

Affiliation: Department of Medicine and Duke Molecular Physiology Institute, Duke University Medical Center, Durham, NC, United States of America

James E. Hixson

Roles Data curation, Methodology, Validation, Writing – review & editing

Affiliation: Human Genetics Center, Department of Epidemiology, Human Genetics and Environmental Sciences, School of Public Health at The University of Texas Health Science Center at Houston (UTHealth), Houston, TX, United States of America

Alanna C. Morrison

Roles Data curation, Supervision, Writing – review & editing

Affiliation: Human Genetics Center, Department of Epidemiology, Human Genetics and Environmental Sciences, School of Public Health at The University of Texas Health Science Center at Houston (UTHealth), Houston, TX, United States of America

Hope Northrup

Roles Conceptualization, Data curation, Funding acquisition, Project administration, Resources, Supervision, Writing – review & editing

Affiliation: Department of Pediatrics, Division of Medical Genetics, McGovern Medical School at The University of Texas Health Science Center at Houston (UTHealth), Houston, TX, United States of America

Kit Sing Au

Roles Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Writing – review & editing

* E-mail: [email protected]

Affiliation: Department of Pediatrics, Division of Medical Genetics, McGovern Medical School at The University of Texas Health Science Center at Houston (UTHealth), Houston, TX, United States of America

Introduction

Myelomeningocele is the most common neural tube defect (NTD), with a prevalence of 3.2 per 10,000 births [1]. Affected people are born with both meninges and spinal cord exposed through a cleft in their vertebral column. People with myelomeningocele ordinarily survive with the appropriate medical care but frequently live with comorbidities such as Chiari malformation type II, sensory and motor issues below the opening, and more [2].

Myelomeningocele is a multifactorial disease, with evidence suggesting genetic susceptibilities play an important contributing role. Although maternal folate deficiency and gestational diabetes are both risk factors for NTDs, not all cases are explained by the environment of the fetus. One study indicates that only 27.6% of myelomeningocele cases can be attributed to known risk factors [3]. In fact, the heritability estimate of myelomeningocele in humans is 0.6 [4]. There is an increasing number of naturally occurring and lab-generated knockout mice with disruption of at least 372 genes exhibiting NTD phenotypes in mouse models [5], illustrating the role that genetic mutation can have on this family of disorders in vertebrates.

NTDs in humans occur when the neural tube fails to close between weeks three and four of gestation. At this crucial time, closure is initiated by convergent extension of neural plate and continued as the neuroepithelium bends to form a tube. Aside from controlling other morphological events across many species [6], the planar cell polarity (PCP) pathway regulates this convergent extension and bending of neural plate during embryogenesis [7]. In humans with NTDs, variants predicted to impair protein function have been found in PCP pathway genes [8]. PCP is one branch of the larger group of WNT signaling pathways [9–11]. Some genes from the WNT signaling pathways outside PCP are also implicated in the development of NTDs in humans [12, 13].

Given previous evidence that genes involved in WNT signaling contribute to NTDs in humans and model organisms, we aimed to comprehensively evaluate rare, likely deleterious, coding variants within all WNT signaling pathway genes. To do so, we leveraged a gene-based mutational burden analysis, which provides the following advantages: it does not require multiplex family data, it lends potentially more power than single-variant approaches, and it has been successfully applied to the publicly-available datasets we chose as references in another study [14]. We hypothesize that genes within the WNT signaling pathways harbor rare deleterious variants (RDVs) that are overrepresented in myelomeningocele subjects.

Materials and methods

Subject population and sample collection

Recruitment of myelomeningocele subjects with written informed consent was in accordance with an institutional review board at the University of Texas Health Science Center at Houston and is described by Au et al. 2008. The protocol, HSC-MS-00-001, for subject data collection was approved by the Committees for the Protection of Human Subjects at The University of Texas Health Science Center at Houston. Only patients with a diagnosis of isolated, non-syndromic myelomeningocele at birth were eligible for the study [15].

We evaluated exome sequence data from both myelomeningocele and publicly available reference populations to look for genomic variation. Genomic DNA samples were used for exome capture with TargetSeq reagents (Life Technologies, Inc.) based on high density oligonucleotide hybridization of GENCODE annotated coding exons, NCBI CCDS, exon flanking sequences (including intron splice sites), small non-coding RNAs and a selection of miRNA binding sites. After capture, libraries were constructed with addition of barcodes (AB Library Builder, Life Technologies, Inc.). Multiplexed sequencing used the Ion Proton platform (Life Technologies, Inc.) based on proton assays for polymerase sequencing of individual DNA molecules in wells of modified semiconductor chips.

Reference population variant data was retrieved from version 2 of the Genome Aggregation Database (gnomAD) [16]. Specifically, we used gnomAD’s 8,556 “control” Admixed American (AMR) exome data that includes Mexican Americans, Puerto Ricans, Medellin Colombians, and Peruvians as well as gnomAD’s 21,384 “control” Non-Finnish European (NFE) exome data. The word “Hispanic” will be used when referring to both Mexican American cases and the AMR references collectively. Likewise, the phrase “European ancestry” will be used when referring to both the European American cases and NFE references together.

WNT signaling gene lists

To analyze all genes within the WNT signaling pathways, the AmiGO2 web-based tool provided by the Gene Ontology (GO) Consortium was used to retrieve 523 unique Homo sapiens genes under WNT signaling GO accession number GO:0198738 as listed in the S1 Table. To evaluate which components of the overarching WNT signaling pathway were affected, we used other GO-specific gene lists including WNT protein secretion (GO:0061355), WNT-related planar cell polarity (GO:0060071), canonical WNT/β-catenin signaling (GO:0060070), and WNT-related calcium modulation (GO:0007223). These gene lists can also be found in the S1 Table [17, 18].

Analysis overview

Primary input files for the analysis of the genes include exome sequencing variant data in the form of variant call format (VCF) files for the myelomeningocele subject samples from Genome Analysis Toolkit (GATK) sequencing and VCF files for the reference population from version 2 of the Genome Aggregation Database (gnomAD) [16]. Variant calls were filtered based on quality control metrics and annotated for genomic function before variant allele burden in the myelomeningocele cases was compared to that of a reference population (Fig 1). The steps of the analysis were largely performed using custom scripts written in Python 2.7 and R 3.5.3.

[Figure omitted. See PDF.]

Fig 1. Variant filtering and annotation workflow.

The analysis pipeline included filtering variants based on sample and variant-level criteria for quality control. The pipeline also included annotating variants with gene names, reference population frequencies, and CADD phred scores so that variant counts could be collapsed over a gene’s coding region and so that the analysis could focus on likely deleterious, rare variants. Abbreviations: VCF = variant call format, GATK = Genome Analysis Toolkit, dbNSFP = Database of Non-synonymous Single-nucleotide Variant Functional Predictions version 4.0b1a, gnomAD = Genome Aggregation Database version 2.

https://doi.org/10.1371/journal.pone.0239083.g001

Annotation with dbNSFP

Annotation of those input VCF variants was performed using functional predictions from the database of non-synonymous single-nucleotide variant functional predictions (dbNSFP) version 4.0b1a [19] and exon start/stop locations were retrieved from the University of California Santa Cruz (UCSC) table browser [20].

Quality control and selection of damaging rare variants

The quality control filters for case exome data were chosen to closely match the criteria published by gnomAD. They can be categorized as being focused on variant-specific parameters and sample-specific parameters.

The Genome Analysis Toolkit (GATK) provides many variant-specific parameters. Variants retained for analysis include those with a GATK variant quality score recalibration value of “PASS”, root mean squared of mapping quality (MQ) ≥ 20, and inbreeding coefficient < -0.3. Of the single nucleotide variants (SNVs) filtered from the WES data, we evaluated only those SNVs that lay within any coding transcript or splice sites of the 523 genes of interest. In order to ensure that the allele counts compared at each locus (single nucleotide position) are representative of each population, we only compared SNVs that have ≥ 89.5% coverage of DP ≥ 10 in both the myelomeningocele cases and the gnomAD population used as references.

The remaining quality control filters target the quality of a variant within an individual sample. We kept samples that met the following at each locus: a genotype with alternate allele depth ≥ 25%, a read depth ≥ 10, and genotype quality score ≥ 20. A variant site is considered “covered” if its position had a DP ≥ 10 in 89.5% of samples.

In addition to evaluating only high quality variants and samples in the case population, the analysis was further focused to include only variants that were rare (defined as having an alternate allele frequency in the gnomAD reference population less than 1%) and predicted to be deleterious (coding for a stop gain, stop loss, splice site missense, or having a combined annotation dependent depletion phred score greater than 20) [21]. This applies to both myelomeningocele cases and gnomAD references. We refer to rare, deleterious variants that met the above quality control standards as “qualifying variants.”

Mutational burden analysis

To find any genes with higher mutational burden in a case population compared to its gnomAD reference, a two-by-two table for each gene was constructed for Fisher’s exact test. The Fisher’s exact test compares the number of affected individuals and unaffected individuals from case and the reference populations, where an “affected individual” is someone containing at least one RDV in the gene being evaluated and an “unaffected individual” is someone containing zero RDVs in that gene.

Individual subject data was not readily available from gnomAD. Therefore, affected and unaffected reference population numbers for each gene were estimated using the Hardy-Weinberg equations which utilize the number of qualifying RDVs and the total allele number within a given gene across the gnomAD population. We tested this estimation approach by first applying it to our case population, generating Q-Q plots that all included a λ value of 2. Any λ value greater than 1 indicates overestimation of affected individuals. For the actual burden analysis, we only applied this estimation to our reference data and not our case data. So, this test indicates that our estimation approach errs on the side of overestimating affected individuals in the reference population. This one-sided overestimation of affected individuals makes our analysis less likely to suggest falsely high mutational burden for genes in the case population.

A Bonferroni correction for multiple comparisons was applied in order to find any genes with significant mutational burden in the case populations. For the Bonferroni correction, the conventional alpha value 0.05 was divided by the number of genes that were compared in a population’s mutational burden analysis. The resulting Bonferroni value is a strict P value cutoff for statistical “significance”. It is important to note that genes were only compared if they harbored qualifying variants, so not all 523 original genes were ultimately compared. Regardless of the Bonferroni correction, the term “nominally significant” refers to genes with Fisher’s exact test P values ≤ 0.05.

Gene ontology enrichment analysis

An enrichment analysis of Gene Ontology (GO) terms was conducted using ToppCluster [22] to compare patterns of known biological function between the genes disrupted in Mexican-American (MA) subjects versus genes disrupted in European ancestry (EA) subjects.

Results

Of the 523 WNT signaling genes, 173 contained qualifying RDVs (GATK PASS, AF < 1%, CADD phred ≥ 20, coverage ≥ 90%) in the MA myelomeningocele population (S2 Table). When comparing RDV mutational burden in the 173 genes between the two Hispanic populations, ten genes associated with a risk for myelomeningocele by yielding Fisher’s exact test P values below 0.05 and odds ratios above 1 (Table 1). This included PORCN, CDH2, PRICKLE2, CPE, FUZ, PTPRU, PSMD3, TNRC6B, PPP2R1A, and FERMT2. Of these, PORCN, CPE, and TNRC6B were detected caudally in the closing neural tube in humans [23]. PORCN, CDH2, and FUZ have been previously associated with NTD phenotypes in mouse models [24–27]. This analysis gave a Bonferroni correction value of 2.9 x 10−4, which none of the genes’ P values fell below.

[Figure omitted. See PDF.]

Table 1. Mutation burden analysis of WNT signaling genes in MA and EA populations.

https://doi.org/10.1371/journal.pone.0239083.t001

In the EA myelomeningocele population, 189 genes contained qualifying variants (S3 Table). After comparing those 189 genes between the case and reference populations of European ancestry, seven genes were associated with a risk for myelomeningocele by giving P values below 0.05 and an odds ratio above 1 (Table 1). These genes included DDB1, SDC1, CSNK1G2, SOSTDC1, PLCB1, DVL2, and TLE3. Of these, DDB1, PLCB1, and DVL2 were detected in the human neural tube during closure and DVL2 is also associated with NTDs in mouse models [23, 28]. This analysis gave a Bonferroni correction value of 2.6 x 10−4, which none of the genes’ P values fell below.

All patient samples that harbored qualifying RDVs in the seventeen associated genes were heterozygous for those variants, except one MA patient who was hemizygous for his variant in PORCN. Subjects from both ethnicities tended to have a similar number of genes that contained qualifying variants, but MA subjects tended to have more nominally significant disrupted genes (Fig 2). Nominally significant disrupted genes were identified in 46 (18%) EA and 66 (26%) MA myelomeningocele subjects. In addition, four combinations of nominally significant disrupted genes, each two genes in length, were found in individual samples from the MA myelomeningocele population. These combinations included CPE and PORCN, PTPRU and PSMD3, PTPRU and FERMT2, FERMT2 and PSMD3, PPP2R1A and TNRC6B (S6 Table). No such combinations were found in individual samples from the EA myelomeningocele population.

[Figure omitted. See PDF.]

Fig 2. Distribution of affected genes per subject.

This figure compares Mexican-American (MA) myelomeningocele subjects to European ancestry (EA) myelomeningocele subjects by illustrating how many disrupted genes each subject possessed. Green circles represent MA myelomeningocele data and blue squares represent EA myelomeningocele data. (A) The horizontal axis measures the number of genes that are disrupted by qualifying variants. The vertical axis depicts the number of subjects who possess that number of disrupted genes. (B) The horizontal axis counts only the number of genes that indicated nominally significant risk in the mutational burden analysis. The vertical axis shows the number of subjects with that number of disrupted genes. The vertical dotted lines from either graph represent the median number of disrupted genes possessed by subjects, with green lines representing MA and blue representing EA medians. Both populations followed a similar distribution for how many genes are disrupted per subject, complete with a similar median. However, the Mexican Americans tended to have more nominally significant disrupted genes.

https://doi.org/10.1371/journal.pone.0239083.g002

The disrupted genes that associated with EA myelomeningocele subjects are different from disrupted genes associated with MA myelomeningocele subjects, i.e. the nominally significant genes from each comparison did not overlap.

Gene Ontology analysis regarding biological processes of the disrupted genes in both populations reveals that both groups have disrupted genes enriched for regulation of WNT signaling and negative regulation of canonical WNT signaling. However, disrupted genes in only the EA subjects were enriched for genes associated inositol phosphate metabolism, disassembly of the destruction complex, and recruitment of AXIN to the membrane (Table 2).

[Figure omitted. See PDF.]

Table 2. Gene Ontology (GO) enrichment analysis.

https://doi.org/10.1371/journal.pone.0239083.t002

Discussion

Exomes of the myelomeningocele cases compared to gnomAD references revealed nominally significant disruption among genes that function in WNT trafficking, PCP, WNT/β-catenin signaling, and WNT/Ca2+ signaling. Our reference value estimation method overestimates the number of RDVs in the gnomAD reference population, so protective variants are likely also overestimated. Therefore, we focus on genes with an odds ratio greater than one, which indicates a risk for myelomeningocele. To remain conservative in our conclusions, we do not evaluate genes whose mutational burden indicates a protective effect, marked by an odds ratio less than one. All nominally significant genes with an odds ratio greater than one are included in Table 1, but our discussion further focuses on ten genes that have either been associated with NTD phenotypes, were expressed during closure of the caudal human neural tube, or harbored individual variants that were much more common in the cases than the reference populations. Nominally significant genes meeting one or more of these criteria are: PORCN, CPE, FUZ, TNRC6B, DVL2, DDB1, PRICKLE2, CDH2, PPP2R1A, and SOSTDC1. For clarity, we organize each of these ten genes based on their known role in WNT trafficking (Fig 3), the PCP pathway (Fig 4), the canonical β-catenin pathway (Fig 5) or noncanonical Ca2+ WNT signaling pathway (S1 Fig). These figures provide visual references for these pathways with special attention drawn to the translational products of the ten disrupted genes mentioned above [29].

[Figure omitted. See PDF.]

Fig 3. Altered WNT trafficking.

The proteins PORCN, SDC1, and CPE are all involved in WNT ligand trafficking. Deleterious variants in PORCN may prevent WNT’s acetylation which is necessary for WNT to leave the endoplasmic reticulum. Loss of PORCN function would stop WNT from leaving the cell. CPE prevents WNT from binding cell surface receptors on the target cell. Loss of CPE function may indirectly increase WNT’s effect on the target cell. Proteoglycans like SDC1 have been implicated as regulators in WNT distribution, though the exact mechanism is not yet known.

https://doi.org/10.1371/journal.pone.0239083.g003

[Figure omitted. See PDF.]

Fig 4. Altered planar cell polarity.

DVL2, PRICKLE2, and FUZ genes were each disrupted in the myelomeningocele populations. The FZD-DVL complex is necessary to establish planar cell polarity. DVL proteins are also needed to translocate FUZ to cilia. FUZ is essential for the development of a cilium. PRICKLE inhibits the formation of the FZD-DVL complex.

https://doi.org/10.1371/journal.pone.0239083.g004

[Figure omitted. See PDF.]

Fig 5. Altered β-catenin.

A visual summary of the WNT/β-catenin cascade including genes disrupted in the myelomeningocele populations. The level of β-catenin, which is coded for by the human homolog CTNNB1, is regulated by many proteins in the WNT ligand’s target cell. The proteins CSNK1G, DVL2, DDB1, PSMD3, PTPRU, Fermitin 2 (coded by FERMT2), CPE, PPP2R1A (not shown), and SOSTDC1 contribute to β-catenin’s regulation and all have higher mutational burdens in one of the myelomeningocele populations compared to gnomAD references.

https://doi.org/10.1371/journal.pone.0239083.g005

Myelomeningocele exomes with WNT trafficking genes disrupted

Disrupting the process of WNT synthesis or secretion has the potential to directly and indirectly inhibit one or more of the three downstream WNT signaling pathways. WNT ligands are expressed, processed, and secreted into extracellular space before finding their target cell membrane, where they can bind several cell surface receptors [30] (Fig 3).

In the Hispanic mutational burden comparison, the X chromosome gene PORCN associated with risk for myelomeningocele. PORCN is necessary for the post-translational modification of the WNT3A ligand. PORCN belongs to an evolutionarily conserved gene family termed “Porcupine,” whose members code for Wnt processing proteins across species [31, 32]. In humans, the PORCN protein catalyzes O-palmitoylation of WNT3A’s Ser209 residue, allowing WNT3A to leave the endoplasmic reticulum [33, 34] with the help of Wntless/WLS [35]. Therefore, the likely deleterious variants found in PORCN in our MA myelomeningocele population may prevent WNT3A from leaving the endoplasmic reticulum, ultimately downregulating all WNT signaling pathways in the target cell because less WNT3A would be secreted. Previously, PORCN’s potential role in the development of NTDs has been suggested in a mouse study where heterozygous constitutive inactivation of its homolog Porcn caused open neural tubes in utero [24]. Myelomeningocele has been documented in a person with focal dermal hypoplasia (FDH), a rare congenital disorder associated with mutations in PORCN [36]. Furthermore, PORCN is differentially expressed in the caudal human neural tube during closure [23]. Three myelomeningocele subjects (two females and one male) carried the three RDVs NC_000023.10:g.48371088G>A p.(Gly223Ser), NC_000023.10:g.48371104G>A p.(Arg228His), and NC_000023.10:g.48374165C>G p.(Ala336Gly), all in highly conserved functional regions of PORCN [37]. The male who carries NC_000023.10:g.48371088G>A may constitute a homozygote loss of PORCN function. Although a p.(Arg228Cys) has been described as benign in FDH subjects, the NC_000023.10:g.48371104G>A p.(Arg228His) variant we report is different because it is predicted to be damaging by multiple functional analysis algorithms [38]. The NC_000023.10:g.48374165C>G p.(Ala336Gly) variant is located within one of the transmembrane domains of PORCN that forms the substrate transportation pore to provide substrate for acylation of WNT in the ER lumen. The NC_000023.10:g.48374165C>G p.(Ala336Gly) variant is located in between several FDH variants that cause amino acid substitutions such as p.(Leu331Arg), p.(Ser337Arg) and p.(Ala338Pro) [39].

The carboxypeptidase E (CPE) gene on chromosome 4, which may affect binding of WNT3A to the target cell’s receptor, associated with myelomeningocele in the MA population. The peptide encoded by full-length CPE interacts with WNT3A and its receptor Frizzled (FZD) to decrease WNT/β-catenin signaling in a proteasome-dependent manner in human cells [40]. Disruption of this protein’s function, therefore, may upregulate the WNT/β-catenin pathway.

Myelomeningocele exomes with PCP genes disrupted

Closure of the developing neural tube occurs via a process of convergent extension that is orchestrated by the PCP pathway [7]. Any disruption of PCP has the potential to prevent proper neural tube closure. When WNT ligand binds only the cell surface receptor Frizzled (FZD), then DVL is recruited and the PCP branch of WNT signaling is initiated (Fig 4). Our results reproduce and expand the findings of several current human NTD studies that show damaging variants in PCP genes play a role in the development of NTDs [8].

DVL2, a human homolog of the Dishevelled gene family, associated with myelomeningocele in the EA population. Dvl2 is a Dvl family protein essential in the PCP pathway [6]. Specifically, Dvl2 is required for endocytosis of the activated Wnt receptor, Frizzled [41]. Frog knockouts for the DVL2 homolog Xdsh lack convergent extension and consequently display open neural tubes [42]. Similarly, mouse knockouts for Dvl2 display a spina bifida phenotype [28]. Also, DVL2 is differentially expressed in the caudal neural tube during neural tube closure [23]. Therefore, loss of correct DVL2 function in some EA population subjects may cause myelomeningocele by preventing convergent extension during neural tube closure through disruption of PCP signaling.

PRICKLE2 on chromosome 3 associated with myelomeningocele in the MA population. PRICKLE2 is one of the two vertebrate homologs of the fruit fly’s Prickle [43]. As summarized in Y. Yang & Mlodzik, 2015, Prickle and Vangl help establish PCP by directly antagonizing the formation of the Frizzled/FZD-Disheveled/DVL complex [44]. Because this polarity is required for convergent extension, loss of PRICKLE2 may prevent proper convergent extension and thus failure of the neural tube to close. Another human myelomeningocele association study revealed more potential single-nucleotide polymorphism associations in PRICKLE2 than any other gene among three of the four ethnicities evaluated; however, similar to our own findings, the study failed to establish any strong associations to myelomeningocele in PCP genes [45]. Similarly, a targeted sequencing study of ninety human cranial NTD cases reported a discovery of likely damaging rare variants in PRICKLE2 as well [46].

CDH2, which codes for a cadherin cell adhesion protein, associated with myelomeningocele in the MA population. Other cadherin genes such as CELSR1 are expressed diversely in the developing neural tube [47]. Also, mutations in cell adhesion protein genes such as Celsr1, EphrinA5, and EphA7 cause NTD phenotypes in mice [48, 49]. Moreover, CDH2 itself is already implicated in the cause of NTDs because homozygous mouse knockouts for the homolog Cdh2 displayed a wavy neural tube phenotype [25].

The gene FUZ on chromosome 19, whose mouse homolog Fuz is a PCP effector protein required for ciliogenesis by transporting DVL to the cilium [26, 27, 50], is also associated with myelomeningocele in the MA population. Five of the eight MA subjects with qualifying variants in FUZ possessed the variant NC_000019.9:g.50315872C>T p.(Ser78Asn), a variant which was 21 times more frequent in the MA cases than the matched gnomAD references. Murdoch and Copp reviewed the complex relationship between cilia and NTDs [51]. Similar to the cilia proteins associated with exencephaly, perhaps FUZ’s transport of DVL to the cilium can also influence neural tube closure in myelomeningocele subjects. Indeed, mice with homozygous loss of Fuz expression display NTD phenotypes such as exencephaly before dying in utero [26, 27]. Also, multiple human myelomeningocele subjects possessed nonsynonymous mutations in FUZ that were not found in reference subjects, and these human variants revealed impaired cilia formation when tested in mouse cell lines [52].

Myelomeningocele exomes with canonical WNT/β-catenin pathway genes disrupted

Out of the seventeen genes that nominally associated with risk for myelomeningocele, eleven were from the canonical WNT pathway. A high proportion in the canonical pathway is surprising, given that most WNT signaling genes previously associated with myelomeningocele are from the noncanonical PCP pathway, whose role in neural tube closure is more established. When a WNT family protein binds FZD and low-density lipoprotein receptor-related protein 5 or 6 (LRP5/6), the canonical WNT signaling pathway is activated (Fig 5). The pathway involves multiple events that influence the level of β-catenin in the target cell. If present in sufficient amount, β-catenin translocates to the nucleus where it displaces transducin-like enhancer/TLE family proteins from binding T-cell factor/TCF proteins and proceeds to upregulate the expression of downstream gene targets [53].

As mentioned in the section on PCP, DVL2 associates with myelomeningocele in the EA population. In addition to their role in PCP, Dishevelled (Dvl/Dsh) proteins are important for canonical Wnt signaling, acting downstream of Wnt, Fzd, and Lrp5/6 [54, 55]. Dvl2 is already implicated in the cause of myelomeningocele, because two to three percent of Dvl2 double knockout mice display spina bifida phenotype [28]. In humans, DVL2 is also differentially expressed in the caudal human neural tube during closure [23].

DDB1, which codes for UV damage-specific DNA-binding protein 1, associated with subjects in the EA population as well. DDB1 is used as an adaptor protein for the ubiquitin ligase protein, CUL4 [56]. Other E3 ubiquitin ligases from this family (CUL1 and CUL3) target β-catenin and Disheveled for degradation [57]. Although more research is needed, it is possible that a nonfunctional DDB1 gene might decrease β-catenin degradation, consequently increasing the β-catenin signal in WNT-targeted cells. DDB1 is particularly interesting, because the Fisher’s exact P value for its mutational burden analysis almost passed the strict Bonferroni correction threshold adjusted for the 189 genes compared in the populations of European ancestry (7.88E-4 versus 2.65E-4) and DDB1 is differentially expressed in the caudal human neural tube during closure [23]. In mice, Ddb1 null embryos were extensively degenerated and died at early embryonic stage. Conditional inactivation of Ddb1 in the brain causes neural progenitor cells to apoptose, leading to neuronal degeneration, brain hemorrhages, and neonatal death [58]. However, DDB1 has not been previously associated with myelomeningocele.

As mentioned in the section on WNT trafficking, CPE associated with myelomeningocele in the MA population with one qualifying variant found in half of the relevant subjects. Whereas full length CPE interacts with the Wnt ligand, CPE’s splice variant CPE-ΔN localizes to the nucleus, increases β-catenin expression, and induces expression of Wnt target genes [40]. Plausibly, loss of CPE-ΔN might indirectly lower gene expression driven by β-catenin. Though CPE has not been previously associated with myelomeningocele, CPE is differentially expressed in the caudal neural tube during human neural tube closure [23].

PPP2R1A on chromosome 19 also associated with myelomeningocele in the MA population. Two of the three variants that created the association were at the same location in two of three MA subjects carrying the qualifying variants. The variant NC_000019.9:g.52725413G>T p.(Arg527Leu) was over 34 times more frequent in the MA population than in the gnomAD AMR reference population and the variant NC_000019.9:g.52725413G>A (p.(Arg527His) was not reported in the gnomAD AMR references. PPP2R1A codes for a conserved subunit of the heterotrimeric protein PP2A, a serine/threonine protein phosphatase [59, 60]. Another subunit, B56, may direct PP2A to down-regulate the expression of Wnt/β-catenin effector genes by decreasing the amount of β-catenin in the cell [61]. It is possible this occurs by directly complexing with Axin [62] and dephosphorylating part of the APC complex [61]. Given these understandings, loss of PPP2R1A function could lead to increased expression of β-catenin effector genes. PPP2R1A is another novel gene association with myelomeningocele.

TNRC6B codes for an Argonaut-associated RNA and shows an association with myelomeningocele in the MA population. The Argonaut molecule is a recognition motif-containing protein that participate in RNA interference [63]. The TNRC6B RNA serves as one component of a RISC complex that inhibits NLK translation [64] and the NLK protein participates in the WNT/β-catenin pathway by causing the dissociation of the β-catenin complex from DNA [65]. Importantly, TNRC6B is differentially expressed in closing caudal human neural tube [23]. So, it is possible that a nonfunctional TNRC6B transcript compromises the RISC complex and indirectly decreases β-catenin role in downstream gene expression in myelomeningocele subjects. Our study is the first to association TNRC6B and myelomeningocele.

SOSTDC1 on chromosome 7, which codes for sclerostin domain-containing protein 1, associated with myelomeningocele in the EA population. The variant NC_000007.13:g.16505280G>A p.(Ala5Val) was found in two of three EA subjects who possessed qualifying variants and both subjects were heterozygous for this variant. The NC_000007.13:g.16505280G>A p.(Ala5Val) variant was not present in any of the gnomAD NFE reference subjects. In mouse tooth development, the homolog Sostdc1 serves as an inhibitor of Lrp5/6-mediated Wnt signaling [66], but SOSTDC1 has not previously been associated with myelomeningocele.

Gene ontology enrichment analysis

None of the nominally significant genes among MA myelomeningocele subjects overlapped with the genes from the EA subjects (Table 1). However, gene ontologies from each gene set did overlap between ethnicities. GO enrichment analysis revealed that the three GO terms titled “cell-cell signaling by wnt”, “Wnt signaling pathway”, and “Signaling by Wnt” were enriched in disrupted genes from both populations. Enrichment of these broad terms were expected because the original 523 genes were retrieved using the overarching “cell-cell signaling by wnt” term. In other words, any subset of the original gene list is likely to be enriched for terms that describe the general WNT signaling pathway.

More interesting are the lower level, more specific GO terms that were enriched in both populations’ signaling pathways. These shared ontologies suggest a shared mechanism behind myelomeningocele in these two ethnicities, despite having high mutational burdens in non-overlapping genes. The clearest example is negative regulation of the canonical WNT/β-catenin component of WNT signaling.

That said, many GO terms were only enriched in one or the other population’s disrupted gene lists. So, while some mechanisms may be shared, others may be unique to each ethnicity. This is one possible explanation behind the two populations’ non-overlapping sets of disrupted genes. The ontology terms enriched in only the MA gene list largely pertain to PCP, whereas the story is less clear for the EA ontology terms.

Another difference between populations can be seen in Fig 2, where MA subjects tended to have more disrupted genes than EA subjects. More research is needed, but the discrepancy could help explain why the prevalence of myelomeningocele among Mexicans is higher than in non-Hispanic whites [67].

Limitations

The current study introduced 17 genes associated with risk for myelomeningocele with nominal significance (Fisher exact P values ≤ 0.05). Bonferroni correction for multiple comparisons was applied to each analysis using the number of compared genes as a denominator (173 in Hispanic comparisons and 189 in European ancestry comparisons) and 0.05 as the numerator alpha to calculate the correction threshold. While some genes came close, none of the P values fell below the correction threshold. A large enough group of myelomeningocele subjects would lend enough statistical power to achieve the strict statistical significance of a Bonferroni correction, but exome data for 511 myelomeningocele subjects is considerable, given the resource-intensive nature of gathering samples. The genes discussed above are suggested candidate genes.

Two of our quality control filters assume that each variant exists in Hardy-Weinberg equilibrium within the myelomeningocele population (inbreeding coefficient < -0.3, alternate allele depth ≥ 25%). However, if a selective pressure acts on a variant locus, that Hardy-Weinberg assumption is not met. Therefore, our analysis may exclude important variant loci that are under intense selective pressure.

A new approach was taken by limiting our analysis to a subset of human genes rather than evaluating the entire exome at once. The high concentration of current NTD candidate genes within the WNT signaling pathway and the relevance of WNT signaling to neural tube closure mechanisms such as convergent extension prompted this focused approach. A more inclusive approach would also be valuable.

Perspective

Authors of a recent review list genes containing deleterious rare variants that have been associated with NTDs in humans. This study proffers ANKRD6, CELSR1, CELSR2, CELSR3, DVL2, DVL3, FZD, PK1, VANGL1, VANGL2, LRP6, PTK7, and SCRIB1 as associated genes [68]. Except for DVL2, our study may appear not to corroborate findings in their review. We offer three explanations for this discrepancy. First, the current study uniquely employs a gene-based mutational burden analysis. Adopting a similar methodology to those studies summarized in the review may yield similar results. Second, our subjects differ in ethnicity from many of the populations discussed in other studies. As the current study suggests, ancestry may influence which genes associate with myelomeningocele. Third, any variants that were not equally covered in both case and reference populations were filtered out before analysis. If, for example, an important variant was discovered in the myelomeningocele subjects while the variant’s location was not covered in the corresponding gnomAD reference population, that variant would not contribute to the mutational burden analysis.

In summary, we report seventeen genes within the known WNT signaling pathways which may play a role in the development of myelomeningocele. As discussed above, the genes PORCN, DVL2, CDH2, FUZ are already suspected to play a role in NTD development from studies in animal models and in some cases humans, however the remaining thirteen genes reported here are new in their possible association with myelomeningocele.

Supporting information

[Figure omitted. See PDF.]

S1 Fig. Calcium WNT pathway.

A depiction of PLCB1’s role in the noncanonical Ca2+ WNT signaling pathway. PLCB1 activates the G-protein coupled receptor’s alpha subunit, which is necessary for the downstream activation of CAMK2A.

https://doi.org/10.1371/journal.pone.0239083.s001

(TIF)

S1 Table. Lists of GO-retrieved WNT signaling genes.

https://doi.org/10.1371/journal.pone.0239083.s002

(XLS)

S2 Table. Hispanic gene-based mutational burden analysis results.

https://doi.org/10.1371/journal.pone.0239083.s003

(XLS)

S3 Table. European Ancestry gene-based mutational burden analysis results.

https://doi.org/10.1371/journal.pone.0239083.s004

(XLS)

S4 Table. Hispanic mutational burden analysis variants.

Variant data for those variants used in the Hispanic mutational burden analysis which fall within nominally significant disrupted genes.

https://doi.org/10.1371/journal.pone.0239083.s005

(XLS)

S5 Table. European ancestry mutational burden analysis variants.

Variant data for those variants used in the European ancestry mutational burden analysis which fall within nominally significant disrupted genes.

https://doi.org/10.1371/journal.pone.0239083.s006

(XLS)

S6 Table. RDV gene combinations.

A list of nominally significant disrupted gene combinations that contained RDVs in the same sample.

https://doi.org/10.1371/journal.pone.0239083.s007

(XLSX)

Acknowledgments

We give special thanks to Lawrence C. Shimmin and Do-Kyun Kim from the Department of Epidemiology, Human Genetics and Environmental Sciences, School of Public Health at The University of Texas Health Science Center at Houston (UTHealth) for their work and expertise in WES. Special thanks to Karen Soldano and Melanie Garett in Department of Medicine at Duke University for their technical support and providing human neural tubes RNA SAGE information to annotate WES results.

Citation: Hebert L, Hillman P, Baker C, Brown M, Ashley-Koch A, Hixson JE, et al. (2020) Burden of rare deleterious variants in WNT signaling genes among 511 myelomeningocele patients. PLoS ONE 15(9): e0239083. https://doi.org/10.1371/journal.pone.0239083

References

1. Zaganjor I, Sekkarie A, Tsang BL, Williams J, Razzaghi H, Mulinare J, et al. Describing the prevalence of neural tube defects worldwide: a systematic literature review. PLoS One. 2016;11(4):1–31.

2. Copp AJ, Adzick NS, Chitty LS, Fletcher JM, Holmbeck GN, Shaw GM. Spina bifida. Nat Publ Gr [Internet]. 2015;1:1–18. Available from: http://dx.doi.org/10.1038/nrdp.2015.7

3. Agopian AJ, Tinker SC, Lupo PJ, Canfield MA, Mitchel LE. Proportion of Neural Tube Defects Attributable to Known Risk Factors. 2013;97(1):42–6.

4. Woolf CM. A genetic study of spina bifida cystica in Utah. Biodemography Soc Biol. 1975;22(3):216–20.

5. Salbaum JM, Kappen C. Neural tube defect genes and maternal diabetes during pregnancy. Birth Defects Res Part A—Clin Mol Teratol. 2010;88(8):601–11. pmid:20564432

6. Henderson DJ, Long DA, Dean CH. Planar cell polarity in organ formation. Curr Opin Cell Biol [Internet]. 2018;55:96–103. Available from: pmid:30015152

7. Nikolopoulou E, Galea GL, Rolo A, Greene NDE, Copp AJ. Neural tube closure: cellular, molecular and biomechanical mechanisms. 2017;2:552–66.

8. 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 Part A—Clin Mol Teratol. 2012;94:824–40. pmid:23024041

9. Chu CW, Sokol SY. Wnt proteins can direct planar cell polarity in vertebrate ectoderm. Elife. 2016;5:1–13.

10. Yang W, Garrett L, Feng D, Elliott G, Liu X, Wang N, et al. Wnt-induced Vangl2 phosphorylation is dose-dependently required for planar cell polarity in mammalian development. Cell Res [Internet]. 2017;27:1–19. Available from: pmid:28057935

11. Wu J, Mlodzik M. Wnt/PCP Instructions for Cilia in Left-Right Asymmetry. Dev Cell [Internet]. 2017;40:423–4. Available from: pmid:28292419

12. Allache R, Wang M, De Marco P, Merello E, Capra V, Kibar Z. Genetic studies of ANKRD6 as a molecular switch between Wnt signaling pathways in human neural tube defects. Birth Defects Res Part A—Clin Mol Teratol. 2015.

13. Lei Y, Fathe K, McCartney D, Zhu H, Yang W, Ross ME, et al. Rare LRP6 variants identified in spina bifida patients Yunping. Hum Mutat. 2015;36(3):342–9. pmid:25546815

14. Guo MH, Plummer L, Chan YM, Hirschhorn JN, Lippincott MF. Burden Testing of Rare Variants Identified through Exome Sequencing via Publicly Available Control Data. Am J Hum Genet [Internet]. 2018;103:1–13. Available from: https://doi.org/10.1016/j.ajhg.2018.08.016

15. Kit SA, Tran PX, Tsai CC, O’Byrne MR, Lin JI, Morrison AC, et al. Characteristics of a spina bifida population including North American Caucasian and Hispanic individuals. Birth Defects Res Part A—Clin Mol Teratol. 2008;82(10):692–700. pmid:18937358

16. Karczewski KJ, Francioli LC, Tiao G, Cummings BB, Alföldi J, Wang Q, et al. Variation across 141,456 human exomes and genomes reveals the spectrum of loss-of-function intolerance across human protein-coding genes. bioRxiv [Internet]. 2019; Available from: https://www.biorxiv.org/content/10.1101/531210v2

17. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene Ontology: tool for the unification of biology. Nat Genet. 2000;25(1):25–9. pmid:10802651

18. Consortium TGO. The Gene Ontology Resource: 20 years and still GOing strong. Nucleic Acids Res [Internet]. 2019;47:D330–8. Available from: https://academic.oup.com/nar/advance-article/doi/10.1093/nar/gky1055/5160994 pmid:30395331

19. Liu X, Wu C, Li C, Boerwinkle E. dbNSFP v3.0: A One-Stop Database of Functional Predictions and Annotations for Human Nonsynonymous and Splice-Site SNVs. Hum Mutat. 2016;37(3):235–41. pmid:26555599

20. Karolchik D, Hinrihs AS, Furey TS, Roskin KM, Sugnet CW, Haussler D, et al. The UCSC Table Browser data retrieval tool. Nucleic Acids Res [Internet]. 2004;32:D493–6. Available from: https://academic.oup.com/nar/article-lookup/doi/10.1093/nar/gkh103 pmid:14681465

21. Rentzsch P, Witten D, Cooper GM, Shendure J, Kircher M. CADD: Predicting the deleteriousness of variants throughout the human genome. Nucleic Acids Res. 2019;47(D1):D886–94. pmid:30371827

22. Kaimal V, Bardes EE, Tabar SC, Jegga AG, Aronow BJ. ToppCluster: A multiple gene list feature analyzer for comparative enrichment clustering and networkbased dissection of biological systems. Nucleic Acids Res. 2010;38(Web Server issue):W96–102. pmid:20484371

23. Krupp DR, Xu PT, Thomas S, Dellinger A, Etchevers HC, Vekemans M, et al. Transcriptome profiling of genes involved in neural tube closure during human embryonic development using long serial analysis of gene expression (long-SAGE). Birth Defects Res Part A—Clin Mol Teratol. 2012;94(9):683–92. pmid:22806986

24. Liu W, Shaver TM, Balasa A, Ljungberg MC, Wang X, Wen S, et al. Deletion of Porcn in mice leads to multiple developmental defects and models human focal dermal hypoplasia (Goltz syndrome). PLoS One. 2012;7(3).

25. Radice GL, Rayburn H, Matsunami H, Knudsen KA, Takeichi M, Hynes RO. Developmental defects in mouse embryos lacking N-cadherin. Dev Biol. 1997;181(1):64–78. pmid:9015265

26. Gray RS, Abitua PB, Wlodarczyk BJ, Szabo-rogers HL, Lee I, Weiss GS, et al. The planar cell polarity effector Fuz is essential for targeted membrane trafficking, ciliogenesis, and mouse embryonic development. 2009;11(10):1225–32. pmid:19767740

27. Heydeck W, Zeng H, Liu A. Planar cell polarity effector gene Fuzzy regulates cilia formation and hedgehog signal transduction in mouse. Dev Dyn. 2009;238(12):3035–42. pmid:19877275

28. Hamblet NS, Lijam N, Ruiz-Lozano P, Wang J, Yang Y, Luo Z, et al. Dishevelled 2 is essential for cardiac outflow tract development, somite segmentation and neural tube closure. Development [Internet]. 2002;129:5827–38. Available from: http://dev.biologists.org/cgi/doi/10.1242/dev.00164 pmid:12421720

29. Nusse R, Clevers H. Wnt/β-Catenin Signaling, Disease, and Emerging Therapeutic Modalities. Cell [Internet]. 2017;169(6):985–99. Available from: pmid:28575679

30. Angers S, Moon RT. Proximal events in Wnt signal transduction. Nat Rev Mol Cell Biol [Internet]. 2009;10(7):468–77. Available from: pmid:19536106

31. Tanaka K, Okabayashi K, Asashima M, Perrimon N, Kadowaki T. The evolutionarily conserved porcupine family is involved in the processing of the Wnt family. Eur J Biochem. 2000;267(13):4300–11. pmid:10866835

32. Caricasole A, Ferraro T, Rimland JM, Terstappen GC. Molecular cloning and initial characterization of the MG61/PORC gene, the human homologue of the Drosophila segment polarity gene Porcupine. Gene. 2002;288(1–2):147–57. pmid:12034504

33. Takada R, Satomi Y, Kurata T, Ueno N, Norioka S, Kondoh H, et al. Monounsaturated Fatty Acid Modification of Wnt Protein: Its Role in Wnt Secretion. Dev Cell. 2006;11(6):791–801. pmid:17141155

34. Gao X, Hannoush RN. Single-cell imaging of Wnt palmitoylation by the acyltransferase porcupine. Nat Chem Biol. 2014;10(1):61–8. pmid:24292069

35. Coombs GS, Yu J, Canning CA, Veltri CA, Covey TM, Cheong JK, et al. WLS-dependent secretion of WNT3A requires Ser209 acylation and vacuolar acidification. J Cell Sci. 2010;123(19):3357–67.

36. Peters T, Perrier R, Haber RM. Focal dermal hypoplasia: Report of a case with myelomeningocele, Arnold-Chiari malformation and hydrocephalus with a review of neurologic manifestations of Goltz syndrome. Pediatr Dermatol. 2014;31(2):220–4. pmid:24387693

37. Rios-Esteves J, Haugen B, Resh MD. Identification of key residues and regions important for porcupine-mediated Wnt acylation. J Biol Chem. 2014;289(24):17009–19. pmid:24798332

38. Leoyklang P, Suphapeetiporn K, Wananukul S, Shotelersuk V. Three novel mutations in the PORCN gene underlying focal dermal hypoplasia. Clin Genet. 2008;73(4):373–9. pmid:18325042

39. Fokkema IFAC, Taschner PEM, Schaafsma GCP, Celli J, Laros JFJ, den Dunnen JT. LOVD v.2.0: The next generation in gene variant databases. Hum Mutat. 2011;32(5):557–63. pmid:21520333

40. Skalka N, Caspi M, Caspi E, Loh Y, Rosin-Arbesfeld R. Carboxypeptidase E: a negative regulator of the canonical Wnt signaling pathway. Oncogene. 2013;32(23):2836–47. pmid:22824791

41. Yu A, Rual JF, Tamai K, Harada Y, Vidal M, He X, et al. Association of Dishevelled with the Clathrin AP-2 Adaptor Is Required for Frizzled Endocytosis and Planar Cell Polarity Signaling. Dev Cell. 2007;12(1):129–41. pmid:17199046

42. Wallingford JB, Harland RM. Neural tube closure requires Dishevelled-dependent convergent extension of the midline. Development [Internet]. 2002;129:5815–25. Available from: http://dev.biologists.org/cgi/doi/10.1242/dev.00123 pmid:12421719

43. Katoh M, Katoh M. Identification and characterization of human PRICKLE1 and PRICKLE2 genes as well as mouse Prickle1 and Prickle2 genes homologous to Drosophila tissue polarity gene prickle. Int J Mol Med. 2003;11(2):249–56. pmid:12525887

44. Yang Y, Mlodzik M. Wnt-Frizzled/Planar Cell Polarity Signaling: Cellular Orientation by Facing the Wind (Wnt). Annu Rev Cell Dev Biol. 2015;31(1):623–46.

45. Wen S, Zhu H, Lu W, Mitchell LE, Shaw GM, Lammer EJ, et al. Planar cell polarity pathway genes and risk for spina bifida. Am J Med Genet Part A. 2010;152(2):299–304.

46. Ishida M, Cullup T, Boustred C, James C, Docker J, English C, et al. A targeted sequencing panel identifies rare damaging variants in multiple genes in the cranial neural tube defect, anencephaly. Clin Genet. 2018;93(4):870–9. pmid:29205322

47. Paulson AF, Prasad MS, Thuringer AH, Manzerra P. Regulation of cadherin expression in nervous system development. Cell Adhes Migr. 2014;8(1):19–28.

48. Curtin JA, Quint E, Tsipouri V, Arkell RM, Cattanach B, Copp AJ, et al. Mutation of Celsr1 Disrupts Planar Polarity of Inner Ear Hair Cells and Causes Severe Neural Tube Defects in the Mouse. Curr Biol. 2003;13:1129–33. pmid:12842012

49. Holmberg J, Clarke DL, Frisen J. Regulation of repulsion versus adhesion by different splice forms of an Eph receptor. Nature. 2000;408:203–6. pmid:11089974

50. Dai D, Zhu H, Wlodarczyk B, Zhang L, Li L, Li AG, et al. Fuz controls the morphogenesis and differentiation of hair follicles through the formation of primary cilia. J Invest Dermatol. 2011;131(2):302–10. pmid:20962855

51. Murdoch JN, Copp AJ. The relationship between sonic hedgehog signaling, cilia, and neural tube defects. Birth Defects Res Part A—Clin Mol Teratol. 2010;88(8):633–52. pmid:20544799

52. Seo JH, Zilber Y, Babayeva S, Liu J, Kyriakopoulos P, de marco P, et al. Mutations in the planar cell polarity gene, fuzzy, are associated with neural tube defects in humans. Hum Mol Genet. 2011;20(22):4324–33. pmid:21840926

53. Steinhart Z, Angers S. Wnt signaling in development and tissue homeostasis. Development. 2018;145.

54. Klingensmith J, Nusse R, Perrimon N. The Drosophila segment polarity gene dishevelled encodes a novel protein required for response to the wingless signal. Genes Dev. 1994;8(1):118–30. pmid:8288125

55. Tamai K, Semenov M, Kato Y, Spokony R, Liu C, Katsuyama Y, et al. LDL-receptor-related proteins in Wnt signal transduction. Nat L. 2000;407:530–5.

56. Angers S, Li T, Yi X, MacCoss MJ, Moon RT, Zheng N. Molecular architecture and assembly of the DDB1CUL4A ubiquitin ligase machinery. Nat Lett. 2006;443(October):3–6.

57. Higa LA, Zhang H. Stealing the spotlight: CUL4-DDB1 ubiquitin ligase docks WD40-repeat proteins to destroy. Cell Div. 2007;2(5).

58. Cang Y, Zhang J, Nicholas SA, Bastien J, Li B, Zhou P, et al. Deletion of DDB1 in Mouse Brain and Lens Leads to p53-Dependent Elimination of Proliferating Cells. Cell. 2006;127(5):929–40. pmid:17129780

59. Mayer-Jaekel RE, Hemmings BA. Protein phosphatase 2A - a “menage a trois.” Trends Cell Biol [Internet]. 1994;4:287–91. Available from: pmid:14731592

60. Mumbly MC, Walter G. Protein Serine/Threonine Phosphatases: Structure, Regulation, and Functions in Cell Growth. Physiol Rev. 1993;73(4):673–99. pmid:8415923

61. Seeling JM, Miller JR, Gil R, Moon RT, White R, Virshup DM. Regulation of β-Catenin Signaling by the B56 Subunit of Protein Phosphatase 2A. Science (80-). 1999;283(March):2089–91.

62. Ikeda S, Kishida M, Matsuura Y, Usui H, Kikuchi A. GSK-3β-dependent phosphorylation of adenomatous polyposis cop gene product can be modulated by β-catenin and protein phosphatase 2A complexed with Axin. Oncogene. 2000;19(4):537–45. pmid:10698523

63. Meister G, Landthaler M, Peters L, Chen PY, Urlaub H, Lührmann R, et al. Identification of novel argonaute-associated proteins. Curr Biol. 2005;15(23):2149–55. pmid:16289642

64. Wang K, Wang X, Zou J, Zhang A, Wan Y, Pu P, et al. MiR-92b controls glioma proliferation and invasion through regulating Wnt/beta-catenin signaling via Nemo-like kinase. Neuro Oncol. 2013;15(5):578–88. pmid:23416699

65. Ishitani T, Kishida S, Hyodo-Miura J, Ueno N, Yasuda J, Waterman M, et al. The TAK1-NLK Mitogen-Activated Protein Kinase Cascade Functions in the Wnt-5a/Ca2+ Pathway To Antagonize Wnt/ -Catenin Signaling. Mol Cell Biol. 2003;23(1):131–9. pmid:12482967

66. Ahn Y, Sanderson BW, Klein OD, Krumlauf R. Inhibition of Wnt signaling by Wise (Sostdc1) and negative feedback from Shh controls tooth number and patterning. Development. 2010;137:3221–31. pmid:20724449

67. Canfield MA, Mai CT, Wang Y, O’Halloran A, Marengo LK, Olney RS, et al. The association between race/ethnicity and major birth defects in the United States, 1999–2007. Am J Public Health. 2014;104(9):1999–2007.

68. Wang M, Marco P, Capra V, Kibar Z. Update on the Role of the Non-Canonical Wnt/Planar Cell Polarity Pathway in Neural Tube Defects. Cells. 2019;8(10):1198.