Developmental and epileptic encephalopathy 89: A novel bi‐allelic variant, molecular dynamics simulation, and a comprehensive clinical and molecular profile

Abstract Objective Gamma‐aminobutyric acid (GABA), the major inhibitory neurotransmitter in the adult central nervous system, plays an important role during embryonic neural network formation. GAD67 is the rate‐limiting enzyme in GABA synthesis, and its deficiency leads to developmental and epileptic encephalopathy 89 (DEE 89). Patients who suffered from this syndrome generally manifested severe to profound neurodevelopmental delay, seizures, and often congenital anomalies such as the cleft palate or/and omphalocele. Up to now, only three papers on this syndrome have been published, and our knowledge about the disease's clinical course and pathophysiology is in its infancy. Methods We used whole‐exome sequencing (WES) and multiple in‐silico tools to detect a potential causative variant in a patient with severe neurodevelopmental delay and refractory epilepsy. Moreover, by molecular docking and molecular dynamics simulation, we investigate the effect of the candidate variant on the GAD67 function and structure. Results WES data analysis revealed a novel deleterious variant (NM_000817.3: c.850C>T; p.Leu284Phe) in the GAD1 gene, which encodes the GAD67 enzyme. Molecular docking and molecular dynamics simulation showed that this variant has deleterious effects on the structure and function of the GAD67. This study's patient, in addition to typical symptoms of the DEE89, showed microcephaly and clonus in the toe, which were novel clinical findings. Significance Our findings expand the mutational and clinical spectrum of DEE 89. Also, by gathering clinical symptoms and genetic findings of previously reported cases, moreover providing a comprehensive clinical picture of the disease, we found that there was no common drug therapy among patients whose epilepsy was controlled. Furthermore, the comparison of clinical symptoms between patients with missense and truncating mutations did not show any significant clinical difference, except that patients with missense mutations did not show cleft palates or omphaloceles.


| INTRODUCTION
Gamma-aminobutyric acid (GABA) is a nonprotein amino acid that enacts as the main inhibitory neurotransmitter in the adult central nervous system (CNS). GABA is involved in different biological functions, including the formation of the human neural network during embryonic development, regulating synaptic transmission, locomotor activity, learning, reproduction, circadian rhythms, and the prevention of insomnia and depression. 1,2 GABA is produced from the decarboxylation of glutamate by the glutamate decarboxylase (GAD) enzyme. 3 There are two isoforms of the GAD enzyme, including GAD67 and GAD65, which are encoded by two distinct GAD1 and GAD2 genes, respectively. 4 GADs for catalysis in their active site demand the pyridoxal phosphate (PLP) cofactor, which is derived from vitamin B6. 3 The two isoforms share nearly 70% similarity in their amino acid sequences, and most of their difference occurs in their N-terminal domain, which determines the distinct sub-cellular location of each isoform. 5,6 GAD67 is a ubiquitously expressed and constitutively active isoform that preferentially synthesizes cytoplasmic GABA and provides more than 90% of GABA in the human CNS. 4,6,7 The other isoform, GAD65, which is localized in nerve terminals and regulates the synthesis of GABA for vesicular release, is transiently activated to meet the need for extra demand for GABA in certain conditions such as stress. 7,8 In the embryonic stages, GAD67 is the major isoform that plays a pivotal role in neuronal development and synaptogenesis, as well as in the normal development of the palate and fetal movements. 4 Previously, several animal model studies have established the role of the Gad1 and Gad2 genes in neural and non-neural development. Animal studies have shown that in Gad1 − / − mice, GABA levels are very reduced, and these mice die due to cleft palate or respiratory dysfunction. 9,10 Inversely, Gad2 − / − mice have a normal level of GABA and do not show any apparent developmental abnormalities but are more prone to seizures, anxiety, and epilepsy. 6,11,12 In 2020, Charton et al. 4 introduced the GAD1 as a causative gene for a new syndrome called developmental and epileptic encephalopathy 89 (DEE 89), but no human disease associated with defects in the GAD2 has been reported.
Patients who suffered from this syndrome generally manifested severe to profound neurodevelopmental delay, seizure, and often congenital anomalies such as the cleft palate or/and omphalocele. 13 Several other diseases associated with nucleotide variations in the GAD1 gene, such as schizophrenia, autism spectrum disorder, and cerebral palsy, have also been reported; however, functional studies are still required to confirm such mutations' pathogenicity. [13][14][15] Up to now, only three papers on DEE 89 have been published, and our knowledge of the disease's clinical course and pathophysiology is in its infancy. 4,13,16 Therefore, investigating more cases gives us a better . clinical and molecular overview of the syndrome and helps us find a genotype-phenotype correlation for GAD1-related diseases.
In this study, we, by whole-exome sequencing (WES) and clinical examination of a patient with early-onset epileptic encephalopathy, found a novel deleterious variant (NM_000817.3: c.850C>T; p.Leu284Phe) in the GAD1 gene, which was evaluated by molecular docking and molecular dynamics (MD) simulation. Also, by collecting the clinical symptoms and genetic findings of previously reported patients manifesting DEE 89, provided a comprehensive overview of the spectrum of clinical symptoms and molecular findings.

| Genetic analysis
WES was performed on the affected sibling (IV-2). Exome-enriched genomic libraries were prepared using SureSelect Human All Exon V7 and sequenced on an Illumina HiSeq 4000 with coverage of 100 × means depth. Read alignment, variant calling, and annotation were performed using the Burrows-Wheeler aligner tool (BWA), the Genome Analysis Toolkit (GATK), and Annovar, respectively. Also, variant filtering and downstream analysis were performed according to our previous study. 17 To confirm the candidate variant, using primer3, two primers (F: GCTTC CTG ACC TCC CTTGT) and (R: CCCTG CCT ACC TTT CATTGC) were designed. Also, these primers were evaluated for SNP absence via BLAT (http://genome.ucsc. edu/cgi-bin/hgBlat). The amplicons were visualized on agarose gel 1.5%, and subsequently, bidirectional Sanger sequencing was performed using an ABI 3130 sequencer (Applied Biosystems). The result was compared with the GAD1 gene reference sequence (NG_021477.2).

| Structure retrieval and analysis
The RCSB Protein Data Bank (PDB) server was used to retrieve the native structure of GAD67. This 3D structure was employed as a template to model the mutated structure of GAD67 using the SWISS-Model web server. 18 The native and mutated structure of GAD67 was comparatively analyzed using PyMol (Version 2.2.3, Schrödinger, LLC.). Also, multiple sequence alignment was performed by PolyPhen-2 web server (http://genet ics.bwh.harva rd.edu/pph2/). The protein structures were visualized using the Chimera software package.

| Molecular docking
The molecular docking study was conducted to predict the best binding mode between the native and mutated proteins and candidate ligands. 19 The study was carried out by using the AutoDock Vina virtual screening tool. Furthermore, the binding interaction of the protein-ligand complexes was analyzed by the BIOVIA Discovery Studio Visualizer v19.1.0.18287 (BIOVIA). The 3D structure of GAD67 was retrieved from the RCSB PDB web server (PDB-ID: 3PV6). The SWISS-Model web server was used to predict the mutated structure of GAD67. Furthermore, the PubChem database was used to retrieve the pyridoxal phosphate structure.

| Molecular dynamics simulation
The dynamic behavior of the native GAD67 protein in comparison to its mutated model was investigated by conducting an MD simulation. The major aim of this study was to give insights into the ability of mutation to produce pathogenic variations or cause disease. The 3D structure of constructed complexes including the mutated model and normal protein structure in interaction with PLP was subjected to energy minimization by the GROMACS 5.1.4 software. 20 The topology of ligands was extracted from Prodrg online server. After generation and subsequently solvation of the water box, counter ions were added to the box. The steepest descent minimization algorithm was run to minimize the system's energy. The equilibration process was used to couple the protein and ligand. During the simulation run, the backbone of the protein was restrained and the solvent molecules and the counter ions were allowed to move. The simulations were continued under periodic boundary conditions for 100 ns in triplicates. Different criteria were measured to analyze the simulation results including root mean square deviation (RMSD), hydrogen bond (Hb), and the radius of gyration (Rg).

| Family description
The proband was a 3-year-old male patient from Iranian parents with consanguineous marriage. No abnormal symptoms were found in prenatal ultrasound scan. The patient was born at term with an uneventful pregnancy but experienced a seizure during delivery. At birth, the head circumference and weight were 35 cm and 2800 g, respectively. After birth due to poor feeding and seizure admitted to neonatal intensive care unit for 20 days. She was also hospitalized for 2 weeks at the ages of one and 1.5 years due to a pulmonary infection. Currently, at the age of three, he manifested severe neurodevelopmental delay with bed confinement, spasticity in limbs, speech absence, and incontinency. He suffers from refractory tonicclonic seizures despite treatment with a combination of vigabatrin and a ketogenic diet. At the age of 6 months, the electroencephalogram (EEG) showed a suppressionburst pattern. The deep tendon reflex (DTR) test showed hyperreflexia and manifested frequent tremors in the toes of both feet. Also, he showed microcephaly, which may be due to cerebral palsy. Ophthalmology assessment revealed horizontal nystagmus and the right eye suffered from strabismus. Moreover, the patient manifested arthrogryposis in his hands, open mouth, high arched palate, large ear, constipation, and dysphasia. He did not manifest cleft palate and omphalocele. Metabolic screening tests were normal. The urogenital examination did not show any abnormalities. MRI showed diffuse cerebral atrophy and severe ventriculomegaly.

| Genetic findings
WES data analysis revealed a single nucleotide variant (SNV) (NM_000817.3: c.850C>T; p.Leu284Phe) in the GAD1 gene. This variant is a missense variant that results in the substitution of a leucine residue with phenylalanine in position 284. The detected variant (rs780519382) has been submitted to dbSNP, but its frequency in the homozygous state in several population databases such as 1000 genome, gnomAD, Exac, and Iranom was 0 (PM2). This variant was located in the well-established functional domain of the GAD67 enzyme (PM1). As best we know, this variant is not reported in LOVD, ClinVar, HGMD databases, or literature. Several pathogenicity predictive tools, including MetaRNN, EIGEN, FATHMM-XF, and MutPred, predicted the variant as pathogenic, and the Combined annotation depletion (CADD) score of this variant was 28.6. Also, alignment among different species revealed the substituted amino acid is a conserved residue Figure 1 (PP3). According to the InterVar tool (https:// winte rvar.wglab.org/) for clinical interpretation of genetic variants, this variant has occurred a gene with a low rate of benign missense variation and in which missense variants are a common mechanism of the disease (PP2). Moreover, the detected variant co-segregated with the phenotype in the family member (PP1). By and large, this variant has (PM1-PM2-PP1-PP2-PP3) criteria and, according to ACMG guidelines categorized as likely pathogenic.

| Structure retrieval and analysis
The native structure of GAD67 was retrieved from the PDB web server (PDB-ID:3VP6). The 3D structure of the mutated sequence of (Leu284Phe) was predicted using the SWISS-Model web server. Figure 2 represents the structure of the native and mutated models. In this figure, the differences between the two structures are colored in red on the structure of the native GAD67. As the result indicated, the helix and sheet structures in the mutated structure are shorter than the native, and these structures are altered to loop. However, increasing the loop structure in comparison with helix and sheet may result in decreased protein stability, and this leads to protein disease-causing features. 21

| Molecular docking
The optimum intermolecular interactions between the ligand and the target macromolecule were determined by conducting the molecular docking study. The intermolecular interactions of the PLP with PubChem CID1051 with the normal and mutated structure of GAD67 were determined using the AutoDock Vina tool. Figure 3 illustrates the active sites and the interactions between the ligand compound and proteins of GAD67. Table 1 shows the Docking Score for the ligand compound and proteins. As expected, the normal protein binding to PLP obtained a higher negative score than the mutated structure, suggesting greater binding efficiency.

| Molecular dynamics simulation analysis
3.5.1 | RMSD analysis of the native and mutated GAD67 The structure of GAD67 (PDB-ID: 3PV6) and its constructed mutated model and their docked structures with PLP (PubChem CID1051) were subjected to 100 ns MD simulations at 300 K to study their stability and mobility. Analyzing the RMSD of the GAD67 normal and mutated protein backbone in Figure 4A indicates their stability during simulation without substantial variation. This observation specifies that the conducted simulation time was long sufficient to reach an equilibrant structure. Furthermore, the RMSD plot of the native GAD generally is lower than its mutated model, which confirms the stability of the native structure. The general behavior of the GAD67 complexes in interaction with PLP is illustrated in terms of RMSD in Figure 4B. In this figure, the RMSD plot of the native GAD67 complex generally is lower than its mutated model with a smooth plot.

| H-bond and radius of Gyration analysis
The steady conformation of a protein can be investigated by counting the number of H-bonds during the simulation. The structural flexibility of the native and mutated proteins can be identified based on the H-bond plot of the complexes. From the number of H-bonds plotted in Figure 4C, the normal complex remarkably established a higher number of H-bonds than the mutated complex structure during the simulation. Therefore, due to the mutation, the structural flexibility of the normal structure is decreased.
The Radius of Gyration is another key factor that identifies the mass-weighted RMSD of atoms from their common center of mass. The factor is used to determine the universal dimension of protein. It can be seen from Figure 4 (D), the complex of the mutated structure with PLP (red color) shows a higher Radius of Gyration than the normal complex. This indicates that a longer variation in the Radius of the Gyration plot of the mutated protein may endure a main structural transition.

| DISCUSSION
Developmental and epileptic encephalopathies (DEEs) are a heterogeneous group of disorders that are characterized by early-onset, often severe epileptic seizures, and EEG abnormalities on a background of developmental impairment. 22 DEEs may be caused by environmental and genetic factors. Genetic DEEs have been associated with mutations in numerous genes involved in diverse processes, such as cell migration, proliferation, and organization, neuronal excitability, synaptic transmission, and plasticity. 22 Next-generation sequencing technology has made a huge revolution in the diagnosis of the underlying genetic factor that causes DEEs. 23,24 Study of large cohorts of affected individuals with next-generation sequencing technology demonstrated that de novo and bi-allelic pathogenic variants are the cause of DEEs in 30%-50% and 11%-38% of patients, respectively. Further genetic mechanisms, including patient's brain mosaicism, oligogenic inheritance, and epigenetic changes, might also contribute to DEE pathophysiology. 22 In this paper, we used WES to investigate an Iranian patient with severe neurodevelopmental delay and refractory seizures. WES data analysis revealed a novel deleterious variant (NM_000817.3: c.850C>T; p.Leu284Phe) in the GAD1 gene, which co-segregated with the phenotype in this family. The GAD1 gene is located on 2q31.1, has 21 exons, and encodes the GAD67 enzyme, which converts L-glutamate to GABA. 3 As a result of mutation, the leucine amino acid has been replaced by phenylalanine, which has a bigger side chain. The mutation is located in the PLP-binding domain, which is important for F I G U R E 2 The visualization of the structure of native (left) and mutated GAD (Leu284Phe; right). The differences between two structures are colored in red on the structure of the native GAD.
the activity of the protein (Figure 1). The bioinformatics study shows that, at the structural level, the mutated protein has shorter helix and sheet elements than the native protein. These regions are replaced with loops in the mutated protein, which may lead to decreased protein stability ( Figure 2). Also, the molecular docking study shows that the normal protein has a higher binding efficiency for the PLP than the mutant protein (Table 1 and Figure 3).
Further investigations using molecular dynamics simulation indicate that the RMSD plot of the native GAD67 in interaction with PLP is generally lower than the mutated model, which confirms the stability of the native structure ( Figure 4B). The result of Figure 4C confirms the reduced number of hydrogen bonds in the mutated protein complex in comparison with the normal protein complex, which confirms the reduced protein stability during the simulation for the mutated complex. Moreover, the complex of the mutated structure in interaction with PLP shows a higher Radius of Gyration than the normal complex ( Figure 4D). This indicates that a longer variation in the Radius of Gyration plot of the mutated protein may endure a main structural transition. Taken together, all these evidences prove the deleterious effect of this mutation on the structure and function of the GAD67 enzyme.
More than a third of the brain's neurons use GABA for synaptic communication, making it the most abundant neurotransmitter. 25 After being synthesized in presynaptic neurons, this molecule is released from axonal terminals into the synaptic cleft and transmits signals through two types of receptors: ionotropic (GABA A and GABA C ) and metabotropic G-protein-coupled receptor (GABA B ). 26 The GABA A and GABAc receptors are hetero-and homo-pentameric ligand-gated ion channels, respectively, which are found only on presynaptic neurons and are responsible for fast synaptic inhibition. 25,27 On the other hand, GABA B receptors are heterodimer metabotropic G-protein-coupled receptors that are located on both pre-and postsynaptic neurons and mediate slow synaptic inhibition. 25,28 GABA B receptors act through phospholipase C and adenylate cyclase and also activate K+ and Ca2+ ion channels via G-coupled proteins. 29 Several studies have revealed that the binding of GABA to these receptors activates signaling pathways in the neurons, which leads to the regulation of important biological functions, including neural cell proliferation, migration, and differentiation. 30,31 Since the GAD67 enzyme is the rate-limiting enzyme in GABA production, its deficiency disrupts GABA production, leading to DEE 89. 4,32 Literature review revealed all reported patients, except for one patient who had moderate ID and ambulation, suffered from severe to profound ID and were bed-confined.
The clinical investigation has shown that spasms and hypotonia are common among reported patients. The EEG was usually abnormal with suppression-burst pattern and hypsarrhythmia. The majority of the patients had normal MRI, although, in some patients, later MRI at older ages revealed defects. Common dysmorphological symptoms in these patients included scoliosis, cleft palate, omphalocele, club foot, plagiocephaly, and arthrogryposis. All of the reported patients manifested seizures. In reported patients, the age-onset of seizures was from the birth to 6 months old, and the patients showed heterogeneous types of seizures such as tonic-clonic, epileptic spasm, myoclonic, and focal motor seizures. Variations in manifested epilepsy were found in both inter-familial and intra-familial, which may be due to genetic background or environmental factors. In some patients, the administration of antiseizure medications led to the control of epilepsy, but in others, epilepsy remained uncontrolled. Von Hardenberg et al. 16 proposed that the combination of vigabatrin and ketogenic diet is an effective epilepsy treatment in patients with DEE 89 syndrome, but this study's patient, as well as the other patient (Family C, Neuray et al.), 13 regardless of the taking of the combination of vigabatrin and ketogenic diet, continued to suffer from epilepsy. Vigabatrin is an irreversible inhibitor of gamma-aminobutyric acid aminotransferase (GABA-AT), the enzyme responsible for GABA catabolism. 28 Considering that more than 90% of the basal level of GABA is provided by the GAD67 enzyme, 6 it seems that prescribing vigabatrin for patients with GAD1 deficiency is not a suitable solution, at least for mutations that result in severe reduction or complete loss of the enzyme function. On the other hand, taking GABA as a food supplement is not an effective method because, according to several studies, this molecule cannot cross the bloodbrain barrier. 28,33 However, the production of GABA derivatives with the ability to cross the blood-brain barrier may be useful for controlling epilepsy in these patients. A detailed investigation of the antiseizure medications used by patients with defects in the GAD1 gene, in whom epilepsy was controlled, revealed no drugs commonalities, indicating that drugs should be prescribed individually for epilepsy treatment in these patients ( Table 2). The patient of this study showed comparable symptoms to the previously reported cases except for the very early-onset seizure (during delivery), clonus in the toe, and microcephaly. The manifested microcephaly in our patient may be due to cerebral palsy. In the GAD1 gene, different types of mutation such as missense, nonsense, deletion, and splice site had been reported, which are distributed along different domains of the protein (Figure 1). The recurrent c.  symptoms of patients with missense mutations and those with truncating mutations were also compared. The manifestations, age of onset, and severity of symptoms did not differ significantly between two groups, except that none of the patients with missense mutations had cleft palates or omphaloceles (Table 2).
In conclusion herein, we reported a patient with a novel deleterious variant (NM_000817.3: c.850C > T; p.Leu284Phe) in the GAD1 gene, which expands the mutational and clinical spectrum of the DEE 89. Also, by gathering clinical symptoms and genetic findings of previously reported cases, in addition to providing an overview on the clinical picture of the disease, we found that there was no common drug therapy among patients whose epilepsy was controlled. Furthermore, the comparison of clinical symptoms between patients with missense and truncating mutations did not show any significant clinical difference, except that patients with missense mutations did not show cleft palates or omphaloceles.

AUTHOR CONTRIBUTIONS
M.A. and M.K. involved in conceptualization and supervision; E.K. involved in formal analysis, investigation, visualization, and writing-original draft. All authors read and approved the final manuscript.