The contribution of de novo coding mutations to meningomyelocele

Abstract

Meningomyelocele (also known as spina bifida) is considered to be a genetically complex disease resulting from a failure of the neural tube to close. Individuals with meningomyelocele display neuromotor disability and frequent hydrocephalus, requiring ventricular shunting. A few genes have been proposed to contribute to disease susceptibility, but beyond that it remains unexplained¹. We postulated that de novo mutations under purifying selection contribute to the risk of developing meningomyelocele². Here we recruited a cohort of 851 meningomyelocele trios who required shunting at birth and 732 control trios, and found that de novo likely gene disruption or damaging missense mutations occurred in approximately 22.3% of subjects, with 28% of such variants estimated to contribute to disease risk. The 187 genes with damaging de novo mutations collectively define networks including actin cytoskeleton and microtubule-based processes, Netrin-1 signalling and chromatin-modifying enzymes. Gene validation demonstrated partial or complete loss of function, impaired signalling and defective closure of the neural tube in Xenopus embryos. Our results indicate that de novo mutations make key contributions to meningomyelocele risk, and highlight critical pathways required for neural tube closure in human embryogenesis.

Extended Data Fig. 1 |. Power calculation of estimating a cohort size for DNM detection.

Power calculation showing potential number of discovered genes compared with cohort size (350 trios), for two different v values (enrichment ratio of loss of function variants in case versus control) and two different k values (assumed number of risk genes). For instance, if there are 50 genes to discover (k=50), a cohort of 400 trios will identify 16 genes if LOF variants are 2.5x more common in affected (v=2.5). All calculations manage a conservative false discovery rate (FDR). Gray dash: FDR.

Extended Data Fig. 2 |. Spatial expression of DNM genes with MERFISH in E9.5 mouse embryos.

Spatial expression of the MM genes with damaging DNMs. a, Gene expression of marker genes for seven selected cell types (neuron, neural progenitor, pre-epithelial to mesenchymal transition neural crest progenitor (NC progenitor), neural crest, mesoderm, dorsal root ganglia, and blood), in two embryonic replicates. b, Six spatial expression pattern of damaging DNM genes with specific (left) and broad (right) expression. Full MERFISH image of 36 genes can be found in the GitHub (https://github.com/Gleeson-Lab/Publications/tree/main/MM_DNM).

Extended Data Fig. 3 |. Cell type expression of the DNM genes in MERFISH.

a, Expression at E9.5 of the 36 damaging DNMs in seven cell types: neuron, neural progenitor, pre-epithelial to mesenchymal transition neural crest progenitor (Pre-EMT-NCP), neural crest, mesoderm, dorsal root ganglia, and blood. Indeterminate refers to the cells that were not specified with the marker genes designed for the seven cell types. b, Expression of marker genes used for specifying the cell types in MERFISH. Marker genes are shown within the cell type category which they represent.

Extended Data Fig. 4 |. A human protein network constructed with damaging DNMs contributing to MM risk.

By using the 187 damaging MM DNM genes, a propagated network was generated with NetColoc with a background protein network PCNet, incorporating 439 nodes and 2,447 edges. Big blue circle: damaging DNM genes, Small purple circle: propagated gene, green border: known mouse NTD genes. The network is visualized with Cytoscape with STRING database.

Extended Data Fig. 5 |. H1149P patient mutation impairs TIAM1 activity and PLCE1 patient mutation E623Q leads to diminished GTP-bound RhoA.

a, H1149P mutation is located within the Dbl homology (DH) domain responsible for GEF activity. TIAM1 contains an N-terminal pleckstrin homology (PH), coiled-coiled (CC), extension (Ex), RAS binding (RBD), PDZ, Dbl-homology (DH) and PH domains, with the patient mutation falling within the DH domain. b, Schematic of PLCE1 protein with domains annotated. Patient E623Q mutation is located in the Ras GEF domain. PLCE1 contains a Guanine nucleotide exchange factor for Ras-like small GTPases (RAS GEF), Pleckstrin Homology (PH), Phospholipase C catalytic domain X (PLCX), Phospholipase C catalytic domain Y (PLCY), Protein Kinase C conserved region 2 (C2), RAS association domain 1 (RA1), and RAS association domain 2 (RA2). c, Construct expression H1149P (n = 76) is equivalent to wildtype (n = 85) in Phalloidin quantification. d, Construct expression H1149P in constitutive active (C.A.) (n = 53). Src Rac1 Förster resonance energy transfer (FRET) is equivalent to wildtype. P value adjusted with Bonferroni. Kruskal-Wallis followed by a two-sided pairwise Wilcoxon test, P value adjusted with Bonferroni. Data shown with Hampel filter. Error bar: standard error of the mean. P values: ns: not significant. e, Active GTP-bound form of RhoA precipitated from HEK293 expressing Myc-tagged PLCE1 using a GST-rhotekin pulldown assay. Overexpression of WT PLCE1 resulted in a substantial decrease in relative RhoA activity compared with mock cells. Compared to WT, cells transfected with variant forms of PLCE1 exhibited marked differences in GTP-bound RhoA.

Extended Data Fig. 6 |. The P168L patient mutation impairs TNK2 activity and intronic mutation at the splice acceptor site before exon 47 leads to alternative splicing of DNAH5.

a, P168L mutation is located within the kinase domain. TNK2 contains sterile alpha motif (SAM), Src homology 3 (SH3), CDS42 and RAC-interactive binding (CRIB), Mig6 homology region (MHR), and ubiquitin-associated domain (UBA). b, Blots for the A156T kinase dead, wildtype, and the patient mutation P168L. Lysates were probed with pY284-Ack1 (top), Ack1-flag (middle), and gamma-tubulin (bottom). Repeated independently with similar results four times. Ack1 refers to TNK2. c, TNK2 P168L patient mutation impaired WASP phosphorylation from immunoprecipitation (IP) kinase assay, compared to WT and kinase dead A156T. d, Location of chr5:13807727 T > G patient mutation in DNAH5 gene. e, Primer design for detecting altered splicing in DNAH5 cDNA - pair (i) spanning exons 46–48, pair (ii) spanning exons 46–49 and control pair (iii) spanning exons 1–4. f, RT-PCR results using primers listed in e showing altered splicing for exon 46–48 and 46–49 in patient cDNA, but not in controls.

Extended Data Fig. 7 |. KDM1A R332C patient mutation and VWA8b patient mutation R230G significantly reduces protein expression levels.

a, Schematic of KDM1A protein domains - R332C patient mutation is located in the amino-oxidase domain of KDM1A protein. b, Protein levels of WT and R332C KDM1A detected by western blot from HEK293T cells transfected with pEGFP-C2-KDM1A WT or R332C plasmids; c, Quantification of GFP-KDM1A band intensities from b normalized with β-actin loading control (n = 3). Bar: median, Error bar: interquartile range. Two-tailed unpaired t test with Welch’s correction, P value * = 0.0316. d, Schematic of protein domains for human VWA8b consisting of NTPase, Walker A (WA), ATP binding, and Walker B (WB) domains and patient mutation R230G in the NTPase domain. e, Protein levels of HA-tagged mVwa8b empty vector (EV), WT and R230 overexpressed in HEK293T cells detected using anti-HA antibody; alpha tubulin used as loading control. f, Quantification of HA band intensity from panel b normalized to loading control (n = 3). Bar: mean, Error bar: standard deviation of mean (SEM). one-way ANOVA ****: P value < 0.0001.

Extended Data Fig. 8 |. Validation of Spen and Mink1 knockdown.

a, Schematics of SPEN and MINK1 protein with domains annotated and patient mutations. RRM, RNA Recognition motif. MINT, Mxs2-interacting protein. SPOC, Spen paralogue and orthologue SPOC. CNH, Citron homology domain. b, Dorsal views of Xenopus laevis embryos subjected to in situ hybridization for Pax3 to visualize the neural folds in Spen morphants. c, RT-PCR confirmed that Spen MO reduced the amount of normally spliced Spen transcript. d. Dorsal views of Xenopus l. embryos subjected to in situ hybridization for Pax3 in Mink1 morphants. e, Validation of Mink1 morpholinos by RT-PCR. c,e, Each experiment was performed independently at least twice with similar results. f, Dorsal views of embryos injected with Spen gRNAs only or gRNAs with Cas9, with the accompanying chromatogram showing Sanger sequencing at the CRISPR target site. Control embryos injected with gRNAs only (#1 and #2) developed normally and exhibited an intact sequence, while embryos injected with Spen gRNAs and Cas9 (#3-#6) displayed neural tube defects and mosaic mutations at the CRISPR target site. g, Dorsal views of embryos injected with Mink1 gRNAs only or Mink1 gRNAs with Cas9, with the accompanying chromatogram showing Sanger sequencing at the CRISPR target site. Mink1 crispants (#2-#4) exhibited neural tube defect phenotypes and mosaic mutations at the CRISPR target site in both the L and S alleles of Mink1.

Extended Data Fig. 9 |. Validation of Whamm knockdown.

a, Schematics of SPEN and MINK1 protein with domains annotated and patient mutations. JMY, Junction-mediating and WASP homolog-associated domain, JMY_N, N-terminal of JMY. WH2, WASP-homology 2 domain b, The neural tube closure defect phenotype induced by Whamm MO (10 ng) was rescued through the injection of Whamm mRNA (700 pg). Embryos injected only with mRNA showed no significant phenotype. c, Dorsal views of Xenopus embryos at Stage 19, quantified with Pax3 for in situ hybridization to visualize the neural folds. d, Quantification of the average distance between neural folds in Whamm MO with rescue Whamm mRNA. The rescue experiment was repeated independently with similar results, with multiple independent experiments; Control (n = 16), Whamm MO (n = 20), Whamm MO + mRNA (n = 19), mRNA (n = 14). Box plot indicates the median (center line), the interquartile range (bounds of the box), and the whiskers represent the minimum and maximum. P-values by one-way ANOVA, followed by Tukey’s multiple comparison test: **** < 0.0001, ns: not significant. e, RT-PCR confirmed that Whamm MO reduced the amount of normally spliced Whamm transcript. f, Schematic showing gRNA regions designed to target Whamm gene and primer sites for genotyping. g, Control embryos (#1-#5) developed normally, while crispants (#6-#10) displayed neural tube defects. h, Genotyping in the target area. PCR products from control embryos (#1-#5) are approximately 631 bp, while those from Whamm crispants (#6-#10, except #8) are around 331, indicating a deletion of approximately 300 bp. e,h, Each experiment was performed independently at least twice with similar results. i, Comparison of sequence between control (#4 in panel g is shown) and crispant (#8 in panel g is shown) at the Whamm CRISPR target site. Although the #8 embryo did not exhibit the 300 bp deletion, Sanger sequencing result shows it has mosaic mutations.

Extended Data Fig. 10 |. Validation of Nostrin knockdown and synergistic effect with Whamm.

a, Schematics of NOSTRIN protein with domains annotated and patient mutations. FCH, Fes/CIP4, and EFC/F-BAR homology domain. F-BAR, Fes/CIP4 homology – Bin-Amphiphysin-Rvs domain. HR1, REM-1 domain. SH3, Src homology 3 domain. b, RT-PCR confirmed that the splice-blocking MO for Nostrin S reduced the amount of normally spliced Nostrin transcript. Experiment was performed independently at least twice with similar results c, Dorsal views of embryos injected with Nostrin gRNAs only or Nostrin gRNAs with Cas9, with the accompanying chromatogram showing Sanger sequencing at the CRISPR target site. Control embryos injected with gRNAs only (#1) developed normally and exhibited intact sequences, while embryos injected with Nostrin gRNAs and Cas9 (#2-#4) displayed neural tube defects and mosaic mutations at the CRISPR target site. d, Neural folds visualized by in situ hybridization for Pax3 in Nostrin splice-blocking MO and/or Wham MO. e, Dorsal views of late neurula embryos injected with Nostrin translation-blocking MO and/or Whamm MO.

Fig. 1 |. Enrichment of damaging DNMs in MM versus control.

a, DNMs categorized by predicted functional impact: LGD, D-mis, D-mis-HC (highly constrained, called from a meta predictor), tolerant mis and synonymous. DNM rates per proband in MM (left sides of the triangle with strips) and control (right sides). Areas represent the ratio of DNM rates between MM and controls in each functional category. b,c, Variant rate (10⁻⁸, left y axis) and theoretical rate (DNM rate per child, right y axis) are shown for different categories. Statistical analysis of the ratio of DNM rates for MM and control cohorts, denoted by the ratio of the two Poisson rates, or rate ratio (RR), calculated for all genes (b, n = 19,658) and constrained genes (c, n = 3,060; pLI ≥ 0.9). The theoretical rate was calculated by normalization of the variant rate with the total size of the hg38 coding region (59,281,518 base pairs). P-values were calculated by a one-sided rate-ratio test. Error bars: 95% confidence interval of two sample t-test; ***P < 0.001, **P < 0.01, *P < 0.05; NS, not significant.

Fig. 2 |. Functional convergence of genes implicated by damaging DNMs.

a, Proportion of singletons (DNM occurred in one trio) and doubletons (DNMs occurred in two independent trios) damaging (LGD or D-mis) genes in controls and MM. Five doubleton genes were annotated with variant functional categories (LGD, D-mis-HC or D-mis). b, A protein–protein interaction network composed of MM damaging DNM genes. Genes connected by at least one edge are shown; unconnected (orphan) genes are shown in Supplementary Fig. 4. Node colours denote the variant functional categories: LGD (purple), D-mis (light pink) and D-mis-HC (dark pink). Doubleton genes have green borders (three doubletons are shown). Edge thickness denotes the confidence score of the protein interaction, defined in the STRING database. c,d, The total number of edges (interactions) was counted in 100,000 networks randomly generated with bootstrapping by selecting gene set size of 108 (that is, 80% of the 135 control DNM genes) in MM and control cohorts. The density of the number of edges (c) is shown in 100,000 bootstrapped networks. Higher numbers of edges denote denser network interconnections. P-values were calculated by two-sided Wilcoxon rank-sum test. ***P < 2.2 × 10⁻¹⁶. The number of nodes with 0–10 edges is shown (d) in controls (green) and MM (yellow) gene sets. Points show the median and error bars denote the first and third quartiles. P-values were calculated by a one-sided Wilcoxon rank-sum test. ****P < 2.2 × 10⁻¹⁶; NS, not significant. e, Gene Ontology (GO) term network visualization for genes overlapping with 187 MM DNM genes with statistical significance (a GO enrichment analysis). GO terms of functional relevance are connected by an edge through commonly involved genes (small circles). Node size indicates the significance of the terms. Degree of connectivity between terms (edges), kappa statistics by ClueGO. NTRK1 is also known as TRKA.

Fig. 3 |. Functional modules that contribute to MM risk.

a, Functional submodules clustered from the propagated network, with the 187 damaging DNM genes from the MM cohort used as seeds. They are clustered with a Leiden algorithm with co-expression values from the STRING database as attributes, five submodules annotated with FDR < 10⁻⁵ by GO biological process, KEGG or Reactome databases shown. Functional terms annotated with FDR < 10⁻⁸ are shown with light purple background. b, Most significantly enriched GO biological processes with FDR < 10⁻⁸ shown for GTPase involved actin cytoskeleton and microtubule-based processes. FDR value: −log₁₀ scale.

Fig. 4 |. Functional validation of damaging DNMs related to actin polymerization.

a, Damaging DNM genes highlight different functional submodules ofnine mutated genes, including the D-mis genes (TNK2, TIAM1, PLCE1, KDM1A and DNAH5) and LGD genes (WHAMM, NOSTRIN, SPEN and MINK1). b, A TIAM1 H1149P patient mutation led to fewer lamellipodia (filamentous actin) than the wild type (WT) imaged with phalloidin. Experiment independently repeated three times.Scale bars,10μm. c, A TIAM1 H1149P patient mutation (n=27) decreased RAC1 activation with constitutively active SRC, observed using Förster resonance energy transfer (FRET) compared with vector (n=37) and WT TIAM1 (n=32). Kruskal–Wallis test followed by a two-sided pairwise Wilcoxon test; P value adjusted with a Bonferroni correction. Data are shown with a Hampel filter. Error bars show s.e.m. d, Quantification of GTP-bound RHOA (n=5) of E623Q PLCE1. PLCE1 contains a RAS-like small GTPase (RAS GEF) domain so there is crosstalk between RAS and RHOA, PLCE1 seems to have an indirect effect on RHOA rather than a direct one. Each dot represents one independent experiment (n=5). One-way analysis of variance (ANOVA) with a posthoc Bonferroni correction Error bar shows s.d.a.u., arbitrary units. e,g,i, In X. laevis, knockdown of spen (e), mink1 (g) and nostrin and whamm (i) cause an open neural tube. Dorsal views of Xenopus embryos injected with MOs at the late neurula stage are shown. f,h,j, Quantification of the average distance between neural folds in Xenopus spen (f; control (n=25), 10ng (n=17), 20ng (n=19)), mink1 (h; control (n=13), 5ng (n=12), 10ng (n=10)), and nostrin (5ng (n=16), 10ng (n=17)) and whamm (5ng (n=16), 10ng (n=17)) (j) embryos injected with MOs with pax3 insitu hybridization (control (n=14), 5ng (n=13) and 10ng (n=14)). Error bars show s.e.m. P-values calculated by one-way ANOVA, followed by Tukey’s multiple comparison test. ****P < 0.0001, *** P < 0.001, **P < 0.01, *P < 0.05; NS, not significant.