常 瑞 罗先锟 何 江 梁美莲 朱秋劲

(1贵州大学酿酒与食品工程学院/贵州省农畜产品贮藏与加工重点实验室，贵州 贵阳 550025；2 贵阳市农产品质量安全监督检验测试中心，贵州 贵阳 550004)

红肉及其制品是供给人体蛋白质、多不饱和脂肪酸、B 族维生素、硒等优质营养素的来源。但近来诸多流行性病学和免疫学研究表明，日常饮食红肉及其加工制品而摄入的 N - 羟乙酰神经氨酸 (Nglycolylneuraminic，Neu5Gc) 可被机体识别产物抗-Neu5Gc 抗体，这种抗体诱导机体处于低度炎症状态，并与结直肠癌等炎症的发生相关[1]。Neu5Gc 属于吡喃酮糖衍生物[2]，是N-乙酰神经氨酸 (Nacetylneuraminic acid，Neu5Ac) 被胞苷单磷酸-N-乙酰神经氨酸羟化酶 (CMP-Neu5Ac hydroxylase，CMAH) 羟化的产物，其在人体中不能正常合成，主要依靠外源饮食猪肉、牛肉等红肉摄入[3]。2018年美国牛肉消费量为1 218 万t，中国猪肉消费量为5 540 万t，且依然持续增涨[4]，因此，Neu5Gc 成为重要的食品风险因子[5]，安全高效降低Neu5Gc 在红肉中的含量关乎红肉制品的安全。

目前已有不同方式改变红肉中Neu5Gc 存在状态和含量的研究，蒋芸等[6]通过对猪肉和牛肉进行水煮、微波加热、有机酸腌制处理使得肉样中Neu5Gc 含量降低70%；β-半乳糖苷酶处理肉样后，虽然可使样品中Neu5Gc 含量降低84%，但这种方式仅破坏了Neu5Gc 与生物分子相连的糖苷键，使结合态变成了游离态。鼠伤寒沙门氏菌LT2 唾液酸酶对4-甲基伞形-Neu5Ac 的效果好于4-甲基伞形-Neu5Gc[7]，但该酶无法用于食品工业。梁美莲等[8]以十多种食品加工常见酶制剂处理肉样和Neu5Gc 标准品的效果研究表明，外切型菊粉酶能够降低牛肉中Neu5Gc 的含量，对纯标准品含量降低率可达50.52%[8]，表明菊粉酶有望成为工业化降低游离态Neu5Gc 的有效工具。

利用分子模拟技术获得底物与酶作用的残基信息，有助于酶的理性设计和作用机制解析。杨倩等[9]为提高米根霉α-淀粉酶(ROAmy) 的热稳定性，基于分子动力学模拟结果，对该酶中的3 个氨基酸残基G128、K269 和G393 进行了突变，获得了热稳定性更好的突变体。目前已有X-衍射晶体法和分子动力学模拟对外切菊粉酶作用果糖催化位点解析[10]、同源建模分析菊粉酶功能域和活性位点[11]、菊粉酶对果糖六磷酸和蔗果三糖相互作用模拟的报道[12]。但有关Neu5Gc 与菊粉酶相互作用机制的研究却鲜有报道。本研究对外切型和内切型菊粉酶与Neu5Gc 复合物进行分子动力学模拟，采用密度泛函理论分析相关残基与Neu5Gc 分子间作用的类型、强弱和分布，以期获得菊粉酶与Neu5Gc 相互作用的机制，为菊粉酶在红肉加工中应用提供参考。

1 材料与方法

1.1 试验材料

来自曲霉的外切型菊粉酶标准品(EC 3.2.1.80)、内切型菊粉酶标准品(EC 3.2.1.7)，购自爱尔兰Megazyme 公司；Neu5Gc 标准品，购自美国Sigma 公司。菊粉酶结构从蛋白质PDB 数据库获取，均来源于曲霉，外切型菊粉酶(PDB 编号:1Y9G)，内切型菊粉酶(PDB 编号: 3SC7)。Neu5Gc 分子结构取自Pubchem 数据库(CID: 440001)。

1.2 主要仪器与设备

Waters Xevo TQ 沃特世三重四级杆液相色谱质谱联用仪，美国Waters 公司；Mill-Q 超纯水仪，美国Milipore 公司。SYBYL-X(2.0)药物设计平台、开源分子动力学程序Gromacs (2 016.4)、Gaussian View6、Gaussian16 (Revision A.03)[13]程序由贵州大学云计算平台支持。Discovery studio 3.5 客户端(蛋白去水、去杂原子和对接空腔参数记录)、对接程序AutoDock-Vina、波函数程序Mulitwfn3.6[14]、酶底复合物二维残基相互作用可视化程序Ligpot＋、氨基酸残基虚拟突变工具Swiss PDB Viewer[15]、等值面图绘制程序VMD(1.9)[16]，均由个人电脑提供计算平台。

1.3 试验方法

1.3.1 菊粉酶作用效果检测 取1 mL 10 mg·mL-1的标准Neu5Gc 溶液于液相瓶中，加入0.8%(w/w)菊粉酶，置于最适温度50℃酶解30 min，酶解完成后用0.22 μm 膜过滤酶解液，备用。采用液相色谱质谱联用仪进行菊粉酶作用效果检测，根据文献[17]的方法略有修改。取1 μL 酶解液进样分析，色谱条件:ACQUITY UPLC C18 色谱柱(1.7 μm，2.1 mm×100 mm)，流动相A:100%乙腈，流动相B:0.1% 甲酸水溶液，梯度洗脱(5% A 0～3 min，10% A 3 ～5 min，100%A 5～6 min，5% A 6～7 min)，流速0.2 mL·min-1，柱温25℃，进样量1 μL。质谱条件: 离子源为电喷雾化离子源(ESI)，源温度100℃，探针温度450℃，毛细管电压3 000 V，锥孔电压19 V，脱溶剂气流速1 000 L·hr-1，负离子方式检测，扫描方式为多重反应监测(multiple reaction monitoring，MRM)。定量分析离子对323.8/116.0 时，碰撞气能量为17 V；定量分析离子对323.8/185.8 时，碰撞气能量为15 V。

1.3.2 Neu5Gc 与菊粉酶相互作用 参照文献[18]的方法，采用分子对接、氨基酸残基虚拟突变、分子动力学模拟研究Neu5Gc 与菊粉酶的分子间相互作用。为更加全面的考察Neu5Gc 与菊粉酶相互作用过程中涉及的蛋白残基，同时采用AutoDock-Vina 和SYBYL(2.0)对接方法进行分子动力学模拟，研究实际水环境下的对接复合物稳定性。

AutoDock-Vina 对接: 对受体菊粉酶和配体Neu5Gc 均以AutodockTool-1.5.6 打开，进行定义原子类型和添加Casteiger 电荷处理，并保存为pdbqt 格式。从PDB 蛋白数据库得到的菊粉酶结构含小分子配体，这些小分子配体构成了对接时的活性空腔，为获得Neu5Gc 与菊粉酶不同配体空腔对接的最佳位置，根据Discovery studio 3.5 客户端中显示的配体空间参数，采用AutoDock-Vina 进行不同配体位置的对接。设置对接时的最大打分构象数目为9，对接盒子大小设置根据空间参数适当调整。对接完成后根据结合能大小和均方位移(root mean square deviation，RMSD)值作为最佳结合构象选择依据。对接后的最佳复合物结构用Ligpot＋进行氨基酸二维相互作用分析。

Surflex-Dock (SFXC)-Geom 对接:受体菊粉酶加氢后删除配体和水分子以及杂原子，加Casteiger 电荷，在Tripos 力场下采用最陡下降法进行1 000 步的能量最小化处理，其余参数均为默认。对接空腔由配体坐标产生，在其周围氨基酸0.5 nm 范围半径内产生活性空腔。最终对接配体可能的构象数目设定为20种，配体的RMSD 最小控制在0.5。取对接后打分最高者进行对接分析。

Gromacs 动力学模拟:依据结合能与打分排序对以上对接后的复合物进行筛选，对最佳Neu5Gc-菊粉酶复合物在水环境中进行动力学模拟以考察其结合稳定性。受体分子菊粉酶的力场选择AMBER99-ILDN，水盒子采用TIP3P 三点水模型，体系大小为盒子边界距离蛋白文件2.5 nm。配体分子Neu5Gc 采用Amber Tools 17 中的GAFF 力场，产生itp 限制势文件和gro拓扑文件，使用Gromacs 中的genrestr 工具生成体系距离限制势文件，使用-DPOSRES 关键词限制配体的位置。对体系加入抗衡离子平衡电荷。在模拟参数文件mdp 中先用最陡下降法进行2 000 步，步长为0.01 ps的能量最小化，邻居列表生成选择cutoff-scheme ＝Verlet，静电作用为库伦力，非键相互作用采用PME 截断方式，LINCS 算法约束成键，每50 步输出一次能量信息。然后进行500 ps 的控温、控压限制性动力学直到平衡温度298.15 K，热浴设定为velocity-rescale，控温的时间常数为0.2 ps，压浴设定为Berendsen，采用各项同性控压，控压时间常数为0.5 ps，只约束氢键，该过程每500 步输入一次能量信息，每1 000 步输出一次坐标信息。限制性动力学平衡后，对配体和受体设置索引和分组，放开限制进行100 ns 的成品模拟，设定每步0.002 ps，每1 000 步输出一次能量信息。成品模拟平衡后，选择平衡末端80～100 ns 处轨迹文件，每隔10帧从轨迹文件中提取平均构象，采用g＿mmpbsa 工具，基于分子力学泊松-波尔兹曼表面积(molecular mechanics poisson-boltzmann surface area，MM-PBSA)法进行蛋白和配体结合自由能△Gbind计算。△Gbind值越低说明受体和配体间的亲和力越高。溶剂中蛋白和配体结合自由能按照公式进行计算[19]:

式中，Gcomplex是总的自由能；Gprotein和Gligand是溶剂中独立的自由能。其中，溶剂中独立的自由能GX按照公式进行计算:

式中，X代表蛋白或配体或复合物；TS 是真空中构象变化导致的熵变对自由能的贡献，由于构象熵对本研究变化影响不大，故不考虑；EMM是真空中的分子力学能；Gsolvation代表溶剂化能；按照公式进行计算:

式中，Ebond代表成键作用，如键长、键角、旋转等；Enobonded代表非键作用力，如范德华力和静电作用力；Gpolar代表极性溶剂化能，由求解波兹曼方程得到；Gnonpolar由基于溶剂可及表面积SASA 模型、溶剂可及体积SAV 和WCA 模型计算。

计算完成后，对动力学过程中重要残基进行能量分解，以获得残基对结合自由能的贡献。为考察模拟过程的稳定性，进行蛋白和复合物的RMSD、残基的RMSF 分析。

1.3.3 残基虚拟突变 当外切菊粉酶1Y9G 的活性残基ASP189 用ASN 取代时，酶活性降低[9]；以外切菊粉酶1Y9G 为模板同源建模构建内切型菊粉酶，当活性残基ASP460 突变为ALA 后，突变体催化活性丧失[20]。内切型菊粉酶3SC7 定点突变表明ASP42 对催化活性有重要作用[21]。采用Swiss PDB Viewer 工具，选择将外切型菊粉酶1Y9G 的活性残基ASP41 突变为ASN41；将内切型菊粉酶3SC7 的活性残基GLU233 突变为ALA233。对突变后的菊粉酶与Neu5Gc 以AutoDock-Vina 对接并产生复合物结构，进行分子动力学模拟考察稳定性和活性位点残基变化。提取菊粉酶1Y9G 和3SC7 与Neu5Gc 动力学模拟平衡后最后一帧，在M062X/6-311G(d，p)水平下对活性残基-Neu5Gc 复合物以分子中原子理论(atom in molecules，AIM) 进行键临界关键点(bond critical point，BCP) 分析，用独立梯度模型(independent gradient model，IGM)考察分子间弱相互作用。

1.4 数据处理与分析

试验数据均用平均值表示，平行3 次。采用SPSS 23. 0 软件进行统计分析，数据显著性差异分析采用t检验，显著性水平为P＜0.05，采用OriginPro 2016 软件制图。

2 结果与分析

2.1 菊粉酶对Neu5Gc 的作用效果

通过液相质谱(liquid chromatography-mass spectrometry，LC-MS) 对2 种酶在最适温度下与Neu5Gc 标准品作用后滤膜液进行Neu5Gc 含量检测。由图1可知，外切型菊粉酶对Neu5Gc 标准品中Neu5Gc 含量的平均降低率为55.33%，内切型菊粉酶对Neu5Gc 标准品中Neu5Gc 含量的平均降低率为43.81%，两组间差异不显著(P>0.05)，但均显著高于对照组(1.73%)。表明外切型菊粉酶对吡喃糖衍生物Neu5Gc 的降低效果略优于内切菊粉酶。这可能与菊粉酶对吡喃糖衍生物，如蔗糖的水解能力差异有关；外切型菊粉酶能将蔗糖转化为单糖中的果糖；而大多数内切酶转化蔗糖的效率较低[22-23]。

2.2 Neu5Gc 与菊粉酶不同位点分子对接

由于2 种菊粉酶对Neu5Gc 作用效果相近，为探究作用过程中二者间活性残基的差异，首先采用分子对接对二者进行初步探究。用Discovery studio 3.5 客户端显示菊粉酶结构，外切型菊粉酶1Y9G 共5 个配体，配体1 由NAG1001、NAG1002 组成，配体2 由NAG2001 组成，配体3 由NAG3001 组成，配体4 由NAG5001 组成，配体5 由FRU801 组成。3SC7 内切型菊粉酶X 链共4 个配体，配体1 由MAN517 组成，配体2 由MAN518 和MAN519 组成，配体3 由NAG750、NAG751 和BMA752 组成，配体4 由EPE520 组成。对菊粉酶不同配体位置与Neu5Gc 进行AutoDock-Vina对接，结果如表1、2 所示，对接最佳复合物处Neu5Gc-菊粉酶二维活性残基作用结果如图2所示。Surflex-Dock 对接的结果如图3所示。

图1 菊粉酶对Neu5Gc 标准品的作用效果Fig.1 Effect of inulinase on Neu5Gc standards

由表1可知，1Y9G 在配体5 处对接结合能最低，为-30.96 kJ·mol-1，因此Neu5Gc 与1Y9G 在配体5 处对接结果最佳。由表2可知，3SC7 在配体2 处对接结合能最低，为-26.77 kJ·mol-1，因此Neu5Gc 与3SC7在其配体2 处对接结果最佳。由图2可知，Neu5Gc 与1Y9G 配体5处的残基ASN40、GLN57、SER103、GLN149、ASN186、ASP189、ASN265 和TRP335 相互作用产生10 个氢键，与残基ASP41、PHE102、ARG188 和GLU241 产生疏水相互作用。由图3中SYBYL 对接结果可知，Neu5Gc 与残基 GLN14、ASN40、GLN57、SER103、GLU241、TRP335 产生8 个氢键。外切型菊粉酶1Y9G 催化果糖的活性位点残基为ASN40、ASP41、TRP65、SER103、ARG188、ASP189、GLU241 和CYS242[10]。对比1Y9G 活性位点可知，Neu5Gc 均可与外切型菊粉酶残基CYS242 外的其余活性残基发生相互作用，且其主要与ASP 型、GLN 型残基相互作用。

Neu5Gc 与3SC7 配体2 周围残基ASN42、GLU43、GLU233、ASN320、GLY323 产生7 个氢键。与残基VAL23、TRP40、PHE99、THR255、ASN265、ASP298 产生疏水相互作用。SYBYL 对接结果显示，与3SC7 的残基GLU258、ASN265、ARG295、ASN320、GLY323、ASP440 产生8 个氢键。内切型菊粉酶催化果糖和戊酮糖时共有的残基为TRP40、GLU43、GLN59、ASN61、TRP67、GLY68、ILE70、PHE99、ARG175、GLU233、ASN320 和GLY323[24]。对比3SC7 的活性位点可知，Neu5Gc 可与内切型菊粉酶活性残基TRP40、GLU43、PHE99、GLU233、ASN320 和GLY323 相互作用，表明其主要与GLU 型残基作用。已知在外切和内切型菊粉酶活性残基中，ASP 和GLU 分别是重要的催化残基[25]。综合以上对接后相互作用残基分析可知，Neu5Gc 主要与外切型菊粉酶活性残基中的ASP 残基和内切型菊粉酶活性残基中的GLU 残基相互作用。这与Singh 等[12]采用分子对接得到外切型菊粉酶和果糖- 6 -磷酸主要结合残基为ASP22、ASP128 和ASP179，内切型菊粉酶与酮糖主要结合残基为GLU34和GLU200 一致；也和动力学模拟果糖-6-磷酸和酮糖分别与外切型菊粉酶ASP285、GLN350 产生氢键作用的结论一致[26]。

表1 1Y9G 与Neu5Gc AutoDock-Vina 对接结果Table1 1Y9G and Neu5Gc AutoDock-Vina docking results

表2 3SC7 与Neu5Gc AutoDock-Vina 对接结果Table2 3SC7 and Neu5Gc AutoDock-Vina docking results

2.3 菊粉酶与Neu5Gc 分子动力学模拟

2.3.1 动力学模拟稳定性评价 为获得2 种菊粉酶与Neu5Gc 对接后最佳酶底复合物在水环境中的稳定结合状态，提取最佳打分复合物进行100 ns 动力学模拟。为验证动力学模拟中菊粉酶催化残基对结合Neu5Gc 并稳定其构象的贡献，对菊粉酶催化活性残基进行虚拟突变。采用蛋白骨架波动、配体与受体蛋白的偏移程度、蛋白残基波动进行模拟稳定性评价，结果分别如图4、5、6 所示。波动数值越大表明稳定性越低。由图4可知，1Y9G 突变后残基蛋白稳定性轻微提高，但突变后复合物稳定性小于突变之前。由图5可知，突变后蛋白残基运动波动性相似，但均大于晶体结构的残基波动，这与蛋白稳定性是一致的。对3SC7突变活性残基后，蛋白稳定性变小，残基波动明显高于晶体试验值，复合物稳定性明显减弱。

图2 Neu5Gc 与菊粉酶1Y9G 和3SC7 AutoDock-Vina 对接结果二维模式图Fig.2 Two-dimensional interaction model of Neu5Gc and inulinase 1Y9G and 3SC7 by AutoDock-Vina docking

图3 Neu5Gc 与菊粉酶1Y9G 和3SC7 Surflex-Dock 对接结果三维模式图Fig.3 Three-dimensional interaction model of Neu5Gc and inulinase 1Y9G and 3SC7 by Surflex-Dock docking

2.3.2 动力学模拟结合模式 为获得稳定对接状态下Neu5Gc 与菊粉酶活性残基相互作用过程中二者间的空间结合模式，对模拟平衡后20 ns 轨迹每隔10 ps进行提取，统计模拟过程中活性残基与Neu5Gc 作用产生的氢键个数和贡献权重，选取模拟平衡后的构象进行结合模式分析。由图7可知，Neu5Gc 与1Y9G 的残基ASP41、GLN57、GLU241、ASN265 和TRP335 产生氢键。当突变1Y9G 的残基ASP41 为ASN41 后，产生氢键的残基为ASN40、GLN57、TRP65 和GLN149，其中活性残基仅ASN40、TRP65 与Neu5Gc 稳定相互作用，这是由于突变后ASN41 不能亲核攻击底物的异构碳原子形成共价中间体，使得其活性丧失。根据外切型菊粉酶活性位点可知，在水溶液环境中Neu5Gc 可与外切型菊粉酶关键活性残基ASP41、GLU241 以及GLN57 稳定相互作用。由于残基ASP41 和GLU241分别在菊粉酶催化底物过程中充当亲核试剂、催化酸碱对[27]，表明外切型菊粉酶可对Neu5Gc 时行催化作用。

图4 菊粉酶1Y9G(A)和1Y9G－Neu5Gc 复合物(B)动力学模拟蛋白骨架均方根位移Fig.4 Dynamics simulation RMSD of inulinase 1Y9G (A) and 1Y9G－Neu5Gc complex(B) protein skeleton

图5 菊粉酶1Y9G (A) 和3SC7 (B)动力学模拟残基均方根涨落Fig.5 Dynamics simulation RMSF of inulinase 1Y9G(A) and 3SC7 (B) protein residues

Neu5Gc 与3SC7 产生氢键的残基为ASN42、GLU43、ASN265、ASP298 和ASN320，其中催化活性残基ASP298 和GLU43 与晶体结构分析结果一致[28]。由于Neu5Gc 是N-乙酰葡萄糖的衍生物，且有学者以3SC7 为模板和N -乙酰葡糖胺对接，表明残基ASP285、PHE352、THR288、VAL284、LEU343、MET289、GLN350、LEU343、HIS351、PRO283、THR271可与其稳定相互作用[29]。对比可知，对接获得的残基相互作用信息与水环境中动力学模拟得到的作用残基结果差异较大，表明与Neu5Gc 产生较稳定氢键作用的残基对结合过程贡献较大。当突变GLU233 为ALA233 后，Neu5Gc 仅与残基ASP176、ARG175 和GLU173 产生氢键作用。根据内切型菊粉酶活性位点可知，Neu5Gc 可与内切型菊粉酶关键残基GLU43、ASN320 稳定相互作用。突变内切型菊粉酶3SC7 的残基GLU233 为ALA233 时，Neu5Gc 仅与3SC7 的活性残基ARG175 产生氢键作用，这是由于GLU 较大的侧链能够造成底物结合过程的崩溃，由此影响了活性空腔的酸性环境。

2.3.3 结合过程残基能量分解 为获得水环境下Neu5Gc 与菊粉酶结合过程中分子间相互作用类型、活性残基稳定Neu5Gc 构象到活性空腔中的能量贡献，对动力学模拟平衡后20 ns 轨迹每隔100 ps 进行提取，采用g＿mmpbsa 工具对酶和底物结合能贡献进行波尔兹曼能量分解，结果如表3所示。并对模拟平衡时菊粉酶与底物Neu5Gc 相互作用贡献较大的残基进行能量贡献分解，结果如图8所示。

图6 Neu5Gc 相对于菊粉酶3SC7(A) 和3SC7-Neu5Gc 复合物动力学模拟均方位移Fig.6 Dynamics simulation RMSD of Neu5Gc on inulinase 3SC7(A) and 3SC7-Neu5Gc(B)

图7 Neu5Gc 与菊粉酶1Y9G(A)、3SC7(B)、1Y9G 突变后(C)和3SC7 突变后(D)动力学模拟平衡结合模式图Fig.7 Binding model of Neu5Gc and inulinase 1Y9G (A)，3SC7 (B)，mutated 1Y9G (C) and mutated 3SC7 (D)under dynamics simulation equilibrium

表3 Neu5Gc 与菊粉酶体系动力学平衡后结合能贡献分解Table3 Binding energy contribution decomposition of Neu5Gc and inulinase system under dynamics simulation equilibrium /(kJ·mol－1)

由表3可知，对Neu5Gc 与菊粉酶结合过程起阻碍作用主要是极性溶剂化能，起促进作用的是范德华力、静电相互作用和溶剂可及表面积。未突变时，外切菊粉酶1Y9G 与Neu5Gc 的结合能明显小于内切型菊粉酶3SC7，表明Neu5Gc 优先与外切型菊粉酶相互作用。当突变1Y9G 的催化残基后，范德华力增强，这可能是由于天冬酰胺使得结合位点的电荷由原来的负电荷变为正电荷，静电相互作用减少。这与Torabizadeh等[30]通过以1Y9M 建模的菊粉酶活性残基经磷酸吡哆醛修饰后，负电荷增加，静电作用增强的结果相似。突变3SC7 的GLU233 为ALA233 后，总结合能增强，但静电相互作用对结合过程的促进作用明显减少，导致复合物稳定性明显降低。由图8可知，结合过程中主要相关残基的能量项主要由对结合过程起促进作用的分子力学能和起阻碍作用的极性溶剂化能组成。从有利于结合过程的分子力学能变化可知，1Y9G 中促进与Neu5Gc 结合的关键残基是ASP41、GLU241 和TRP335，而当突变残基ASP41 后主要为TRP65 残基。3SC7 中促进与Neu5Gc 结合的关键残基是GLU43、ASN265 和ASP298，当突变GLU233 后主要为ASP176残基。表明突变关键活性残基后，菊粉酶对于Neu5Gc稳定到活性空腔的静电能将下降，进一步证明了Neu5Gc 与菊粉酶主要以催化活性残基相互稳定作用。

2.3.4 AIM 和IGM 分析 为考察动力学平衡后Neu5Gc 与活性残基形成的氢键强度和分布，提取1Y9G 和3SC7 与Neu5Gc 动力学模拟平衡后关键结合残基，选择对弱相互作用计算较好的M062X 泛函，在6-311G(d，p)基组水平下产生波函数文件对临界键关键点BCP 进行AIM 分析，结果如表4所示。采用独立梯度约化泛函IGM 对氢键投影进行等值面分析，结果如图9所示；分子间和分子内相互作用散点图如图10所示。

图8 Neu5Gc 与菊粉酶1Y9G(A)、3SC7(B)、1Y9G 突变后(C)和3SC7 突变后(D)动力学平衡时残基能量分解Fig.8 Residual energy decomposition of Neu5Gc and inulinase 1Y9G (A)，3SC7 (B)，mutated 1Y9G (C)and mutated 3SC7 (D) under dynamics simulation equilibrium

表4 Neu5Gc 与菊粉酶1Y9G、3SC7 活性残基分子间氢键AIM 分析(能量项ρ、▽2ρ、V、G、HBCP、EHB)Table4 AIM analysis of intermolecular hydrogen bonding between Neu5Gc and inulinase 1Y9G，3SC7 active residues (energy items ρ、▽2ρ、V、G、HBCP、EH B)

在AIM 理论中，体系坐标函数的二阶鞍点称为键临界点。键临界点处的电子密度、拉普拉斯密度与键的强度和类型有关[31]。从键临界点搜索结果可以看出，1Y9G 和3SC7 活性残基与Neu5Gc 均产生4 个氢键临界点。由Rozas 等[32]的氢键强度划分标准:弱氢键(EHB＜50.20 kJ·mol-1，▽2(rBCP) > 0，G(rBCP) ＋V(rBCP) >0)，中等氢键(50.20 ＜EHB＜100.4 kJ·mol-1，▽2(rBCP) > 0，G(rBCP) ＋V(rBCP) ＜0)，强氢键(EHB>100.4 kJ·mol-1，▽2(rBCP) ＜0，G(rBCP) ＋V(rBCP) ＜0)。由表4可知，外切型菊粉酶1Y9G 活性残基与Neu5Gc 作用的最大氢键强度为-9.65 kJ·mol-1，内切型菊粉酶3SC7 活性残基与Neu5Gc 作用的最大氢键强度为-11.53 kJ·mol-1，可知2 种菊粉酶与Neu5Gc均以弱氢键相互稳定结合，其中，就氢键强度总和而言，外切型菊粉酶稍高于内切型菊粉酶的作用效果，这与结合自由能分解得到Neu5Gc 更易于和内切型菊粉酶在水环境下稳定结合的结果是一致的。

约化密度梯度(reduced density gradient，RDG)函数是基于电子密度拓扑的非共价相互作用分析方法，电子密度梯度降低的区域与非共价相互作用区域相关，该区域用等值面图可以将非共价作用的位置分布和类型进行可视化[33]，以S 作为弱相互作用的等值面函数。

图9 Neu5Gc 与菊粉酶1Y9G(A)、3SC7(B)活性残基分子间作用IGM 等值面图Fig.9 Intramolecular IGM isosurface map of Neu5Gc with inulinase 1Y9G (A) and 3SC7 (B) active residues

图10 Neu5Gc 与菊粉酶1Y9G、3SC7 活性残基分子间作用IGM 散点图Fig.10 Intramolecular IGM scatter plot of Neu5Gc with inulinase 1Y9G and 3SC7 active residues

式中，ρ 代表电子密度，▽ρ 代表电子密度梯度的模量。在IGM 理论中，当SIGM基于｜▽ρ(r)IGM｜时，SIGM和S 最大的差值代表键临界点处的RDG 值，定义δg表示相互作用区域，划分为片段内相互作用的δgintra和片段间相互作用的δginter。其表达式如下:δg(r)＝｜▽ρ(r)IGM｜-｜▽ρ(r)｜，δg(r)＝δg(r)inter＋δg(r)intra。将δg 描述符以符号函数sign(λ2)ρ 投影。当sign(λ2)ρ值为-0.04 时，δginter值显示明显的峰，代表了分子间相互作用的氢键；当sign(λ2)ρ 值为正值时，代表了空间排斥作用；当sign(λ2)ρ 值在-0.3 附近，δgintra有个较大峰，代表了分子内相互作用。默认的sign(λ2)ρ 值在-0.05～0.05 范围内。等值面图中蓝色越强表示吸引作用越强，绿色的代表弱的范德华力作用[34]。

由图9可知，与Neu5Gc 稳定结合的菊粉酶活性残基和Neu5Gc 产生氢键的氮氧原子与受体之间有明显的等值面分布。从等值面分布上看，1Y9G 活性残基与Neu5Gc 产生的蓝色区域明显少于3SC7 的活性残基，表明分子间静电作用力弱于3SC7，这与结合能分解得到其静电作用能较低相一致。由颜色分布可知范德华力是主要作用力，氮原子为供体产生的较强氢键仅为1Y9G 的GLN57 处，这与AIM 分析的结论一致。由3SC7 活性残基和Neu5Gc 弱相互作用的等值面分布图可知，产生氢键的位置主要为供体羟基氧原子和受体氧之间，与AIM 分析结果相一致。同样，散点图10中观察到菊粉酶与Neu5Gc 分子间存在的弱氢键相互作用峰(sign(λ2)ρ -0.04 处)。以上结果表明Neu5Gc 羟基氧和酰胺氮与菊粉酶活性残基提供的亲核性质是导致该结合过程稳定存在的原因之一。

3 讨论

菊粉酶广泛用于功能性低聚糖的制备，但由于菊粉酶催化活性域具有一定保守性，不同底物与不同类型菊粉酶之间相互作用的过程和机制有一定共性。菊粉酶可以作用于单糖衍生物[12]，对酮糖也有一定的降低效果[26]，推测菊粉酶具有消减吡喃糖衍生物Neu5Gc 的潜在功能。因此，基于前期对多种食品酶制剂作用于Neu5Gc 的筛选结果进行深入研究，采用液相质谱检测以克服Neu5Gc 液相衍生化检测干扰，其中，外切菊粉酶结果(55.33%)较液相法检测结果(50.52%)提高了4.81%[7]，这是由于衍生化试剂1，2-二氨基-4，5-亚甲二氧基苯(DMB) 会与菊粉酶活性中心竞争性结合，干扰检测。外切型和内切型菊粉酶活性域的差异性导致其对相同底物的作用效果不同，因此本研究观察到内切型菊粉酶对Neu5Gc 含量的降低率较外切型菊粉酶低6.71%左右。Neu5Gc 虽然不是菊粉酶的高效底物，但高达50%以上的降低率表明菊粉酶对Neu5Gc 具有一定的消减能力。分子动力学模拟广泛用于酶底物结合过程细节的讨论，故采用该方法研究2 种菊粉酶作用于Neu5Gc 的效果和差异机制。

目前对菊粉酶有关的分子动力学模拟研究多集中在理性设计和底物识别方面[35]，但对其于分子间弱相互作用的讨论不够深入。国内武希箭[36]采用分子动力学模拟研究内切型菊粉酶对菊粉的分子间相互作用，考察了结合能、残基突变情况，发现主导菊粉酶发挥作用的关键分子间作用力是静电力和范德华力，且活性位点残基侧链减少后，对底物稳定到活性口袋的作用减弱，这与本研究结论一致，但其未定量研究活性残基与底物结合的氢键强度和分布情况[36]。目前对于酶和底物分子对接的研究多采用单一平台进行分析，形成复合物后采用动力学模拟进行验证，但由于不同分子对接软件的打分函数和评价规则不同，导致一些潜在的酶活性残基不能被识别[37]。本研究综合利用AutoDock-Vina 的高并行性和SYBYL 平台的高精度分子力场优势进行对接，残基与底物分子互作结果也表明了不同对接平台在识别活性残基方面的差异，采用多平台对接可获取更多酶底物结合过程信息。本研究在水环境下对于Neu5Gc 和菊粉酶对接复合物的动力学模拟表明，Neu5Gc 能够与菊粉酶活性残基以静电作用力、范德华力和疏水作用力稳定结合，氨基酸虚拟突变也证实了菊粉酶催化活性残基ASP 和GLU 的重要贡献，这与大多数研究结论一致[27]。但本研究通过分子中原子理论对活性残基和Neu5Gc 复合物的分子间氢键定量分析发现，外切型菊粉酶1Y9G 与Neu5Gc的分子间氢键强度较内切型菊粉酶3SC7 高18.7 kJ·mol-1左右，这可能是二者对Neu5Gc 含量降低效果不同的原因之一。氢键强度差异对于酶催化过程中过渡态的形成和反应势垒具有一定影响，随着量子力学/分子力学(quantum mechanics/molecular mechanics，QM/MM) 法在研究糖苷酶水解、糖基转移酶的反应机制方面的应用[38]，活性残基对底物作用过程中热力学和动力学变化逐渐被解析，但该过程涉及复杂力场参数的新建和优化，较为困难。对于开发高效消减Neu5Gc 及其他食源性风险物质的食品加工酶类，在获知酶底物互作信息基础上，利用分子生物学技术即可为实现酶的理性设计提供参考和支持。

4 结论

体外菊粉酶降低Neu5Gc 标准品含量的试验结果表明，外切型菊粉酶和内切型菊粉酶分别对Neu5Gc降低效果最佳为55.33%和43.81%。选择两类菊粉酶(PDB ID: 1Y9G 和3SC7)与Neu5Gc 稳定结合的分子对接复合物进行分子动力学模拟和弱相互作用，结果表明，Neu5Gc 可以与外切型菊粉酶1Y9G 的活性残基ASP41 和GLU241 通过分子间氢键产生稳定的相互作用，内切型菊粉酶3SC7 的活性残基ASP298 和GLU43 对结合Neu5Gc 到其活性口袋过程有重要贡献。突变1Y9G 中的关键活性残基ASP41 为ASN41，3SC7 的活性残基GLU233 为ALA233 后的结果显示，活性残基中ASP 的亲核性和GLU 的酸碱对性质对菊粉酶结合Neu5Gc 并导致其含量降低过程具有重要作用。残基能量分解、AIM 分析和IGM 等值面分析表明，Neu5Gc 通过羟基氧原子和酰胺氮原子与菊粉酶活性残基发生静电作用主导的分子间弱相互作用，稳定占据到菊粉酶的活性位点。红肉是人体多种营养物质的重要来源，本研究通过解析菊粉酶作用Neu5Gc 的活性残基信息，为菊粉酶的理性设计和定向进化提供了新的理论参考。