巩建业 ， 吴喆瑜 ， 刘嘉男 ， 廖 辉 ， 李利君 *，2， 倪 辉 ，2

（1.集美大学 食品与生物工程学院，福建 厦门361021；2.福建省食品微生物与酶工程重点实验室，福建 厦门361021）

柑橘是全球第一大种植和消费水果，是仅次于小麦和玉米的第三大国际贸易农产品[1]。目前对于柑橘类水果的主要加工方式为外果皮生产蜜柚香精、皮囊渣等提取果胶和鲜切柚果等，但对比橙子、柠檬等柑橘类水果，产量最大、附加值最高的蜜柚果汁类产品始终无法深入开发，主要原因是蜜柚压榨果汁具有明显的苦味，这大大制约了蜜柚的高附加值利用[2]。柚皮苷是柑橘类果汁加工中后呈苦味的重要物质之一，而α-L-鼠李糖苷酶能高效降解柚皮苷，生成有医用价值的普鲁宁和L-鼠李糖[3-4]。但工业生产中往往需要性质更优良的α-L-鼠李糖苷酶，这就需要进行蛋白质工程改造。研究酶与底物相互作用方式，可为工程改造奠定基础[5]。

随着核磁共振、X射线衍射、冷冻电镜及同源建模、从头计算等解析蛋白质结构技术的发展，越来越多的蛋白质三维结构被测出来。这为研究蛋白质与底物小分子奠定了基础。分子对接是一种静态的研究受体大分子与配体小分子结合作用的分子模拟方法，并通过对接相互作用能、关键残基等方面评价受体与配体的相互作用[6]。分子动力模拟按照给定温度下的麦克斯韦速率分布给原子赋予一定的初始速度，在这一定的初始条件下求解牛顿第二定律的微分方程，得到蛋白质中所有原子的运动轨迹[7]，MM-PBSA计算结合自由能可以在一定程度上反映受体与配体之间结合的主要推动力[8]。因此，在基于黑曲霉α-L-鼠李糖苷酶和柚皮苷结构的基础上进行分子对接和分子动力模拟，进而利用MMPBSA法计算结合自由能，分析黑曲霉α-L-鼠李糖苷酶和柚皮苷结合主要推动力，并揭示结合过程中起关键作用的氨基酸，探讨α-L-鼠李糖苷酶和柚皮苷的结合机理，可为开发更高效降解柚皮苷的α-L-鼠李糖苷酶提供理论依据。

1 材料与方法

1.1 模型构建及分子对接

黑曲霉 α-L-鼠李糖苷酶 （NCBI登录号:AGN92963.1）三维空间结构由实验室在前期研究中使用Modeller 9.15构建[9]，柚皮苷三维空间结构由Chembio office 2015构建并进行能量最小化。分子对接使用Autodock 4.0软件，对接算法采用Autodock 4.0软件自带的遗传算法，对接中心为（69.602，7.157，22.046）， 对接盒子大小为 （82×62×66）Å3，共对接100次，选择构象最佳的对接结果，使用Pymol 1.6.5软件进行分析并作图。

图1 柚皮苷二级结构和三维模型Fig.1 Secondary structures and three-dimensional model of naringin

1.2 柚皮苷与黑曲霉α-L-鼠李糖苷酶复合物分子动力模拟

柚皮苷-黑曲霉α-L-鼠李糖苷酶复合物结构取自分子对接产生的最佳对接构象，分子动力模拟采用Gromacs5.14[10]软件，向体系中加入SPC模型的水分子，为保证整个体系的整体电中性，加入了12个Na+离子。黑曲霉α-L-鼠李糖苷酶使用力场为GROMOS96 53a6，柚皮苷采用全原子OPLS力场。复合体系远程静电相互作用使用PME[11]方法计算，其中栅格宽度为1.0 Å；采用LINCS[12]算法约束所有键长。范德华和库伦相互作用的截断半径均为1.4 nm，短程邻近截断距离为1.0 nm。先对复合体系进行了能量最小化以移除体系的空间立体冲突，再对体系进行了1 ns的 NVT（the constant Number of particles，Volume and Temperature） 和51 ns 的 NPT（the constantNumberofparticles，Pressure and Temperature）模拟，之后对体系进行了20 ns的非限制性分子动力模拟，并分别分析体系中α-L-鼠李糖苷酶柚皮苷均方根偏差（RMSD）。

1.3 柚皮苷与黑曲霉α-L-鼠李糖苷酶复合物MM-PBSA结合自由能计算

分子力学泊松-玻尔兹曼表面积（molecular mechanics Poisson-Boltzmann surface area，MMPBSA）法[13-14]常被用于计算受体与配体分子动力模拟后结合自由能。采用G_mmpbsa工具计算柚皮苷与黑曲霉α-L-鼠李糖苷酶分子动力模拟平衡后结合自由能，结合自由能越低说明受体与配体之间的亲和力越高[15-16]。将总结合自由能分解到各个残基上，可以直观得看到各种相互作用的贡献。

1.4 柚皮苷与黑曲霉α-L-鼠李糖苷酶复合物疏水作用及氢键分析

取柚皮苷-黑曲霉α-L-鼠李糖苷酶复合物分子动力模拟最后的构象，使用Ligplot软件分析疏水作用残基。对分子动力模拟平衡后的运动轨迹，使用VMD软件进行氢键分析，并记录供体、受体、键长。

2 结果与讨论

2.1 柚皮苷与黑曲霉α-L-鼠李糖苷酶分子对接

许多报道[17-19]已经证实（α/α）6桶装结构域是α-L-鼠李糖苷酶的催化结构域。在100次对接结果中，柚皮苷与α-L-鼠李糖苷酶的对接位置均在（α/α）6的桶装结构域底部的loop区，这与Grandits等人[19]的分子对接结果完全吻合。图2（a）显示了柚皮苷与α-L-鼠李糖苷酶最佳对接结果得相对位置，图2（b）显示了柚皮苷范围内可能有结合作用的的氨基酸。Autodoc4.2依据半经验的自由能计算方法来评价受体和配体之间的能量匹配，其值越小表明配体与受体之间的亲和力越好[20]。在100次对接实验结果中结合自由能均为负值，最佳对接结果的结合作用能为-6.81 kcal/mol，这表明柚皮苷与黑曲霉α-L-鼠李糖苷酶有良好的结合作用。

2.2 柚皮苷与α-L-鼠李糖苷酶复合物分子动力模拟

我们对柚皮苷与α-L-鼠李糖苷酶的模型复合物进行了20 ns的分子动力模拟，均方根偏差（RMSD）反应出2个分子结构相同原子间的距离，可以用来评价体系稳定性，若均方根偏差随时间的变化较小，则表明体系达到平衡状态[21]。图3显示了在20 ns模拟过程中RMSD随时间变化的结果，α-L-鼠李糖苷酶在8.0 ns左右达到了平衡状态，柚皮苷在2.5 ns达到了平衡状态。分子动力模拟平衡后的轨迹用于进一步分析。

图2 柚皮苷与黑曲霉α-L-鼠李糖苷酶分子对接结果Fig.2 Molecular docking of naringin and Aspergillus niger α-L-rhamnosidase

图3 柚皮苷与α-L-鼠李糖苷酶分子动力模拟RMSD Fig.3 Molecular dynamics simulation RMSD of naringin and Aspergillus niger α-L-rhamnosidase

2.3 柚皮苷与黑曲霉α-L-鼠李糖苷酶MMPBSA总结合自由能

分子力学泊松-玻尔兹曼表面积 （MM-PBSA）是用于受体与配体分子动力模拟后计算结合自由能的一种方法，结合自由能由分子力学能量（MM部分），即真空中范德华力和库仑静电相互作用能，极化溶剂化能（PB部分），非极性溶剂化能（SA部分）三部分组成，计算所得的结合自由能可用来反映配体与受体结合的稳定性[22]。表1列出了柚皮苷与α-L-鼠李糖苷酶对结合自由能贡献的各项能量值。柚皮苷与α-L-鼠李糖苷酶之间的范德华力、库伦静电相互作用，极性溶剂化能和非极性溶剂化能分别为-285.931，-24.353，119.325 kJ/mol和-27.068 kJ/mol，这表明柚皮苷与α-L-鼠李糖苷酶相互作用以范德华作用力为主，静电相互作用及非极性溶剂化能力对结合自由能贡献较小，极性溶剂化能阻碍柚皮苷与α-L-鼠李糖苷酶的结合。此外，在开发酶抑制剂时，往往也需要考虑抑制剂与酶分子的相互作用力，计算柚皮苷与α-L-鼠李糖苷酶相互作用力为开发酶抑制剂提供了思路。

表1 柚皮苷与黑曲霉α-L-鼠李糖苷酶复合物结合自由能Table 1 Binding free energy of complex of naringin andAspergillus niger α-L-rhamnosidase

2.4 柚皮苷与黑曲霉α-L-鼠李糖苷酶 MMPBSA结合自由能分解

为了揭示黑曲霉α-L-鼠李糖苷酶中与柚皮苷结合的关键氨基酸，将MM-PBSA总结合自由能分解到各个氨基酸残基上（见图4），发现Trp236、Asp239、Pro276、Ser339、Ala340、Asp341、Tyr359、Phe461、Phe465、Phe501、Glu503、Tyr516、Pro520、Val522 这 14 个 残基对总结合自由能贡献值较大。表2列出了这些关键氨基酸的MM-PBSA结合能分解，除Glu503与柚皮苷的主要作用方式为极性溶剂化能外，其它残基与柚皮苷的主要作用方式均为范德华力和库伦电荷作用，另外Asp239也发挥了有利于柚皮苷结合的极性溶剂化能。这为蛋白质工程设计更高效α-L-鼠李糖苷酶分子提供了理论数据。

图4 柚皮苷与α-L-鼠李糖苷酶复合体系结合自由能分解图Fig.4 Binding free energy of complex of naringin and Aspergillus niger α-L-rhamnosidase

2.5 柚皮苷与黑曲霉α-L-鼠李糖苷酶疏水作用分析

非极性相互作用能对结合自由能有促进作用，柚皮苷可与黑曲霉α-L-鼠李糖苷酶形成较多的疏水相互作用。图5为通过Ligplot软件统计显示的柚皮苷与α-L-鼠李糖苷酶结合状态下的疏水氨基酸残基。

由图5（a）可知在柚皮苷与α-L-鼠李糖苷酶结合 过 程 中 ，Trp236、Ala340、Ile462、Phe461、Tyr516、Val522、Trp528等残基均与柚皮苷形成了疏水相互作用。将这些残基对应于三维空间模型，可以发现这些氨基酸残基在α-L-鼠李糖苷酶（α/α）6桶装催化结构域底部形成了一个疏水性口袋，而底物柚皮苷被包裹在这个疏水口袋里（图5（b））。由此可见疏水作用可能是α-L-鼠李糖苷酶（α/α）6桶能够结合柚皮苷的重要原因。

2.6 柚皮苷与黑曲霉α-L-鼠李糖苷酶氢键作用分析

通过对柚皮苷与α-L-鼠李糖苷酶结合自由能和疏水作用的分析，发现极性相互作用虽然很大程度上表现为阻碍柚皮苷与α-L-鼠李糖苷酶，但并不是所有形成极性相互作用的氨基酸都不利于结合的形成，例如Ser286在结合过程中也表现出了有利于结合自由能的极性作用，因为在结合过程中它与柚皮苷形成了氢键。利用VMD软件我们分析了分子动力模拟中平衡后的氢键，将体系中稳定的氢键统计在表3中。

表2 黑曲霉α-L-鼠李糖苷酶中与柚皮苷复合体系结合自由能分解表Table 2 Binding free energy of complex of naringin and Aspergillus niger α-L-rhamnosidase

图5 柚皮苷与α-L-鼠李糖苷酶间疏水作用残基及疏水作用口袋Fig.5 Hydrophobic residue and hydrophobic pockets of naringin and Aspergillus niger α-L-rhamnosidase

在分子动力模拟平衡后的分析中，柚皮苷与黑曲霉α-L-鼠李糖苷酶形成3条氢键，这3个关键的氨基酸残基分别为Ser286、Phe465、Pro520，且氢键的占有率均为100%。氢键分析帮助我们进一步确定了柚皮苷与α-L-鼠李糖苷酶之间结合的关键氨基酸位点。Rafanl等[22]利用分子动力模拟方法设计提高甘露聚糖酶活力时，提出与底物形成氢键的关键氨基酸可以作为改变酶活力的候选点。本研究分析得到的Ser286、Phe465、Pro520可以作为进一步α-L-鼠李糖苷酶活力的候选点，为蛋白质工程改造黑曲霉α-L-鼠李糖苷酶活力提供新思路。

表3 柚皮苷与黑曲霉α-L-鼠李糖苷酶分子动力模拟过程中氢键分析Table 3 Hydrogen bond of molecular dynamic simulation of naringin and Aspergillus niger α-L-rhamnosidase

3 结 语

越来越多的人采用分子动力学研究生物分子相互作用时的动态变化过程。作为一个被广泛应用的重要手段，分子动力学模拟在揭示配体-受体相互作用关系方面变得日益重要[23-24]。本研究表明范德华力是柚皮苷与黑曲霉α-L-鼠李糖苷酶结合的主要驱动力，库伦静电力和非极性溶剂化能对于结合作用较小。疏水作用分析发现柚皮苷结合在Trp236、Ala340、Ile462、Phe461、Tyr516、Val522、Trp528等残基形成在一个疏水性口袋里面，而Ser286、Phe465、Pro520与柚皮苷形成了稳定的3条氢键。这些位点都是黑曲霉α-L-鼠李糖苷酶与柚皮苷相互作用的关键位点，为深入研究α-L-鼠李糖苷酶与柚皮苷的构效关系确定了实验靶点。因此本研究为了解α-L-鼠李糖苷酶与柚皮苷之间的相互作用提供了重要信息，同时也为α-L-鼠李糖苷酶的抑制剂的开发和α-L-鼠李糖苷酶的定向改造提供了理论支持。