Thriving in the heat – Lysine acetylation stabilizes the quaternary structure of a Mega-Dalton hyperthermoactive PEP-synthase

Over time structural adaptations enabled proteins and enzymes to have sufficient stability and flexibility to perform the basic functions of life under various environmental conditions. The catalytic cores of key metabolic enzymes of hyperthermophilic archaea work at a temperature range of 80-120 °C, similar to the conditions wher the earliest life forms may have thrived. Here we characterize a key enzyme of the central carbon metabolism of Pyrococcus furious, through an integrative approach combining structural mass spectrometry, cryo-electron microscopy, mass photometry and molecular modelling with molecular dynamics simulations. From our investigation, we unveil the structural organization of phosphoenolpyruvate synthase (PPSA). Its 24-meric assembly - weighing over 2 MDa - harbors flexible distal domains, whose proper functioning and coordination depends on widespread chemical acetylation of lysine residues. This non-enzymatic post-translational modification, along with other types of lysine modifications, also occurs on most other major protein complexes of P. furiosus. These modifications likely originated in the chemically favorable primordial conditions and gradually became highly specialized and enzyme-driven in more distantly related mesophiles and Eukaryotes.

Molecular dynamics simulations and analysis – The all-atom structures of the full c1(X4) 24-mer PPSA models, carrying highly acetylated sites at 14 positions (106, 120,185, 187, 427, 466, 492, 496, 557, 574, 641, 726, 737, 805) or unmodified lysines were coarse grained (CG), mapping their atoms to the SIRAH force field (ff), that uses a classical Hamiltonian common to most all-atom potentials to describe particle–particle interactions, and recently extended to support the PTMs most commonly found on proteins(Garay et al., 2020). Both starting structures contained a disulfide bond between Cys42-Cys189 as well as phosphorylation of Thr440. All starting structures, simulation boxes and parameters files used for the minimization, equilibration and production are provided (Supplementary Data 5). Simulation was conducted in GROMACS 2020.4(Hess et al., 2008), for which code can be found here: https://doi.org/10.5281/zenodo.5636522. The 24-aly system consisted of 94128 CG atoms for the protein representation, solvated with 179271 WT4 CG water beads (Garay et al., 2020). A neutral charge was achieved adding 6617 NaW (Na⁺) and 5633 ClW (Cl^-), corresponding to a concentration of approximately 150 mM NaCl concentration. The 24-lys system was also prepared accordingly and consisted of 93456 CG atoms for the protein, 186690 WT4 CG water beads and neutralized with 6503 NaW (Na⁺) and 5855 ClW (Cl^-). Solvation was done using the default radii of 0.105 nm for atoms not present in the VdW database (vdwradii.dat) and then removing the WT4 molecules within 0.3 nm from the solute. In all cases, eventual clashes were relaxed during the solute-restrained energy minimization. Due to the length of the loops connecting the 3 domains, an alternative configuration of the tetramers where the CD-NBD domains are sitting on top of the neighboring PPSA subunit is also possible (see Supplementary Note 1). This was named “alternative” (a) configuration, as opposed to the “original” (o), thus producing 4 starting models: a24-aly, o24-aly, a24-lys and o24-lys (Supplementary Data 5). All these configurations were subjected to a first round of production consisting of 125 ns in duplicate to define the stable configuration for further extension of the simulation time to 500 ns.

The main stages of the simulation can be summarized as follows:1) solvent and side-chain relaxation by 2 stages of 20’000 steps of energy minimization, imposing positional restraints of 1000 kJ mol^-1 nm^-2 on the whole protein (stage 1) and only on backbone beads (GN and GO, stage 2); 2) solvent NVT ensemble equilibration with a first stage where the temperature was slowly increased from 303K to 363K in 7 steps of 4 ns each, and a second stage to equilibrate the protein by gradually releasing the positional restraints from 1000 kJ mol^-1 nm^-2  on the backbone beads (GN and GO) to 100 kJ mol^-1 nm^-2 on the C-terminus to compensate for the missing stabilization effect of the Met799-Fe cluster; 3) production simulation of an additional 80 ns in NPT ensemble at 363 K and 1 bar imposing positional restraints of 100 kJ mol^-1 nm^-2 on the C-terminus. Non-bonded interactions were treated with a 1.2 nm cutoff and PME for long-range electrostatics. An integration time-step of 15 fs was used during MD production runs. The system pressure was controlled by the Parrinello-Rahman barostat (Parrinello and Rahman, 1981) with a coupling time of 4 ps. The positional restraints during the production simulation were necessary to maintain the overall system at the high energy of the particles at 363 K (90 ˚C), the temperature mimicking the near-optimum temperature for PPSA catalytic activity. All simulations were run in duplicate. Although the simulations do not accurately reflect the possible dynamics of the protein in native conditions, the computational challenges of simulating >350000 atoms for the protein only drove us to use a CG representation of the system.  To analyze the trajectories, every functional module (tetramer) was extracted from the simulation of the full 24-mer for each replicate, resulting in two sets of 12 independent trajectories (6 for each 24-mer) either with and without acetylated Lysine residues. To produce morphed movies and analyze secondary structure and interfaces, the CG tetramers trajectories were back mapped to atomistic detail using SIRAH tools via VMD. A back mapped atomistic model was produced every 35.7 ns, and subjected to 100 cycles of energy minimization in AMBER, resulting in 14 atomistic models describing each trajectory of a given tetramer over the 500 ns of simulation, further interpolated in ChimeraX using the morph command to produce the atomistic representation of the dynamics (Supplementary Data 5). RMSD calculation and trajectory analysis were done with the MDanalysis suite(Michaud-Agrawal et al., 2011).

Each probability density in Figure 5B and 5C is calculated over two independent coarse-grained molecular dynamics simulations of the last 250 ns of production of the PPSA 24-mer in the acetylated and non-acetylated form, using a snapshot frequency of 5 ns. In Figure 5B, the RMSD of the Cα atoms (equivalent to backbone GC beads in CG MD) of the NBD region (residues 1-365) in the tetramer form was calculated with respect to their conformation in the tetramer PPSA resting state. The six trajectories of all six tetramers included in the simulated PPSA complex, resulting in 3 μs ([250ns*6]*2) of accumulated tetramer simulation, were used in the probability density calculation. In Figure 5C, the RMSD of the Cα atoms (equivalent to backbone GC beads in CG MD) was calculated over the CD (residues 379-481) and PBD (residues 510-790) regions with respect to the modelled CD-PBD state of PPSA in the monomer form. All 24 PPSA monomers of which the simulated PPSA complex was composed were included in the probability density calculation, resulting in 12 μs ([250ns*24]*2) of accumulated PPSA monomer simulation.