In silico prediction of structural and functional effects of the E71G mutation on GARS and its role in CMT2D and HMN5A pathogenesis

Background The E71G mutation in GARS results in disease type 2D (CMT2D) with overlapping distal hereditary motor neuropathy type 5A (HMN5A) phenotype, but the pathogenesis is unclear. This study aimed to predict the structural and functional effects of this mutation on GARS’ interaction with some of its putative or actual ligands in relation to CMT2D and HMN5A manifestation. Methods Putative or actual ligands of GARS were identified from 10 protein databases using Cytoscape 3.4.0. Thirty-one (31) ligands met our criteria and were selected for the study. GRAMM-X was utilized for docking simulations. The RMSD value difference between wild-type and docked ligands, change in binding free energy, and changes in interface between wild type and native GARS were analyzed. Results Our in silico analysis showed that GARS interaction with 19 distinct putative or actual ligands increased in stability, 10 decreased in stability, and two (2) were unchanged energetically. Data mining and network analysis showed that the ligands are involved in apoptosis, autophagy, and immune response, as well as the ubiquitin proteasome system and NGF/TrkA signaling. Conclusions These findings indicate that several biological processes (apoptosis, autophagy, and immune response) and pathways (ubiquitin proteasome system and NGF/TrkA signaling) may be affected by the E71G mutation. This suggests that GARS is a key molecule in CMT2D and HMN5A pathogenesis and shows promise as a target for novel drugs. Further experimentation is needed to verify these findings.

is to catalyze the attachment of glycine to glycyl-tRNA during translation [4]. It is also involved in ATP binding, where it synthesizes diadenosine tetraphosphate, a signaling molecule important in cell regulatory pathways [5]. Its gene products may be in the cytoplasm, mitochondria [6] and axons [7,8] and are associated with neuronal granules [8].
Mutations in GARS are known to cause two main neurodegenerative diseases: Charcot-Marie-Tooth Disease 2D (CMT2D) [9,10] and Distal Hereditary Motor Neuropathy Type 5A (HMN5A) [9,11]. In CMT2D, there is a distal sensory loss [9] and both the motor nerves and the sensory nerves are affected. In HMN5A, there is no sensory loss and the motor nerves can be severely damaged, while the sensory nerves remain normal [12]. There are 14 known mutations related to CMT2D and HMN5A. This study focused on the E71G point mutation (glutamic acid → glycine at position 125). It is known to result in CMT2D with overlapping HMN5A phenotype [7][8][9]13]. Both CMT2D and the HMN5A show an autosomal dominant mode of inheritance [14]. These diseases are very rare. For the past ten years, only about a dozen families have been identified to be affected [14]. There is no known mechanism as to how these diseases manifest.
GARS was the first ARS protein to be characterized in relation to CMT [8]. Recent studies suggest that the mechanism is largely by gain of function effects consistent with autosomal dominant mutations rather than loss of function, as previously suspected [8,[15][16][17]. Neomorphic interactions of GARS with other ligands are believed to play a role in mechanisms for CMT2D manifestation. Some mutations resulted in the same neomorphic conformational opening, exposing a region hidden in the folding of the wild type GARS protein and providing a possible site for neomorphic interactions [18]. There is strong molecular support on aberrant dimerization resulting in neomorphic interactions as an underlying cause, but other mechanisms cannot be ruled out completely [18]. The observed differences in localization of wild type and mutant proteins as well as granule formation is still unexplained [8,19]. The disease manifestation also seems to involve defects in developmental pathways involving the VEGF/Nrp1 signaling pathway and neuromuscular junction machineries, which may or may not be related [20][21][22]. Although it can be explained why degeneration in motor neurons result from CMT2D mutation, it is still not known why the disease mostly affects distal motor neurons specifically [23]. In addition, why the upper limbs display more severe symptoms than the lower limbs in CMT2D remains a mystery. The connection to sensory neurons, the second type of neurons affected in the disease, is also unresolved.
Similar to CMT2D, the molecular mechanism behind the occurrence of HMN5A is not widely understood. In vitro cell-based and murine models suggest that GARS and YARS (Tyrosyl tRNA synthetase) mutations are the cause of the diseases through a toxic gain of function [24]. Experiments suggest that the mutations do not cause the disease because of the loss of amino-acylation. In addition, it is also not clear whether the mutations have a neuronspecific effect on glycine. It was hypothesized that the axonal transport of mutant GARS may be damaged, thus leading to the degeneration of the distal axon [25]. This is supported by the immunohistochemistry of normal human tissue where the GARS containing puncta in spinal cord, nerve root, and peripheral nerve sections suggest that the GARS protein may be transported to the distal axon [8]. With these, the impaired axonal transport of GARS, mitochondrial dysfunctions, and loss of a non-canonical function have been the suggested mechanisms for the manifestation of HMN5A [25].
Mutations that change the protein sequence may affect its interaction with other proteins by changing its structural conformation and stability [26]. Changes in protein-protein interactions may lead to aberrant cellular functions [27]. These interactions usually occur through complex formation which can be modeled through protein-protein docking [28,29].
Understanding protein interactions and computationally predicting the functional and structural effects of a mutation involved can provide insight on the pathogenesis of CMT2D and HMN5A.

Methods
Identification of proteins that interact with GARS Cytoscape 3.4.0 [30] was used to retrieve the network of all proteins that are known to interact with GARS. Information in this network came from ten (10) [36], InnateDB-All [37], Reactome-FIs [38,39], IntAct [40], mentha [41], and Database of Interacting Proteins (DIP) [42,43]. The resulting network was filtered for proteins that only occur in humans. Proteins that are validated across at least two different databases were included in the study. In addition, only the proteins with available experimentally determined structures that exist freely or in complex were considered. Most of these proteins are smaller than GARS and from here on are referred to as "ligands". These ligands were used as basis for the elucidation of gene function and structural changes in protein-protein interactions.

Docking of proteins
The respective structures of GARS and of the ligands included in the initial network were obtained from the RCSB Protein Data Bank (PDB). The rigid body docking of the ligands with wild type GARS and mutant GARS were done using GRAMM-X [44,45]. GRAMM-X was selected after systematic comparison of its performance with three other docking software (manuscript in prep). The results from the docking procedure of wild type GARS are referred to as "GARS wt complexes" and that of mutant GARS as "GARS m ut complexes." All the top scoring solutions from GRAMM-X were analyzed.

Analysis of Docking Results
The root mean square deviation (RMSD) values of the resulting solutions from the docking of the wild type GARS were compared with that of the solutions from the docking of the mutant GARS. The binding free energy (ΔG) in kcal/mol of each solution was computed using PRODIGY (Protein binding energy prediction) [46]. The change in the binding free energy (ΔΔG) was then computed using the equation: ΔΔG = ΔG Mutant-ΔG Normal [47].
The interfaces of the native and mutant complexes were compared by analyzing DIMPLOT [48] results for the interfaces. For all complexes, the number of residues involved in hydrophobic contact and the number of H-bonding interactions were determined, as well as the average length of the H-bonds.

Results
The network generated using Cytoscape 3.4.0 is shown in Figure 1. It has 398 nodes and 451 edges connecting to GARS. Of these, 85 proteins interact with GARS and are validated across at least two of the ten (10) databases that were searched. The experimentally determined structures of 31 of the 85 proteins are available in the Protein Databank (PDB). They were the ones investigated in this study and are in Table 1.
The results from the protein-protein docking of wild type GARS and mutant GARS with the ligands are in Figure 2. The best solutions for GARS wt complexes, in blue, are superimposed with the best solutions for GARS m ut complexes, in red.
Analysis with PRODIGY showed that upon mutation, 19 protein-ligand complexes became more stable, ten (10) decreased in stability and two were unchanged energetically based on the change in binding free energy (Table 2). Most stabilized protein complexes had an increased number of residues involved in hydrophobic contact, number of H-bonding interactions, and decreased average length of H-bonds (Table 3). On the other hand, most destabilized complexes had a decreased number of residues involved in hydrophobic contact and number of H-bonding interactions, and increased average length of H-bonds (Table 4). The complexes with catenin alpha 1 (CTNNA1) and myotrophin (MTPN) remained unchanged energetically and showed properties that were similar to stabilized and destabilized complexes, respectively (Table 5).
These results are expected, as stronger protein-ligand binding leads to a more stable complex, while weaker binding results in a less stable complex. However, not all proteinligand interactions showed this. These peculiarities in binding properties and ΔΔG values can be attributed to the limitations of the programs used. These limitations are discussed in a separate manuscript (in prep)..

Apoptosis
Neuronal cell death is characteristic of many neurodegenerative diseases such as Amyotrophic Lateral Sclerosis (ALS), Huntington's disease (HD), Alzheimer's disease (AD) and Parkinson's disease (PD) [49]. It is known that neurons have minimal regenerative potential [50]. Since most of the neurodegenerative diseases including CMT2D and HMN5A progress slowly, it is likely that although the individual neurons may be dysfunctional for extended periods, rapid death of neurons may occur once the apoptotic cascade is fully activated [51]. Thus, the manifestation/development of the neurodegenerative diseases may be due to the progressive death of the individual neurons as well as its limited regeneration capacity.
The ligands Proteasome 26S Subunit, ATPase 3 (PSMC3) and Mitogen-Activated Protein Kinase 8 (MAPK8) have functions relating to apoptosis [48]. PSMC3 is associated with apoptosis through its role in the degradation of ubiquinated proteins [48]. However, the role of the ubiquitin-proteasome system in relation to apoptosis in neurons remains to be established.
MAPK8 belongs to the c-JUN-N-terminal kinases (JNK) family [52]. The JNKs are crucial for neuronal apoptosis as well as neurite growth, and the level of JNK activity is correlated with neurite length [53]. Thus, the inhibition of JNK will result in the inhibition of neurite growth and decreased neurite length. Studies using Jnk-knockout mice have revealed additional roles for the JNKs in brain morphogenesis, neuronal pathfinding, axodendritic architecture maintenance, and neuronal death after excitotoxic insults [54,55].
Autophagy is important in removing accumulations of aggregated and ubiquinated proteins. These often result in neurodegenerative diseases such as AD, PD, HD, ALS, static encephalopathy of childhood with neurodegeneration in adulthood (SENDA), hereditary spastic paraplegia (HSP) and Lafora disease. Therefore, impairment of autophagy can lead to disease [57,58]. Puncta formation and localization in axons has been observed but it is expected in wild type GARS. The puncta formation is even inhibited by mutations in GARS [8]. CMT2D and HMN5A may not be due to aggregations of misfolded proteins. However, autophagy has also been found to be required in the maintenance of homeostasis in axons and its loss may cause axonal dystrophy [59]. Both CMT2D and HMN5A are characterized by axonal peripheral neuropathies [9,12]. GARS is also expressed in the mitochondrial matrix [9]. As mentioned, GABARAPL2 also play a role in mitophagy. Drosophila studies have shown that peripheral neuropathies can result from aberrant mitochondrial dynamics such as the organelle's enlargement and mislocalization [60]. It was also suggested that oxidative stress and mitochondrial dysfunction in relation to autophagy may be contributing factors to the development of neurodegenerative diseases [61].

Ubiquitin Proteasome System
It is notable that many of the ligands that were considered have functions relating to The ubiquitin-proteasome system (UPS) has been of particular interest in neurodegenerative disorders because of its role in the removal of misfolded proteins [65].
Since the accumulation of misfolded proteins is absent in CMT2D and HMN5A, it is unlikely that the disease is caused by defects in the UPS that prevent it from removing such proteins. The UPS has been shown to have a role in axonal growth and axonal guidance [66]. It has also been linked to axonal degeneration in mice and Drosophila studies [67].
The studies suggest that the inhibition of the UPS is generally harmful to the cell, but the molecules involved in its role in the axons are yet to be determined [67]. A deficiency of the E3 ubiquitin ligase TRIM2 has been shown to cause axonal neuropathy following accumulation of neuronal filaments that failed to be ubiquinated [68]. If the UPS is involved in the manifestation of CMT2D and HMN5A, it is likely that the defects in the system affect the maintenance or development of axons, which result in axonal neuropathy.
Nerve Growth Factor/Tropomyosin receptor kinase A (NGF/TrkA) Signaling Pathway One interesting interaction shown to have been affected by mutation in GARS in this study is that with Neurotrophic Receptor Tyrosine Kinase 1 (NTRK1). NTRK1 belongs to the Trk family of receptor tyrosine kinases that binds with the neurotrophin, NGF (neuronal growth factor). Neurotrophins can bind to two types of receptors, Trk tyrosine kinase receptors or to p75 neurotrophin receptor (p75NTR). These two types of receptor can enhance or inhibit the effects of the other to facilitate the actions of neurotrophins [69]. The binding of neurotrophins to Trk receptors can activate the phosphatidylinositol 3-kinase/Akt protein kinase signaling pathway. It is essential for the survival of neurons. It can also activate the extracellular-signal-regulated kinase (ERK) pathway, which is a neuroprotective mechanism [70]. In summary, neurotrophins that bind to Trk-receptors activate signal transduction pathways that ultimately lead to the activation of gene expression, neuronal survival and neural outgrowth [71].
The importance of the NGF/TrkA signaling in both the peripheral and central nervous systems was known when mice carrying a defective receptor gene developed severe sensory and sympathetic neuropathies and died within one month of birth [72]. This was confirmed by another study which showed that the NGF/TrkA signaling is important in the development of correct peripheral target innervation and in full biochemical differentiation of sensory neurons [73]. If NGF/TrkA signaling is affected upon mutation in GARS, the sensory deficits experienced by patients of CMT2D and HMN5A could be explained by atypical differentiation of peripheral sensory neurons during embryonic development. However, symptoms of CMT2D and HMN5A manifest at ages 16 to 30 and 5 to 25, respectively, and not at birth. This suggests that their sensory neurons are normal until they experience axonal degeneration at the said ages. It has been shown that in adult neurons, the regeneration of axons of sensory neurons is mediated by cytokine signaling pathways rather than NGF/TrkA signaling pathways [74].
NGF is derived from its precursor, proNGF, which plays a role in apoptosis upon binding with a complex of p75NTR and sortilin. The levels of proNGF are shown to be upregulated in the central and peripheral nervous system in old age [75]. Moreover, a study showed an increase in proNGF levels in the peripheral nervous system of old mice where it plays a role in cell death of sympathetic neurons [75]. Furthermore, brains of Alzheimer disease patients show elevated levels of proNGF, which, upon extraction, caused the death of sympathetic neurons [75,76]. It was also found that proNGF is more abundant in sympathetic and peripheral sensory ganglia as well as in their target tissues than mature NGF [77]. TrkA and p57NTR have different affinities for NGF and proNGF. TrkA has a higher affinity for NGF, while p75NTR has a higher affinity for proNGF. In the presence of both receptors and even with sortilin, NGF and proNGF promote cell survival. However, a reduction in TrkA results in proNGF signaling for cell death and loss of activity for NGF.
When there is no TrkA, proNGF still signals for cell death even when p75NTR and sortilin levels are reduced [78]. If the interaction between mutant GARS and NTRK1 results in a reduction or loss in function of NTRK1 or its failure to be transported to neurites, the proNGF present in peripheral neurons and their target tissues can induce apoptosis in the neurons. This could help explain the neuropathy observed in CMT2D and HMN5A patients.
However, this cannot explain the specificity of the neuropathy to peripheral nerves.
Changes in interaction between GARS and NTRK1 can possibly affect the NGF/TrkA signaling pathway which, in turn, can affect its functions in the maintenance and development of neurons. However, the biological processes affected by the pathway are numerous and diverse, therefore, much more experimental evidence is needed before its implication on the manifestation of CMT2D and HMN5A can be inferred.

Immune Response
It is notable that many of the protein-ligand interactions affected by the mutation in GARS KARS signals for inflammatory response by activating monocytes or macrophages. Its secretion from uninjured cells is induced by tumor necrosis factor alpha (TNF-α). KARS then binds to the macrophages and induces the increased production and migration of TNF-α [79]. Macrophages, in turn, play a role in Wallerian degeneration and axonal regeneration in peripheral nerves upon injury. They promote nerve regeneration possibly through several mechanisms such as stimulating the extracellular matrix for regeneration or regulating Schwann cell activities [80]. KARS is also the fourth protein belonging to the ARS family to be associated with CMT, providing further evidence of the importance of the ARS family in axon function [81]. This suggests that disease manifestation upon mutation of GARS and KARS share a mechanism. The inhibition of the macrophage-activating function of KARS could possibly prevent axonal regeneration which, when acting in conjunction with other factors that destroy axons in these mutations, may result in peripheral neuropathies.

Conclusions
We found that the E71G mutation in GARS affected the stability and binding properties of the complexes it formed with other proteins. This mutation can possibly lead to disease manifestation by affecting the biological processes of apoptosis and autophagy, the pathways of the ubiquitin/proteasome system and NGF/TrkA signaling, as well as immune responses found to be significant in axonal regeneration. The apoptotic process is related to CMT2D and HMN5A through MAPK8, to autophagy through GABARAPL2 and MAPK8, to the ubiquitin/proteases system through UFD1L, PSMC3, CANAD1, CUL5, OTUB1, UBE2N, PSMC3 and TRAF6, to NGF/TrkA signaling through NTRK1, and to the immune response through KARS. It is suspected that disease manifestation is caused by an interplay of the disruptions among these.   Table 3. Summary of DIMPLOT results for interfaces of stabilized complexes.  Figure 1 Network of all proteins interacting with GARS generated using Cytoscape 3.4.0.

List Of Abbreviations
The protein interactions that are validated across at least two of the ten databases that were searched are in yellow.