Structural Evolution of Primate Glutamate Dehydrogenase 2 as Revealed by in Silico Predictions and Experimentally Determined Structures

1. Introduction

Glutamate dehydrogenase (GDH) reversibly interconverts glutamate to a-ketoglutarate and ammonia using NAD(P)+ as cofactors [1]. The enzyme interconnects carbon and nitrogen metabolism and is found in almost all living organisms[2]. In eukaryotes, GDH is abundantly expressed in mitochondrial matrix where it contributes to energy homeostasis. Specifically, α-ketoglutarate, produced via oxidative deamination of glutamate, feeds the Krebs Cycle, serving anaplerotic functions and leading to ATP synthesis[2].

In addition to the GLUD1 gene (encoding for hGDH1), humans possess GLUD2 (encoding for hGDH2), an intronless X-linked gene thought to have evolved through retroposition of a spliced GLUD1 mRNA (retroposon) [3]. Subsequent phylogenetic studies revealed that the retrotransposition of the GLUD1 gene to the X chromosome occurred during primate evolution more than 23 million years ago. [4]. After emerging in the common ancestor of humans and other modern apes, GLUD2 underwent rapid evolutionary adaptation concurrently with brain evolution [4]. This adaptation involved 15 amino acid substitutions in the mature hGDH2 that provided unique functional properties [5].

Both GLUD genes encode for a 558-amino-acid long polypeptide sequence. The first 53 amino acids located on the N-terminus domain correspond to the leader peptide, which is responsible for the transportation of the enzyme inside the mitochondrial matrix [6,7,8]. The mature hGDH1 and hGDH2 isoenzymes, resulting from cleavage of the leader peptide inside the mitochondrion, share all but 15 of their 505 amino acids [3]. Despite its sequence similarity to hGDH1, hGDH2 has unique enzymatic and regulatory properties, including GTP resistance, relatively low basal activity markedly responsive to activation by ADP and/or L-leucine, lower optimal pH and relative sensitivity to thermal inactivation [5]). As shown by enzymatic studies, these highly divergent properties are to a large extent related to two (Arg443Ala and Gly456Ala) of the 15 amino acid substitutions that occurred during hGDH2 evolution [9,10,11].

In addition to distinct enzymatic and regulatory properties, hGDH2 displays a unique expression pattern. hGDH1 is encoded by the housekeeping GLUD1 gene and is expressed in all human tissues, with the highest levels found in the liver. Gain of function amino acid changes lead to the hyperinsulinism hyperammonemia syndrome [12], a serious metabolic disorder with childhood-onset. On the other hand, hGDH2 is expressed mainly in the human brain, kidney, testis, steroidogenic organs and shows low expression levels in the human liver [13]. Recently, the possibility has emerged that hGDH2 is involved in the pathogenesis of neurodegenerative and neoplastic disorders[14,15].

The 3D-structure of hGDH1 [16], and the structures of several other mammalian and non-mammalian GDH1s [17] have been determined by X-ray crystallography. The mammalian GDH1 structure is a symmetric homo-hexamer, with each subunit consisting of the N-terminal glutamate-binding domain, the NAD+-binding domain, the antenna, the pivot helix and the C-terminal helices. Recently, we have determined the crystal structure of the hGDH2 protein at 2.9 Å resolution, showcasing important differences compared to hGDH1 [18]. However, no information is available about the structure of GDH2 in modern apes, other than humans, or in their, now extinct, common ancestors (nodes B, C, D, and E in Figure 1). Importantly, the modern hGDH1 corresponds to the original GDH sequence present 23 million years ago, from which the line that led to the modern primate GDH2 emerged. This evolutionary conservation shows that GDH1 is a crucial metabolic enzyme with little tolerance for changes.

Advanced protein structure prediction algorithms have been recently developed to supplant the experimentally determined protein structures. One such algorithm is AphaFold, developed by DeepMind, which uses artificial intelligence to accurately predict protein structures from their amino acid sequence [19,20]. Here, we examine the accuracy of the predicted hGDH1 and hGDH2 AphaFold Colab models by comparing them with the experimentally determined human enzyme structures. Furthermore, AlphaFold Colab is used to predict the structures of GDH2 of modern apes and their ancestors, going back 23 million evolutionary years. Finally, we present the effects of the amino acid substitutions that occurred at each evolutionary step.

2. Materials and Methods

2.1. Phylogenetic tree analysis

The phylogenetic tree, based on the GLUD2 sequences encoding the mature polypeptide, was constructed by the Molecular Evolutionary Genetics Analysis (MEGA) program [21] using the neighbor-joining method (Figure 1). On the branches of this tree, the amino acid substitutions that led to the emergence of current GDH2 proteins in great apes are depicted.

2.2. Protein structural prediction and analysis

The experimental crystallographic structure of hGDH1 and hGDH2 was retrieved from the Protein Data Bank (pdb code “1L1F” and “6G2U”, respectively). AlphaFold colab (https://colab.research.google.com/github/sokrypton/ColabFold/blob/main/AlphaFold2.ipynb ) was used to predict the structures of GDH2 proteins in modern day great apes and in extinct primates. This server predicts protein structure from their amino-acid sequence, using a simplified version of AlphaFold v2.0 that does not require homologous structures (templates).

The best five models were selected according to the ranking by the predicted local-distance difference test (pLDDT) confidence values (higher=better to lower=worse). The AlphaFold pLDDT scores for the proteins studied are shown in Table 1. The resulting models were examined, aligned, and compared to each other and to the experimentally determined structures using PyMOL. The command “super” was used comparing protein backbones. To evaluate the differences between predicted or experimentally determined structures, we used the root-mean-square deviation (RMSD) values resulting from the alignments. An RMSD value below <1.8 Å was considered as suggestive of high accuracy.

2.3. Mutational analysis

Mutant GDH2 stability for each evolutionary step was estimated by changes in free energies, ΔΔG (kcal/mol). The predicted structures of mutant GDH2s generated from AphaFold were used to perform the analysis. Five different web servers were used: the sequence-based iSTABLE (http://predictor.nchu.edu.tw/iStable/about.php) [22] and the structure-based PremPS (https://lilab.jysw.suda.edu.cn/research/PremPS/) [23], MaestroWEB (https://pbwww.services.came.sbg.ac.at/maestro/web/) [24], SDM (http://marid.bioc.cam.ac.uk/sdm2) [25] and DynaMut (http://biosig.unimelb.edu.au/dynamut/) [26].

The evaluation of structural stability of various GDH2s considered the results of the five methodologies to reach a majority consensus. In the framework of these methodologies, the application of iSTABLE, PremPS, MaetroWEB, SDM and DynaMut provided the estimate of the unfolding and total free energy as well as the vibrational entropy (Table 2). Differences in the results obtained by these servers when calculating the stabilizing/destabilizing impact of each amino acid change, are due to the use of different algorithms.

iSTABLE combines the results from different predictors such as iMUTANT and MUpro to determine the effect of point mutations on protein stability [22], by calculating the difference in folding energy change (ΔΔG in Kcal/mol) between the wild type and the mutant protein. DynaMut predicts changes in protein stability (ΔΔG in Kcal/mol), variation in entropy energy, changes in protein flexibility and allows visualization of non-covalent molecular interactions [26]. ΔΔG > 0 corresponds to stabilizing effect whereas ΔΔG < 0 to destabilizing effect. PremPS predicts changes in protein stability (ΔΔG in Kcal/mol) as well as the location of the mutation, either in the hydrophobic core or on the protein surface [23]. Positive ΔΔG values indicate a destabilizing effect on protein stability whereas negative ΔΔG a stabilizing effect. The Multi AgEnt STability pRedictiOn (MAESTRO) webserver estimates the changes in unfolding free energy upon point mutation through a machine learning system [24]. ΔΔG > 0 corresponds to destabilizing effect whereas ΔΔG < 0 to stabilizing effect. Site Directed Mutator (SDM) uses environment-specific amino-acid substitution frequencies within homologous protein families to calculate a stability score, which is analogous to the free energy difference between the wild-type and mutant protein[25]. Positive ΔΔG values indicate a stabilizing effect whereas negative ΔΔG values a destabilizing effect.

3. Results

3.1. AlphaFold predicted Versus experimentally determined hGDH1 and hGDH2 structure

The AlphaFold provided a satisfactory prediction of the experimental 3D-structures of the hGDH1 and hGDH2 protein (Figure 2). The predicted protein structures presented in this study, show all the important domains found in each subunit of the hexameric glutamate dehydrogenases. These domains include a glutamate binding region towards the N terminus, a NAD binding domain, and a regulatory domain consisting of the antenna and the pivot helix.

We initially explored whether the hGDH1 and hGDH2 structures predicted from their sequences using AlphaFold Colab were accurate. To answer this, the predicted hGDH1 structure derived from AlphaFold Colab and the experimentally determined hGDH1 structure (PDB entry 1L1F) were superimposed (at a total of 3,418 atoms) using PyMOL. The RMSD value between the two superimposed structures was estimated to be 1.745 Å (Figure 2a). Similarly, the predicted hGDH2 structure derived from AlphaFold Colab and the experimental hGDH2 structure (PDB entry 6G2U) were superimposed using PyMOL, at a total of 3,278 atoms. The RMSD between the AlphaFold Colab structure and the experimental template was 0.895Å (Figure 2b). Thus, the comparisons of the AlphaFold predicted structures with the experimentally determined, reported in PDB, structures highlight the ability of this approach to adequately predict the structures of the individual domains.

Also, comparison of the AlphaFold derived structures corresponding to the nodes B, C, D and E in Figure 1 with the experimentally determined hGDH1 and hGDH2 structures gave comparable results with these described above (Figure 3 and Figure 4). Specifically, the comparison of experimental hGDH1 with proteins predicted for nodes B, C, D and E gave RMSD values of 1.695 Å, 1.714 Å, 1.766 Å and 1.680 Å, respectively (Figure 3). For hGDH2, these values were calculated to be 0.905Å, 0.896Å, 0.943Å and 0.895Å, respectively (Figure 4).

3.2. hGDH2 AlphaFold Colab predicted structures during evolution that led to humans

Based on the phylogenetic tree of the primates and the gene sequence in different modern-day species, it was found that the retrotransposition event, that led to the emerge of the GLUD2 gene, occurred after the separation of the phylogenetic branches of the great apes of the Old World and the African green monkey, almost 23 billion years ago [4]. The human GLUD1 gene that encodes for hGDH1 has remained unchanged for the last 23 million years. This indicates that it is an ortholog of and essentially identical to the original GDH gene (node A, Figure 1) in the common ancestor of modern great apes which gave rise to the GLUD2 gene through retrotransposition. Thus, we have good reason to support that the experimentally determined hGDH1 structure corresponds to that of the common ancestral enzyme.

In the common ancestor of humans and modern apes, seven amino acid substitutions occurred during the first evolutionary step following the retrotransposition event (node B, Figure 1). These were Ala3Val, Glu34Lys, Asp142Glu, Ser174Asn, Arg443Ser, Gly456Ala and Asn498Ser. During the second evolutionary step, after the separation of the gibbon branch, six amino acid substitutions (Val3Leu, Arg39Gln, Lys299Arg, Ser331Thr, Met370Leu, Arg470His) appeared (node C, Figure 1). Finally, on the last two steps (nodes D and E, Figure 1) three (Ile166Val, Gly247Arg, Ala321Val) and two (Arg299Lys, Met415Leu) substitutions, respectively, led to the current hGDH2 protein in humans.

The GDH predicted structures corresponding to node A and node B were superimposed (at a total of 3,171 atoms) using PyMOL and the RMSD value between the two models was 0.112Å (Figure 1, Figure 5). Similarly, the GDH2 node B and node C predicted structures were superimposed (at 3,212 atoms) and the RMSD value was 0.102Å (Figure 1, Figure 5). The RMSD value between the superimposed GDH2 node C and node D structures (at 3,088 atoms) as well as the GDH2 D and node E structures (at 3,1997 atoms) were 0.095Å and 0.080Å, respectively (Figure 1, Figure 5). Finally, GDH2 node A and node E predicted structures were superimposed (at 3,088 atoms) and the RMSD value was 0.122Å (Figure 1, Figure 5).

3.3. Great ape GDH2 AlphaFold Colab predicted structures and comparison with predicted hGDH2

The predicted structure models for each great ape (chimpanzee, gorilla, gibbon, orangutan) and the predicted structure model for hGHD2 were superimposed using PyMOL. The RMSD value between the chimpanzee predicted structure and the hGDH2 predicted structure was 0.078Å (3,100 atoms), whereas the RMSD value between the gorilla predicted structure and the hGDH2 predicted structure was 0.083Å (3,190 atoms). Correspondingly, the RMSD values from the superimposition of the gibbon predicted structures and the orangutan predicted structures with hGDH2 predicted structure were 0.136Å (3,078 atoms) and 0.129Å (3,155 atoms), respectively (Figure 1, Figure 6).

Figure 6. Superposition of the model structure corresponding to hGHD2 (blue) with every predicted structure model for each ape(green). (a) chimpanzee GDH2- hGHD2. The RMSD value was estimated to be 0.078Å. (b) gorilla- hGHD2. The RMSD value was estimated to be 0.083Å. (c) orangutan- hGHD2. The RMSD value was estimated to be 0.147Å. (d) Gibbon- hGHD2. The RMSD value was estimated to be 0.136Å.

Figure 6. Superposition of the model structure corresponding to hGHD2 (blue) with every predicted structure model for each ape(green). (a) chimpanzee GDH2- hGHD2. The RMSD value was estimated to be 0.078Å. (b) gorilla- hGHD2. The RMSD value was estimated to be 0.083Å. (c) orangutan- hGHD2. The RMSD value was estimated to be 0.147Å. (d) Gibbon- hGHD2. The RMSD value was estimated to be 0.136Å.

3.4. Great ape GDH2 AlphaFold Colab predicted structures during evolution

During great apes’ evolution, after the separation of the gibbon branch (node B, Figure 1) three substitutions (Thr101Ala, Leu377Val, Glu8Lys) emerged and led to the establishment of the current gibbon GDH2 protein. Similarly, five amino acid substitutions (Ile239Asn, Leu240Val, Ile275Val, Leu375Val, Gln441Arg) appeared after the separation of the orangutan branch (Node C, Figure 1) that led to emerge of the current orangutan GDH2 enzyme. The establishment of the gorilla and chimpanzee protein was due to the emergence of four (Ser66Cys, Lys362Arg, Leu365Gln, Glu439Asp) and two (Ile305Leu, Val321Ile) amino acid substitutions, respectively, after the separation of their phylogenetic branches (node D and E, respectively, Figure 1).

The model structure corresponding to the common ape ancestor was superimposed, using PyMOL, with every predicted structure model for each ape (Figure 7). The RMSD value between the chimpanzee protein and the common ancestor protein was 0.151Å (3,106 atoms), while the RMSD value between the gorilla predicted structure and the ancestor predicted structure was 0.139Å (3,180 atoms). Similarly, RMSD values from the superposition of the gibbon protein and the orangutan protein with the common ancestor predicted structure were 0.109Å (3,289 atoms) and 0.147Å (3,280 atoms), respectively. These results are comparable to the same calculations for modern day hGDH2 (0.122Å).

Figure 7. Superposition of the model structure corresponding to the common ape ancestor (blue) with every predicted structure model for each ape(green). (a) chimpanzee GDH2-common ape ancestor. The RMSD value was estimated to be 0.151Å. (b) gorilla-common ape ancestor. ). The RMSD value was estimated to be 0.139Å. (c) orangutan-common ape ancestor. The RMSD value was estimated to be 0.129Å. (d) Gibbon-common ape ancestor. The RMSD value was estimated to be 0.109Å.

Figure 7. Superposition of the model structure corresponding to the common ape ancestor (blue) with every predicted structure model for each ape(green). (a) chimpanzee GDH2-common ape ancestor. The RMSD value was estimated to be 0.151Å. (b) gorilla-common ape ancestor. ). The RMSD value was estimated to be 0.139Å. (c) orangutan-common ape ancestor. The RMSD value was estimated to be 0.129Å. (d) Gibbon-common ape ancestor. The RMSD value was estimated to be 0.109Å.

3.4. Mutational and intramolecular interactions analysis

Totally 18 evolutionary amino acid substitutions (with 15 of them still present in modern humans) were analyzed to predict the result of each amino acid substitution during hGDH2 evolution (nodes A to E, Figure 1). The evaluation of the effect of the amino acid substitutions on protein stability by the consensus indicated that 50% of the mutated sites generated a stabilizing effect and 50% a destabilizing effect (Table 2 and Table S1, Figure S1, Figure 8, 9 and 10). Since our findings revealed that the amino acid substitutions occurring during great apes’ evolution are altering the free energy and the dynamicity of the enzyme, we aimed to investigate the impact of these amino acid replacements on the intramolecular interactions (Table 3). Structure-based analysis by DynaMut, using the hGDH1 structure as template, revealed that the amino acid substitutions were significantly affecting these intramolecular interactions (Table 3).

In specific, during the separation of the phylogenetic branches of the Old-World apes and the African green monkey seven amino acid changes emerged (Node A, Figure 1). Ala3Val, Asp142Glu, Ser174Asn and Gly456Ala increased protein stability based on the consensus of methods (Table 2, Figure 9). On the other hand, Glu34Lys, Arg443Ser and Asn498Ser decreased protein stability (Table 2, Figure 9). The Ala3Val substitution led to loss of a bond with Ser1 and an interaction with Ala5. The substitution of the negatively charged Glu34 by Lys led to the establishment of new interatomic interactions. The Asp142Glu substitution led to significant changes in the interatomic interactions as old were lost, and new interactions were observed. New interactions with Tyr99 and Pro137 were observed in the Ser174Asn substitution. Interactions with Glu439 and Phe440, and Ala447 were lost in the Arg443Ser substitution, while new ones were observed with Phe440 and Ala447. The Gly456Ala substitution led to significant changes in the interatomic interactions, as two bonds with His454 and Tyr459 and one with Thr460 were lost. On the other hand, 2 new hydrophobic interactions with Phe387, 3 bonds with Val453, His454 and Tys459, and one with Ile452 were observed. Finally, the Asn498Ser substitution led to the loss of the hydrophobic interactions with Phe494 and Ile52. On the other hand, new interactions were observed with Gly501, Ala500 and Phe494, respectively.

Five of the six amino acid substitutions that occurred after the separation of the gibbon phylogenetic branch (node B, Figure 1; Table 2) had a destabilizing effect on protein structure, with only Met370Leu increasing protein stability. No significant changes were observed in the interatomic interactions during the Val3Leu and Ser331Thr substitutions as only a bond with Ser1 was gained and an interaction with Gln334 was lost, respectively. Lysine 299 from α1 helix is able to make H-bonds and electrostatic interactions with residues form α2 helix and β1 strand (as it is depicted on Figure 8), therefore connecting all these elements together. The Lys299Arg substitution leads to even more interactions (Figure 8D) and a higher intraconnection of these secondary structure elements (Figure 8). Several hydrophobic interactions with Ile347, Phe230, Met237 were lost in the Met370Leu substitution, while new hydrophobic interactions were observed with Tyr236, Leu479 and Leu481. Also, new bonds were observed with Ile347. The Arg470His substitutions led to the loss of a bond with Met473 and Ala472. Finally, no changes in the intramolecular interactions were noted in the Arg39Gln mutant.

During the separation of the orangutan phylogenetic branch (node C, Figure1) three amino acid substitutions emerged. Ile166Val and Gly247Arg were found to destabilize the protein structure whereas Ala321Val was found to have an opposite stabilizing effect (Table 2; Figure 9; Figure 10). The Ile166Val substitution led to significant changes in the hydrophobic interactions as new bonds with Gly160 and Ile162 emerged. Also, bonds with Gly163 were observed. A bond with Ile318 was lost in the Ala321Val mutant, while several new hydrophobic interactions were observed with Tyr314, Ile318, Val252, Cys323. Additionally, two new bonds were observed with Cys323 and Lys344. Significant changes were observed in the interatomic interactions during the Gly247Arg substitution as new r bonds with Lys249 were gained.

Arg299Lys and Met415Leu that emerged during the separation of the Homo branch from the chimpanzee branch (node E, Figure 1), decreased and increased protein stability, respectively. On the other hand, new hydrophobic interactions were observed with Gln301 and Phe256. Finally, several interatomic interactions were lost in the Met415Leu mutant. On the other hand, new hydrophobic interactions were observed with Gln301 and Phe256. Figure 1. Table 2). Results showed that most of these substitutions were destabilizing.

Figure 1. Effect of Lysine 299 substitution by Arginine. Intramolecular interactions from position 299 (green GDH1/PDB ID 1L1F and cyan GDH2/PDB ID6G2U) when it is occupied by a lysine (panels A, B and C) and an arginine (modelled in 6G2U PDB file). Lys 299 lies on α1 helix and makes H-bonds with the carbonyl group of Phe 256 (α2 helix) and electrostatic interactions/H-bonds with Aso 279 ( β1 strand). When this position is occupied by an Arginine the number of possible interactions with the same elements increases.

Figure 1. Effect of Lysine 299 substitution by Arginine. Intramolecular interactions from position 299 (green GDH1/PDB ID 1L1F and cyan GDH2/PDB ID6G2U) when it is occupied by a lysine (panels A, B and C) and an arginine (modelled in 6G2U PDB file). Lys 299 lies on α1 helix and makes H-bonds with the carbonyl group of Phe 256 (α2 helix) and electrostatic interactions/H-bonds with Aso 279 ( β1 strand). When this position is occupied by an Arginine the number of possible interactions with the same elements increases.

4. Discussion

GLUD2 is a novel human gene that emerged though duplication in the hominoid ancestor (approximately 23 million years ago) [4] and underwent rapid evolutionary adaptation concurrently with primate brain evolution. The encoded human GDH2 (hGDH2) diverged substantially from its ancestor, the conserved hGDH1, in its functional, expressional and structural profile [9,10,11].

Although the 3D-structures of modern hGDH1 and hGDH2 have been experimentally determined using X-ray crystallography [16,18], the structural and functional properties of ancestral GDH2 enzymes that appeared during evolution are currently unknown. In this respect, we do not know if the primate GDH2 enzyme acquired its modern-day structural characteristics upon its emergence more than 23 million years ago or during subsequent evolutionary steps. In addition, due to the lack of experimental structure, is presently unclear whether hGDH2 differs from that of other modern primates.

To approach this question in terms of predicted structures, we used AlphaFold, a server that provides high quality 3D-structure predictions based on amino-acid sequence. It is widely accepted that AlphaFold predictions are accurate and often comparable to the experimentally determined structures, even though this is not absolute [19,27]. The AlphaFold algorithm uses Al-ML to predict 3D model structures across 21 proteomes of human (98.5% of the human proteins) and non-human (model organisms, agricultural crops and pathogens) organisms [27,28]. However, it has not been widely used for the delineation of structures of proteins of extinct species, as done here. Initially, we used the amino acid sequence of hGDH1 and hGDH2 as templates to obtain the AlphaFold predicted structures of these proteins and to compare them with the true experimentally determined hGDH1 and hGDH2 structures. These comparisons revealed that AlphaFold predictions were fairly accurate, thus highlighting the ability of this approach to adequately predict the structures of the individual domains.

Then, to gain insight into the evolutionary emergence of hGDH2, we compared the AlphaFold predicted structures of GDH2 of modern time apes and their now extinct ancestor. Our most important result using AlphaFold is that the predicted structure of GDH2 of the common ancestor of humans and other extant apes (chimpanzee, gorilla, orangutan, and gibbon) was the steppingstone for the structural and functional evolution of GDH2s in primates, with the first evolutionary step being associated with a higher RMSD value than subsequent steps. Indeed, judging by the RMSD values (Figure 1), the first evolutionary step was the most crucial one for the evolution of GDH2 in primates. In addition, we find that the gibbon GDH2 structure was more divergent from hGDH2 than those of other modern primates that are more closely related to humans.

It is common for proteins to acquire their most important properties upon emergence, otherwise they are doomed to become non-functional pseudogenes [29]. This initial evolutionary step before the separation of the human and gibbon lineages lasted about 5 million years and coincided with increasing functional properties of the primate brain [4]. Similar evolutionary processes to those of hGDH2 were in action for several other proteins that evolved during this period, such as opsins [30].

Also, given that we found differences in GDH2 structure between non-human primates and humans, this could form part of the diversification of brain function between primates. Humans are very close to chimpanzees at the genetic level (nucleotide difference of 1-2% at the level of the genome), making chimpanzees their closest living relatives [31] . However, despite the great nucleotide similarity, only 20% of the proteins are identical between the two species, even though research on this is ongoing [31]. An example of a protein that differs between humans and other apes is the digestive enzyme amylase. All vertebrates, including primates, express the enzyme in their pancreas. However, Old-World monkeys and humans, but not New World monkeys, express a-amylase additionally in their saliva. The ability to express a-amylase in saliva in Old-World monkeys, apes, and humans evolved after several duplications of the pancreatic amylase gene AMY2 within the primate lineage [32]

In addition, we examined separately the effect of the amino-acid substitutions emerging during apes evolution, using 5 different webservers, as done for other proteins, including spike glycoprotein of SARS-COV2 [33]. Of note, we cannot make comparisons between different amino acid substitutions since results vary across different stability prediction webservers. However, these properties are indicative of the general properties of the enzyme, given that our results are compatible with previous enzymatic studies [9,10,11].

It is worth mentioning that two important evolutionary changes, Arg443Ser and Gly456Ala, had a destabilizing and stabilizing effect, respectively. These findings corroborate previous enzymatic studies by us and others that Arg443Ser and Gly456Ala, which occurred in this first step, before the separation of the gibbon lineage, gave GDH2 its most important functional properties (low basal activity markedly activated by ADP and resistance to GTP inhibition, respectively [8,9,10]. There is evidence that these regulatory properties provided a novel role for hGDH2 in primate biology by enabling enzyme recruitment (through an ADP-dependent mechanism) under conditions of high energy utilization (increased conversion of ATP to ADP).

Even though these two changes (Arg443Ser and Gly456Ala) conferred most of the properties of the modern day hGDH2 enzyme they were not adequate to fully convert the ancestral enzyme to hGDH2. Indeed, similarly to the Arg443Ala single mutant, the double mutant of hGDH1 (Arg443Ala/Gly456Ala) was found to be essentially inactive, (basal activity <1%; little activation by physiologically relevant ADP concentrations), suggesting that additional evolutionary substitutions substantially modified the drastic effect of the Arg443Ala mutation thus providing the unique properties of hGDH2 [5]. It is of interest in this respect that three additional amino acid changes that occurred in the first step (Ala3Val, Asp142Glu and Ser174Asn) had a stabilizing effect (by consensus, Figure 9), probably rendering the enzyme active and ADP responsive (functional retroposon) and therefore contributing to its survival and subsequent evolution.

The results of this study contributed to the elucidation of the structure-function relationships of hGDH2 through the evolutionary lens. This is important as hGDH2 is involved in human physiology and pathophysiology. Specifically, there is accumulating evidence concerning hGDH2’s putative role in neurodegenerative processes, including Alzheimer’s disease and Parkinson’s disease, and tumorigenesis [13,14,15]. In this respect, we have known that tumors occur in primates in different organs and glands [34]. Also, it has been proven that Alzheimer’s disease is not human-specific but also affects non-human primates [35]. Thus, given its pathophysiological importance for severe human diseases, hGDH2 is becoming an attractive drug target and study of the structural evolution presented here could assist in rational drug design strategies.

Strengths and limitations

Of note, AlphaFold Colab can predict the structure only at the single subunit level and not of the functional hGDH hexamer. However, the predicted subunit structure includes all the important domains needed for glutamate dehydrogenase function: a glutamate binding region towards the N terminus, a NAD binding domain, and a regulatory domain consisting of the antenna and the pivot helix. It is known that the interactions between different subunits are present in the experimentally determined hexamer and these interactions influence the structure of each individual subunit [36]. However, when we compared the experimentally determined and the Alpha Fold predicted structure of hGDH1 and hGDH2 we did not detect significant deviations. In addition, since in our studies we compare the in silico predicted structures at the individual subunit level during evolution, the hexameric influences are not present in all predicted structures included here. Finally, given that the structural predictions of AlphaFold are similar but not identical to the experimentally determined structures, there is need to for experimental structural data that will verify, complement, and expand in silico AlphaFold produced data.

5. Conclusions

In summary, our most important results from AlphaFold structure predictions were that 1) GDH2 of modern-day apes is different from hGDH2 and 2) GDH2 in the common ancestor of humans and modern apes (node B in Figure 1) was the steppingstone for the structural and functional evolution of GDH2s in primates. Following this, primate GDHs underwent minor modifications that fine-tuned its enzymatic properties to adapt to the functional needs of modern-day primate nervous and other tissues. These results shed light on the structural/functional relationships of an enzyme that is important for human physiology and disease pathogenesis.