Changes in hydration structure are necessary for collective motions of a multi-domain protein

Conformational motions of proteins are necessary for their functions. To date, experimental studies measuring conformational fluctuations of a whole protein structure have revealed that water molecules hydrating proteins are necessary to induce protein functional motions. However, the underlying microscopic mechanism behind such regulation remains unsolved. To clarify the mechanism, multi-domain proteins are good targets because it is obvious that water molecules between domains play an important role in domain motions. Here, we show how changes in hydration structure microscopically correlate with large-amplitude motions of a multi-domain protein, through molecular dynamics simulation supported by structural analyses and biochemical experiments. We first identified collective domain motions of the protein, which open/close an active-site cleft between domains. The analyses on changes in hydration structure revealed that changes in local hydration in the depth of the cleft are necessary for the domain motion and vice versa. In particular, ‘wetting’/‘drying’ at a hydrophobic pocket and ‘adsorption’/‘dissociation’ of a few water molecules at a hydrophilic crevice in the cleft were induced by dynamic rearrangements of hydrogen-bond networks, and worked as a switch for the domain motions. Our results microscopically demonstrated the importance of hydrogen-bond networks of water molecules in understanding energy landscapes of protein motions.

The most open -110.0 a ± 26.8 b 171.4 ± 18.4 The most closed -167.9 ± 12.4 -140.4 ± 14.5 a The shown dihedral angles were the averages calculated from a 50-ps time window around each conformation.
b The fluctuations of these dihedral angles were also calculated from the same time window.

Supplementary note 1: Analyses of surface topography of GDH crystal
Here we describe the procedures for analysing the measured surface topography in the following four subsections. We selected images composed of the unit cell arrays displaying almost identical structural features with a resolution to resolve individual subunits. Images were first corrected regarding the tilt of crystal surface against the plane movement of the piezo-scanner of the AFM equipment. Then, we carried out the assignment of unit cells followed by the calculation regarding the height distribution of N-domains in specified subunits. In addition, we describe the procedure to calculate the simulated height distribution of N-domain from the MD trajectory.

Supplementary note 1.1: Assignment of unit cells and N-domains
Because the N-domains of subunits D were found to be free from crystal contacts and protruded from the surface with the largest heights, these domains were used as an experimental reference to examine the magnitude of domain motions in the MD simulation. Several pairs of images of subunit D were selected ( Supplementary Fig. S2a) to calculate the positions and dimensions of unit cells from topography images. In each selected subunit image, we detected peak region of 29×29 Å 2 (a box of 6×3 pixels) and calculated the centroids of this region. Then, the vectors connecting the centroids were defined for each pair of images ( Supplementary Fig. S2a). From the lengths of the vectors, the unit cell dimensions of the surface lattice were calculated to be a = 129 Å, b = 170 Å, and γ = 90°. They were slightly larger than those determined by X-ray diffraction experiment for GDH crystal at ambient temperature (a = 112.48 Å and b=163.12 Å) prior to the AFM measurement.

Supplementary note 1.2: Measured height distribution of N-domain motion
We used the height distribution of N-domain in subunit D to describe the conformational variety of GDH on crystal surface. An image of subunit D in each unit cell was masked by using a rectangular box of 15×5 pixels corresponding to 73×49 Å 2 ( Supplementary Fig. S2b). We scanned each masked region using a box of 5×1 pixels (24×10 Å 2 ) to find which area had the largest height. Finally, we obtained the height distribution for the selected areas (Fig. 3e), the total number of which was 133 from three topographies.
To estimate the influences of the Brownian motion of cantilever probe at the tip 53,54 on the surface profile, we also measured surface topography of haemoglobin crystal 30 . We obtained a height distribution of haemoglobin subunits exposed to solvent from 99 images of the unit cells with the dimensions of 98×82 Å 2 (Fig. 2f) The obtained standard deviation of the distribution was used as a reference of the Brownian motion of a cantilever probe for analysing the height distribution of N-domains in D subunits of GDH. It should be noted that the viscosity of the stabilization buffer for haemoglobin crystal (1.5 M ammonium sulphate, 0.35 M ammonium phosphate, and 50 mM sodium cacodylate (pH 6.4)) was similar to that of the stabilization buffer for GDH crystal.

Supplementary note 2: Time window to calculate the solvent density map
Here we describe why we used a time window of 50 ps to calculate the solvent density map from the MD trajectory. Supplementary Fig. S4 illustrates the variety of solvent density maps calculated by using the different sizes of time windows. Peaks in the maps calculated using time windows of 50 and 100 ps are consistent with the crystal water sites 22,23 . In contrast, due to the noise in maps of 10 and 20 ps, it was difficult to identify the hydration sites in which water molecules resided stably. Thus, we used a time-window of 50 ps as the minimum period to reduce the noisy densities in solvent density maps.

run and the crystal structure
To examine whether the solvent density maps calculated from the MD trajectory reproduce hydration structures found in crystal structure 22,23 , we compared the solvent density maps calculated from the 1.2-ns solvent equilibration run with the crystal water sites. Because we could not identify water molecules from the disordered electron density maps at HS1 in the crystal structure, the comparison was conducted only for HS2. The calculated solvent density maps at HS2 for the six subunits were consistent with the water positions in the crystal structure ( Supplementary Fig. S5). The consistency strongly suggests that the current MD simulation, using the MARBLE software with the CHARMM27 force-field 42 and the TIP3P water model 40 , is good for the investigation of the hydration structure in the cleft region.

Supplementary note 4: Conformational variety of the sidechain of W89 at HS1
In In the trajectory of 200-ns run, we identified three rotamers designated as I, II and III with respect to χ 1 and χ 2 angles of W89 ( Supplementary Fig. S6b,c). In rotamer I, hydration water molecules freely penetrate into HS1, while rotamers II and III prevent the penetration. The populations of rotamers I, II, and III were approximately 92, 7 and 0.3 %, respectively ( Supplementary Fig. S6d). Therefore, we could analyse a huge number of snapshots of rotamer I in the analysis of hydration state of HS1 after excluding the snapshots of minor and unnecessary rotamers II and III. The exclusion of these rotamer states had little influences on our present results because of their small populations.
Supplementary note 5: Four parameters, d HS1 , Q HS1 , d HS2 and Q HS2 , to monitor hydration and structure changes at HS1 and HS2 To monitor the hydration structure changes along with the conformational changes at HS1 and HS2, we calculated the time courses of four parameters, d HS1 , Q HS1 , d HS2 and Q HS2 . Here we describe the details of the calculations.
The size of the hydrophobic pocket of HS1 along the direction of the N-domain motion was represented by parameter d HS1 , the distance between the H δ1 or H ε1 atom of F340 in the upper jaw of the active-site cleft and the midpoint of the C δ1 atom of W89 and the C γ atom of W92 in the lower jaw.
Parameter Q HS1 monitored the amount of hydration water molecules residing in the hydrophobic pocket 23 as the sum of the solvent densities of only 1-Å cube voxels inside the hydrophobic pocket of HS1. To judge whether a voxel is inside the pocket, we used line segments drawn between atoms of the upper jaw (all atoms of F340) and of the lower jaw (C α , C β , C γ , C δ1 atoms of W89 and all side-chain atoms of W92) ( Supplementary Fig. S7a). When a voxel was located within 1.4 Å from any line segment, that voxel was judged to be inside of the pocket (Supplementary Fig. S7b).
The size of the crevice in HS2 was measured by parameter d HS2 , the distance between the A190 C α atom and the E354 O ε1 atom. Parameter Q HS2 reported the amount of hydration water molecules inside the crevice of HS2 as the sum of the solvent densities calculated through the following procedure. We first approximated the shape of the crevice by a cylinder with the radius of 2.5 Å (Supplementary Fig.   S8a). The length of the cylinder is variable depending on the N-domain motion. The cylinder was set to the crevice so that its principal axis coincided with the line segment connecting the C β atom of T191 and the O ε atom of E354. We detected high solvent densities within the cylinder. Then, Q HS2 was calculated as the sum of the solvent densities in 1-Å cube voxels located within 1.7 Å from any centre positions of high solvent densities.
The time-courses of the four parameters were calculated for the six subunits throughout the 200-ns MD trajectory. Supplementary Fig. S9 compiles mutual dependences among the parameters. The results for the six subunits are almost consistent with those for subunit A shown in Figs. 6-8 in the main text.

Supplementary note 6: Enzymatic activities of wild type and W89F mutant of GDH
To investigate the effects of the hydrophobicity of HS1 on the function of GDH, the enzymatic activities of the wild-type and W89F-mutated GDH at 293 K were determined by measuring the reaction to reduce NADP + to NADPH 55 (Supplementary Fig. S11). The reaction mixture contained 100 mM Na phosphate buffer (pH 7.0), 10 mM sodium L-glutamate, 6.25-300 μM NADP + , and 4 μg/ml GDH (49).
The amount of the product, NADPH, was measured by monitoring the absorption at 340 nm with a spectrophotometer U-2900 (Hitachi High Technologies, Japan). The initial velocity of catalysis, V 0 , was calculated from linear least-square fitting of the progress curves at the initial stages within 50 and 300 s for wild type and W89F, respectively. To estimate a catalytic rate, k cat , and dissociation constant, K m , the dependence of V 0 on NADP + concentration was analysed by the Michaelis-Menten equation. The resultant values of these parameters are listed in Table S1. In calculating the contact maps, we searched pairs of residues in contact. When the distance between any two non-hydrogen atoms belonging to different residues was less than 4.5 Å, the residues were defined to be in contact. There was little significant inter-domain interactions except those at HS1 in both the most open and closed conformations (Supplementary Fig S13b,c). At HS1, the increase in the number of inter-domain interactions was observed upon the domain closure, reflecting the structural packing of the hydrophobic pocket described in the main text. These results indicate that the domain motions of subunits are free from reorganization of inter-domain interactions, which causes an energy barrier for the motion.
We also analysed the dihedral angles of the residues belonging to the hinge-bending regions in the most open and closed conformations. Significant changes in dihedral angles between the two conformations were found only in G183. G183 is in the loop connecting C-and N-domains and near