Recombinant protein production and purification

The GDAP1Δ303–358 and GDAP1Δ319–358 constructs, with an N-terminal His₆ tag and a Tobacco Etch Virus (TEV) protease digestion site, for producing soluble recombinant human GDAP1 in E. coli, have been described [11]. The point mutations R120Q, R161H, A247V, H256R, and R282H were generated in GDAP1Δ303–358, and the mutations R310Q and R310W in GDAP1Δ319–358, by a site-directed mutagenesis protocol with Pfu polymerase. All constructs were verified by DNA sequencing.

Recombinant GDAP1 variants were expressed in E. coli BL21(DE3) using autoinduction [28], and purified as described [11]. Briefly, GDAP1 was separated from the lysate by Ni²⁺-NTA chromatography, and the affinity tag was cleaved using TEV protease. Another Ni²⁺-NTA affinity step removed the tag and TEV protease. Size exclusion chromatography (SEC) was performed on a Superdex 75 10/300 GL increase column (Cytiva) using 25 mM HEPES (pH 7.5) and 300 mM NaCl (SEC buffer) as mobile phase. SEC peak fractions were analysed with SDS-PAGE and concentrated by centrifugal ultrafiltration.

X-ray crystallography

Mutant GDAP1Δ303–358 crystals were obtained using sitting-drop vapour diffusion at +4°C. Proteins were mixed with mother liquor on crystallisation plates using a Mosquito LCP nano-dispenser (TTP Labtech). The protein concentration was 10–30 mg/ml in 75 nl, and 150 nl of reservoir solution were added. R120Q crystals were obtained in 0.1 M succinic acid, 15% (w/v) PEG 3350. A247V crystals were obtained in 0.15 M DL-malic acid (pH 7.3), 20% (w/v) PEG3500. R282H crystals were obtained in 0.1 M succinic acid, 20% (w/v) PEG 3350. For cryoprotection, crystals were briefly soaked in a mixture containing 10% PEG200, 10% PEG400, and 30% glycerol, before flash cooling in liquid N₂.

X-ray diffraction data were collected at the PETRA III synchrotron (DESY, Hamburg, Germany) on the P11 beamline [29, 30] and the EMBL/DESY P13 beamline [31] at 100 K and processed using XDS [32]. The structure of wild-type GDAP1Δ303–358, PDB entry 7ALM [11], was used as the search model for molecular replacement in Phaser [33]. The models were refined using Phenix.Refine [34] and rebuilt using COOT [35]. The structures were validated with MolProbity [36]. The data processing and structure refinement statistics are in Table 1, and the refined coordinates and structure factors were deposited at the Protein Data Bank with entry codes 7B2G (R120Q), 8A4J (A247V), and 8A4K (R282H).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Data processing and refinement statistics.

Data in parentheses refer to the highest-resolution shell.

https://doi.org/10.1371/journal.pone.0284532.t001

Modelling

A model for full-length human GDAP1 was obtained from AlphaFold2 [37]. In addition, for an alternative model missing loops of the human wild-type GDAP1 crystal structure were built with CHARMM-GUI [38, 39]. Earlier structure-based bioinformatics results [12] were analysed further with respect to the mutational spectrum of GDAP1.

Synchrotron small-angle X-ray scattering

The structure and oligomeric state of the GDAP1 mutants were analysed with SEC-coupled small-angle X-ray scattering (SAXS). SEC-SAXS experiments were performed on the SWING beamline [40] (SOLEIL synchrotron, Saint Aubin, France). Samples were dialyzed against SEC buffer and centrifuged at >20000 g for 10 min at +4°C to remove aggregates. 100 μl of each protein sample at 1.7–36 mg/ml were injected onto a BioSEC3-300 column (Agilent), run at a 0.3 ml/min flow rate. SAXS data were collected at +15°C, over a q-range of 0.003–0.5 Å⁻¹. SAXS data analysis, processing, and modelling were done in ATSAS 3.0 [41]. Scattering curves were analysed and particle dimensions determined using PRIMUS [42] and GNOM [43]. Chain-like ab initio models were generated using GASBOR [44], dummy atom models were built with DAMMIN [45], and model fitting to data was analysed with CRYSOL [46].

To complement the SEC-SAXS, batch mode SAXS experiments were carried out for wild-type GDAP1 and four variants at the concentration range 1–4 mg/ml, to observe possible clear differences in oligomeric state at similar concentrations. These measurements were carried out on the P12 beamline at EMBL/DESY (Hamburg, Germany) [47].

Synchrotron radiation circular dichroism spectroscopy

Synchrotron radiation circular dichroism (SRCD) spectra were collected from 0.5 mg/ml samples on the AU-SRCD beamline at the ASTRID2 synchrotron (ISA, Aarhus, Denmark). The samples were prepared in 10 mM HEPES pH 7.5 and 100 mM NaF buffer, equilibrated to room temperature in 0.1 mm closed circular quartz cuvettes (Suprasil, Hellma Analytics). SRCD spectra were recorded from 170 nm to 280 nm at +25°C. Three scans per measurement were repeated and averaged. The spectra were processed using CDToolX [48].

Thermal stability

Thermal stability of GDAP1 variants was studied by nanoDSF using a Prometheus NT.48 instrument (NanoTemper), in SEC buffer. NanoDSF experiments were done at 0.5 mg/ml, and 3 aliquots of each sample were run simultaneously (technical replicates). Tryptophan fluorescence was excited at 280 nm, and emission was recorded at 330 and 350 nm, while the samples were heated from +20 to +90°C at a rate of 1°C/min. Changes in the fluorescence ratio (F₃₅₀/F₃₃₀) were used to determine apparent melting points. The sample analysis was done using Nanotemper PR.ThermControl analysis package utilizing Boltzmann fit analysis and average of three F₃₅₀/F₃₃₀ inflection points was taken. The curves presented in the figures are the averaged curves. The data were analyzed using Origin (OriginLab Corporation, Northampton, MA, USA).

Sequence entropy

The NCBI non-redundant protein database was used for all sequence searches. Starting from the human GDAP1 reference sequence (NP_061845.2), iterative PSI-BLAST [49] searches were performed, initially accepting sequences with an e-value of 10⁻⁹⁹ or less. In a second iteration, sequences with e-value < 10^-7 were accepted, resulting in 5986 sequences. This initial set was divided into three subsets according to their annotations: (a) GDAP1, (b) GDAP1L1 and (c) GST-labelled. For many calculations, GDAP1 and GDAP1L1 were combined to improve the statistics of conservation calculations, as they form a clear subgroup within GST-like sequences. Since the separation was based on annotations, sequences labelled as hypothetical or putative were removed. Sequence calculations were done either on subgroups or across the full set of 5065 sequences.

All multiple sequence alignments were calculated with MAFFT [50]. Results are always considered with respect to a reference sequence. This was done by removing any columns corresponding to a gap in the reference. This does lose information but is effectively essential for interpretation. For the GDAP1L1 group, the reference was the human sequence. For the larger GST group, the sequence of PDB entry 1PKZ was used to allow interpretation in structural terms.

Per-site entropy at site i in a multiple sequence alignment was given by where p_j is the probability of seeing amino acid type j at the site and the summation runs over the 20 amino acid types. Gaps were treated as missing data. The use of base 20 in the logarithm ensures that S ranges from 0 (fully conserved) to 1 (random). In the interpretation below, S<0.2 was considered conserved and S<0.1 as highly conserved.

Kullback-Leibler divergence

The most interesting sequence sites are those which are conserved within groups of proteins, but different between groups. The Kullbach-Leibler divergence (D^KL) captures the difference between two distributions. With discrete distributions (types of amino acids), the value at position i is given by where the summation runs over the 20 amino acid types. p_i,j and q_i,j are the frequencies (probabilities) of seeing residue type j at position i in the first and second protein groups [51, 52].

Phylogenetic analyses

Maximum likelihood phylogenetic trees were calculated with MrBayes [53, 54] with four independent chains of length 2 × 10⁶ steps, with 25% of the steps discarded for burn-in and fixed amino acid exchange rate matrices. The trees in results are all consensus/average trees displayed with Interactive Tree of Life (iTOL) [55] using a bacterial sequence (WP_173192278.1) to root the tree.

Sequence data sets for phylogeny were different to those used for conservation calculations, the aim being to cover sequences from GDAP1 to common GSTs and archaea. The first sequence searches and guide trees had identified the archaeal and bacterial sequences closest to GDAP1. The similarity matrix suggested that bacterial GSTs were remotely related to both the eukaryotic and archaeal sequences. The first 799 database hits of GDAP1 formed one group, and the first 1000 hits from human GDAP1L1 (NP_001243666.1) the next. The first 1000 hits from the bacterial GST (TNE50161.1) and the human GST (sequence from PDB entry 1PKZ) were the third group. Finally, the first 1000 hits for the archaeal sequence (MAE98075.1) formed a fourth group.

Each of these four sets were then reduced to 100 representatives based on the similarity matrix calculated during a sequence alignment [56]. The procedure sorts the distance matrix to get the closest sequence pairs and removes one member at a time, ensuring the most even possible spread of sequences. Combining the 4 × 100 sets and removing duplicates left 397 sequences used for phylogeny.

Helices α3, α6, α7, and α8 form a scaffold for GDAP1 intramolecular networks

The majority of CMT-linked missense mutations in GDAP1 are located within the vicinity of the hydrophobic clusters of the GST-like domains and the dimer interface [11], and the variants may induce changes in intramolecular hydrogen bonding networks [12]. In addition to the R120W and H123R studied earlier [12], we determined three new mutant crystal structures: R120Q, A247V, and R282H. These mutations reside in helices α3 (R120Q), α7 (A247V), and α8 (R282H), which are core elements of the GDAP1 fold. We shall first look at the central helices regarding GDAP1 folding.

The GST-like core fold of GDAP1 is supported by the α7 helix, surrounded by helices α3, α6, and α8. The helix α3 is connected to α6 via the α6-α7 loop, and Cys240 in this loop–itself being a CMT mutation site–is central to many interactions. The α6 helix can either extend or turn back towards the dimer interface, as seen in earlier crystal structures and models [12]. A new model built here indicates that the extended conformation is predictable (Fig 1C). Open/close movements of α6 can be functionally relevant for GDAP1 interactions with other proteins, such as cofilin and tubulin [9]. The α8 helix is positioned perpendicular to the others, and its orientation could be related to the transmembrane helix position, as it is expected to face the mitochondrial outer membrane surface. Together, these helices form an intramolecular network of polar (Fig 1D) and non-polar contacts, and many CMT-linked missense mutations are found on these helices (Fig 1E).

There are 15 designated hydrogen bond contacts or ionic interactions between helices α3, α6, α7, and α8. Helix α7, which is central to the GDAP1 fold, makes only a few hydrogen bonds to the surrounding helices, and hydrophobic residue clustering houses the α7 helix in the core of the fold. This is in line with earlier centrality analyses of the GDAP1 fold, showing helix α7 to be the most central part of the 3D structure [12]. Some of the disease mutations on helices α3, α6, and α8 correspond to solvent-exposed residues, and these helices show more polar interactions and higher flexibility compared to α7.

Structural effects of the individual CMT mutations

We have previously analysed the effects of R120W and H123R at the protein level, both residing on helix α3 [12]. Here, we extend the human GDAP1 crystal structure analyses to three more CMT mutations: R120Q, A247V, and R282H (Fig 2, Table 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Structural details for each point mutation studied by X-ray crystallography in this study.

A. Immediate environment of the R120Q mutation. Arg120 participates in a hydrogen bond network between helices α3 and α6. B. Effects of the A247V mutation. Ala247 on helix α7 is a central residue of the GDAP1 hydrophobic core. C. The R282H mutation. Arg282 interacts with the α6-α7 loop, and the mutation causes loss of these interactions. The structure of wild-type GDAP1 is shown in grey in all panels.

https://doi.org/10.1371/journal.pone.0284532.g002

The R120Q mutation is located close to the N-terminal end of the α3 helix. Due to the mutation, residue 120 loses contact with the α6-α7 loop tip at Cys240, whereby the hydrogen bond between the Arg120 side chain and the backbone carbonyl of Cys240 is lost (Fig 2A). This is further linked to alterations in surrounding side chain conformations.

The A247V mutation site resides near the N terminus of helix α7, tightly surrounded by helix α3 and the α6-α7 loop, and having van der Waals contacts to Val121, Tyr124, Cys240, and Thr245. Comparing to the wild-type protein, changes in the crystal structure are small, but an overall movement of surrounding protein segments is caused by the presence of Val in this position, due to the increased volume of the side chain. Both helices α3 and α6 move slightly away, without altering hydrogen-bonding patterns (Fig 2B).

R282H is located on helix α8, and the side chain of Arg282 points inwards in the wild-type GDAP1 structure; it is stacked against Trp238 and makes hydrogen bonds to the backbone carbonyl groups of residues 236 and 237 in the α6-α7 loop (Fig 2C). All these interactions are lost upon the mutation, and the His residue in the mutant is observed in a double conformation, making no side-chain hydrogen bond contacts.

Conformation and stability in solution

Earlier computational predictions suggested an overall destabilising effect of CMT-linked mutations in the GDAP1 protein [12]. Similarly, protein destabilisation was observed experimentally for the myelin protein P2 in the context of all identified CMT mutations [57, 58]. A comparative analysis of seven GDAP1 variants in solution was therefore carried out here, to support the crystal structures and other complementary data from the current and earlier studies, and to identify general trends linking CMT mutations and GDAP1 stability in vitro. While the 3D shape and dimensions were studied using SAXS, SRCD was employed to follow secondary structure content and nanoDSF to compare thermal stability.

The SAXS analysis of GDAP1Δ303–358 showed that the protein particle dimensions in solution correspond to a dimer, and the scattering profiles showed only minor shape differences, ruling out large-scale conformational differences or aggregation (Fig 3A, Table 2). All variants in this construct background eluted mainly as dimers, and the top of the peak was selected for the analyses. The largest deviation was observed for R161H, which had a larger R_g than the other variants, suggesting a more open structure (Fig 3B). Also R120Q had a slightly different conformation, most clearly seen in the distance distribution. The pair distance distribution function showed that the maximum particle dimension in all samples was ~90–110 Å, which corresponds to a dimer (Fig 3B). Ab initio modelling, here carried out on the R282H mutant, which was measured at the highest concentration, supported the presence of a dimeric GDAP1, corresponding to the conformation of the dimer we previously showed to fit the solution SAXS data for wild-type GDAP1 [12] (Fig 3C). As the SEC-SAXS experiments were carried out at the maximum concentration for each variant, we complemented these experiments with a set of batch-mode SAXS experiments with wild-type GDAP1 and four mutant variants (S1 Fig). All variants in this supplementary experiment are predominantly dimeric in the concentration range 1–4 mg/ml.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. SAXS analysis.

A. SAXS curves for the mutants in the GDAP1Δ303–358 construct. The curves have been displaced along the y axis for clarity. B. Distance distributions for the curves in A. R161H shows a more open conformation than the other variants. C. Top: Dummy atom model of R282H superimposed with a structure based on the collapsed conformation of wild-type GDAP1Δ295–358 [12]. Bottom: Dummy atom model of the R310Q mutant in GDAP1Δ319–358, superimposed on the same structure, indicating additional volume for the C-terminal extension. D. SAXS curves for the mutants in the GDAP1Δ319–358 construct. The dominant monomer peak from SEC was used for R310W and the dimer peak for R310Q. E. Distance distributions for the curves in D. F. Dimensionless Kratky plots for all constructs show similar levels of rigidity and globularity. The cross marks the location of the peak in a perfect globular particle.

https://doi.org/10.1371/journal.pone.0284532.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. SAXS parameters.

Data for wild-type GDAP1 as well as the monomeric mutation Y29EC88E are taken from [11]. The MW estimate corresponds to the Bayesian estimate from PRIMUS.

https://doi.org/10.1371/journal.pone.0284532.t002

For the two mutations affecting Arg310, we used the longer construct GDAP1Δ319–358, and both R310Q and R310W eluted as two peaks in SEC-SAXS, corresponding to a dimer and monomer (Fig 3D, 3E and Table 2). This behaviour has previously been linked to protein concentration, in that GDAP1 monomers can be observed at low concentration [11], and the Arg310 mutants were here studied at lower concentration than all the other mutants. R310Q mainly eluted as a dimer, while the highest peak for R310W was monomer; since the protein concentrations were similar, this result indicates that R310W may affect GDAP1 dimerisation and that the dimerisation equilibrium is rather slow. The SAXS data for the dominant peak for each mutant were analysed further. The R310Q dimer was more elongated than R282H of the shorter construct (Fig 3C), reflecting the presence of the C-terminal amphipathic domain in the construct. Comparing the Kratky plot for all variants studied here, it is evident that they all are rigidly folded, with little flexibility (Fig 3F).

For thermal stability and secondary structure analysis, nanoDSF and SRCD were performed. The SRCD spectra showed similar secondary structure composition for all variants and wild-type GDAP1, indicated by the fact that the spectral shapes were nearly identical, despite some fluctuation in amplitude. The latter could be caused by errors in concentration determination or aggregation status in the strong synchrotron UV beam, or by differences in dynamics of the secondary structure elements. The SRCD spectrum for R161H showed a slight change in the spectral shape between 205–225 nm, which could be linked to its more open conformation observed in SAXS above. It can be concluded that none of the studied mutations interfere strongly with the overall folding of GDAP1 (Fig 4A and 4B). The nanoDSF experiment (Table 3) revealed a ~1–12°C decrease in apparent melting temperature for all core domain mutants, compared to the wild-type GDAP1Δ303–358 (Fig 4C). This is similar to the previously studied R120W and H123R mutants, which had melting temperatures ~52–54°C, i.e. a drop of ~10°C compared to the wild-type protein [12]. Notably, the F₃₅₀/F₃₃₀ ratio varies between the variants at the starting point, especially so for H256R, which may be an indication of subtle differences in structure. All in all, the biophysical data indicate that the structural effects of all studied CMT variants in the core domain are local, and the mutations do not cause large-scale conformational changes at a level detectable with SRCD or SAXS. However, all studied mutations in the core domain caused a decrease in thermal stability, suggesting a breakdown of stabilising intramolecular interactions within the GDAP1 molecule. As the most dramatic example, the mutation A247V inside the folded core decreased the stability by 13°C, indicating the importance of a correctly packed hydrophobic core for GDAP1 stability.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Folding and stability of GDAP1 as affected by CMT mutations.

A. SRCD spectra for the mutants in the GDAP1Δ303–358 construct. The spectral shape is most different for R161H, and for R282H, the spectral amplitude is increased, but the shape does not change compared to wild-type GDAP1. B. SRCD spectra for the mutants in GDAP1Δ319–358. C-D. nanoDSF analysis for the mutants in GDAP1Δ303–358 (C) and GDAP1Δ319–358 (D). The curves shown are each average of three independent nanoDSF curves run in parallel. E. Non-reducing SDS-PAGE analysis of all studied variants. Top: variants in GDAP1Δ303–358; wt1Δ302 refers to the wild-type construct. Bottom: variants in GDAP1Δ319–358; wt1Δ319 refers to the wild-type construct. For the longer constructs (panel below), the added segment is most likely a membrane-binding motif, and the mutations may affect membrane interactions; this could explain the difference in electrophoretic mobility between the wild-type and mutant variants. The segment in question is likely to bind to SDS and form a helical structure instead of getting denatured. The uncropped gel images are in S2 Fig.

https://doi.org/10.1371/journal.pone.0284532.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. nanoDSF apparent melting points.

All values are average ± standard deviation from 3 replicates.

https://doi.org/10.1371/journal.pone.0284532.t003

For the mutations at position 310, it was observed that both R310Q and R310W in fact increased the stability of GDAP1Δ319–358 by ~1°C (Fig 4D). This can be explained by the fact that the segment carrying these mutations is no longer part of the core domain, but rather likely to represent an ⍺ helix attached to the membrane surface, possibly linking the membrane to the core domain. Overall, GDAP1 stability correlates with the dimer-monomer ratio observed in solution; the most destabilised mutants also show a higher fraction of monomeric protein on non-reducing SDS-PAGE (Fig 4E), even though in the 3D structure, they are far from the dimer interface. The exact mechanism for this observation is currently unknown, but it could mean that altered protein dynamics and local conformations affect the dimer formation and/or that dimerisation in turn stabilises the protein fold.

Implications of the point mutations in the context of the whole GDAP1 protein

The studied mutations cause subtle variations in hydrogen bond and van der Waals distances in nearby residues in the crystal structures, while there are no drastic structural differences, when the mutant structures are superposed onto the wild-type crystal structure. However, when comparing the variation of the residues in helices α3, α6, α7, and α8, in comparison to the wild-type Cα-atoms, a specific pattern arises. We analysed the residues participating a hydrogen-binding network in the mutant structures (Fig 5A), and the mapping shows that the residues pointing towards the dimer interface have only minor structural variations compared to the wild-type protein. In contrast, at the C-terminal end of helix α6, in all mutant crystal structures, the variation of the Cα atoms is high (Fig 5B). This suggests that flexibility of helix α6 could arise from altered intramolecular contacts in the vicinity. The C-terminal end of the α6 helix, the conformation of which apparently is affected by the mutations studied here, is itself a target for multiple CMT mutations, affecting Gln218, Va219, Glu222, Arg226, and Glu227 [18, 59–62].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Structural bioinformatics analysis of the crystallised CMT variants.

A. The C⍺ deviation of each mutant vs. wild-type GDAP1 structure. The highest deviations can be found in the residue range 220–230. Shown are only the residues participating in the hydrogen bonding network. B. Mapping the results onto the structure, it becomes evident that the C-terminal end of helix ⍺6 deviates the most from wild-type GDAP1 on average. C. Average predicted ΔΔG effect of GDAP1 mutations at CMT sites, as defined by CUPSAT. D. CUPSAT predictions for R120Q, A247V, and R282H being mutated into all possible amino acids. Note that a negative ΔΔG in CUPSAT means destabilisation, and that A247V is falsely predicted as stabilising. The mutations studied here are marked with spheres. Ala247, red; Arg120, blue; Arg282, green.

https://doi.org/10.1371/journal.pone.0284532.g005

We additionally investigated the more global mutational effects using bioinformatics tools. In our previous study [12], we analysed the wild-type GDAP1 core domain structure and the full-length AlphaFold2 coordinates with CUPSAT and MAESTRO [63, 64]. This provided predicted ΔΔG values and geometric properties; the results from CUPSAT analyses are further depicted in Fig 5C. On the average, most CMT mutation sites are predicted to cause destabilisation and in general, the few mutations predicted to be stabilising lie on the protein surface.

While the above analyses indicate an average effect that may destabilise GDAP1 structure, we looked at predictions for each mutation crystallised here in more detail (Fig 5D). In essence, all substitutions to Arg120 are unfavourable; this shows that the interactions made by Arg120 are important for folding and stability. This is in line with our experimental data for both R120Q and R120W, which indicate minor changes in structure, but destabilisation of the fold. For Ala247, a mutation into Val is predicted to be slightly stabilising, indicating that such a replacement, in a tightly confined pocket and with potential long-range effects, is difficult for the prediction algorithm. In the case of Arg282, mutation to His is indeed one of the most destabilising variants in the prediction, linked to a loss of several key interactions by the Arg282 side chain.

Phylogenetic analysis of GDAP1 and the GST superfamily

The ubiquitous GST superfamily is inherently very large and diverse. There are massive differences between families and even within the subfamilies. Starting a BLAST search with the reference sequence will quickly end up in GSTs. This means that GSTs will swamp other proteins in any naïve analysis. To avoid this, GDAP1, GDAP1L1 and GSTs were considered separately.

The first calculations focussed on GDAP1 (Fig 6). Intriguingly, based on large-scale sequence alignments, GDAP1 turns out to be more closely related to prokaryotic GSTs than eukaryotic proteins. In this situation, we must remember that sequence conservation across the GST superfamily is overall low. The phylogenetic calculations do not resolve to a clear bifurcating consensus tree, and the uncertainty is very high around the divergence of the GDAP1 group. One can say that GDAP1 is at least as far from generic eukaryotic GSTs as it is from prokaryotic proteins. The analysis provides little clues as to any enzymatic function of GDAP1, but it supports the findings that it has no conventional GST activity.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Phylogenetic analyses.

A. Colour scheme for phylogenetic trees. Note that not all organisms seen in the colour legend are present in every picture. This is because the trees are very large and to increase visibility parts of the tree will be collapsed (denoted by a grey circle), to hide underlying branches. B. The first branch of the phylogenetic tree. The tree was rooted at the root sequence. On the right side are eukaryotic GST sequences (mammals, fish, birds, reptiles). The branch includes 48 sequences. The rest of the 398 sequences that are represented are at hierarchically lower levels of the phylogenetic tree. The corresponding branches were collapsed (grey circle). C. In contrast to panel B, the left branch is partially extended. The right-sided branch shows eukaryotic GST sequences, while the left side shows a total of 110 bacterial GST sequences. The remaining sequences are at hierarchically lower levels of the phylogenetic tree. The corresponding branches were collapsed (grey circle). The red line indicates that the time could not be resolved at this level. D. The last hierarchical level of the phylogenetic tree. The tree was rooted at the root sequence. Purple sequences on the right side correspond to bacterial GST sequences. Sequence for eukaryotic GDAP1 and GDAP1L1 are in the left side branches. The grey circle denotes the collapsed branch for the eukaryotic GSTs (see panel B) and the magenta circle denotes the second group of bacterial GST sequences (see panel C).

https://doi.org/10.1371/journal.pone.0284532.g006

Fig 6B and 6C show that eukaryotic GSTs split off from other GSTs very early in evolutionary history. This was not to be expected. In addition, Fig 6C shows that the evolutionary history at this level could not be resolved (represented by a red line) despite sampling over millions of MCMC steps. Fig 6D shows the last level of the phylogeny, where the consensus does converge to a proper bifurcating tree. Towards the right side, there is a branch of bacterial GST sequences (purple), and in the left branch, the GDAP1/GDAP1L1 sequences are found again. The figure also shows that orthologs of GDAP1 and GDAP1L1 are found in mammals, birds, amphibians, and fishes and likely orthologs of those genes in invertebrates as well.

One may well want to interpret the results in terms of the endosymbiotic theory and remember that archaea are closer to eukaryotes than bacteria [65]. Since GDAP1 is an outer mitochondrial membrane protein, one might expect it to be more closely related to archaeal proteins. The tree, however, does not confirm or refute this. Calculating for longer simply generates more tree topologies with similar probabilities. Ultimately, one might admit that the history of the protein family involves too many duplication or gene transfer events to fit neatly into a simple timeline.

Sequence entropy in GDAP1

Sequence conservation (entropy) calculations were carried out on the different groups to identify their conserved residues and sites which differ between the groups. We first consider the GDAP1 subset of sequences (Fig 7A). 44 residues were found with an entropy score ≤ 0.10 (highly conserved) and a fraction of non-gap characters in the alignment > 70% (Table 4); the entropy is mapped onto the GDAP1 structure in Fig 7B. The data highlight several interesting residues, which are known to be sites of CMT mutations, relevant for folding and structure, or with possible functions in ligand binding. These aspects are discussed more below.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Sequence entropy analysis.

A. Entropy plot for GDAP1. The positions of the mutations studied here are shown with dots, and these positions indicate a high level of conservation (low entropy). B. Mapping of entropy onto the GDAP1 crystal structure monomer. Blue indicates high entropy and red low.

https://doi.org/10.1371/journal.pone.0284532.g007

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Entropy analysis of the GDAP1 subfamily.

Shown are the positions (human GDAP1 reference sequence numbering) with the lowest entropy (S < 0.1).

https://doi.org/10.1371/journal.pone.0284532.t004

One intuitive approach to search for functionally important sites is to focus on residues involved in ligand binding, i.e. residues that are located in the active centre. However, the question of whether GDAP1 is an active GST, or even lacks glutathione binding completely, has not been clearly resolved. Hexadecanedioic acid has been identified as possible ligand and subsequently co-crystallised with GDAP1 [11]. Gln235, Trp238, Arg282, Arg286, and Lys287 contact the ligand; Trp238, Arg282 and Arg286 are highly conserved (Table 4). This suggests that they could be important for ligand binding and that mutations in these positions could affect GDAP1 function; indeed, Trp238 and Arg282 are CMT mutation sites. These 3 residues cluster on the side of GDAP1, which must face the MOM surface. One option is that the bound fatty acid, in fact, mimics the membrane surface, and that these residues are directly involved in membrane surface binding at the outer surface of the MOM.

Table 5 shows the entropy calculated for a combined data set (GDAP1, GDAP1L1, and GST). Since the GST superfamily is so diverse, we identified sites with entropy < 0.2 (less conserved than the calculation within a family). The 14 sites that are conserved and present in at least 70% of the sequences are listed in Table 5. These include Leu239 and Gly241 in the ⍺6-⍺7 loop and Ala247 and Asp 248 on helix ⍺7. The effect of the A247V mutation on protein stability may have been therefore expected based on its conservation across this very broad set of GSTs.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Entropy analysis for the large data set including the GST dataset.

Shown are residue positions with the lowest entropy (S < 0.2).

https://doi.org/10.1371/journal.pone.0284532.t005

Kullback-Leibler divergence

Table 6 shows the sites with a high KL divergence. These are those which are conserved in at least one group and different in the second group. The results have not been filtered by ignoring gaps since an insertion or deletion is of interest. Of the GDAP1 residues Cys51, Gly83, Tyr124, and Glu228 that are not assigned to a gap symbol, Cys51 and Gly83 are located close to the dimerisation site of GDAP1. Gly83 is a target of a non-pathological polymorphism, G83A [66]. The β sheet is the dimerisation site of GDAP1, and the larger amino acids Phe and Asp found in canonical GSTs may lead to steric hindrance during dimer formation. On the other hand, Phe and Asp are predestined for stacking and polar interactions respectively. Since the binding site for glutathione in GST takes place near this segment, these residues could be important for ligand binding in the GSTs, which may not be relevant for GDAP1, due to the apparent lack of glutathione binding. The observations support the assumption that the binding site for a small-molecule ligand, if any, in GDAP1 may be located elsewhere.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 6. KL divergence.

KL₁ and KL₂ are the Kullback-Leibler divergence using the GDAP1 or GST group as reference distribution. S₁ and S₂ are the sequence entropy within the two groups. GDAP1 and GST show the most common residue within each group.

https://doi.org/10.1371/journal.pone.0284532.t006

On the other hand, Tyr124 and Glu228 in GDAP1 are close to each other in 3D space, being located on the core helices ⍺3 and ⍺6, respectively. The conversion of tyrosine to tryptophan (Tyr124 → Trp) and glutamic acid to glutamine (Glu228 → Gln) between GDAP1 and canonical GSTs causes only small changes in physicochemical properties. However, the central location of these residues in the GDAP1 fold suggests this finding may reflect a structural or functional aspect specific to GDAP1.

Of additional interest are Lys193 and Asn166, which are both close to the GDAP1-specific insertion. In contrast to the crystal structure, where Asn166 is not visible, the AlphaFold2 model predicts that the helix ⍺6 is bent, which brings Lys193 and Asn166 into close proximity. This could enable intramolecular interactions or provide a specific site for protein-protein interactions. Overall, the results indicate that certain subdomains in GDAP1 may have evolved to fulfill a different function compared to canonical GSTs.

CMT mutations at the GDAP1 protein level

CMT-linked missense mutations are relatively common in the GDAP1 gene compared to other CMT target genes, especially when the size of the protein is taken into consideration [22, 71]. With careful examination of structural models and using them as inputs for further bioinformatics analyses, the mutations are observed to cause subtle changes in intramolecular residue interaction networks, in line with experimental data. A comprehensive understanding of the effects of single mutations requires observation of the local structure coupled to experimental data at the protein and cellular levels.

Of the residues highlighted by the entropy analysis either within the GDAP1 sequence set or the entire GST set, several are targets for missense mutations. Here, we shall briefly compare selected CMT mutations in GDAP1. The immediate environment of each residue is considered to shed light on local effects of each mutation.

Tyr29 is conserved at the symmetry axis of the GDAP1 dimer, forming a H bond between OH groups of Tyr29 from the two protomers. Y29S [72] would both remove this polar interaction and make the dimer interface much less hydrophobic.

Leu239 is at the tip of the ⍺6-⍺7 loop, inserted into the structural core, and its mutation to Phe has been reported in CMT [73]. The side chain is close to those of Cys240 and Ala247, and hence, a larger residue at this position could similarly affect protein stability as A247V. Similar effects could be foreseen for the C240Y mutation [74]. Pro111 is in the ⍺2-⍺3 loop, close in 3D space to the N terminus of helix ⍺7 and the ⍺6-⍺7 loop. Its mutation to His [75] would cause steric hindrance and alter the conformation of the ⍺2-⍺3 loop, possibly destabilising the fold.

Ala247, studied in the current work, lies centrally in the GDAP1 fold on helix ⍺7. A247V is linked to CMT, showing that no larger residue fits into this tight space and A247V destabilised the GDAP1 fold. Ala247 is one of the most conserved residues in the GST superfamily, indicating a role in the GST fold. Other residues conserved in the GDAP1 set include Leu255 and Gly262, which lie in the middle and at the end of helix ⍺7, respectively. L255F [76] would be expected to cause similar steric hindrance as A247V, and G262E [72] will likely disturb the tight turn right after helix ⍺7 and cause steric clashes. Both Leu255 and Gly262 are highly conserved in GDAP1 sequences.

Arg282, studied here in the form of the CMT mutation R282H, is one of the most conserved residues in the GDAP1/GDAP1L1 subfamily. Its strong hydrogen bonding to the carbonyl groups in the ⍺6-⍺7 loop indicates a central role in GDAP1 structure. In addition to R282H, the mutation R282C has been reported [77], and we can expect it to have the same kind of effect as R282H, losing the side-chain interactions of Arg282 to the backbone of the ⍺6-⍺7 loop. Arg282 is, in addition, optimally stacked against the aromatic side chain of Trp238, which is both a conserved residue and a CMT mutation target itself [78]. This Arg-Trp-(⍺6-⍺7) unit is likely important for the stable folding of GDAP1.

Arg310 is not present in the constructs we used for crystallisation, and hence, we have no high-resolution data on its conformation. However, the AlphaFold2 model of GDAP1 extends our crystal structures and predicts that Arg310 points away from the MOM surface and forms two salt bridges to acidic residues in the GDAP1 fold (Fig 8). Interestingly, Arg310 resides in the segment originally coined the hydrophobic domain; in our view, this fragment corresponds to an amphipathic helix, which could bind onto a phospholipid membrane surface. The AlphaFold2 prediction supports this view. Arg310 is likely to link this membrane-bound helix to the folded core of GDAP1 and, therefore, directly affect the conformation of full-length, membrane-bound GDAP1. This is supported by the lack of a destabilising effect on recombinant GDAP1 structure in our experiments; Arg310 is located outside the folded core and has a more subtle role.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Location of Arg310 outside the GST-like core.

Arg310 within the amphipathic helix preceding the transmembrane domain is predicted to make salt bridges back towards the core domain.

https://doi.org/10.1371/journal.pone.0284532.g008

The clinical profiling of CMT patients and experiments in cell-based models have shown that mutation of residues located within or near the vicinity of the transmembrane helix are severe [79, 80]. These mutations likely affect the proper localisation into the MOM, leading to impaired GDAP1 folding and function. So far, studies in vitro on the GDAP1 protein have been done only on soluble constructs, lacking the transmembrane domain. The stability of full-length GDAP1 in vivo will also involve interactions with the lipid bilayer, which is a topic of ongoing work.

The central role of the ⍺6-⍺7 loop

Despite the broad spectrum of CMT disease mutations affecting GDAP1, certain general conclusions hold. The intramolecular interaction network described here shows that the critical helices in the C-terminal GST-like core domain are linked by residues that correspond to CMT target locations and/or are conserved in evolution. Therefore, these residues are likely essential for the structural integrity of GDAP1 and its function.

The ⍺6-⍺7 loop is a central feature in the structure of GDAP1, but also in all canonical GSTs. This loop inserts itself back into the protein core, being a central interaction hub between helices ⍺3, ⍺6, ⍺7, and ⍺8 (Fig 9). In GDAP1, the 3 residues at the tip of the loop (238–240) are all targets for CMT disease mutations, highlighting the crucial importance of this segment. Our entropy analysis confirms the strong conservation of this segment. We can then take a broader scope and further look at residues interacting with the ⍺6-⍺7 loop in GDAP1. Of note, several of the residues in direct contact with the ⍺6-⍺7 loop are known CMT mutation target sites and/or highly conserved (Table 4). Hence, the inward-bending ⍺6-⍺7 loop appears to be a central structural feature in GDAP1, and its alteration either directly or via disturbance of intramolecular interactions may be a common mechanism for CMT.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. An overview of the ⍺6-⍺7 loop and its surroundings in light of known CMT mutations and conserved positions.

In this stereo view, the ⍺6-⍺7 loop is coloured orange, green shows positions for CMT mutations, yellow positions that are highly conserved, and blue the positions that are both targeted by CMT mutations and highly conserved. Selected residues are labelled for clarity.

https://doi.org/10.1371/journal.pone.0284532.g009

GDAP1 as a member of the GST superfamily

The predictions that GDAP1 would be structurally related to GSTs have been confirmed by recent structural studies by us and others [10–12]. Important differences to canonical GSTs exist, however. The sequence identity between GDAP1 and any enzyme with known GST activity is very low, and contradictory results have been obtained as to the GST activity of GDAP1. In our hands, GDAP1 does not bind glutathione or act as a GST [11]. The latter is logical, since the GDAP1 dimers are formed completely differently from canonical GST, in which the active site in fact lies at the dimer interface [81]. We hypothesise that the residues conserved across the complete GDAP1/GST set are crucial for correct folding. On the other hand, residues additionally highlighted in the GDAP1/GDAP1L1 dataset may relate to GDAP1-specific functions and/or unique structural aspects of this subfamily.

Our conservation analyses highlight the low conservation of GDAP1 (and GDAP1L1) in the GST superfamily. While the entire evolutionary pathway cannot be traced based on the analyses, it is intriguing that GDAP1 is closer to bacterial than eukaryotic GSTs. It is possible that during evolution, GDAP1 has lost the characteristic GST activity, while becoming an integral membrane protein of the MOM. Its functions could, therefore, be mediated through protein-protein interactions instead of enzymatic activity. The functions can be redox-regulated, which could explain the observation that the human GDAP1 dimer is mediated by a disulphide bridge via Cys88 [11].

To complement the above analyses, we superposed crystal structures of human GDAP1 and a canonical GST, that from S. japonicum [81], and analysed the current GDAP1 crystal structures with respect to GST. This is the GST widely used in molecular biology applications as a fusion tag for affinity purification. We were interested in the residues affected by CMT mutations, especially those crystallised here. Hence, of specific interest were the apparent non-conservation of Arg120 and the conservation of Ala247 and Arg282 (Fig 10A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Comparisons of the studied mutations to a canonical GST from S. japonicum.

GDAP1 is in gray and GST in blue. A. Overall view of the GDAP1-GST superposition, with the CMT mutation sites crystallised here highlighted in pink. B. Arg120 (pink) in GDAP1 and the corresponding interaction in GST, made by an Arg from a nearby helix. C. Tight packing of Ala247 (green) in the GDAP1 structure; Ala at this position is highly conserved across the whole GST family. Note Tyr124, which was highlighted in the KL divergence analysis, making direct van der Waals contact with Ala247. D. The interactions of Arg282 (pink) towards the ⍺6-⍺7 loop and Trp238 are conserved in GST.

https://doi.org/10.1371/journal.pone.0284532.g010

Arg120 in GDAP1 is a target of several CMT mutations [82–84], and its interactions with the backbone carbonyl of Cys240 in the ⍺6-⍺7 loop appear central. Structures of the mutants R120Q and R120W indicate loosening of the structure locally, as well as loss of the hydrogen bond, which is accompanied by a decrease in observed heat stability in R120Q and R120W [12]. Arg120 is not conserved in GSTs at the sequence level; however, an Arg residue from the neighbouring helix in GST reaches the same position and makes similar interactions to the same carbonyl group in the ⍺6-⍺7 loop (Fig 10B). These observations highlight the importance of the central ⍺6-⍺7 loop in intramolecular interaction networks and fold stability of GDAP1.

Ala247 of GDAP1, although perhaps a mundane residue per se, appears surprisingly conserved in the GST superfamily in our dataset (Table 4). In the superposed individual GST structure (Fig 10A), this residue is Pro, which fits well into the structure due to slightly different conformations of the main helices in GST. In GDAP1, Ala247 is so snugly packed (Fig 10C) that even the addition of two methyl groups in the A247V variant causes protein instability and disease, despite only minor effects in the crystal state. A247V is an example of a mutation introduced into the hydrophobic core that may have long-range effects on the entire protein structure stability. One of the residues in GDAP1 highlighted by the KL divergence analysis, Tyr124, closely interacts with Ala247, indicating a GDAP1-specific arrangement at this site. Additionally, one of the most conserved GDAP1 residues, Asp248, is central in a hydrogen bonding network linking the conserved CMT target Tyr279 into the picture (Fig 10C).

GST has an Arg corresponding to Arg282, and this residue is highly conserved. Arg282 in GDAP1 interacts directly with the backbone of the ⍺6-⍺7 loop, and a similar interaction is observed in GST (Fig 10D). This is another example of an Arg–loop backbone interaction in GDAP1 that is both conserved and relevant for human disease mutations.

Biological implications

The function of GDAP1 at the molecular level remains enigmatic. While indications exist from functional studies [2, 3, 80], evidence is still incomplete for any enzymatic activity as well as direct protein-protein interactions. Redox regulation seems to play a role [9, 74, 85], but is this related to an enzymatic GST-like activity, or regulation of oligomeric state and/or protein-protein interactions?

At the molecular level, we believe to have identified important residue interaction networks between the core helices in the GST-like domain of GDAP1, strongly interacting with the inward ⍺6-⍺7 loop. These networks could be important for both GDAP1 stability and its interactions with other molecules, such as the membrane or the cytoskeleton. The correct conformation of GDAP1 on the MOM, as well as its interactions with other proteins, will then directly or indirectly affect mitochondrial dynamics to promote correct development and function of the nervous system. The disease mutations may–due to their involvement in the same intramolecular networks–cause similar overall effects on GDAP1 stability and properties, which then leads to the CMT disease phenotype in patients. In line with this, all the missense mutations we have studied at the protein level allow for GDAP1 folding, while at the same time decreasing the stability of the fold.

It should be noted that due to the large pool of GDAP1 mutants causing CMT, our experimental sample size is still relatively small, and the hypothesis may not be correct in all cases. However, we have picked mutations from different core secondary structure elements for experimental analyses to account for an incomplete dataset, and predictions have been done for all mutations [12], indicating a general trend of structural destabilisation upon CMT missense mutations in GDAP1.

Concluding remarks

Considering mitochondrial dynamics and interactions with other organelles, the homology to GSTs brings attractive prospects for GDAP1 function. GDAP1 takes part not only in mitochondrial fission and fusion, but also in interactions with the endoplasmic reticulum, peroxisomes, Golgi, and the cytoskeleton. These findings coupled with structural and biophysical data shall aid the understanding of the pathophysiological mechanism of GDAP1-linked CMT and may affect future treatment approaches. Future studies are needed to identify proteins and small molecules directly interacting with full-length GDAP1 in the physiological setting, allowing further structural investigations on the related molecular processes in nervous system function and disease. In a wider setting, we hypothesise that decrease in overall protein stability upon missense mutations is one common mechanism for CMT at the molecular level.