Ion mobility mass spectrometry and molecular dynamics simulations unravel the conformational stability of zinc metallothionein-2 species

Ion mobility-mass spectrometry (IM-MS) unraveled different conformational stability in Zn4–7-metallothionein-2. We introduced a new molecular dynamics simulation approach that permitted the exploration of all of the conformational space confirming the experimental data, and revealed that not only the Zn–S bonds but also the α–β domain interactions modulate protein unfolding.


Optimization of salt concentration for ESI-MS studies.
Our calibration was based on protein calibrants dissolved in similar solutions conditions as for the subsequent analysis (50 mM AmAc, 1 mM TCEP). In contrast, Russell and coworkers reported a TW calibration based on denaturing calibrants. 6 Moreover, no differences in the CCSD were obtained when spraying the proteins at higher salt content, 200 mM AmAc ( Figure S1).
In addition, calibration of the TW cell with different wave heights and velocities yielded similar CCS values. To note, similar wave height and velocity as it was used in their studies were used here. Therefore, the CCS differences could be more likely attributed to these calibration differences.
Avoiding drastic pH drop is especially relevant for metallothioneins, where 20 thiolates coordinating Zn 2+ ions can be easily protonated, resulting in Zn 2+ dissociation.
However, ammonium acetate solution does not constitute a buffer at neutral pH since the buffering properties of ammonium acetate are at pH 4.75 ± 1 and 9.25 ± 1, which corresponds to the acetic acid and ammonium pKas. 7 During the desolvation process in the ESI plume, protonation of acetate generates acetic acid (in ESI positive mode) and favors the formation of a [M + zH] z+ ions. While this process leads to inevitable acidification, likely close to the pKa values of acetic acid, the higher the AmAc concentration, the lower the pH is shifted. For instance, a pH drop to 6.5 can be estimated when using 100 mM ammonium acetate as a solution. 7 Lowering the pH would, in turn, protonate Cys residues directly affecting the stability of the Zn 2+ -binding sites. 8 Therefore, to mitigate the acidification process, the following experiments were mostly performed with 200 mM AmAc. deconvoluted by means CIUSuite 2 software. 9 The following parameters were used: maximum number of components = 10; peak amplitude = 0.05; peak overlap penalty mode = relaxed; expected fwhm were obtained similarly to Deslignière et al. 10 Briefly, a portion of ions from the convoluted ATD was selected and ejected to the pre-store while remaining ions ejected to the TOF. Then, those ions were reinjected from the pre-store to the array, separated, and the fwhm was calculated for the ATD recorded (1.50 ms). Slicing out the convoluted ATD to calculate fwhm provides accurate initial values for the peak modeling process. After one pass, a fwhm of 1.50 ms was calculated for the isolated slice.
As the fwhm scales as the √n, where n is the number of passes for a peak with a single conformer, we could estimate the initial fwhm for 2 (2.12 ms) and 3 passes (2.59 ms).
The resolving power was calculated as CCS/∆CCS, where ∆CCS is the extracted fwhm of the mobility peak.

IMS-CA-IMS.
To study the presence of different ion populations in the mass-selected Zn7MT2 5+ , we used a multistage IMS 2 . As above, we utilized the same parameters in the trap (5 V), post-trap bias (25 V) and helium cell (20 V) that minimize ion activation prior to the separation in the cIM device. A mobility-selected ion population or slice was S6 ejected and trapped in the pre-store array, while the remaining ions were ejected to the TOF. Then, the isolated ions were reinjected from the pre-store into the array. Unfolding of the mobility-selected ions was done by increasing both the pre-array gradient and the pre-array bias while keeping the voltage difference between these two parameters the same. Activated ions were subjected to one pass around the cyclic array and then ejected to the TOF. Worth comment is the need to perform two control experiments, named background ion signal and ion aging experiments. 11

Computational studies
Model building. The X-ray structure PDB ID 4MT2, which contained four Cd 2+ and two Zn 2+ ions, was selected as the initial structure. The initial Zn7MT2 was obtained by replacing the Cd 2+ with Zn 2+ ions, and point mutations were done to match with the human MT2 sequence with the VMD mutator plugin. 14 All of the simulations were performed using the GROMACS 2018.4 software. 15 The AMBER FF19SB force field was used to model the protein, and derived parameters were used to describe cysteine-Zn 2+ interactions. 16 Because of the lack of structural X-ray of NMR models for the partially Here, we incorporated three alternative solutions to tackle this issue, namely: (i) a mobile Na + charge scheme in which all titratable residues were set up to their default pH 7.0 protonation state, except for the Cys residues, and Na + ions carried out the positive charges. Here, the Cys residues that had bound Zn 2+ were modeled in the deprotonated form as this is how they are found experimentally in the X-ray structure and determined in this study by MS. In the partially Zn 2+ -depleted and metal-free MT2 species, the free Cys residues were modeled as protonated in agreement with our MS data; (ii) a mobile Na + /static H + charge scheme in which the negative charges (Asp and Glu residues) were neutralized while protonating the positive residues (Lys residues). 18 We should consider that Zn 2+ carries a 2+ charge, each Cys-binding residue has a -1 charge, and a free Cys residue is neutral. In our case, metallothionein does not contain enough residues to distribute the z excess protons on the protein. In this approach, it is not fully correct to consider that the acidic sites are generally neutral since carboxylates R-COOsites involved in salt bridges have been found. 18 Then, in all of the systems, several Na + were added to obtain a total system charge of 5+ and; (iii) a mobile H + approach that considers that the protons are highly mobile and the preferred proton-binding residues can change as the simulation, and therefore the protein structure evolves. [19][20] We used a recently released charge placement algorithm (ChargePlacer), 21

albeit modified to incorporate
Zn 2+ -binding residues in determining of the protonation pattern. The MD simulation was divided into multiple 50 ps NVT runs, and at the beginning of each simulation, the charge placement algorithm redistributed the protons along all of the titratable residues and all of the 20 Cys residues, independently if they had bound a Zn 2+ ion. Briefly, the ChargePlacer algorithm finds the protonation pattern that minimizes the total energy of the system, which accounts for both Coulomb repulsion and proton affinity, resembling the approach described by Konermann. [19][20] Gas-phase desolvation MD. Each protein system was solvated in a rhombic dodecahedron box with ~ 2500 TIP4P/2005 water molecules, since they provide more realistic results of the ESI droplet evaporation. 18 The three-site TIP3P water molecule is treated as a nonpolarizable molecule and exhibits a lower ~30% surface tension than the real water molecules. Here we employed the Na + mobile approach by which the aqueous droplet was charged by randomly replacing water molecules with excess Na + to obtain a system total charge of 16+. According to the charge residue mechanism (CRM) the maximum charge of positive charges that an ion can obtain is the Rayleigh limit (zR). 22 For globular spherical protein, zR = 0.0778 m 1/2 , where m is the molecular weight of the protein in Da units. 23 The value of zR can deviate for proteins that are not structurally globular. 24 To account for any deviation, we employed a 2.5zR excess of Na + ions in the initial droplets. MD runs with z ~ zR yielded similar final [M+ zNa] z+ ions and z/zR.
We carried out a pseudo-PBC approach in which the Coulomb and Lennard-Jones cutoffs were set to 300 nm, and the PBC box dimensions were set up to 900 nm 3 . Although this approach uses periodic boundary conditions, the box size and cutoffs exclude interactions between PBC images. Each system was subjected to energy minimization by 10 000 steps of steepest descent minimization, followed by 10 ps NVT equilibration to 350 K by using the Nosé-Hoover thermostat with a coupling constant of 0.1 ps. The LINCS algorithm was used to constraint bonds involving hydrogen atoms to be able to use a 2 fs time step, and neighbor list was updated every 100 steps using the Verlet method. To simulate droplet desolvation, the MD simulations were split into multiple consecutive 250 ps length NVT runs. At the end of each window, water molecules and Na + and/or Zn 2+ ions further than 30 Å from the center of mass of the protein were removed. The system was then recentered in the box, and the velocities reassigned from a Maxwell-Boltzmann distribution. This approach avoids evaporative cooling problems and speeds up the MD simulations as the number of particles is reduced. After 250 runs or 62.5 ns of NVT production at 350 K, the system was equilibrated to 500 K, and run for 5 ns to remove the last water molecules, also called "sticky" waters. Two independent runs were performed for each system, amounting to 0.9 µs of dynamics.
Collision cross section values were calculated every 250 ps using the trajectory method implemented in IMPACT software. 25 The simulations were analyzed using MDAnalysis 2.0 26-27 , MDTraj 1.98 28 , pytraj 2.0.5, 29 and in-house Python 3.5 scripts.
Gas-phase MD simulation. In contrast to the desolvation protocol, here, each protein system was first equilibrated in the presence of solvent molecules and then directly placed in pseudo-vacuum conditions. To compare it with our experimental data, the (i) mobile Na + charge scheme, (ii) mobile Na + /static H + charge scheme, and a (iii) mobile H + approach were considered to capture the charge state distribution for the [M + 5H] 5+ ions.
As above, we used a pseudo-PBC approach and ran each system for 100 ns at 298 K, and 798 K on three replicate runs for the (i) and (ii) approaches. While approaches (i) and (ii) were run for the seven protein systems, approach (iii) was run twice for 1500 ns at 298 K exclusively for the Zn7MT2 structure. The cumulative simulation time approached 11.4 µs.
Simulated Annealing. To simulate protein unfolding as obtained by CIU experiments, structures obtained after gas-phase desolvation underwent a simulated annealing (SA) protocol with the pseudo-PBC approach. In the SA protocol the temperature was rise linearly from 298 to 798 K during 10 ns, we also compared results when applying a 100 ns run. Each system was run on three replicate runs. A total of 210 ns of SA were run.
Data from three temperature ranges (300-350, 450-550, 750-850 K) was then extracted, and modeled using Gaussian kernel density estimation (KDE). Results are shown in kcal·mol -1 ·nm -2 ) with a pulling speed of 10 Å·ns -1 were used to produce extensions up to 50 Å so that Zn 2+ dissociation does not occur, as observed during collision-induced unfolding experiments. The relationship between the pulling speed and rupture or unfolding force was also considered by using several pulling rates: 1, 10, and 100 Å·ns -1 . The simulation time needed for each pulling speed was calculated as (rF-r0)/k, where rF is the final distance, r0 is the initial distance, and k is the pulling speed. Thus, 50, 5, or 0.5 ns were run for each pulling rate assayed, respectively. A pulling speed of 10 Å·ns -1 S10 was chosen as it gave comparable unfolding forces as the lowest speed considered (1 Å·ns -1 ) but saved computational time. In another set of SMD simulations, Rg was employed as a CV. Initials trials established the maximum Rg value needed to promote a conformational transition from a compact to an unfolded structure with intact Zn 2+ binding sites. Thus, Rg moved from an initial value of 11 Å to 18 Å. As above, three different force constants of 10, 25, and 50 kcal·mol -1 ·nm -2 were used. Gaussian kernel density estimation (KDE) was used to model each force-CCS dataset and presented in Figure 4C. In total, 225 ns of SMD simulation time was considered for analysis. The

Well-Tempered Metadynamics (WT-MetaD) simulations. Metadynamics 31 simulations
was used to estimate the free energy of unfolding dynamics of the protein systems obtained after gas-phase desolvation. To do so, a CV that described the distance between the center of mass between both protein domains was built. The temperature was set to a KBT of 0.6 kcal/mol and a CV width of 0.5 Å was selected. Gaussian hills deposited every 500 time steps with an initial height of 0.6 kcal/mol and rescaled with a bias factor of 80.
Such "aggressive" acceleration of the sampling ensures that all energy barriers are easily overcome at room temperature and provide a fast exploration of the conformational space, at expenses of convergence. 32 We used upper walls to avoid sampling unphysical states, that is when Zn 2+ dissociation occurs. The free energy was reconstructed after the recrossing event, and three independent runs were considered for the analysis.