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 <em>Pyrococcus furious</em>, 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 <em>P. furiosus</em>. 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. <strong>Molecular dynamics simulations and analysis – </strong>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^-2on 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^-2on the backbone beads (GN and GO) to 100 kJ mol^-1 nm^-2on 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^-2on 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 &gt;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.随着演化进程，结构适应性使得蛋白质与酶得以获得足够的稳定性与柔性，从而在多样环境条件下执行生命的基本功能。超嗜热古菌核心代谢酶的催化结构域工作温度范围为80~120℃，与早期生命形式可能繁盛的环境条件相仿。本研究通过整合结构质谱（structural mass spectrometry）、冷冻电镜（cryo-electron microscopy）、质量光度法（mass photometry）以及结合分子动力学（molecular dynamics, MD）模拟的分子建模手段，对激烈热球菌（Pyrococcus furiosus，原文疑似笔误为furious）中央碳代谢的关键酶进行了表征，该酶为磷酸烯醇式丙酮酸合酶（phosphoenolpyruvate synthase, PPSA）。其24聚体组装体分子量超过2 MDa，包含柔性远端结构域，这些结构域的正常功能与协调依赖于赖氨酸残基广泛的化学乙酰化修饰。 这种非酶促翻译后修饰（post-translational modification, PTM），连同其他类型的赖氨酸修饰，同样存在于激烈热球菌的绝大多数其他主要蛋白质复合物中。这类修饰可能起源于化学环境适宜的原始生命条件，并在亲缘关系更远的嗜温生物与真核生物中逐渐演变为高度特化的酶促修饰。 <strong>分子动力学模拟与分析——</strong> 携带14个高度乙酰化位点（106、120、185、187、427、466、492、496、557、574、641、726、737、805位）或未修饰赖氨酸的完整c1(X4)24聚体PPSA模型的全原子结构，被粗粒化（coarse grained, CG）处理，其原子映射至SIRAH力场（force field, ff）——该力场采用与多数全原子势能通用的经典哈密顿量描述粒子间相互作用，且近期已扩展支持蛋白质上最常见的翻译后修饰（Garay等，2020）。所有初始结构均包含Cys42-Cys189之间的二硫键，以及Thr440的磷酸化修饰。本研究提供了所有用于能量最小化、平衡及生产模拟的初始结构、模拟盒子与参数文件（补充数据5）。模拟在GROMACS 2020.4（Hess等，2008）中完成，其代码可在https://doi.org/10.5281/zenodo.5636522获取。 24-乙酰化（24-aly）系统的蛋白质表示包含94128个粗粒化原子，使用179271个WT4粗粒化水分子珠（Garay等，2020）进行溶剂化。通过添加6617个NaW（Na⁺）与5633个ClW（Cl⁻）实现电荷中性，对应约150 mM的NaCl浓度。24-未修饰赖氨酸（24-lys）系统也按相同方式制备：蛋白质部分含93456个粗粒化原子，186690个WT4粗粒化水分子珠，并用6503个NaW与5855个ClW完成电荷中和。溶剂化采用范德华数据库（vdwradii.dat）中未包含原子的默认半径0.105 nm，随后移除溶质0.3 nm范围内的WT4分子。所有体系均在溶质约束的能量最小化过程中松弛了潜在的空间冲突。 由于连接3个结构域的环区长度较长，四聚体存在另一种构象：CD-NBD结构域位于相邻PPSA亚基之上（详见补充说明1）。该构象被命名为“替代型（a）”，以区别于“原始型（o）”，由此得到4种初始模型：a24-aly、o24-aly、a24-lys及o24-lys（补充数据5）。所有构型均先进行两轮重复的125 ns生产模拟，以确定稳定构象，随后将模拟时长延长至500 ns。模拟的主要步骤可总结如下： 1) 溶剂与侧链松弛：分两个阶段进行20000步能量最小化，第一阶段对整个蛋白质施加1000 kJ·mol⁻¹·nm⁻²的位置约束，第二阶段仅对主链珠（GN与GO）施加约束； 2) 溶剂NVT系综平衡：第一阶段通过7个时长为4 ns的步骤将温度从303 K缓慢升至363 K，第二阶段通过逐步释放主链珠（GN与GO）的位置约束（从1000 kJ·mol⁻¹·nm⁻²降至C端的100 kJ·mol⁻¹·nm⁻²）来平衡蛋白质，以弥补Met799-Fe簇缺失带来的稳定效应； 3) NPT系综下额外80 ns的生产模拟，温度维持363 K、压强1 bar，对C端施加100 kJ·mol⁻¹·nm⁻²的位置约束。非键相互作用采用1.2 nm截断值处理，长程静电相互作用采用PME方法。分子动力学生产模拟的积分时间步长为15 fs。体系压强通过Parrinello-Rahman恒压器（Parrinello与Rahman，1981）控制，耦合时间为4 ps。生产模拟中的位置约束是必要的，以维持体系在363 K（90℃）下的高粒子能量——该温度模拟了PPSA催化活性的近最适温度。所有模拟均重复进行两次。 尽管本模拟未能准确反映蛋白质在天然条件下的动态变化，但由于模拟仅蛋白质部分就需处理超过350000个原子的计算挑战，我们不得不采用体系的粗粒化表示。为分析轨迹，我们从完整24聚体的每一次重复模拟中提取每个功能模块（四聚体），由此得到两组各12条独立轨迹（两种PPSA聚体各6条），分别对应带有乙酰化赖氨酸残基与未修饰赖氨酸残基的体系。为生成动态动画并分析二级结构与相互作用界面，我们使用SIRAH工具通过VMD将粗粒化四聚体轨迹映射回全原子细节。每35.7 ns生成一个反向映射的全原子模型，并在AMBER中进行100轮能量最小化，由此得到14个全原子模型，用于描述给定四聚体在500 ns模拟过程中的轨迹，随后使用ChimeraX的morph命令进行插值，生成动态过程的全原子可视化结果（补充数据5）。均方根偏差（root mean square deviation, RMSD）计算与轨迹分析使用MDanalysis工具包完成（Michaud-Agrawal等，2011）。 图5B与5C中的每个概率密度均基于两次独立的粗粒化分子动力学模拟计算得到，模拟对象为乙酰化与非乙酰化形式的PPSA 24聚体生产模拟的最后250 ns，快照采样频率为5 ns。在图5B中，我们以四聚体形式PPSA的静息态构象为参照，计算了四聚体形式NBD区域（残基1~365）的Cα原子（对应粗粒化分子动力学中的主链GC珠）的RMSD。本研究纳入了模拟PPSA复合物中所有6个四聚体的6条轨迹，累计得到3 μs（[250 ns×6]×2）的四聚体模拟时长，用于概率密度计算。在图5C中，我们以单体形式PPSA的建模CD-PBD状态为参照，计算了CD区域（残基379~481）与PBD区域（残基510~790）的Cα原子（对应粗粒化分子动力学中的主链GC珠）的RMSD。本研究纳入了组成模拟PPSA复合物的全部24个PPSA单体，累计得到12 μs（[250 ns×24]×2）的PPSA单体模拟时长。