Thermodynamic Stability of Human γD-Crystallin Mutants Using

Article pubs.acs.org/JPCB

Cite This: J. Phys. Chem. B 2019, 123, 5671−5677

Thermodynamic Stability of Human γD-Crystallin Mutants Using Alchemical Free-Energy Calculations Rodrigo Aguayo-Ortiz,† Augusto Gonzaĺ ez-Navejas,† Giovanni Palomino-Vizcaino,§ Oscar Rodriguez-Meza,‡ Miguel Costas,‡ Liliana Quintanar,§ and Laura Dominguez*,† Facultad de Química, Departamento de Fisicoquímica and ‡Facultad de Química, Laboratorio de Bio-fisicoquímica, Departamento de Fisicoquímica, Universidad Nacional Autónoma de México, Mexico City 04510, Mexico § Departamento de Química, Centro de Investigación y de Estudios Avanzados (Cinvestav), Mexico City 07360, Mexico Downloaded via UNIV OF SOUTHERN INDIANA on July 17, 2019 at 15:12:06 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

S Supporting Information *

ABSTRACT: γD-Crystallin (HγDC) is a key structural protein in the human lens, whose aggregation has been associated with the development of cataracts. Single-point mutations and post-translational modifications destabilize HγDC interactions, forming partially folded intermediates, where hydrophobic residues are exposed and thus triggering its aggregation. In this work, we used alchemical free-energy calculations to predict changes in thermodynamic stability (ΔΔG) of 10 alanine-scanning variants and 12 HγDC mutations associated with the development of congenital cataract. Our results show that W42R is the most destabilizing mutation in HγDC. This has been corroborated through experimental determination of ΔΔG employing differential scanning calorimetry. Calculations of hydration free energies from the HγDC WT and the W42R mutant suggested that the mutant has a higher aggregation propensity. Our combined theoretical and experimental results contribute to understand HγDC destabilization and aggregation mechanisms in age-onset cataracts.

1. INTRODUCTION Crystallins (α, β, and γ) constitute nearly 90% of the total soluble protein content in the human lens.1 α-Crystallins are small heat shock proteins (sHSP) with a chaperone function that prevents aggregation of partially unfolded β- and γcrystallins.2 Oligomeric β- and monomeric γ-crystallins are the major structural proteins in the soft human lens.3 γ-Crystallins are proteins of approximately 175 amino acids in length that exhibit four Greek key motifs (eight intercalated β-strands) organized at the N-terminal (N-td) and C-terminal domains (C-td) (Figures 1A and S1). The Greek key motifs in γcrystallins are linked to higher refractive index increments and more stabilizing salt-bridges than other proteins containing these supersecondary structures.4 Seven different types of γcrystallins have been reported (γA-γF, γS), from which γC, γD, and γS are located in the human eye lens.5 Human γDcrystallin (HγDC) is one of the most abundant crystallins in the central lens nucleus, and its aggregation has been associated with the development of age-onset cataracts.1 Cataract is the pathological opacification of the lens caused by either the breakdown of the lens microarchitecture or upon aggregation of crystallins.3,6 It has been proposed that the mechanism of HγDC aggregation involves the formation of partially folded intermediates (N-td unfolded and C-td folded) that expose hydrophobic regions and favor crystallin polymerization and aggregation (Figure 1B).7 Single-point mutations and post-translational modifications in crystallins can cause protein unfolding due to destabilization of native interac© 2019 American Chemical Society

Figure 1. (A) HγDC structure showing the N-terminal (red) and Cterminal (teal) domains. (B) Schematic representation of the proposed HγDC unfolding model.7 (C) Depiction of 12 congenital missense mutations and three deamidation modifications (*) reported in HγDC.

Received: February 25, 2019 Revised: June 13, 2019 Published: June 14, 2019 5671

DOI: 10.1021/acs.jpcb.9b01818 J. Phys. Chem. B 2019, 123, 5671−5677

Article

The Journal of Physical Chemistry B

AMBER99SB-ILDN29 force field implemented in GROMACS 4.6.7.30 The model was solvated by using a three-point water model (TIP3P) in a cubic box and brought to a concentration of 0.15 M with sodium and chloride ions. The system was energy minimized using the steepest-descent algorithm and equilibrated for 5 ns under the NVT and NPT ensembles. Equilibration was followed by a 200 ns production run with a 2 fs time step recorded every 5 ps. The temperature was set to 300 K using the v-rescale coupling thermostat algorithm, and the pressure was set to 1.0 bar using the isotropic Parrinello− Rahman barostat algorithm. The Lennard-Jones potential was calculated within a cut-off radius of 1.0 nm. Electrostatic interactions were calculated between 0 and 1.2 nm, after which the electrostatic interactions were calculated using the particle mesh Ewald (PME) approach. Bonds involving hydrogens were constrained using the Linear Constraint Solver algorithm. 2.1.3. Free-Energy Calculations. The double free-energy differences (ΔΔG) of the 22 HγDC mutations were calculated using the alchemical molecular dynamics methodology implemented by Gapsys and co-workers.31−33 Hybrid structures and topologies of the HγDC mutant models were prepared with the PMX software,31 and the tripeptide hybrids were retrieved from the PMX web-server.33 For this study, we employed a double-system/single-box setup to keep a neutral simulation box during charge changing modifications and to directly obtain the ΔΔG value during the transitions.25 The tripeptide was positioned approximately 3 nm away from the HγDC structure, applying position restraints to the Cα atom of the middle residue. The states A (HγDCWT + tripeptidemutant) and B (HγDCmutant + tripeptideWT) were energy minimized and equilibrated (NVT and NPT ensembles), followed by 50 ns of production runs for the two states using the AMBER99SB-ILDN force field implemented in GROMACS 4.6.7. The MD parameters employed during the simulations are the same as those described above. A total of 100 snapshots from the last 10 ns of each production simulations were extracted to perform 100 independent alchemical transitions of 200 ps each. The work values were retrieved to estimate the free-energy differences employing the Crooks Gaussian Intersection method.34 2.1.4. 3D-RISM Calculations. Three 300 ns MD simulations for each wildtype (WT) and W42R mutant HγDC models were performed using AmberTools 16.35 The topology file of the system was prepared in the tleap module with the ff14SB36 force field using a 8 M urea aqueous urea solution box.23 A chlorine ion was added to the W42R variant systems for charge neutrality. Long-range electrostatic forces were treated using the particle mesh Ewald (PME) method,37 and van der Waals forces were truncated at 1.0 nm. Bonds containing hydrogens were constrained using the SHAKE algorithm.38,39 Energy minimized systems were heated at constant volume, varying the target temperature linearly from 0 to 100 K for 5 ps, subsequently heated with an isotropic pressure scaling and a 2 ps pressure relaxation time for the Berendsen barostat,40 varying linearly from 100 to 425 K over 1 ns. A 1 ns equilibration at 425 K with isotropic pressure scaling was performed for each system, followed by three independent 300 ns NPT simulations. A total of 300 snapshots were extracted from each replica to calculate the hydration free energy using the three-dimensional reference interaction site model (3DRISM) with the Kovalenko−Hirata closure approximation.41−45 The C-td tilt angle of each snapshot was calculated using an in-house python script employing the inbuild libraries

tions.8−10 Twenty mutations in HγDC have been associated with congenital cataracts: 12 missense, two frameshift (premature termination), and six nonsense mutations.5,11 Moreover, deamidation, oxidation, methylation, carbamylation, and the formation of disulfide bonds are the most common post-translational HγDC modifications identified in the aging lens.9,12,13 Truncation of the polypeptide chain due to the premature termination of frameshift (E103fs and G164fs) and nonsense (Y55X, Q100X, Y133X, E134X, R139X and W156X) mutations leads to the partially unfolded HγDC structure with high aggregation propensity.14 On the other hand, most of the missense mutations do not generate structural changes that destabilize HγDC, but they have been related to the formation of specific types of cataracts.15−17 Although the mechanism of destabilization of some mutants is not known, analysis of structural changes and intradomain interactions could facilitate their understanding. Several experimental and computational studies have shown HγDC instability when exposed to elevated temperatures and/ or high concentrations of denaturing agents (i.e., urea and guanidine hydrochloride).7,18−23 However, single-point mutations in the HγDC gene promote destabilization of its Greek motifs at lower temperatures, leading to its early aggregation. To understand the effect of point mutations on the native to intermediate transition of the HγDC structure, King and coworkers conducted experimental alanine-scanning mutagenesis studies of four key amino acids comprising interdomain sidechain interactions and six N-td aromatic residues.18,19 The authors concluded that aromatic to alanine single mutations highly destabilize the N-td Greek-motif, whereas only interdomain double mutations affected the stability of HγDC. These studies have promoted the interest to characterize and predict variations in the thermodynamic stability (ΔΔG) of native and mutant HγDC structures.24 In this work, we aim to understand the effect of selected mutations in the stability of HγDC. Therefore, we performed alchemical molecular dynamics (MD)-based free-energy calculations to predict the thermostability of 22 mutations in HγDC: 10 alanine-scanning variants, 9 missense mutations (of 12 reported), and 3 deamidation modifications (Figure 1C). Since the alchemical method does not currently support proline transitions, the P23T, A35P, and R36P missense mutations were discarded from this study.25 The reliability of the method was verified by comparing calculated ΔΔG values of the 10 alanine-scanning mutations with available experimentally measured ΔΔG energies,18,19 and only for the case of W42R, the ΔΔG value obtained here using differential scanning calorimetry (DSC). Furthermore, we calculated hydration free energies of the WT and the most destabilizing variant (i.e., W42R) using the three-dimensional reference interaction site model (3D-RISM) methodology,26 advocating that higher hydration energies should be related to protein aggregation susceptibility and protein thermodynamic stability.

2. MATERIALS AND METHODS 2.1. Computational Methods. 2.1.1. Model Preparation. The X-ray crystal structure of the human γD-crystallin (HγDC) structure at 1.25 Å resolution (PDB ID: 1HK027) was retrieved from the Protein Data Bank (PDB).28 We added the missing side chain of E17 and fixed the alternate positions of the residues with the highest average occupancy values. 2.1.2. Molecular Dynamics Simulations. An initial MD simulation of the HγDC structure was performed using the 5672

DOI: 10.1021/acs.jpcb.9b01818 J. Phys. Chem. B 2019, 123, 5671−5677

Article

The Journal of Physical Chemistry B of MD Analysis software.46 Figures were generated with PyMOL v0.9,47 and plots were constructed with Gnuplot 5.0.48 2.2. Experimental Methods. 2.2.1. Cloning, Expression, and Purification. Missing His tag pET16b plasmids, encoding tagless HγDC WT, and W42R mutant were used.21 As formerly described,49 recombinant proteins were expressed in BL21-CodonPlus (DE3)-RIL Escherichia coli and purified by precipitation with (NH4)2SO4, followed by gel filtration chromatography. Briefly, the cells were grown at 37 °C in super broth to an OD of ∼5 at 600 nm. Recombinant protein production was triggered using 1 mM IPTG at 18 °C for 16 hrs. Cells from half-liter cultures were harvested and resuspended in 25 mL of 10 mM CH3COONH4 buffer pH 7 and 50 mM NaCl. Cellular pellets were enzymatically lysed with lysozyme (1 mg mL−1) and DNase (20 ng mL−1) at 4 °C for 30 min in the presence of ethylenediaminetetraacetic acidfree general protease inhibitor (Roche, cOmplete), followed by 10−15 sonication cycles (30 s each) in ice. Total lysate was clarified by centrifugation at 14 000g for 45 min. The protein fraction salted out with 30% (NH4)2SO4 was discarded, whereas the protein precipitated with 50% (NH4)2SO4 was resuspended in 10 mM CH3COONH4 buffer pH 7 and 50 mM NaCl. Protein purification was performed using a gel filtration column (Hi-Prep Sephacryl S-100), using a fast protein liquid chromatography instrument (GE Life Sciences). Protein purity was verified by sodium dodecyl sulfate-polyacrylamide gel electrophoresis, and the fractions were pooled and stored at 4 °C until further use. 2.2.2. Differential Scanning Calorimetry. Differential scanning calorimetry (DSC) experiments were performed using a VP-DSC capillary calorimeter (Microcal, Malvern, U.K.). Protein solutions were prepared by exhaustive dialysis against either 5 mM potassium phosphate buffer pH 7.0 with 5 mM NaCl or 20 mM MOPS buffer pH 7.0 with 50 mM NaCl. The protein concentration was 0.4 mg mL−1 (20 μM), as determined by electronic absorption spectroscopy, using the extinction coefficients at 280 nm: ε = 41 cm−1 mM−1 for WT and ε = 36 cm−1 mM−1 for W42R, and a Cary-50 spectrophotometer (Varian, Australia). The scanning rate was 90 K h−1. Buffer−buffer traces (physical base line) were subtracted from sample endotherms, followed by the subtraction of the chemical base line joining (using a spline function) the native with the unfolded regions of the trace. Five and six independent determinations were made for HγDC and W42R, respectively.

in the loop connecting both domains. Based on this, states A (HγDCWT + tripeptidemutant) and B (HγDCmutant + tripeptideWT) were subjected to 50 ns of production runs using the AMBER99SB-ILDN29 force field implemented in GROMACS 4.6.7.30 The last 10 ns of each production simulation was extracted to perform 100 independent alchemical transitions. Finally, we employed the experimental ΔΔG = ΔGmutant − ΔGWT values at 298 K, i.e., the difference between the unfolding free energies for the mutant (final state) and the WT (initial state) to validate our results (Figure S3). Figure 2A shows that there is a strong correlation between calculated and experimental ΔΔG values (r = 0.97),

Figure 2. (A) Experimental against calculated ΔΔG values obtained with the Amber99sb-ILDN force field, for the 10 HγDC alaninescanning mutations. The experimental ΔΔG values are at 298 K and were obtained from King and co-workers.18,19 The dashed gray line represents ΔΔGexp = ΔΔGcalc. The HγDC structure on the right shows the mutated residues color coded as in the plot. (B) Calculated HγDC ΔΔG values for the 10 alanine-scanning variants (red), nine congenital mutations (violet), and three deamidation modifications (green).

confirming that the alchemical method predicts with great accuracy the experimentally observed stability. Similar to King and co-workers’ observations, aromatic single mutations (mainly for Y6A, F11A, Y45A, and Y50A) caused high destabilization of the HγDC structure (ΔΔGcalc < −1.5 kcal mol−1). On the other hand, interdomain single mutations showed a minor or no destabilizing effect on HγDC (ΔΔGcalc > −1.5 kcal mol−1). It is worth mentioning that the alchemical method produced the best correlation coefficient for the ΔΔGexp vs ΔΔGcalc plot when compared to 20 web-server predictors (Table S1), followed by ELASPIC,50 DUET,51 mCSM,52 and PoPMuSIC v2.1.53 These results strongly suggest that alchemical free-energy calculations are reliable to predict ΔΔGs of other HγDC mutations accurately. 3.2. ΔΔG Energies of Congenital Mutations. Figure 2B shows our calculated HγDC ΔΔG values for nine congenital missense mutations and three deamidation modifications (Table S2). Experimental reports have shown that R14C/S mutations do not produce significant protein structural modifications or thermal destabilization but they promote aggregation due to post-translational modifications such as disulfide-linked oligomers and O-phosphorylation.15,16 However, our calculations indicate that R14C/S mutations cause

3. RESULTS AND DISCUSSION 3.1. Comparison between Alchemical ΔΔG Calculations and Experimental Values. Free-energy differences of the 10 HγDC alanine-scanning mutations were calculated using a rigorous alchemical MD methodology implemented by Gapsys and co-workers, as follows: ΔΔGcalculated = ΔGmutant unfolding − 31−33 ΔGWT . Hybrid structures and topologies of folded and unfolding unfolded (tripeptide) models were prepared with the PMX software and PMX web-server.31,33 A double-system/singlebox setup was employed to keep a neutral simulation box during charge changing modifications and directly obtain the ΔΔG value during the transitions.25 To evaluate the stabilization period of the WT structure, we first performed a conventional 200 ns MD simulation of the HγDC structure (Figure S2). This analysis showed that the system reached a steady state in the first 25 ns of simulation, with small fluctuations in the N-td and C-td motifs and a higher flexibility 5673

DOI: 10.1021/acs.jpcb.9b01818 J. Phys. Chem. B 2019, 123, 5671−5677

Article

The Journal of Physical Chemistry B slight destabilization effects on HγDC (−3.0 to −2.8 kcal mol−1) due to the disruption of the E7-R14 salt-bridge, as previously described by Das et al.23 Similarly, the R36S mutant exhibits a marginally destabilizing behavior (−1.8 kcal mol−1) in the HγDC structure due to the disruption of the R36-D61 salt-bridge, which is associated with the nuclear cataract.54,55 Interestingly, R58H, G60C, and R76S mutations (related to nuclear, coralliform, and polar coronary cataracts, respectively) showed positive ΔΔG values (a stabilizing effect), which are in agreement with experimental reports that demonstrate that these mutations have neither structural nor significant thermodynamic destabilization effects in HγDC.17,27,54,56,57 The presence of the M43V (blue dot opacities) mutation in healthy controls indicates that this variant does not promote HγDC aggregation, confirming its predicted null destabilizing effect shown in Figure 2B.32 Meanwhile, modification of the glutamic acid at position 106 by an alanine (E106A, nuclear cataract) was only associated with protein solubility reduction.58 Similar to E106A and the alanine-scanning mutations made in the C-td, the N160D deamidation did not affect the thermodynamic stabilization of HγDC, unlike Q12E and N49D N-td post-translational modifications. 3.3. W42R Mutant. 3.3.1. Destabilizing Effect. Among the 22 mutations studied here, W42R (nuclear cataract) displayed the highest predicted destabilizing effect in HγDC (ΔΔGcalc = −12.6 kcal mol−1). This mutation was first reported by Wang and co-workers,59 who described that the W42R mutant is less stable and more prone to aggregate than the WT protein. This was further supported by a very large decrease in the denaturation midpoint temperature (ΔTm = −28.9 K) when comparing the WT to the mutant protein (Tm associated with unfolding of W42R was considered), and the W42R ΔGunfolding calculated by Gronenborn and co-workers.21,60 Given this situation, it appeared important to corroborate experimentally the very large ΔΔGcalc value for W42R. Using differential scanning calorimetry (see Figures S7 and S8), we obtained the following denaturation midpoint temperatures Tm (average from the replicas): 356.6 ± 0.21 K (83.4 ± 0.21 °C) for HγDC WT and 330.75 ± 0.5 K (57.6 ± 0.5 °C) for the low temperature peak and 346.4 ± 0.4 K (73.3 ± 0.4 °C) for the high temperature peak of HγDC W42R. Similarly, high temperature peaks of HγDC WT (84.3 ± 0.8 and 84.9 ± 0.3 °C) and W42R (77.8 ± 0.7 and 70.4 ± 0.5 °C) were reported by Wang et al.59 and Serebryany et al.,21 the latter also finding a low temperature peak for the mutant (56.1 ± 0.7 °C). The two transitions observed for W42R have been attributed to the unfolding of the N-terminal (low T peak) and C-terminal domains (high T peak).21 Taking as a reference point ΔG(Tm) = 0, ΔS = ΔH/Tm and ΔGunfold = ΔHunfold − T(ΔHunfold/Tm) = ΔHunfold(1 − T/Tm), where ΔHunfold is the area under the calorimetric trace. Extrapolation to 298.15 K provided the ΔGunfold of both proteins (see Table S3), obtaining a ΔΔGexp = −8.4 ± 2.0 kcal mol−1 for HγDC W42R, in excellent agreement with the calculated value (Figure 3A). With this experimental support, it is then clear that the large ΔΔG value for W42R, as compared with all other mutations (Figure 2B), is not due to inadequacies of the calculation scheme or errors in the experimental methodologies. In addition, eight of the 10 best web-server predictors identified during Ala-scanning validation also predicted the destabilizing effect of W42R mutation (see Table S4). The W42 residue is located in the inner part of the N-td Greek key motif (Figures 3A and 1C),60 a region where any modification could generate large

Figure 3. (A) Structure of HγDC highlighting the W42R mutation, with the experimental and calculated ΔΔG values determined in this work. (B) Porcupine representation of the domain displacement motion between the HγDC N-td and C-td. (C) Depiction of HγDC WT (PDB: 1HK0) C-td tilt angle. (D) Distributions from simulated ensembles of HγDC WT and W42R MD replicas (R1 to R3) projected onto (1) the approximation of the hydration free energy (ΔG*hyd) calculated with the 3D-RISM theory and (2) the C-td tilt angle. The black square depicts the ΔG*hyd and C-td tilt values obtained with the crystallographic WT structure.

alterations in its structure, as demonstrated by Serebryany et al.20 with the W42Q mutation. 3.3.2. Aggregation Propensity. In addition to the prediction and experimental corroboration of the destabilization effect of the W42R mutation, we compared the unfolding mechanism of this HγDC variant with that of the WT protein. This was achieved by performing three atomistic MD replicas (R1 to R3) using an Amber14 force field with the AMBER16 program35 (Figures S4 and S5). Based on the denaturing conditions previously described by Das et al.,22,23 each replica was simulated during 300 ns at 425 K in a 8 M aqueous urea solution. The principal component analysis of the simulated systems showed that HγDC’s unfolding process involves an initial interface breakdown that leads to a gradual distancing between domains (Figure 3B). The latter agrees with the first 5674

DOI: 10.1021/acs.jpcb.9b01818 J. Phys. Chem. B 2019, 123, 5671−5677

Article

The Journal of Physical Chemistry B

that penetration of water molecules into the dehydrated hydrophobic core induced the domain/protein destabilization in different crystallins.67

denaturing intermediate configuration (Int1) prior to N-td unfolding observed during the MD simulations of Das et al.23 and Xi et al.61 Superposition of crystal structures of HγDC WT (PDB: 1HK0) and its W42R variant (PDB: 4GR7), as performed by Gronenborn and co-workers,60 demonstrated that W42R C-td has a greater tilt angle between its two domains than the WT structure. This suggests that the W42R mutation induces an initial rupture of its interdomain interactions that could facilitate its dissociation and denaturation. To analyze the C-td tilt differences during our MD simulations, we defined a tilt angle at the HγDC WT crystal structure (PDB: 1HK0) using the center of mass of (1) the complete protein, (2) S39 (N-td) and E127 (C-td) Cα atoms, and (3) the C-td (H87 to R167) (Figure 3C). For C-td tilt calculation, the N-td of the snapshots extracted from each simulation was superimposed on the crystal structure, keeping the first two points of the reference and assigning the third to the snapshot C-td center of mass. In addition to the tilt angle, we employed the three-dimensional reference interaction site model (3D-RISM) to calculate the hydration free energy (ΔGhyd) of the systems through the unfolding simulations.26,62 Since we are employing nonequilibrium samples to calculate an equilibrium property, we considered the energy values obtained with the 3D-RISM method only as an approximation of the ΔGhyd (ΔG*hyd). Previous studies demonstrated that ΔG*hyd distribution calculated with the 3D-RISM theory have been useful to study protein misfolding (unfolded and partially folded intermediates) and aggregation propensity.41−44,63−66 Figure 3D displays the simulated distributions of the calculated ΔG*hyd values and C-td tilt angles for the three atomistic MD replicas performed with the HγDC WT and W42R proteins. The large C-td tilt angles confirm the formation of the unfolding Int1 intermediate of HγDC during the simulated replicas in both systems. However, a greater number of snapshots with C-td tilt angles higher than 80° were observed during the simulations with the W42R model, whereas the WT simulations displayed a population density close to the native structure. Moreover, the domain displacement in the W42R model simulations promoted an increase in the ΔG*hyd values. The observed energy distribution increase in the mutant during the simulated time suggests a decrease in HγDC solvation when compared to the WT model. Interestingly, ΔG*hyd increase in the W42R model is not only related to C-td tilt angle alterations, since we observed that the energy increase during the denaturation process was also due to the exposure of a cavity that contains the W42 residue (Figure S6). Additionally, when studying the solventaccessible surface area (SASA) of the W42 cavity (N-td core), we observed that the presence of the mutated arginine residue promoted an increment of the SASA and, therefore, an increase in the hydration energy. For the WT model, no solvent accessibility increase was observed in any of the simulated replicas. This suggests that the change of Trp to Arg in position 42 promotes the protein misfolding and solvent exposure of the hydrophobic N-td core, resulting in HγDC denaturation and loss of solubility. It is worth noting that the observed N-td core misfolding agrees with the N-td unfolding intermediate observed in previous7 and our experimental studies. The ΔG*hyd and C-td tilt angle increase, along with the exposure of N-td core internal hydrophobic residues, could explain the higher aggregation propensity of the W42R mutant, as compared to the WT protein, during cataract formation. A recent study performed by Yan and co-workers also reported

4. CONCLUSIONS Alchemical free-energy calculations could predict the thermostability of 10 Ala-scanning HγDC mutations, in excellent agreement with the experimentally measured ΔΔG energies reported by King and co-workers.18,19 The precision of the method allowed the correlation of the destabilizing effect of 12 congenital mutations with available qualitative and quantitative experimental information. Our results showed that only four of these mutations (R14C, R14S, R36S, and W42R) destabilize the structure of HγDC, with W42R being the most destabilizing variant due to its position and nature in the inner part of the N-td Greek key motif. Furthermore, our MD simulations on denaturing conditions showed the formation of partially folded intermediate configurations, in both WT and W42R proteins. Nevertheless, the ΔG*hyd calculations with the 3D-RISM method suggest that W42R denaturing intermediates have a higher aggregation propensity than those of WT. Therefore, our combined theoretical and experimental study shows that the W42R mutation favors the destabilization and N-td misfolding of HγDC, associated with congenital cataracts. This study highlights the potential of the methods employed here to gain important insights into the mechanisms of destabilization and partial unfolding of human lens crystallins.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.9b01818. Secondary structure depiction of the human γDcrystallin (Figure S1); RMSD, RMSF, radius of gyration, and number of contacts of γD-crystallin during MD (Figure S2); free-energy results from nonequilibrium transitions of the γD-crystallin mutations (Figure S3); MD analysis of the WT γD-crystallin simulated replicas (Figure S4); MD analysis of the W42R γD-crystallin simulated replicas (Figure S5); simulated distribution of HγDC replicas projected onto the ΔG*hyd and SASA (Figure S6); DSC results for HγDC and W42R (Figure S7); and deconvolution of the calorimetric trace for W42R (Figure S8); ΔΔG values calculated with 22 different protein stability predictors (Table S1); ΔΔG values calculated with the PMX software for the 22 γDcrystallin mutations (Table S2); ΔGunfold at 298.15 K for HγDC and W42R (Table S3); and W42R ΔΔG values calculated with 10 different protein stability predictors (Table S4) (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Rodrigo Aguayo-Ortiz: 0000-0001-9455-5397 Oscar Rodriguez-Meza: 0000-0003-1183-8821 Liliana Quintanar: 0000-0003-3090-7175 Laura Dominguez: 0000-0003-3666-2789 Notes

The authors declare no competing financial interest. 5675

DOI: 10.1021/acs.jpcb.9b01818 J. Phys. Chem. B 2019, 123, 5671−5677

Article

The Journal of Physical Chemistry B

■

the CRYBB2 and GJA3 Genes and Rare Polymorphisms. Mol. Vis. 2010, 16, 1837−1847. (15) Pande, A.; Pande, J.; Asherie, N.; Lomakin, A.; Ogun, O.; King, J.; Lubsen, N. H.; Walton, D.; Benedek, G. B. Molecular Basis of a Progressive Juvenile-Onset Hereditary Cataract. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 1993−1998. (16) Zhang, L.-Y.; Gong, B.; Tong, J.-P.; Fan, D. S.-P.; Chiang, S. W.-Y.; Lou, D.; Lam, D. S.-C.; Yam, G. H.-F.; Pang, C.-P. A Novel γD-Crystallin Mutation Causes Mild Changes in Protein Properties but Leads to Congenital Coralliform Cataract. Mol. Vis. 2009, 15, 1521−1529. (17) Zhang, W.; Cai, H. C.; Li, F. F.; Xi, Y. B.; Ma, X.; Yan, Y. Bin. The Congenital Cataract-Linked G61C Mutation Destabilizes γDCrystallin and Promotes Non-Native Aggregation. PLoS One 2011, 6, 1−9. (18) Flaugh, S. L.; Kosinski-Collins, M. S.; King, J. Interdomain Side-Chain Interactions in Human γD Crystallin Influencing Folding and Stability. Protein Sci. 2005, 14, 2030−2043. (19) Kong, F.; King, J. Contributions of Aromatic Pairs to the Folding and Stability of Long-Lived Human γD-Crystallin. Protein Sci. 2011, 20, 513−518. (20) Serebryany, E.; King, J. A. Wild-Type Human γD-Crystallin Promotes Aggregation of Its Oxidation-Mimicking Misfolding-Prone W42Q Mutant. J. Biol. Chem. 2015, 290, 11491−11503. (21) Serebryany, E.; Woodard, J. C.; Adkar, B. V.; Shabab, M.; King, J. A.; Shakhnovich, E. I. An Internal Disulfide Locks a Misfolded Aggregation-Prone Intermediate in Cataract-Linked Mutants of Human γD-Crystallin. J. Biol. Chem. 2016, 291, 19172−19183. (22) Das, P.; King, J. A.; Zhou, R. β-Strand Interactions at the Domain Interface Critical for the Stability of Human Lens γDCrystallin. Protein Sci. 2010, 19, 131−140. (23) Das, P.; King, J. A.; Zhou, R. Aggregation of γ-Crystallins Associated with Human Cataracts via Domain Swapping at the CTerminal β-Strands. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 10514− 10519. (24) Sahin, E.; Jordan, J. L.; Spatara, M. L.; Naranjo, A.; Costanzo, J. A.; Weiss, W. F., IV; Robinson, A. S.; Fernandez, E. J.; Roberts, C. J. Computational Design and Biophysical Characterization of Aggregation-Resistant Point Mutations for γD Crystallin Illustrate a Balance of Conformational Stability and Intrinsic Aggregation Propensity. Biochemistry 2011, 50, 628−639. (25) Gapsys, V.; Michielssens, S.; Peters, J. H.; de Groot, B. L.; Leonov, H. Calculation of Binding Free Energies. Methods Mol. Biol. 2015, 1215, 173−209. (26) Luchko, T.; Gusarov, S.; Roe, D. R.; Simmerling, C.; David, A.; Tuszynski, J.; Kovalenko, A. Three-Dimensional Molecular Theory of Solvation Coupled with Molecular Dynamics in Amber. J. Chem. Theory Comput. 2010, 6, 607−624. (27) Basak, A.; Bateman, O.; Slingsby, C.; Pande, A.; Asherie, N.; Ogun, O.; Benedek, G. B.; Pande, J. High-Resolution X-Ray Crystal Structures of Human γD Crystallin (1.25 Å) and the R58H Mutant (1.15 Å) Associated with Aculeiform Cataract. J. Mol. Biol. 2003, 328, 1137−1147. (28) Berman, H. M.; Wstbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissing, H.; Shindyalov, I. N.; Bourne, P. E. The Protein Data Bank. Nucleic Acids Res. 2000, 28, 235−242. (29) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved Side-Chain Torsion Potentials for the Amber Ff99SB Protein Force Field. Proteins: Struct., Funct., Bioinf. 2010, 78, 1950−1958. (30) Abraham, M. J.; Murtola, T.; Schulz, R.; Páall, S.; Smith, J. C.; Hess, B.; Lindah, E. GROMACS: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1−2, 19−25. (31) Gapsys, V.; Michielssens, S.; Seeliger, D.; De Groot, B. L. Pmx: Automated Protein Structure and Topology Generation for Alchemical Perturbations. J. Comput. Chem. 2015, 36, 348−354. (32) Gapsys, V.; Michielssens, S.; Seeliger, D.; de Groot, B. L. Accurate and Rigorous Prediction of the Changes in Protein Free

ACKNOWLEDGMENTS R.A.-O. (No. 510728/288862), A.G.-N. (No. 794315/ 620044), and O.R.-M. (No. 894945) are very grateful to CONACyT for the fellowships granted. L.D. thanks the ́ de Dirección General de Cómputo y de Tecnologias Información for the support received in the use of the HP Cluster Platform 3000SL supercomputer (LANCAD-UNAMDGTIC-306). M.C. thanks financial support from Dirección General de Asuntos del Personal Académico (PAPIITUNAM), grant number IN220519 and from Programa de Apoyo a la Investigación y el Posgrado (FQ-UNAM), grant number: 5000−9018. This research was supported by CONACyT (National Council for Science and Technology) by grant PN2076 and fellowship to G.P.-V. The authors would like to thank Dr. Eugene Serebryany and Prof. Jonathan A. King (MIT) for useful discussions and providing the plasmidencoding HγDC W42R.

■

REFERENCES

(1) Lampi, K. J.; Ma, Z.; Shih, M.; Shearer, T. R.; Smith, J. B.; Smith, D. L.; David, L. L. Sequence Analysis of BetaA3, BetaB3, and BetaA4 Crystallins Completes the Identification of the Major Proteins in Young Human Lens. J. Biol. Chem. 1997, 272, 2268−2275. (2) Acosta-Sampson, L.; King, J. Partially Folded Aggregation Intermediates of Human γD-, γC- and γS-Crystallin Are Recognized and Bound by Human AB- Crystallin Chaperone. J. Mol. Biol. 2010, 401, 134−152. (3) Bloemendal, H.; De Jong, W.; Jaenicke, R.; Lubsen, N. H.; Slingsby, C.; Tardieu, A. Ageing and Vision: Structure, Stability and Function of Lens Crystallins. Prog. Biophys. Mol. Biol. 2004, 86, 407− 485. (4) Mahendiran, K.; Elie, C.; Nebel, J. C.; Ryan, A.; Pierscionek, B. K. Primary Sequence Contribution to the Optical Function of the Eye Lens. Sci. Rep. 2014, 4, No. 5195. (5) Vendra, V. P. R.; Khan, I.; Chandani, S.; Muniyandi, A.; Balasubramanian, D. Gamma Crystallins of the Human Eye Lens. Biochim. Biophys. Acta, Gen. Subj. 2016, 1860, 333−343. (6) Hejtmancik, J. F. Congenital Cataracts and Their Molecular Genetics. Semin. Cell Dev. Biol. 2008, 19, 134−149. (7) Kosinski-Collins, M. S.; King, J. In Vitro Unfolding, Refolding, and Polymerization of Human γD Crystallin, a Protein Involved in Cataract Formation. Protein Sci. 2003, 12, 480−490. (8) Moreau, K. L.; King, J. A. Protein Misfolding and Aggregation in Cataract Disease and Prospects for Prevention. Trends Mol. Med. 2012, 18, 273−282. (9) Sharma, K. K.; Santhoshkumar, P. Lens Aging: Effects of Crystallins. Biochim. Biophys. Acta, Gen. Subj. 2009, 1790, 1095−1108. (10) Hains, P. G.; Truscott, R. J. W. Age-Dependent Deamidation of Lifelong Proteins in the Human Lens. Invest. Ophthalmol. Visual Sci. 2010, 51, 3107−3114. (11) Yang, G.; Chen, Z.; Zhang, W.; Liu, Z.; Zhao, J. Novel Mutations in CRYGD Are Associated with Congenital Cataracts in Chinese Families. Sci. Rep. 2016, 6, No. 18912. (12) Searle, B. C.; Dasari, S.; Wilmarth, P. A.; Turner, M.; Reddy, A. P.; David, L. L.; Nagalla, S. R. Identification of Protein Modifications Using MS/MS de Novo Sequencing and the OpenSea Alignment Algorithm. J. Proteome Res. 2005, 4, 546−554. (13) Wilmarth, P. A.; Tanner, S.; Dasari, S.; Nagalla, S. R.; Riviere, M. A.; Bafna, V.; Pevzner, P. A.; David, L. L. Age-Related Changes in Human Crystallins Determined from Comparative Analysis of PostTranslational Modifications in Young and Aged Lens: Does Deamidation Contribute to Crystallin Insolubility? J. Proteome Res. 2006, 2554−2566. (14) Santhiya, S. T.; Kumar, G. S.; Sudhakar, P.; Gupta, N.; Klopp, N.; Illig, T.; Söker, T.; Groth, M.; Platzer, M.; Gopinath, P. M.; et al. Molecular Analysis of Cataract Families in India: New Mutations in 5676

DOI: 10.1021/acs.jpcb.9b01818 J. Phys. Chem. B 2019, 123, 5671−5677

Article

The Journal of Physical Chemistry B Energies in a Large-Scale Mutation Scan. Angew. Chem., Int. Ed. 2016, 7364−7368. (33) Gapsys, V.; De Groot, B. L. Pmx Webserver: A User Friendly Interface for Alchemistry. J. Chem. Inf. Model. 2017, 57, 109−114. (34) Goette, M.; Grubmüller, H. Accuarcy and Convergence of Free Energy Differences Calculated from Nonequilibrium Switching Processes. J. Comput. Chem. 2009, 30, 447−456. (35) Case, D. A.; Betz, R. M.; Botello-Smith, W.; Cerutti, D. S.; Cheatham, T. E., III; Darden, T. A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W. et al. AMBER 2016; University of California: San Francisco, 2016. (36) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. Ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from Ff99SB. J. Chem. Theory Comput. 2015, 11, 3696−3713. (37) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N· log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089. (38) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (39) Miyamoto, S.; Kollman, P. A. Settle: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 952−962. (40) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (41) Blinov, N.; Dorosh, L.; Wishart, D.; Kovalenko, A. 3D-RISMKH Approach for Biomolecular Modelling at Nanoscale: Thermodynamics of Fibril Formation and Beyond. Mol. Simul. 2011, 37, 718− 728. (42) Maruyama, Y.; Mitsutake, A. Stability of Unfolded and Folded Protein Structures Using a 3D-RISM with the RMDFT. J. Phys. Chem. B 2017, 121, 9881−9885. (43) Chong, S. H.; Ham, S. Interaction with the Surrounding Water Plays a Key Role in Determining the Aggregation Propensity of Proteins. Angew. Chem., Int. Ed. 2014, 53, 3961−3964. (44) Sosnin, S.; Misin, M.; Palmer, D. S.; Fedorov, M. V. 3D Matters! 3D-RISM and 3D Convolutional Neural Network for Accurate Bioaccumulation Prediction. J. Phys.: Condens. Matter 2018, 30, No. 32LT03. (45) Sindhikara, D. J.; Yoshida, N.; Hirata, F. Placevent: An Algorithm for Prediction of Explicit Solvent Atom Distribution  Application to HIV-1 Protease and F-ATP Synthase. J. Comput. Chem. 2012, 33, 1536−1543. (46) Gowers, R. J.; Linke, M.; Barnoud, J.; Reddy, T. J. E.; Melo, M. N.; Seyler, S. L.; Domański, J.; Dotson, D. L.; Buchoux, S.; Kenney, I. M. et al. In MDAnalysis: A Python Package for the Rapid Analysis of Molecular Dynamics Simulations, Proceedings of the 15th Python in Science Conference, 2016; pp 98−105. (47) DeLano, W. L. The PyMOL Molecular Graphics System; DeLano Scientific LLC: Palo Alto, CA., 2007. http://www.pymol.org. (48) Williams, T.; Kelley, C. Gnuplot: An Interactive Plotting Program, 2016. (49) Domínguez-Calva, J. A.; Haase-Pettingell, C.; Serebryany, E.; King, J. A.; Quintanar, L. A Histidine Switch for Zn-Induced Aggregation of γ-Crystallins Reveals a Metal-Bridging Mechanism That Is Relevant to Cataract Disease. Biochemistry 2018, 57, 4959− 4962. (50) Witvliet, D. K.; Strokach, A.; Giraldo-Forero, A. F.; Teyra, J.; Colak, R.; Kim, P. M. ELASPIC Web-Server: Proteome-Wide Structure-Based Prediction of Mutation Effects on Protein Stability and Binding Affinity. Bioinformatics 2016, 32, 1589−1591. (51) Pires, D. E. V.; Ascher, D. B.; Blundell, T. L. DUET: A Server for Predicting Effects of Mutations on Protein Stability Using an Integrated Computational Approach. Nucleic Acids Res. 2014, 42, W314−W319.

(52) Pires, D. E. V.; Ascher, D. B.; Blundell, T. L. MCSM: Predicting the Effects of Mutations in Proteins Using Graph-Based Signatures. Bioinformatics 2014, 30, 335−342. (53) Dehouck, Y.; Kwasigroch, J. M.; Gilis, D.; Rooman, M. PoPMuSiC 2.1: A Web Server for the Estimation of Protein Stability Changes upon Mutation and Sequence Optimality. BMC Bioinf. 2011, 12, No. 151. (54) Pande, A.; Pande, J.; Asherie, N.; Lomakin, A.; Ogun, O.; King, J.; Benedek, G. B. Crystal Cataracts: Human Genetic Cataract Caused by Protein Crystallization. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 6116−6120. (55) VanderVeen, D. K.; Andrews, C.; Nihalani, B. R.; Engle, E. C. Crystalline Cataract Caused by a Heterozygous Missense Mutation in γD-Crystallin (CRYGD). Mol. Vis. 2011, 17, 3333−3338. (56) Li, F.; Wang, S.; Gao, C.; Liu, S.; Zhao, B.; Zhang, M.; Huang, S.; Zhu, S.; Ma, X. Mutation G61C in the CRYGD Gene Causing Autosomal Dominant Congenital Coralliform Cataracts. Mol. Vis. 2008, 14, 378−386. (57) Ji, F.; Jung, J.; Gronenborn, A. M. Structural and Biochemical Characterization of the Childhood Cataract-Associated R76S Mutant of Human γD-Crystallin. Biochemistry 2012, 51, 2588−2596. (58) Messina-Baas, O. M.; Gonzalez-Huerta, L. M.; CuevasCovarrubias, S. A. Two Affected Siblings with Nuclear Cataract Associated with a Novel Missense Mutation in the CRYGD Gene. Mol. Vis. 2006, 12, 995−1000. (59) Wang, B.; Yu, C.; Xi, Y. B.; Cai, H. C.; Wang, J.; Zhou, S.; Zhou, S.; Wu, Y.; Yan, B.; Ma, X.; et al. A Novel CRYGD Mutation (p.Trp43Arg) Causing Autosomal Dominant Congenital Cataract in a Chinese Family. Hum. Mutat. 2011, 32, E1939−E1947. (60) Ji, F.; Jung, J.; Koharudin, L. M. I.; Gronenborn, A. M. The Human W42R γD-Crystallin Mutant Structure Provides a Link between Congenital and Age-Related Cataracts. J. Biol. Chem. 2013, 288, 99−109. (61) Xi, Y. B.; Chen, X. J.; Zhao, W. J.; Yan, Y. Bin. Congenital Cataract-Causing Mutation G129C in γC-Crystallin Promotes the Accumulation of Two Distinct Unfolding Intermediates That Form Highly Toxic Aggregates. J. Mol. Biol. 2015, 427, 2765−2781. (62) Beglov, D.; Roux, B. An Integral Equation To Describe the Solvation of Polar Molecules in Liquid Water. J. Phys. Chem. B 1997, 101, 7821−7826. (63) Chong, S. H.; Lee, C.; Kang, G.; Park, M.; Ham, S. Structural and Thermodynamic Investigations on the Aggregation and Folding of Acylphosphatase by Molecular Dynamics Simulations and Solvation Free Energy Analysis. J. Am. Chem. Soc. 2011, 133, 7075−7083. (64) Chong, S. H.; Ham, S. Atomic Decomposition of the Protein Solvation Free Energy and Its Application to Amyloid-Beta Protein in Water. J. Chem. Phys. 2011, 135, No. 034506. (65) Chong, S.-H.; Ham, S. Impact of Chemical Heterogeneity on Protein Self-Assembly in Water. Proc. Natl. Acad. Sci. U.S.A. 2012, 109, 7636−7641. (66) Chong, S. H.; Ham, S. Distinct Role of Hydration Ater in Protein Misfolding and Aggregation Revealed by Fluctuating Thermodynamics Analysis. Acc. Chem. Res. 2015, 48, 956−965. (67) Zhang, K.; Zhao, W.-J.; Yao, K.; Yan, Y.-B. Dissimilarity in the Contributions of the N-Terminal Domain Hydrophobic Core to the Structural Stability of Lens β/γ-Crystallins. Biochemistry 2019, 58, 2499−2508.

5677

DOI: 10.1021/acs.jpcb.9b01818 J. Phys. Chem. B 2019, 123, 5671−5677