Next Article in Journal
Construction and Validation of a Tumor Microenvironment-Based Scoring System to Evaluate Prognosis and Response to Immune Checkpoint Inhibitor Therapy in Lung Adenocarcinoma Patients
Previous Article in Journal
Genetic Expression between Ageing and Exercise: Secreted Protein Acidic and Rich in Cysteine as a Potential “Exercise Substitute” Antiageing Therapy
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Genes
Volume 13
Issue 6
10.3390/genes13060952
genes-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Eph and Ephrin Variants in Malaysian Neural Tube Defect Families

by
Siti Waheeda Mohd-Zin
1,
Amelia Cheng Wei Tan
1,
Wahib M. Atroosh
2,
Meow-Keong Thong
3,
Abu Bakar Azizi
4,
Nicholas D. E. Greene
5 and
Noraishah Mydin Abdul-Aziz
1,*
1
Invertebrate and Vertebrate Neurobiology Laboratory, Department of Parasitology, Faculty of Medicine, Universiti Malaya, Kuala Lumpur 50603, Malaysia
2
Department of Parasitology, Faculty of Medicine, Universiti Malaya, Kuala Lumpur 50603, Malaysia
3
Department of Paediatrics, Faculty of Medicine, Universiti Malaya, Kuala Lumpur 50603, Malaysia
4
Neurosurgery Unit, Department of Surgery, Faculty of Medicine, Universiti Kebangsaan Malaysia, Cheras 56000, Malaysia
5
Developmental Biology and Cancer Department, Great Ormond Street Institute of Child Health, University College London, London WC1N 1EH, UK
*
Author to whom correspondence should be addressed.
Genes 2022, 13(6), 952; https://doi.org/10.3390/genes13060952
Submission received: 7 April 2022 / Revised: 20 May 2022 / Accepted: 23 May 2022 / Published: 26 May 2022
(This article belongs to the Section Molecular Genetics and Genomics)
Browse Figure
Review Reports Versions Notes

Abstract

:
Neural tube defects (NTDs) are common birth defects with a complex genetic etiology. Mouse genetic models have indicated a number of candidate genes, of which functional mutations in some have been found in human NTDs, usually in a heterozygous state. This study focuses on Ephs-ephrins as candidate genes of interest owing to growing evidence of the role of this gene family during neural tube closure in mouse models. Eph-ephrin genes were analyzed in 31 Malaysian individuals comprising seven individuals with sporadic spina bifida, 13 parents, one twin-sibling and 10 unrelated controls. Whole exome sequencing analysis and bioinformatic analysis were performed to identify variants in 22 known Eph-ephrin genes. We reported that three out of seven spina bifida probands and three out of thirteen family members carried a variant in either EPHA2 (rs147977279), EPHB6 (rs780569137) or EFNB1 (rs772228172). Analysis of public databases shows that these variants are rare. In exome datasets of the probands and parents of the probands with Eph-ephrin variants, the genotypes of spina bifida-related genes were compared to investigate the probability of the gene–gene interaction in relation to environmental risk factors. We report the presence of Eph-ephrin gene variants that are prevalent in a small cohort of spina bifida patients in Malaysian families.

1. Introduction

Neural Tube Defects (NTDs) are a group of congenital malformations in which the formation of the brain and/or spinal cord are compromised as a result of failed closure of the embryonic neural tube, during the fourth week of post-fertilization development [1]. NTDs are among the most common congenital malformations and occur in 0.5 to 10 per 1000 pregnancies globally [2]. The clinical severity of NTDs varies [3,4], with open lesions affecting the brain (anencephaly) and craniorachischisis being the most severe phenotypes. Spina bifida describes lesions affecting the spinal region and affected individuals often exhibit motor and sensory neurological deficits below the level of the lesion [5,6]. In view of the life-changing health and economic consequences of NTDs, considerable effort should be invested in exploring the pathophysiological mechanisms that govern the etiology of NTDs in lieu of the ultimate goal of primary prevention [4].
Spina bifida is etiologically known for both genetic and environmental factors [2]. For instance, several nongenetic risk factors have been identified as having possibilities for prevention by way of maternal folic acid supplementation. On the other hand, genetic alterations have been found to cause NTDs in mice, with more than 250 genes having been shown to cause NTDs when mutated [7,8]. Planar cell polarity (PCP) genes are among the many genes implicated in the mechanism underlying neural tube closure [2,7]. So far, our knowledge of the causative genes in humans is less complete and around 59 genes have been reported to be potentially associated with spina bifida in humans [8,9,10] (Supplementary Table S1). Each gene has been implicated in only a small proportion of NTD patients or in specific populations [11,12,13,14], suggesting that there is considerable heterogeneity underlying the genetic basis of NTDs.
Ephs and ephrins are attractive candidates to be involved in NTDs in view of the involvement of Eph-ephrin interactions and signaling during adhesion and fusion processes [15,16,17,18,19,20,21,22]. Evidence of the possible involvement of Ephs and ephrins in neural tube development comes from neural tube defects models in mice [15,17,22], Xenopus [18,19] and zebrafish [21]. The first murine Eph and ephrin knockouts exhibited an open neural tube defect (ephrin-A5 null mice) [15], the second was a spina bifida occulta model (double heterozygous EphA2tmjrui/+ EphA4rb−2J/+) [22] and the most recent one (Efnb1-deficient embryo) was undefined as the phenotype was not assessed during the closure of primary neurulation [17]. In Xenopus embryo models, Eph and ephrin knockdown (EphA7-morpholino oligonucleotide (MO) and ephrinB2-MO) disrupt cranial neural tube closure [18,19]. Although neurulation does not occur by closure of neural folds, mosaic Eph and ephrin morpholino (loss-of-function) treatment of zebrafish embryos suggest that the Eph and ephrin are specifically and individually required to facilitate integration of progenitor cells during the cross-midline cell division that occurs at the neural keel [21]. Although teleost neural keel formation and mammalian neural tube closure are developmentally distinct events, defects in either process result in severe neural tube defects [1,23,24].
Ephs are the largest group of receptor tyrosine kinases (RTKs) and are known to not only bind to their ligand ephrins but Eph-ephrin complexes are also known to interact or couple with co-receptors such as the TrkB neurotrophin receptor [25,26,27], p75 neurotrophin receptor [28], and Ret tyrosine kinase receptor [29,30]. Eph-ephrin complexes can also activate or inhibit signaling effectors such as protein tyrosine phosphatase (PTPase) SHP2 protein [31,32,33], Ras/Rho family GTPases [34], ADP-ribosylation factor 6 (Arf6) [17], and focal adhesion kinase (FAK) [18]. Ephs are integral membrane receptors, whereas ephrinA ligands are linked to the plasma membrane via a cell surface glycosyl phosphatidylinositol (GPI)-anchor, and ephrinBs are transmembrane ligands. EphrinA5 is implicated in interactions with different splice forms of EphA7 to mediate cellular adhesion or repulsion during neural fold fusion in mouse embryos [15]. EphrinB1 is associated with the maintenance of the structural integrity of apical cell and extracellular matrix (ECM) adhesions in mouse neuroepithelial development [17]. Moreover, a functional role for Eph–ephrin interactions during neural tube closure was suggested by the finding that closure was delayed in mouse embryos cultured with EphA1 and EphA3 fusion proteins used to specifically disrupt ephrinA–EphA receptor interactions [16]. In mouse embryos, EphA2 is shown to be expressed in a lamellipodium-like protrusion structure which extends towards the opposite neural fold [16].
In mouse knockout embryos with perturbation of both alleles of one gene and a single allele of the second gene simultaneously in EphA2 and EphA4 crosses (Epha2tm1Jrui/+Epha4rb−2J/rb−2J and Epha2tm1Jrui/tm1JruiEpha4rb−2J/+), a large number of rounded cells were seen in the open cranial and open spinal neuropores [22]. Furthermore, the double heterozygous embryos carrying the loss of function alleles of EphA2 and EphA4 (Epha2tm1Jrui/+Epha4rb−2J/+) exhibit spina bifida occulta and exencephaly at a penetrance of more than 50%. These findings suggest a dual compensatory role of EphA2 and EphA4 during murine spinal neural tube closure [22]. To date there has been no systematic study to implicate Ephs and ephrins in a human spina bifida cohort although a multitude of Ephs and ephrins have been implicated in human cancers [35,36,37,38,39].
It is important to assess the genetic basis of NTDs in diverse populations. Progress has been made in identifying ‘risk’ variants for NTDs in a number of genes and some of these studies have identified risk factors that may show differing genetic predisposition among ethnic groups [2,40]. Potential variation in genetic predisposition among ethnic groups is suggested by differences in the NTD prevalence between ethnic groups, which in some cases persists after migration to other geographical locations. In the current study, we focused on a cohort from Malaysia. The Malaysian population is multi-ethnic, and this study included individuals of Malay, Chinese and Indian origin.
In particular, the Malaysian population seems biased towards the closed type of Neural Tube Defect seen in preliminary datasets ([41], unpublished data from the Invertebrate and Vertebrate Neurobiology Laboratory, Universiti Malaya). In the preliminary datasets, based on MRI and surgical reports [4], lipomyelomeningocele, lipomeningocele, and meningocele were classified as spina bifida occulta (closed type NTD) whilst myelomeningocele and meningocele were classified as spina bifida aperta (open type NTD). The potential relationship of open and closed spina bifida is further supported in genetic models of mouse knockouts which exhibit both open and closed neural tube defects in the same family of molecules [15,16,22]. Furthermore, there exists a number of mouse models which display spinal dysraphism, whose etiology remains unresolved as the phenotype was not ascertained during primary neurulation [42,43,44]. Therefore, in this study, both types of human spina bifida—4 occulta and 3 aperta—were recruited. The objective for this study was to screen pathogenic variants in Eph and ephrin genes through whole exome sequencing in a spina bifida cohort.

2. Materials and Methods

2.1. Proband Selection for Whole Exome Sequencing and Candidate Gene Validation

Seven unrelated individuals with sporadic spina bifida (5 Malays, 1 Chinese, and 1 Indian) and their healthy family members, where available, were recruited for this study. The 7 patients presented with spina bifida with neurological deficits, encompassing 4 spina bifida occulta and 3 spina bifida aperta. In total, 31 individuals participated in this study. The 7 patients (5 Malays, 1 Chinese and 1 Indian) with sporadic spina bifida were grouped into distinct families, comprising 5 complete trios (mother–father–proband), 1 quartet family (mother–father–twin sibling–proband), and 1 single-parent family (mother–proband) (Supplementary Figure S1). The remaining 10 subjects (6 Malays, 2 Chinese, and 2 Indians) were healthy individuals unrelated to the probands and were included as comparable controls.
There were no descriptions of NTD cases in other family members and none of the patients were the result of consanguineous parents. An assumption was made based on the possibility that NTD phenotypes may be due to mutations in shared genes, thus NTD datasets were grouped based on type of spina bifida and similarity of neurological deficits (Supplementary Table S2).
Informed consent was obtained from all study participants. Collection of samples was performed in accordance with the ethical approval given by the Medical Ethics Committee (MEC) at the University of Malaya Medical Centre (UMMC) (MEC reference number: 914.5).

2.2. Whole Exome Sequencing

DNA was isolated from peripheral blood, saliva samples (Oragene®-DNA OG-500 manufactured by DNA Genotek®, Ottawa, ON, Canada) or buccal mouthwash of spina bifida individuals. Samples were prepared in Illumina TruSeq kit and Agilent SureSelect Target Enrichment kit for whole exome sequencing (WES). The samples were sequenced using the Illumina HiSeq 2000 and HiSeq 4000 platforms. The platforms were able to generate an average read length of 100 bp and had a median depth of 50x coverage per sample. Data alignment and variant analysis were performed [45,46]. The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.

2.3. Exome Datasets Analysis

The variants from the 7 probands of spina bifida were filtered for non-synonymous (missense), frameshift variants, and non-frameshift variants in the exonic region. The variants within the splicing region located adjacent to exons were also included. Variants with minor allele frequency (MAF < 0.01) based on 1000 Genome Phase 1 (annotated based on database build dbSNP135 with November 2010 and October 2011 allele frequency data) were first short-listed. The variant location and nucleotide changes were compared in probands and parents using the platform Galaxy biomedical data tool. Further, the variants were filtered based on a mapping quality of above 50 and a read depth higher than 15 to ensure that the variants were mapped to the reference genome with a high degree of confidence and reduced error probabilities [47,48]. The homozygous variants with alternate allele to total read depth ratio (AD/TD) above 85% and heterozygous variants range between 30% and 70% will be included to get the most accurate candidate variants [49].
Subsequently, the MAF of candidate variants were checked against 1000 Genome Phase 3 (annotated based on dbSNP 142 April 2020), gnomAD version 2.1.1 [50], ExAC, ESP and TOPMED. Candidate variants that are rare (MAF < 0.01), low frequency (MAF 0.01–0.05), and not reported in the 5 global databases were short-listed [51,52]. The MAF of variants in the specific population database, which were the Asian sub-database (gnomAD and ExAC), the East Asian (EAS) sub-database (1000 Genome Phase 3 and gnomAD), the South Asian (SAS) sub-database (1000 Genome Phase 3) and the Singaporean (Singapore Genome Variation Project (SGVP)) variant calling files (VCF) (genotype data build on human genome 18 project assembly, Hg18, March 2006) [53] were included for the inclusion of ancestry-matched controls [54].
Whole exome sequencing analysis was performed to identify variants in 22 known EphAs, EphBs, ephrinAs and ephrinBs genes. All candidate variants in 22 known EphA (EPHA1–EPHA8, and EPHA10), EphB (EPHB1–EPHB4, and EPHB6), ephrinA (EFNA1–EFNA5) and ephrinB (EFNB1–EFNB3) genes were subjected to bioinformatics analysis to determine possible pathogenicity. The loss-of-function (LOF) effect was checked against in silico protein function and conservation tools: Polyphen-2 (Polymorphism Phenotype-2) [55], SIFT (Sorting Tolerance from Intolerance) [56], PROVEAN (Protein Variant Effect Analyser) [57], CADD [58], FATHMM-MKL [59], MutationTaster [60], GERP [61], and PhyloP [62,63]. The position of the variants and the predicted structure of the protein were analyzed using SMART (Simple Modular Architecture Research Tool).

2.4. PCR and Sanger Sequencing for Validation

PCR amplification and Sanger sequencing were used to verify candidate variants from the WES data and to analyze the variants in the family members (where available) and 10 controls (without NTD). A total of 3 sets of primers were designed (Supplementary Table S3).

3. Results

3.1. Whole Exome Sequencing Analysis of Eph and Ephrin Identifies 3 Variants

The possible contributions of the Eph-ephrins were investigated. There were no shared exonic variants in any known mammalian Eph and ephrin genes (EPHA1–EPHA8, EPHA10, EPHB1–EPHB4, EPHB6, EFNA1–EFNA5, and EFNB1–EFNB3) that fulfilled the MAF < 0.01 criteria in all seven spina bifida probands. Each exome dataset was screened individually in all seven probands to look for potential variants. We compared whole exome analysis on five complete trios (SB1A, SB2A, SB5A, SB7A, and SB13A), one quartet-families (SB27A) and one single-parent family (SB3A) and evaluated all possible means of transmission for the Eph and ephrin variant: de novo variants and homozygous and/or compound heterozygous variants for a recessive transmission. From 13 Eph and ephrin variants, there were no variants that represented de novo variants or emulated recessive transmission and appear to segregate exclusively with the phenotype (Supplementary Table S4).
Nevertheless, from the thirteen candidate variants, three variants in EPHA2 (rs147977279), EPHB6 (rs780569137) and EFNB1 (rs772228172) were unreported or MAF < 0.01 and were predicted as pathogenic and conserved in at least seven in silico prediction tools. The EPHA2 variant (rs147977279) was found in a heterozygous state in proband SB2A and the father of SB2A. The EPHB6 variant (rs780569137) was found to be heterozygous in the exome of proband SB5A and the mother of SB5A. The EFNB1 (rs772228172) was found to be heterozygous in the exome of proband SB1A and was genotyped as hemizygous in the father of SB1A (Table 1). Each of the candidate variants found in one of SB2A (rs147977279 in EPHA2), SB5A (rs780569137 in EPHB6) or SB1A (rs772228172 in EFNB1) were detected in single families but not in other spina bifida families or the control group.
The rs147977279 variant in the EphA-type receptor, EPHA2, is a G to C base change at 16,477,423 on chromosome 1, in an individual with spina bifida occulta, SB2A (Table 2). This alters the coding sequence, resulting in a leucine to valine substitution in the ligand binding domain of EPHA2 (Figure 1). The rs780569137 variant in EphB-type receptor EPHB6 is an A to G base change at 142,562,247 on chromosome 7, in an individual with spina bifida aperta, SB5A (Table 2). This alters the coding sequence resulting in a tyrosine to cysteine substitution in the ligand binding domain of EPHB6 (Figure 1). The rs772228172 variant in ephrinB-type transmembrane ligand, EFNB1 is a C to T base change at 68,049,626 on chromosome X, in an individual with spina bifida occulta, SB1A (Table 2). This alters the coding sequence, resulting in an arginine to tryptophan substitution in the signal peptide domain of EFNB1 (Figure 1).

3.2. Unreported and Rare Eph and Ephrin Variants

Among the three variants, MAF of rs780569137 (EPHB6) were not listed in the five global population databases with total number of alleles between 5008 to 282,912 alleles, suggesting that the variant is novel (Table 2). The rs147977279 (EPHA2) was considered rare (MAF < 0.01) in the five global databases of total alleles between 5008 and 251,420 alleles. Similarly, the MAF of rs772228172 (EFNB1) was considered very rare (MAF < 0.01) in gnomAD (MAF = 0.000014; number of alleles (2n) = 4 alleles), and TOPMED (MAF = 0.000008; 2n = 1 allele) in a total of 124,568 to 282,912 alleles. The MAF of rs772228172 (EFNB1) was not reported in the ExAC or ESP global population databases which include a total of 13,006 and 121,250 alleles, suggesting that the variant is novel in the population covered by the ExAC and ESP databases (Table 3).
From the databases surveyed, rs147977279 (EPHA2) was found at allele frequency MAF < 0.01 (rare) in the Asian population in gnomAD exome samples (total number of alleles = 49,008 alleles) or ExAC (total number of alleles = 25,142 alleles). Similarly, rs147977279 (EPHA2) was not reported in the East Asian sub-population (EAS) in 1000 Genome Phase 3 (total number of alleles = 1008), or gnomAD (total number of alleles = 1556 alleles). However, rs147977279 (EPHA2) was found at a low frequency of 0.011 among a total of 978 alleles in the South Asian sub-population (SAS) (Table 3). Likewise, rs772228172 (EFNB1) was found at allele frequency MAF < 0.01 in the Asian population in gnomAD exome samples (total number of alleles = 26,746 alleles). Interestingly, in the SGVP database, whose population group (98 Malays, 99 Chinese, and 95 Indians) is a closer representation of the ethnic origins of Malaysians, none of the three reported SNPs were listed in the genotype data built on the human genome 18 project assembly (Hg18; March 2006), suggesting the variants were principally absent in 292 Singapore individuals (Table 3).

3.3. In Silico Prediction of the Effect of the Variants on the Protein

The deleterious impact scores relating to the protein function of each variant were predicted as probably damaging or possibly damaging in Polyphen2 HumDiv with scores between 0.709 and 1.00, where 1.00 is predicted to be the most damaging. In SIFT, the variants were predicted as damaging with scores between 0.000 and 0.001 with 0.05 as the cut-off value for variants with a deleterious impact. In CADD prediction analysis, the variants were ranked among the 1% most deleterious in the genome with scores above 20 as a cut-off value (Table 4). Similarly, in Polyphen2 HumVar and FATHMM-MKL, rs147977279 (EPHA2) and rs780569137 (EPHB6) were predicted to be damaging with high deleterious impact scores (between 0.94 and 0.997) compared to rs772228172 (EFNB1)’s scores (between 0.61563 and 0.709), whereby 1.00 is predicted to be the most damaging or pathogenic.
In MutationTaster, rs147977279 (EPHA2) was predicted as probably deleterious (disease causing) whereas rs780569137 (EPHB6) and rs772228172 (EFNB1) were predicted to be probably harmless (polymorphism). In Provean, all variants were predicted as neutral with scores ranging between −0.51 and −2.27 with a cut-off value of −2.5. In sequence evolutionary conservation analysis, the wild-type nucleotide at c.G121 for rs147977279 (EPHA2), c.A689 for rs780569137 (EPHB6) and c.C7 for rs772228172 (EFNB1) are highly constrained with GERP scores of 4.17 to 5.21 (cut-off value is −12.36 to 6.16), and PhyloP scores of 1.538 to 3.28894 (cut-off value is −14 to 6) (Table 4).

3.4. Spina Bifida-Related Genes in Probands with Eph and Ephrin Variants

This study extends the investigation into the possible involvement of spina bifida-related genes (Supplementary Table S1, [4,10]) amongst the datasets with Eph and ephrin variants which we validated in SB1A, SB1C, SB2A, SB2C, SB5A and SB5B. Thirty selected variants in reported spina bifida-related genes were short-listed in Table 5. The lists were filtered based on variants that segregate with spina bifida in probands (Table 5, Row A) and variants that had different genotype annotation in probands and parent of probands with the Eph and eph variants (Table 5, Row B). In the latter list, the pattern of the Eph-ephrin variants from this study and the spina bifida-related gene were compared in probands (SB1A, SB2A, and SB5A) and parents of probands with the Eph-ephrin variants (SB1C, SB2C, and SB5B).
In the exome dataset of proband SB2A with the EPHA2 variant (rs147977279), four variants in ALDH1L1 (rs2886059C>A), CUBN (rs62619939C>G), GRHL3 (rs2486668C>G) and PARD3 (rs118153230C>T) were de novo variants with genotypes in both father (SB2C) and mother (SB2B) of proband SB2A and were wildtype (Table 5). Specifically, the allele frequency of the PARD3 variant (rs118153230C>T) was considered rare (MAF < 0.01) and two other variants in ALDH1L1 (rs2886059C>A) and CUBN (rs62619939C>G) were predicted with a damaging effect on the protein function (based on PolyPhen2 HumDiv, PolyPhen 2 HumVar, SIFT, Provean, and CADD; Supplementary Table S5). On the other hand, three variants in BHMT (rs3733890G>A), MTHFD1 (rs2236225G>A) and TNIP1 (rs2233311C>A) that did not comply with the de novo pattern were predicted to have a damaging effect on the protein function.
Likewise, in the exome dataset of proband SB5A with EPHB6 variant (rs780569137), five variants in spina bifida-related genes were likely pathogenic based on MAF and in silico predictions (Table 5, Supplementary Table S5). One of the five variants was found in MTRR (rs1801394A>G) and was annotated as homozygous in proband SB5A and by contrast as heterozygous in SB5B (mother of proband SB5A with EPHB6 variant). The MTRR variant (rs1801394A>G) was predicted to have a damaging effect on the protein function. Another four variants were found in SCRIB (rs781978489G>A), MTHFR (rs200947520G>T, and rs1801133G>A) and XPD (rs1799793C>T) and were annotated as heterozygous in proband SB5A, and as wildtype in SB5B. Notably, the SCRIB variant (rs781978489G>A) was considered a rare variant and was predicted to have a damaging effect on the protein function (Supplementary Table S5). Other than the SCRIB (rs781978489G>A) and MTRR (rs1801394A>G) variants, the allele frequency of MTHFR (rs200947520G>T) was considered rare and two other variants in MTHFR (rs1801133G>A) and XPD (rs1799793C>T) were predicted to have a damaging effect on the protein function. Remarkably, de novo variants were not found in the exome dataset of proband SB5A (Table 5).
As in the exome dataset of proband SB1A with the EFNB1 variant (rs772228172), four variants were annotated as homozygous (CUBN, COMT, PTCH1 and TRDMT1) and eight variants as heterozygous (CUBN, MTRR, and VANGL1) in the reported spina bifida-related genes (Table 5). However, no comparison could be made as the exome dataset was not available for the father of proband SB1A (SB1C) with the EFNB1 variant (rs772228172). Among the twelve variants listed as potential de novo variants (if the variants for SB1C were validated), one variant in VANGL1 (Chr01:116206438C>A) was unreported (novel) as well as predicted to have a damaging effect on protein function (Table 5, Supplementary Table S5). Other than the variant in VANGL1, two other variants in CUBN (rs143400113C>T, and rs369981313G>T) were considered rare (MAF < 0.01) and five other variants in PTCH1 (rs357564G>A), CUBN (rs1801224G>T, and rs2271460A>C), and MTRR (rs1801394, and rs2287780C>T) were predicted to have a damaging effect on protein function.

4. Discussion

Based on whole exome sequencing and subsequent validation by PCR, nearly half of the Malaysian spina bifida families (three out of seven families) analyzed in this study were found to carry variants in Eph or ephrin encoding genes. In each of the families, other than the proband, at least one unaffected member was detected as a heterozygous or hemizygous carrier of the variant allele, suggesting an uncaptured additive gene–gene or gene–environmental interaction that potentially disrupted neural tube closure in the probands.
Based on the interrogation of the 1000 Genome Phase 3, gnomAD, ExAC, ESP and TOPMED databases, these variants are novel or rare with an unreported allele (novel) or MAF of 0.0024 and below (Table 3), based on the accepted convention that any variants recorded in less than 1% or 0.01 of the general population are considered to be rare [52]. Furthermore, EPHB6 (rs780569137) and EFNB1 (rs772228172) variants were not reported, not listed (based on variant calling files, vcf) or rare (MAF < 0.01) in the population specific databases (East Asian, South Asian and Singaporean). The Singaporean database (SGVP) is especially relevant, as the population group genotyped consisted of 292 Singaporean Malay, Chinese and Indian individuals who are a closer representation in genetic ancestries to the Malaysian individuals in our study. The EPHA2 variant (rs147977279) is reported to have low frequency (0.01 < MAF < 0.1) in the South Asian population (1000 Genome Phase 3 database) but the variant was not listed in the Singaporean database (SGVP). Thus, the reported rs147977279 (EPHA2) allele number (2n = 978 alleles) in the South Asian subpopulation (SAS) was perhaps unlikely to have been identified in Singaporean or Malaysian individuals (Table 2).
By considering the potential contributing functions of Ephs and ephrins in neural tube closure [15,16,22], and the in silico prediction of the effect of the variants on the protein (in seven to eight out of the nine prediction tools), we suggest that these variants may play a causative role in the phenotypes. Moreover, the variants scored high in terms of the deleterious impact of its protein function and are remarkably close to the cut-off values (scores 0.94351 to 1.00 for cut-off value 1.00; scores 0.000 to 0.001 for cut-off value 0.05) indicating the high probability for pathogenicity by amino acid substitutions (Table 4). However, the prediction that the EPHB6 and EFNB1 variants are probably harmless in FATHMM-MKL might be due to differences in the categorization and the algorithm’s predefined threshold scores compared to other prediction tools [64]. In PROVEAN, all three variants were predicted to have a neutral effect on the protein function, which might be due to the implicit assumption in the PROVEAN method that large changes to proteins are deleterious (the more dissimilar to the query sequence the worse the score) [57,65]. In this way the PROVEAN scoring algorithm is more generic in evaluating the amino acid sequence (multiple amino acids) compared to others, such as PolyPhen2 and SIFT, which are more specific at the position of interest (single amino acid) [55,56,57,63]. Ultimately, functional studies will be required to confirm the variant pathogenicity on the neural tube development and spina bifida susceptibility.
Eph-ephrin proteins such as EphA7 [15] and ephrinB1 [17], and Eph-ephrin mRNAs such as ephrinA5 [15], ephrinA1, ephrinA3, EphA1, EphA2, EphA4, EphA5 [16] and EphA7 (Wang et al., 2016) have been reported to be expressed in the neural tubes of mouse and Xenopus embryos. Specifically, EphA2 was found to be expressed in surface ectoderm and lamellipodium-like protrusions at the neural fold tips of E9.5 mouse embryos [16]. EphrinB1 was detected in the neuroepithelial cells of E10.5 mouse embryos [17]. However, an association of EphB6 function with neural tube development has not been reported in human or animal model embryos, to date.
Based on the exome datasets, variants in a set of genes relating to environmental factors, and planar cell polarity genes were also found in the three probands with Eph and ephrin variants (Table 5). It would be interesting to further investigate potential gene–gene interactions of variants in Eph and ephrin genes and spina bifida-related genes. For instance, there were three de novo variants with allele frequency that were considered rare (PARD3 rs118153230C>T) or were predicted to have a damaging effect on protein function (ALDH1L1 rs2886059C>A, and CUBN rs62619939C>G) found in proband SB2A with the EPHA2 variant (rs147977279) (Table 5 and Table S5). On the other hand, it is interesting that in the proband with a severe type of spina bifida, such as spina bifida aperta (SB5A), de novo variants were not found in comparison to proband SB2A that were diagnosed as spina bifida occulta. However, larger cohorts of individuals with spina bifida occulta and aperta will be required for a conclusive comparison.
Although the exome dataset of the father of proband SB1A (SB1C) was not available to screen for a de novo or recessive inheritance pattern, it would be interesting to further test if the novel variant in VANGL1 (Chr01:116206438C>A; planar cell polarity gene) that was found in proband SB1A and predicted to have a damaging effect on the protein, could play an additive or combinatorial role with the EFNB1 variant (rs772228172) in spina bifida susceptibility. Another variant in VANGL1 (rs4839469) was also detected in proband SB2A and proband SB5A. This variant had a different genotype in the probands with the Eph variants (EPHA2, and EPHB6) compared to their parents. Furthermore, another variant in SCRIB (rs781978489G>A), which is also a planar cell polarity gene, was also found in proband SB5A. The SCRIB (rs781978489G>A) variant had a rare allele frequency and predicted a damaging effect on the protein as well.
Among the three probands with Ephs and ephrins variants, probands SB1A and SB5A both have variants in Ephs and ephrins, and MTRR (rs1801394) (Table 5). Moreover, in exome datasets of proband SB2A and the father of proband SB2A (SB2C), who are both heterozygous for EPHA2 variant (rs147977279), two different variants in MTRR (rs162036, and rs10380) were also annotated as heterozygous in SB2A but not in SB2C (wildtype). Our previous study on the same cohort revealed that 57% of patients and 83% of parents carried the rs1801394 variant in their MTRR gene, based on either homozygous (G/G) or heterozygous (A/G) genotypes [66]. In that study, we concluded that the MTRR rs1801394 variant may be an NTD risk factor in the Malaysian population based on the prevalence of this variant in other studies and that its association with NTDs differed across populations worldwide [66]. Hereto, the polymorphisms in Eph and ephrin genes, single one-carbon metabolism genes or planar cell polarity genes could play a limited role in overall NTD risk determination [10,67,68,69].

5. Conclusions

In summary, we report the existence of variants in Eph and ephrin genes in three out of seven NTD families in Malaysia. Even though the sample size in this study is small, it provides initial evidence of the need to screen for these predicted pathogenic variants in a larger Malaysian cohort. In each of the families, other than the proband, at least one unaffected member was detected as a heterozygous or hemizygous carrier of the variant allele, suggesting Eph and ephrin variants as a potential screening tool for NTD development among Malaysians. Further validation experiments at the transcript level and a gene perturbation study would help to assess the gene–gene interactions and gene function to support the pathophysiology of NTD etiology amongst Eph and ephrin candidate genes.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/genes13060952/s1, Table S1: List of human spina bifida genes in addition to lists that were reviewed in Greene et al. (2009) Hum. Mol. Genet. 18(R2):p.R113-29 (40 genes) and Mohd-Zin et al. (2017) Scientifica (Cairo) 2017, p. 5364827 (8 genes); Figure S1: Pedigree of the 7 unrelated individuals with sporadic spina bifida and their unaffected family members, where available; Table S2: Clinical features of the 7 probands sent for whole exome sequencing; Table S3: Primers used in this study for screening and validation of the human variants; Table S4: Eph and ephrin variants screened in 7 spina bifida probands (Proband SB1A to Father of SB7A), (Proband SB5A to Father of SB13A); Table S5: Minor allele frequency (MAF) and protein function prediction (PolyPhen 2 HumDiv, PolyPhen 2 HumVar, SIFT, PROVEAN and CADD) of thirty selected variants of reported spina bifida-related genes.

Author Contributions

Conceptualization, N.M.A.-A. and N.D.E.G.; methods, N.M.A.-A., N.D.E.G., A.C.W.T., W.M.A. and S.W.M.-Z.; validation, N.M.A.-A., N.D.E.G., M.-K.T. and S.W.M.-Z.; formal analysis, N.M.A.-A., N.D.E.G., A.C.W.T. and S.W.M.-Z.; investigation, N.M.A.-A., M.-K.T., A.B.A. and S.W.M.-Z.; resources, N.M.A.-A., M.-K.T. and A.B.A.; data curation, N.M.A.-A., N.D.E.G. and S.W.M.-Z.; writing—original draft preparation, N.M.A.-A., N.D.E.G. and S.W.M.-Z.; writing—review and editing, N.M.A.-A., N.D.E.G., M.-K.T. and S.W.M.-Z.; visualization, N.M.A.-A., N.D.E.G., M.-K.T. and S.W.M.-Z.; supervision, N.M.A.-A., N.D.E.G. and A.B.A.; project administration, N.M.A.-A., A.B.A., W.M.A. and S.W.M.-Z.; funding acquisition, N.M.A.-A. and A.B.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Ministry of Higher Education (MOHE) Malaysia under Fundamental Research Grant Scheme, grant number FRGS/1/2019/SKK08/UM/02/17 (FP090-2019A) (N.M.A.-A., and A.B.A.) and FRGS/1/2019/SKK06/UKM/02/06 (A.B.A. and N.M.A.-A.)

Institutional Review Board Statement

The study was conducted in accordance with the Declaration of Helsinki, and approved by the Medical Ethics Committee of University of Malaya Medical Centre (protocol code 914.5 and date of approval 23 August 2016).

Informed Consent Statement

Informed consent was obtained from all subjects involved in the study.

Data Availability Statement

The data presented in this study are available on request from the corresponding author. The data are not publicly available due to concerns regarding participant/patient anonymity.

Acknowledgments

We would like to thank all the patients and families for their support. The authors would like to thank Azlina Ahmad-Annuar from Department of Biomedical Science, Faculty of Medicine, Universiti Malaya for the consultation on the exome annotation management and research data interpretation for initial used. The authors would also like to thank Nivrenjeet Singh and Nisheljeet Singh for copy-editing of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Nikolopoulou, E.; Galea, G.L.; Rolo, A.; Greene, N.D.E.; Copp, A.J. Neural tube closure: Cellular, molecular and biomechanical mechanisms. Development 2017, 144, 552–566. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Greene, N.D.; Copp, A.J. Neural tube defects. Ann. Rev. Neurosci. 2014, 37, 221–242. [Google Scholar] [CrossRef] [Green Version]
  3. Copp, A.J.; Stanier, P.; DE Greene, N. Neural tube defects: Recent advances, unsolved questions, and controversies. Lancet Neurol. 2013, 12, 799–810. [Google Scholar] [CrossRef] [Green Version]
  4. Mohd-Zin, S.W.; Marwan, A.I.; Chaar, M.K.A.; Ahmad-Annuar, A.; Abdul-Aziz, N.M. Spina Bifida: Pathogenesis, Mechanisms, and Genes in Mice and Humans. Scientifica 2017, 2017, 5364827. [Google Scholar] [CrossRef] [PubMed]
  5. Copp, A.J.; Adzick, N.S.; Chitty, L.S.; Fletcher, J.M.; Holmbeck, G.N.; Shaw, G.M. Spina bifida. Nat. Rev. Dis. Primers 2015, 1, 15007. [Google Scholar] [CrossRef] [PubMed]
  6. Oakeshott, P.; Hunt, G.M.; Poulton, A.; Reid, F. Open spina bifida: Birth findings predict long-term outcome. Arch. Dis. Child. 2011, 97, 474–476. [Google Scholar] [CrossRef] [PubMed]
  7. Harris, M.J.; Juriloff, D.M. An update to the list of mouse mutants with neural tube closure defects and advances toward a complete genetic perspective of neural tube closure. Birth Defects Res. Part A Clin. Mol. Teratol. 2010, 88, 653–669. [Google Scholar] [CrossRef] [PubMed]
  8. Wilde, J.J.; Petersen, J.R.; Niswander, L. Genetic, Epigenetic, and Environmental Contributions to Neural Tube Closure. Annu. Rev. Genet. 2014, 48, 583–611. [Google Scholar] [CrossRef] [Green Version]
  9. Pangilinan, F.; Molloy, A.M.; Mills, J.L.; Troendle, J.F.; Parle-McDermott, A.; Signore, C.; O’Leary, V.B.; Chines, P.; Seay, J.M.; Geiler-Samerotte, K.; et al. Evaluation of common genetic variants in 82 candidate genes as risk factors for neural tube defects. BMC Med. Genet. 2012, 13, 62. [Google Scholar] [CrossRef] [Green Version]
  10. Greene, N.; Stanier, P.; Copp, A.J. Genetics of human neural tube defects. Hum. Mol. Genet. 2009, 18, R113–R129. [Google Scholar] [CrossRef]
  11. Lemay, P.; Guyot, M.-C.; Tremblay, E.; Dionne-Laporte, A.; Spiegelman, D.; Henrion, E.; Diallo, O.; De Marco, P.; Merello, E.; Massicotte, C.; et al. Loss-of-function de novo mutations play an important role in severe human neural tube defects. J. Med. Genet. 2015, 52, 493–497. [Google Scholar] [CrossRef] [Green Version]
  12. Lemay, P.; De Marco, P.; Emond, A.; Spiegelman, D.; Dionne-Laporte, A.; Laurent, S.; Merello, E.; Accogli, A.; Rouleau, G.A.; Capra, V.; et al. Rare deleterious variants in GRHL3 are associated with human spina bifida. Hum. Mutat. 2017, 38, 716–724. [Google Scholar] [CrossRef]
  13. Chen, Z.; Lei, Y.; Cao, X.; Zheng, Y.; Wang, F.; Bao, Y.; Peng, R.; Finnell, R.H.; Zhang, T.; Wang, H. Genetic analysis of Wnt/PCP genes in neural tube defects. BMC Med. Genom. 2018, 11, 38. [Google Scholar] [CrossRef] [Green Version]
  14. Wang, M.; de 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, 1198. [Google Scholar] [CrossRef] [Green Version]
  15. Holmberg, J.; Clarke, D.L.; Frisén, J. Regulation of repulsion versus adhesion by different splice forms of an Eph receptor. Nature 2000, 408, 203–206. [Google Scholar] [CrossRef]
  16. Abdul-Aziz, N.M.; Turmaine, M.; Greene, N.D.; Copp, A.J. EphrinA-EphA receptor interactions in mouse spinal neurulation: Implications for neural fold fusion. Int. J. Dev. Biol. 2009, 53, 559–568. [Google Scholar] [CrossRef]
  17. Arvanitis, D.; Béhar, A.; Tóth, P.; Bush, J.O.; Jungas, T.; Vitale, N.; Davy, A. Ephrin B1 maintains apical adhesion of neural progenitors. Development 2013, 140, 2082–2092. [Google Scholar] [CrossRef] [Green Version]
  18. Wang, X.; Sun, J.; Li, C.; Mao, B. EphA7 modulates apical constriction of hindbrain neuroepithelium during neurulation in Xenopus. Biochem. Biophys. Res. Commun. 2016, 479, 759–765. [Google Scholar] [CrossRef]
  19. Ji, Y.J.; Hwang, Y.-S.; Mood, K.; Cho, H.-J.; Lee, H.-S.; Winterbottom, E.; Cousin, H.; Daar, I.O. EphrinB2 affects apical constriction in Xenopus embryos and is regulated by ADAM10 and flotillin-1. Nat. Commun. 2014, 5, 3516. [Google Scholar] [CrossRef] [Green Version]
  20. Laussu, J.; Audouard, C.; Kischel, A.; Assis-Nascimento, P.; Escalas, N.; Liebl, D.J.; Soula, C.; Davy, A. Eph/Ephrin Signaling Controls Progenitor Identities In The Ventral Spinal Cord. Neural Dev. 2017, 12, 10. [Google Scholar] [CrossRef] [Green Version]
  21. Kemp, H.A.; Cooke, J.E.; Moens, C.B. EphA4 and EfnB2a maintain rhombomere coherence by independently regulating intercalation of progenitor cells in the zebrafish neural keel. Dev. Biol. 2009, 327, 313–326. [Google Scholar] [CrossRef] [Green Version]
  22. Abdullah, N.L.; Mohd-Zin, S.W.; Ahmad-Annuar, A.; Abdul-Aziz, N.M. A Novel Occulta-Type Spina Bifida Mediated by Murine Double Heterozygotes EphA2 and EphA4 Receptor Tyrosine Kinases. Front. Cell Dev. Biol. 2017, 5, 105. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  23. Ciruna, B.; Jenny, A.; Lee, D.; Mlodzik, M.; Schier, A.F. Planar cell polarity signalling couples cell division and morphogenesis during neurulation. Nature 2006, 439, 220–224. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Tawk, M.; Araya, C.; Lyons, D.; Reugels, A.M.; Girdler, G.C.; Bayley, P.R.; Hyde, D.R.; Tada, M.; Clarke, J.D.W. A mirror-symmetric cell division that orchestrates neuroepithelial morphogenesis. Nature 2007, 446, 797–800. [Google Scholar] [CrossRef]
  25. Meier, C.; Anastasiadou, S.; Knöll, B. Ephrin-A5 Suppresses Neurotrophin Evoked Neuronal Motility, ERK Activation and Gene Expression. PLoS ONE 2011, 6, e26089. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Marler, K.J.M.; Becker-Barroso, E.; Martínez, A.; Llovera, M.; Wentzel, C.; Poopalasundaram, S.; Hindges, R.; Soriano, E.; Comella, J.; Drescher, U. A TrkB/EphrinA Interaction Controls Retinal Axon Branching and Synaptogenesis. J. Neurosci. 2008, 28, 12700–12712. [Google Scholar] [CrossRef] [Green Version]
  27. Marler, K.J.; Poopalasundaram, S.; Broom, E.R.; Wentzel, C.; Drescher, U. Pro-neurotrophins secreted from retinal ganglion cell axons are necessary for ephrinA-p75NTR-mediated axon guidance. Neural Dev. 2010, 5, 30. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  28. Lim, Y.-S.; McLaughlin, T.; Sung, T.-C.; Santiago, A.; Lee, K.-F.; O’Leary, D.D. p75NTR Mediates Ephrin-A Reverse Signaling Required for Axon Repulsion and Mapping. Neuron 2008, 59, 746–758. [Google Scholar] [CrossRef] [Green Version]
  29. Bonanomi, D.; Chivatakarn, O.; Bai, G.; Abdesselem, H.; Lettieri, K.; Marquardt, T.; Pierchala, B.A.; Pfaff, S.L. Ret Is a Multifunctional Coreceptor that Integrates Diffusible- and Contact-Axon Guidance Signals. Cell 2012, 148, 568–582. [Google Scholar] [CrossRef] [Green Version]
  30. Lisabeth, E.M.; Falivelli, G.; Pasquale, E.B. Eph Receptor Signaling and Ephrins. Cold Spring Harb. Perspect. Biol. 2013, 5, a009159. [Google Scholar] [CrossRef] [Green Version]
  31. Miao, H.; Burnett, E.; Kinch, M.; Simon, E.; Wang, B. Activation of EphA2 kinase suppresses integrin function and causes focal-adhesion-kinase dephosphorylation. Nat. Cell Biol. 1999, 2, 62–69. [Google Scholar] [CrossRef]
  32. Miao, H.; Nickel, C.H.; Cantley, L.G.; Bruggeman, L.A.; Bennardo, L.N.; Wang, B. EphA kinase activation regulates HGF-induced epithelial branching morphogenesis. J. Cell Biol. 2003, 162, 1281–1292. [Google Scholar] [CrossRef] [Green Version]
  33. Saxton, T.M.; Pawson, T. Morphogenetic movements at gastrulation require the SH2 tyrosine phosphatase Shp2. Proc. Natl. Acad. Sci. USA 1999, 96, 3790–3795. [Google Scholar] [CrossRef] [Green Version]
  34. Pasquale, E.B. Eph-Ephrin Bidirectional Signaling in Physiology and Disease. Cell 2008, 133, 38–52. [Google Scholar] [CrossRef] [Green Version]
  35. Lee, H.-H.; Wang, Y.-N.; Yang, W.-H.; Xia, W.; Wei, Y.; Chan, L.-C.; Wang, Y.-H.; Jiang, Z.; Xu, S.; Yao, J.; et al. Human ribonuclease 1 serves as a secretory ligand of ephrin A4 receptor and induces breast tumor initiation. Nat. Commun. 2021, 12, 2788. [Google Scholar] [CrossRef]
  36. Gao, W.; Zhang, Q.; Wang, Y. EphB3 protein is associated with histological grade and FIGO stage in ovarian serous carcinomas. APMIS 2017, 125, 122–127. [Google Scholar] [CrossRef]
  37. Efazat, G.; Novak, M.; Kaminskyy, V.O.; De Petris, L.; Kanter, L.; Juntti, T.; Bergman, P.; Zhivotovsky, B.; Lewensohn, R.; Hååg, P.; et al. Ephrin B3 interacts with multiple EphA receptors and drives migration and invasion in non-small cell lung cancer. Oncotarget 2016, 7, 60332–60347. [Google Scholar] [CrossRef] [Green Version]
  38. Karidis, N.P.; Giaginis, C.; Tsourouflis, G.; Alexandrou, P.; Delladetsima, I.; Theocharis, S. Eph-A2 and Eph-A4 expression in human benign and malignant thyroid lesions: An immunohistochemical study. Med. Sci. Monit. 2011, 17, BR257–BR265. [Google Scholar] [CrossRef]
  39. Barquilla, A.; Pasquale, E.B. Eph Receptors and Ephrins: Therapeutic Opportunities. Annu. Rev. Pharmacol. Toxicol. 2015, 55, 465–487. [Google Scholar] [CrossRef] [Green Version]
  40. Agopian, A.J.; Tinker, S.C.; Lupo, P.J.; Canfield, M.A.; Mitchell, L.E.; The National Birth Defects Prevention Study. Proportion of neural tube defects attributable to known risk factors. Birth Defects Res. Part A Clin. Mol. Teratol. 2013, 97, 42–46. [Google Scholar] [CrossRef] [Green Version]
  41. Sahmat, A.; Gunasekaran, R.; Mohd-Zin, S.W.; Balachandran, L.; Thong, M.-K.; Engkasan, J.P.; Ganesan, D.; Omar, Z.; Azizi, A.B.; Ahmad-Annuar, A.; et al. The Prevalence and Distribution of Spina Bifida in a Single Major Referral Center in Malaysia. Front. Pediatr. 2017, 5, 237. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Tian, T.; Cao, X.; Chen, Y.; Jin, L.; Li, Z.; Han, X.; Lin, Y.; Wlodarczyk, B.J.; Finnell, R.H.; Yuan, Z.; et al. Somatic and de novo Germline Variants of MEDs in Human Neural Tube Defects. Front. Cell Dev. Biol. 2021, 9, 641831. [Google Scholar] [CrossRef]
  43. Rocha, P.P.; Bleiss, W.; Schrewe, H. Mosaic expression of Med12 in female mice leads to exencephaly, spina bifida, and craniorachischisis. Birth Defects Res. Part A Clin. Mol. Teratol. 2010, 88, 626–632. [Google Scholar] [CrossRef] [Green Version]
  44. Pallerla, S.R.; Pan, Y.; Zhang, X.; Esko, J.D.; Grobe, K. Heparan sulfate Ndst1 gene function variably regulates multiple signaling pathways during mouse development. Dev. Dyn. Off. Publ. Am. Assoc. Anat. 2007, 236, 556–563. [Google Scholar] [CrossRef] [PubMed]
  45. Li, Y.; Chen, W.; Liu, E.Y.; Zhou, Y.-H. Single Nucleotide Polymorphism (SNP) Detection and Genotype Calling from Massively Parallel Sequencing (MPS) Data. Stat. Biosci. 2012, 5, 3–25. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  46. Nielsen, R.; Paul, J.S.; Albrechtsen, A.; Song, Y.S. Genotype and SNP calling from next-generation sequencing data. Nat. Rev. Genet. 2011, 12, 443–451. [Google Scholar] [CrossRef] [PubMed]
  47. Li, H.; Ruan, J.; Durbin, R. Mapping short DNA sequencing reads and calling variants using mapping quality scores. Genome Res. 2008, 18, 1851–1858. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  48. Sims, D.; Sudbery, I.; Ilott, N.E.; Heger, A.; Ponting, C.P. Sequencing depth and coverage: Key considerations in genomic analyses. Nat. Rev. Genet. 2014, 15, 121–132. [Google Scholar] [CrossRef]
  49. Patel, Z.H.; Kottyan, L.; Lazaro, S.; Williams, M.S.; Ledbetter, D.H.; Tromp, H.; Rupert, A.; Kohram, M.; Wagner, M.; Husami, A.; et al. The struggle to find reliable results in exome sequencing data: Filtering out Mendelian errors. Front. Genet. 2014, 5, 16. [Google Scholar] [CrossRef] [Green Version]
  50. Karczewski, K.J.; Francioli, L.C.; Tiao, G.; Cummings, B.B.; Alföldi, J.; Wang, Q.; Collins, R.L.; Laricchia, K.M.; Ganna, A.; Birnbaum, D.P.; et al. The mutational constraint spectrum quantified from variation in 141,456 humans. Nature 2020, 581, 434–443. [Google Scholar] [CrossRef]
  51. Mani, A. Pathogenicity of De Novo Rare Variants: Challenges and Opportunities. Circ. Cardiovasc. Genet. 2017, 10, e002013. [Google Scholar] [CrossRef]
  52. Chen, H.; Hendricks, A.E.; Cheng, Y.; Cupples, A.L.; Dupuis, J.; Liu, C.-T. Comparison of statistical approaches to rare variant analysis for quantitative traits. BMC Proc. 2011, 5, S113. [Google Scholar] [CrossRef] [Green Version]
  53. Teo, Y.-Y.; Sim, X.; Ong, R.T.; Tan, A.K.; Chen, J.; Tantoso, E.; Small, K.S.; Ku, C.-S.; Lee, E.J.; Seielstad, M.; et al. Singapore Genome Variation Project: A haplotype map of three Southeast Asian populations. Genome Res. 2009, 19, 2154–2162. [Google Scholar] [CrossRef] [Green Version]
  54. Manrai, A.; Funke, B.H.; Rehm, H.L.; Olesen, M.S.; Maron, B.A.; Szolovits, P.; Margulies, D.M.; Loscalzo, J.; Kohane, I.S. Genetic Misdiagnoses and the Potential for Health Disparities. N. Engl. J. Med. 2016, 375, 655–665. [Google Scholar] [CrossRef]
  55. Adzhubei, I.A.; Schmidt, S.; Peshkin, L.; Ramensky, V.E.; Gerasimova, A.; Bork, P.; Kondrashov, A.S.; Sunyaev, S.R. A method and server for predicting damaging missense mutations. Nat. Methods 2010, 7, 248–249. [Google Scholar] [CrossRef] [Green Version]
  56. Kumar, P.; Henikoff, S.; Ng, P.C. Predicting the effects of coding non-synonymous variants on protein function using the SIFT algorithm. Nat. Protoc. 2009, 4, 1073–1081. [Google Scholar] [CrossRef]
  57. Choi, Y.; Sims, G.E.; Murphy, S.; Miller, J.R.; Chan, A.P. Predicting the Functional Effect of Amino Acid Substitutions and Indels. PLoS ONE 2012, 7, e46688. [Google Scholar] [CrossRef] [Green Version]
  58. Rentzsch, P.; Witten, D.; Cooper, G.M.; Shendure, J.; Kircher, M. CADD: Predicting the deleteriousness of variants throughout the human genome. Nucleic Acids Res. 2019, 47, D886–D894. [Google Scholar] [CrossRef]
  59. Shihab, H.A.; Rogers, M.F.; Gough, J.; Mort, M.; Cooper, D.N.; Day, I.N.M.; Gaunt, T.R.; Campbell, C. An integrative approach to predicting the functional effects of non-coding and coding sequence variation. Bioinformatics 2015, 31, 1536–1543. [Google Scholar] [CrossRef] [Green Version]
  60. Schwarz, J.M.; Cooper, D.N.; Schuelke, M.; Seelow, D. MutationTaster2: Mutation prediction for the deep-sequencing age. Nat. Methods 2014, 11, 361–362. [Google Scholar] [CrossRef]
  61. Davydov, E.V.; Goode, D.; Sirota, M.; Cooper, G.M.; Sidow, A.; Batzoglou, S. Identifying a High Fraction of the Human Genome to be under Selective Constraint Using GERP++. PLoS Comput. Biol. 2010, 6, e1001025. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  62. Pollard, K.S.; Hubisz, M.J.; Rosenbloom, K.R.; Siepel, A. Detection of nonneutral substitution rates on mammalian phylogenies. Genome Res. 2009, 20, 110–121. [Google Scholar] [CrossRef] [Green Version]
  63. Krupp, D.R.; Soldano, K.L.; Garrett, M.E.; Cope, H.; Ashley-Koch, A.E.; Gregory, S.G. Missing genetic risk in neural tube defects: Can exome sequencing yield an insight? Birth Defects Res. Part A Clin. Mol. Teratol. 2014, 100, 642–646. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  64. Raychaudhuri, S. Mapping Rare and Common Causal Alleles for Complex Human Diseases. Cell 2011, 147, 57–69. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  65. Sandell, L.; Sharp, N.P. Fitness Effects of Mutations: An Assessment of PROVEAN Predictions Using Mutation Accumulation Data. Genome Biol. Evol. 2022, 14, evac004. [Google Scholar] [CrossRef] [PubMed]
  66. Tan, A.C.W.; Mohd-Zin, S.W.; Ishak, N.; Thong, M.-K.; Ahmad-Annuar, A.; Azizi, A.B.; Abdul-Aziz, N.M. MTRR gene variant rs1801394 found in Malaysian patients with neural tube defects. Neurosci. Res. Notes 2020, 3, 24–31. [Google Scholar] [CrossRef]
  67. Relton, C.; Wilding, C.; Pearce, M.; Laffling, A.J.; Jonas, P.A.; Lynch, S.A.; Tawn, E.; Burn, J. Gene-gene interaction in folate-related genes and risk of neural tube defects in a UK population. J. Med. Genet. 2004, 41, 256–260. [Google Scholar] [CrossRef]
  68. Kibar, Z.; Salem, S.; Bosoi, C.; Pauwels, E.; De Marco, P.; Merello, E.; Bassuk, A.; Capra, V.; Gros, P. Contribution of VANGL2 mutations to isolated neural tube defects. Clin. Genet. 2010, 80, 76–82. [Google Scholar] [CrossRef] [Green Version]
  69. Kibar, Z.; Torban, E.; McDearmid, J.R.; Reynolds, A.; Berghout, J.; Mathieu, M.; Kirillova, I.; De Marco, P.; Merello, E.; Hayes, J.M.; et al. Mutations in VANGL1 Associated with Neural-Tube Defects. N. Engl. J. Med. 2007, 356, 1432–1437. [Google Scholar] [CrossRef]
Genes 13 00952 g001 550
Figure 1. The EPHA2, EPHB6 and EFNB1 variants in their protein form. The position of the non-synonymous variants (EPHA2:rs147977279, EPHB6:rs780569137 and EFNB1:rs772228172) were represented by red dashed bar (L = Leucine; V = Valine; Y = Tyrosine; C = Cysteine; R = Arginine; W = Tryptophan; EPH_lbd = Eph ligand binding domain; EGF like = epidermal growth factor-like; FN3 = fibronectin-3; Pfam EphA2_TM = protein family of EphA2 transmembrane; TyrKc = tyrosine kinase catalytic domain; SAM = sterile alpha domain; Pfam Ephrin_rec like = protein domain tyrosine-protein kinase ephrin type A/B receptor like; STYKc = Protein Ser/Thr/Tyr kinase or phosphotransferases; Pfam ephrin = protein family ephrin; red box = signal peptide; blue box = transmembrane region; pink box = low complexity region).
Figure 1. The EPHA2, EPHB6 and EFNB1 variants in their protein form. The position of the non-synonymous variants (EPHA2:rs147977279, EPHB6:rs780569137 and EFNB1:rs772228172) were represented by red dashed bar (L = Leucine; V = Valine; Y = Tyrosine; C = Cysteine; R = Arginine; W = Tryptophan; EPH_lbd = Eph ligand binding domain; EGF like = epidermal growth factor-like; FN3 = fibronectin-3; Pfam EphA2_TM = protein family of EphA2 transmembrane; TyrKc = tyrosine kinase catalytic domain; SAM = sterile alpha domain; Pfam Ephrin_rec like = protein domain tyrosine-protein kinase ephrin type A/B receptor like; STYKc = Protein Ser/Thr/Tyr kinase or phosphotransferases; Pfam ephrin = protein family ephrin; red box = signal peptide; blue box = transmembrane region; pink box = low complexity region).
Genes 13 00952 g001
Table
Table 1. Genotypes of spina bifida probands from exome datasets.
Table 1. Genotypes of spina bifida probands from exome datasets.
ProbandsGenotypes
EPHA2
(rs147977279)		EPHB6
(rs780569137)		EFNB1
(rs772228172)
G/G
 		G/C
(Heterozygous)		A/A
 		A/G
(Heterozygous)		C/C
 		C/T
(Heterozygous)		T
(Hemizygous)
Proband-1
(SB2A)		-1
(SB5A)		-1
(SB1A)		-
Parent of probands1
(Mother of SB2A)		1
(Father of SB2A)		1
(Father of SB5A)		1
(Mother of SB5A)		1
(Mother of SB1A)		-1
(Father of SB1A)
Other NTD Probands6-6-6--
Parents and unaffected twin-sibling of other NTD probands12-12-12--
Controls10-10-10--
Total2922922911
Table
Table 2. Position and changes induced by the 3 gene variants.
Table 2. Position and changes induced by the 3 gene variants.
GeneLocationTranscriptcDNA ChangeAmino Acid ChangeNucleotide Change
EPHA2
rs147977279		Chr1: 16477423NM_004431
(exon 2)		c.G121Cp.L41VG to C
EPHB6
rs780569137		Chr7: 142562247NM_004445
(exon 7)		c.A689Gp.Y230CA to G
EFNB1
rs772228172		ChrX: 68049626NM_004429
(exon 1)		c.C7Tp.R3WC to T
Table
Table 3. Minor Allele Frequency (MAF) databases information specific for EPHA2 (rs147977279), EPHB6 (rs780569137), and EFNB1 (rs772228172).
Table 3. Minor Allele Frequency (MAF) databases information specific for EPHA2 (rs147977279), EPHB6 (rs780569137), and EFNB1 (rs772228172).
MAF DatabasesEPHA2
(rs147977279)		EPHB6
(rs780569137)		EFNB1
(rs772228172)
1000 Genome Phase 3
(2N = 5008 alleles)		Global: 0.0022 (2n = 11)
EAS = Not reported (total number of alleles = 1008)
SAS = 0.011 (2n = 978)		Not listedNot listed
GnomAD v2.1.1
(2N = 282,912 alleles)		GnomAD exomes (Global) = 0.002013 (2n = 506/251,420)
GnomAD exomes (Asian) = 0.00853 (total number of alleles = 49,008)
GnomAD (Global) = 0.000019 (2n = 6/31,384)
GnomAD (EAS) = Not reported (2n = 0/1556)		Not listedGnomAD exomes (Global) = 0.000014 (2n = 2/142,298)
GnomAD exomes (Asian) = 0.00004 (total number of alleles = 26,746)
ExAC
(2N = 121,250 alleles)		Global = 0.0024 (2n = 291)
Asian = 0.00927 (total number of alleles = 25,142)		Not listedNot reported
ESP
(2N = 13,006 alleles)		Global = 0.0000308 (2n = 4/13,006)Not listedNot reported
TOPMED
(2N = 124,568 alleles)		Global = 0.000334 (2n = 42/125,568)Not listedGlobal = 0.000008 (2n = 1/125,568)
Singapore Genome Variation Project (SGVP)Not listedNot listedNot listed
2n is allele count of variants in 1 individual (n); EAS = East Asian population; SAS = South Asian population; Not listed = the reported SNP (rs) were not called in variant call database (VCF); Not reported = The reported SNP (rs) is not found (2n = 0) in the database.
Table
Table 4. Predictions of pathogenicity potential and sequence conservation analysis of Eph and ephrin candidate variants.
Table 4. Predictions of pathogenicity potential and sequence conservation analysis of Eph and ephrin candidate variants.
FunctionToolScore Cut-Off/RangeEPHA2
(rs147977279)		EPHB6
(rs780569137)		EFNB1
(rs772228172)
ScorePredictionScorePredictionScorePrediction
PathogenicityPolyphen2
HumDiv		0–10.998Probably Damaging1.00Probably Damaging0.999Probably Damaging
Polyphen2 HumVar0–10.997Probably Damaging0.962Probably Damaging0.709Possibly Damaging
SIFTCut-off = 0.050.001Damaging0.000Damaging0.001Damaging
ProveanCut-off = −2.5−2.27Neutral−1.92Neutral−0.51Neutral
CADD0–10 = Bottom 90%
10–20 = Top 10%
>20 = Top 1%		25.2Top 1% most deleterious
in the genome		24.0Top 1% most deleterious
in the genome		32Top 1% most deleterious
in the genome
FATHMM-MKL0–10.94351Pathogenic0.94876Pathogenic0.61563Pathogenic
MutationTaster-Disease CausingProbably deleteriousPolymorphismProbably harmlessPolymorphismProbably harmless
Sequence
conservation		GERP−12.36 to +6.164.78Evolutionary constrained5.21Evolutionary constrained4.17Evolutionary constrained
PhyloP−14 to +63.25157Conserved3.28894Conserved1.538Conserved
Table
Table 5. The predicted impact of the spina bifida-related variants on the potential interaction with the Eph and ephrin candidates in probands.
Table 5. The predicted impact of the spina bifida-related variants on the potential interaction with the Eph and ephrin candidates in probands.
SB1ASB1C
(Father of Proband SB1A)		SB1B
(Mother of Proband SB1A)		SB2ASB2C
(Father of Proband SB2A)		SB2B
(Mother of Proband SB2A)		SB5ASB5B
(Mother of Proband SB5A)		SB5C
(Father of Proband SB5A)
Ephs and ephrinsEFNB1 (het)EFNB1 (hemi)EFNB1 (wt)EPHA2 (het)EPHA2 (het)EPHA2 (wt)EPHB6 (het)EPHB6 (het)EPHB6 (wt)
(A) Reported spina bifida-related genes with variants segregated with spina bifida**/***PTCH1 (hom)
* CUBN.1 (het)
** CUBN.2 (het)
CUBN.3 (het)
** CUBN.4 (het)
* CUBN.5 (het)
** MTRR.1 (het)
** MTRR.4 (het)
*/** VANGL1.1 (het)
*** COMT (hom)
*** CUBN.6 (hom)
*** TRDMT1 (hom)		Exome dataset not availablePTCH1 (wt)
CUBN.1(wt)
CUBN.2(wt)
CUBN.3 (wt)
CUBN.4 (wt)
CUBN.5 (wt)
MTRR.1 (wt)
MTRR.4(wt)
VANGL1.1 (wt)
COMT (het)
CUBN.6 (het)
TRDMT1 (het)		** ALDH1L1 (het)
** CUBN.7 (het)
GRHL3 (het)
* PARD3 (het)		ALDH1L1(wt)
CUBN.7 (wt)
GRHL3(wt)
PARD3(wt)		ALDH1L1(wt)
CUBN.7 (wt)
GRHL3(wt)
PARD3(wt)		No variants foundNo variants foundNo variants found
# (B) Reported spina bifida-related genes in probands and parent of probands with Ephs and ephrins variantNot relevant due to exome dataset SB1C was not availableExome dataset not availableNot relevant due to exome dataset SB1C was not available** BHMT (het)
** MTHFD1 (het)
MTRR.2 (het)
MTRR.3 (het)
** TNIP1 (het)
VANGL1.2 (het)
*** CUBN.6 (hom)
*** SOD2 (hom)		BHMT(wt)
MTHFD1(wt)
MTRR.2 (wt)
MTRR.3 (wt)
TNIP1(wt) VANGL1.2 (wt)
CUBN.6 (het)
SOD2 (het)		** BHMT (het)
** MTHFD1 (het)
MTRR.2 (het)
MTRR.3 (het)
** TNIP1 (het)
*** VANGL1.2(hom)
*** CUBN.6 (hom)
*** SOD2 (hom)		APEX1 (het)
* MTHFR.2 (het)
** MTHFR.1 (het)
PCMT1 (het)
*/** SCRIB (het)
VANGL1.2 (het)
** XPD (het)
ALDH1A2 (het)
**/*** MTRR.1 (hom)		APEX1(wt)
MTHFR.2(wt)
MTHFR.1(wt)
PCMT1(wt)
SCRIB(wt)
VANGL1.2(wt)
XPD(wt)
ALDH1A(wt)
MTRR.1 (het)		APEX1 (het)
* MTHFR.2 (het)
** MTHFR.1 (het)
PCMT1 (het)
*/** SCRIB (het)
VANGL1.2 (het)** XPD (het)
*** ALDH1A2 (hom)
**/*** MTRR.1 (hom)
The rs number of the variants were EFNB1 = rs772228172C>T, EPHA2 = rs147977279G>C, and EPHB6 = rs780569137A>G, whereas the rs numbers for variants found in the reported spina bifida-related genes were represented in the following alphabetical order as: ALDH1A2 = rs4646626C>T, ALDH1L1 = rs2886059C>A, APEX1 = rs1130409T>G, BHMT = rs3733890G>A, COMT= rs4680G>A, CUBN.1 = rs143400113C>T; CUBN.2 = rs1801224G>T, CUBN.3 = rs3740168G>C, CUBN.4 = rs2271460A>C, CUBN.5 = rs369981313G>T, CUBN.6 = rs1801231G>A, CUBN.7: rs62619939C>G, GRHL3 = rs2486668C>G, MTHFD1 = rs2236225G>A, MTHFR.1 = rs1801133G>A, MTHFR.2 = rs200947520G>T, MTRR.1 = rs1801394A>G, MTRR.2 = rs162036A>G, MTRR.3 = rs10380C>T, MTRR.4 = rs2287780C>T, PARD3 = rs118153230C>T, PCMT1 = rs4816G>A, PTCH1 = rs357564G>A, SCRIB = rs781978489G>A, SOD2 = rs4880A>G, TNIP1 = rs2233311C>A, TRDMT1 = rs11254413G>A; VANGL1.2 = rs4839469G>A, VANGL1.1 = Chr01:116206438C>A, and XPD = rs1799793C>T; # Although, the variants were not segregated with NTD, however, the genotype of the parents with Ephs and ephrins was wildtype and heterozygous compared to the proband; * The predicted impact of the variants on the likely interaction with the Ephs and ephrins candidates from this study was annotated for MAF < 0.01 (gnomAD public database); ** Predicted as deleterious in one of the five in silico predictions. An extended list of MAF and five protein function predictions are available in Supplementary Table S5). wt/bold font = wildtype/reference genotype; *** hom/bold font = homozygous genotype; het = heterozygous genotype).
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Mohd-Zin, S.W.; Tan, A.C.W.; Atroosh, W.M.; Thong, M.-K.; Azizi, A.B.; Greene, N.D.E.; Abdul-Aziz, N.M. Eph and Ephrin Variants in Malaysian Neural Tube Defect Families. Genes 2022, 13, 952. https://doi.org/10.3390/genes13060952

AMA Style

Mohd-Zin SW, Tan ACW, Atroosh WM, Thong M-K, Azizi AB, Greene NDE, Abdul-Aziz NM. Eph and Ephrin Variants in Malaysian Neural Tube Defect Families. Genes. 2022; 13(6):952. https://doi.org/10.3390/genes13060952

Chicago/Turabian Style

Mohd-Zin, Siti Waheeda, Amelia Cheng Wei Tan, Wahib M. Atroosh, Meow-Keong Thong, Abu Bakar Azizi, Nicholas D. E. Greene, and Noraishah Mydin Abdul-Aziz. 2022. "Eph and Ephrin Variants in Malaysian Neural Tube Defect Families" Genes 13, no. 6: 952. https://doi.org/10.3390/genes13060952

APA Style

Mohd-Zin, S. W., Tan, A. C. W., Atroosh, W. M., Thong, M. -K., Azizi, A. B., Greene, N. D. E., & Abdul-Aziz, N. M. (2022). Eph and Ephrin Variants in Malaysian Neural Tube Defect Families. Genes, 13(6), 952. https://doi.org/10.3390/genes13060952

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop