基于分子动力学模拟的扩展青霉棒曲霉素MFS 蛋白转运机制研究

杨 琪，王艳玲

（兰州理工大学生命科学与工程学院，甘肃兰州 730050）

青霉病是全球最常见且经济上最重要的采后水果腐烂病之一[1]，其中扩展青霉（Penicillium expansum，P.expansum） 是引起果实青霉病或软腐病的主要病原菌[2]。P. expansum可在水果作物（苹果、梨、柑橘和樱桃等）及其衍生品、加工产品的采后、贮藏、运输和加工期间侵袭[3-7]。同时，P. expansum还可产生具有危害的次级代谢产物——棒曲霉素（Patulin，PAT）[8-10]，通过诱导DNA 链断裂、氧化损伤以及DNA 链之间相互交联、产生核质桥等机制影响微核形成，从而引起DNA 结构性损伤[11]，产生遗传毒性、细胞毒性和免疫毒性作用等[12]。

20 世纪80 年代，Puel 等提出PAT 的生物合成途径是一个酶级联反应[13-14]。近几年，Li 等[15]鉴定得到参与PAT 生物合成的基因簇，包括15 个基因（PatA～PatO），具体合成途径见图1。其中PatC基因编码MFS 超家族转运蛋白（Major Facilitator Superfamily，MFS）。MFS 的主要功能是将PAT 的前体E-ascladiol 转运至胞外，最终合成PAT[16]，但具体的转运机制尚不清楚。研究表明，E-ascladiol 对来自人类肝、肾、肠和免疫系统的细胞无毒性，因此，增加E-ascladiol 累积以减少棒曲霉素含量是控制棒曲霉素风险的有效手段[17]。目前，MFS 转运蛋白包括105 个家族，是最大的次级跨膜转运蛋白超家族[18]。MFS 转运蛋白大多数都含有400～600 个氨基酸残基[19]，具有12、14 或24 个不规则跨膜α-螺旋（TMSs）包围形成深中央亲水腔[20-21]，催化底物在两个方向上进行跨膜运输，可以是运输单个底物的单转运子，可以是将底物和偶联离子运输在一起的同转运子，也可以是按相反方向依次运输两个不同分子的反转运子（图2a）。根据交替开放模型，MFS 能够在向胞内开放，闭塞和向外开放构象之间变化（图2b），转运蛋白上的底物结合位点暴露在细胞膜的不同侧，底物从细胞膜的一侧结合到蛋白上，而从另一侧被释放（图2c）。此外，MFS 转运蛋白还可以间接调节真菌的内部pH和应激反应机制[22]，参与多种真菌毒素的输出。如拟枝孢镰刀菌（Fusarium sporotrichioides）的TRI12和大豆紫斑菌（Cercospora kikuchii）的CFP[23-24]，两者均编码MFS 转运蛋白，分别参与真菌毒素-单端孢霉烯（Trichothecenes）和尾孢菌素（Cercosporin）的分泌。植物病原体，玉米圆斑真菌（Cochliobolus carbonum）的TOXA也与HC 毒素的宿主特异性毒素的分泌有关[25-26]。TOXA、Tri12及CFP缺失突变体都显示出毒素产生的急剧减少，Tri12缺失突变体也观察到生长迟缓，但CFP缺失突变体与野生型菌株一样正常生长，只是对大豆的毒力降低[26-28]。此外，在携带毒素合成基因的基因簇中也发现了一些编码MFS 的基因[29]。然而，它们在运输特定代谢产物中的作用仍待阐明。可惜的是，对MFS 转运机制的研究甚少，目前只解析得到几个细菌蛋白的完整拓扑结构。相关PAT 毒素输出的研究更是没有。因此致力于探索PatC 的转运机制，有利于苹果青霉病病原菌的防治的研究。

图1 PAT 的生物合成途径[30]Fig.1 PAT biosynthetic pathway[30]

图2 MFS 蛋白的运输机制[31]Fig.2 Mechanism of transport of MFS[31]

随着计算生物学的发展，越来越多的研究者采用分子动力学模拟生物大分子在三维空间的运动状态，从而分析其发挥生物功能的作用机制。本研究采用生物信息学及分子动力学模拟方法探索扩展青霉棒曲霉素PatC 的转运机制，研究其与底物分子Eascladiol 的结合位点，并针对PRO188 定点突变，分析复合体及突变体的RMSD、Rg、RMSF 等变化，以期为研究该蛋白的结构和转运机制奠定基础。

1 材料与方法

1.1 生物信息学分析

通过美国国家生物技术信息中心（National Center for Biotechnology Information，NCBI）搜索Penicillium expansumT01（登录号：AYHP00000000.1），下载Penicillium expansumT01 中全部527 段鸟枪法测序结果。利用DNAMAN 生物信息软件BLAST其他亚种Penicillium expansum strainNRRL 35695（登录号：KF899892.1）中PatC基因及其编码的氨基酸信息，结果与Penicillium expansumT01 中PatC基因序列完全一致。采用Mega 5.0 邻近法（Neighbor Joining，NJ）进行氨基酸序列同源性分析及构建系统进化树；通过在线网站分析PatC 的理化性质、跨膜结构域、蛋白质功能结构域、二级结构及三级结构等，具体所用网站见表1。

表1 生物信息学相关分析工具Table 1 Bioinformatics correlation analysis tools

1.2 分子对接

配体分子和蛋白三维结构的优化：分子对接用于验证活性成分与关键靶点之间的结合活性。通过AlphaFold 蛋白质三维结构数据库在线搜索PatC 的三维结构，再对其三维模型导入AutoDocktools（v1.5.6）进行加氢、计算电荷、分配电荷、指定原子类型等优化处理；从Pubchem 数据库下载小分子或蛋白的3D 结构，导入ChemBio3D Ultra 14.0 进行能量最小化。将优化好的小分子导入AutodockTools-1.5.6 进行加氢、计算电荷、分配电荷等处理。

分子对接：AutoDock Vina1.1.2 是一个采用半柔性对接方法运行的程序，对接精度高达78%，被用来作为本研究的分子对接程序[32]，分子对接流程如图3。

图3 分子对接流程示意图Fig.3 Schematic diagram of the molecular docking process

1.3 分子动力学模拟

先将Wild PatC 和E-ascladiol 对接的最佳结合构象，作为初始复合物结构进行MD 模拟。随后使用蛋白定点突变软件，构建新型PatC 突变体结构（P188A）观察突变后的蛋白与底物对接后复合物的活性。

具体操作如下：使用分子动力学软件Gromacs2020.6 进行模拟。复合物小分子在GAFF 力场下生成拓扑文件，蛋白在AMBER99SB-ILDN 力场生成相应参数化文件。以复合物为中心，添加SPC/E 水模型分子，生成立方体水盒子，向复合物水盒子里添加4 个CL-离子保持模拟体系的电中性。在体系完成后，使用最陡能量下降法进行50000 步能量最小化，使体系充分稳定。采用正则系综（NVT），Berendsen 热浴法控制温度为298.15 K，模拟50000 步，步长2 fs。采用等温等压系综（NPT），Parrinello-Rahman 控压，1 个大气压，模拟50000 步，步长2 fs。进行mdrun，步长2 fs，模拟10000000，总计20 ns 的模拟，对模拟后的数据计算MMPBSA。

1.4 数据处理

利用Discovery studio 2016 对分子对接结果进行可视化，且利用PyMOL（4.3.0 版）软件绘制蛋白与分子之间相互作用的氨基酸残基。分子动力学模拟使用软件gromacs 分析各组参数变化，最终通过程序统计结果文件计算出体系各种宏观性质和微观性质，并使用xmgrace 绘制相应的数据图，如RMSD、RMSF、SASA 及氢键变化等。PAT 生物合成图（图1）使用在线网站FigDraw（https://www.figdraw.com）绘制，图2 使用Adobe Illustrator 软件绘制。

2 结果与分析

2.1 构建系统进化树

选取与PatC 同源性较高的多条已知菌株的氨基酸序列（表2），利用MEGA5.0 软件构建NJ 系统进化树（图4）。所选的11 个物种序列与PatC 的氨基酸序列一致性较高。系统进化树结果表明，Penicillium expansum（扩展青霉）与Penicillium griseofulvum（灰黄青霉）处于同一分支，亲缘关系较近。此现象可能是因为两者的PAT 合成基因簇均包含相同的15 个基因，且基因排布顺序相同。此外，与Penicillium digitatum（指状青霉）和Penicillium italicum（意大利青霉）关系较远，推测是因为两者均只含有棒曲霉素生物合成基因簇的部分基因且不产生PAT[33]。

表2 11 个MFS 转运蛋白登录号Table 2 11 MFS transporters accession numbers

图4 PatC 的系统发育树Fig.4 Phylogenetic tree of PatC

2.2 蛋白结构功能域预测

利用NCBI Domains & Structures 对PatC 进行功能结构域预测分析。如图5 结果表明，该蛋白氨基酸序列53aa～480aa 为MFS 结构域，验证PatC 属于MFS 转运蛋白超家族。Paulsen 等[34]研究MFS家族时发现，大多数MFS 蛋白的保守特征序列是基序A（GrLaDrfGrRRv）、基序B（LvaaRvlQGxGA）和基序C（gxxxGPxxxGGxl），其中位于TMS4 内的基序B 参与大多数MFS 成员的质子转移。Jiang 等[35]提出基序A 是MFS 成员中高度保守的基序，位于TMS 2 和TMS 3 之间的环中，在构象调节中起着至关重要的作用，特别是对外构象的稳定性。此外，Motif C 被称为反转运基序，在TMS5 内形成纽扣，可能与底物结合位点的取向有关[36]。

图5 PatC 的功能结构域预测Fig.5 Prediction of functional domains of PatC

2.3 理化性质分析

利用在线软件ProtParam 对PatC 的理化性质进行预测，结果显示（表3），PatC 的氨基酸数量是546 个，分子量为58733.15；理论等电点PI：5.57，为酸性蛋白；该蛋白氨基酸序列中亮氨酸（Leu）、甘氨酸（Gly）及丙氨酸（Ala）三种氨基酸含量较高，所占比值分别为12.1%、10.1%和9.7%；蛋白不稳定指数：40.3，为不稳定蛋白；从脂肪系数数据分析表明该蛋白脂肪系数：113.77，为脂溶性蛋白。亲水性平均系数（GRAVY）值为0.712，表明该蛋白为疏水性蛋白。

表3 PatC 的理化性质分析Table 3 Physicochemical property analysis of PatC

2.4 跨膜结构预测

采用在线软件TMHMM 对PatC 的跨膜结构进行预测，由图6 和上述理化性质结果表明该蛋白具有强疏水性氨基酸，具有能够行使跨膜转运的功能，符合转运蛋白的特点。

图6 PatC 跨膜区预测Fig.6 Transmembrane structure prediction of PatC

2.5 蛋白质二级结构预测

利用在线SOPMA 软件对PatC 二级结构进行预测和分析（图7），PatC 中二级结构主要存在4 种：α-螺旋、β-转角、β-折叠和无规则卷曲。二级结构元件中α-螺旋所占比例最高，此外无规则卷曲、β-折叠和β-转角所占比例依次减少分别为30.95%、24.73%和6.23%。

图7 PatC 的二级结构预测Fig.7 Secondary structure prediction of PatC

2.6 蛋白质三级结构预测

PatC 同源性蛋白较少，无法利用SWISS MODEL同源建模方法得到其三级结构。因此通过AlphaFold预测PatC 的三级结构。由图8 说明该蛋白是14 次α-螺旋结构（与跨膜结构预测相吻合），可能是进化过程中以12 次跨膜α-螺旋为基础产生的[37]，分为2 个结构域：N 端结构域和C 端结构域。每个结构域都由6 个α-螺旋组成，2 个结构域呈现二次赝对称（twofold psudosymmetry）[20]。与大肠杆菌中的二肽转运蛋白（Escherichia coli YbgH，EcYbgH）结构类似，为经典的12 次跨膜α-螺旋构象加两个TM 螺旋的插入（HA 和HB）[38]。

图8 PatC 的三级结构预测Fig.8 Tertiary structure prediction of PatC

2.7 分子对接

采用AutoDock Vina1.1.2 将PatC 与E-ascladiol进行半柔性对接，分子对接结果如图9，其之间的平均结合能为-5.0 kcal/mol（＜0），说明配体和受体蛋白可以自发结合。PatC 与E-ascladiol 之间主要以传统氢键（Conventional hydrogen bond）、堆积键（Pi-Pi Tshaped）及疏水作用力（Pi-Alkyl）相互作用。具体相互结合模式：E-ascladiol 与PatC 结合至蛋白空腔，且与该蛋白受体的氨基酸残基SER353、TYR336、PRO339、及PRO188 分别形成2.74Å氢键，2.96ÅPi-Pi 堆积作用力，5.21Å和529Å疏水作用力，这些相互作用使得蛋白质与小分子紧密结合。

图9 PatC 与E-ascladiol 的分子对接Fig.9 Molecular Docking of PatC and E-ascladiol

2.8 分子动力学模拟

2.8.1 RMSF 分析 RMSF（Root Mean Square Fluctuation，RMSF）是用来衡量一段时间内分子结构中原子位置的平均移动，是目的结构特定原子相对于参照结构在整个时间范围内的距离均方根，是判断蛋白质结构结构变化的重要参数之一[39]。

通过Protein（纯蛋白）与Wild（野生型复合体）的RMSF 数据分析，如图10 所示，PatC 在多个氨基酸位点发生柔性变化。与野生型复合物对比，纯PatC 分别在第127、181、188～197、231～241、273、296、333～341、457～463 位残基区域即（Ala127、Gly181、Pro188～Ser197、Gly231～Val241、Asp273、Phe296、Val333～Try341、Gln457～Ala463）发生了形变拉伸或缩短。其中位于Pro188～Ser197aa、Gly231～Val241aa 氨基酸区域内，结合底物分子后的PatC 发生柔性变化，在一定程度上增加其灵活度，改变其空间构象。因此可推测出上述氨基酸区域内可能存在PatC 与E-ascladiol 之间的结合位点。

图10 各模拟体系的RMSF 结果分析Fig.10 Analysis of RMSF results for each simulation system

2.8.2 点突变分析 通过上述分子对接和复合物分子动力学模拟两种方法得出不同的结合位点。由图11可知，MD 计算得到6 个结合位点，其中第188 位氨基酸与分子对接结果重叠，且与RMSF 数据相吻合，因此使用Pymol 软件对PatC 中第188 位氨基酸残基位点突变，将原有的脯氨酸（P）突变为丙氨酸（A），即产生新的P188A 的PatC 突变体。通过将PatC 与P188A 的突变位点进行构象对比，如图12，发现在第188 位氨基酸位置上的点突变并未使该位点的二级结构出现明显差别，且不影响整体蛋白质构型，因此继续使用此结果进行后续实验研究。

图11 各个结合位点的总结合能Fig.11 Total binding energy of each binding site

图12 PatC 及其突变体结构Fig.12 Structure of PatC and its mutant

2.8.3 RMSD 分析 RMSD（Root Mean Square Distance/Deviation，RMSD）是用来衡量同一分子结构不同构象之间的差异，常用来分析体系分子整体结构的稳定性，RMSD 变小表明结构趋于稳定，反之趋于不稳定[39]。

将纯Protein（纯蛋白）、Wild（野生型复合体）、P188A（突变体复合体）三个复合体系体系分别进行了20 ns 的分子动力学模拟实验。在对其数据进行采集和分析中，三个模拟体系结构中的原子骨架的根均方偏差RMSD 值如图13 所示。当反应进行至15 ns时，三个反应体系曲线趋于平稳，虽有一定的波动，但范围在0.05 nm 之内。说明三个体系均达到了平衡状态，且三个反应体系的RMSD 终值都停留在0.3 nm以下，其中Protein、Wild 复合体系和Mutant 复合体系的RMSD 平均值分别为0.17、0.23 和0.18 nm 左右，可看出三者在模拟过程中结构均保持稳定。且相较于Wild，Mutant 稳定性更佳。

图13 各模拟体系的RMSD 结果分析Fig.13 Analysis of RMSD results for each simulation system

2.8.4 Rg 分析 Rg（Radius of Gyratio，Rg）回旋半径用来衡量分子结构的紧密度或折叠程度，Rg 考察主要选择平衡时状态进行[39]，因此主要选择后10 ns进行分析，由图14 可得，后10 ns 两个反应体系中Rg 变化均很小，表明这三个反应体系稳定且在反应过程中均为平衡状态且复合物结构致密，且三个反应体系的平均值均在2.5 nm 以下，并且空白蛋白的Rg值最高，两种结合底物分子的复合物体系中Rg 值均有明显减少，由此可见在PatC 结合了底物分子后，整体构象进行了收缩。

图14 各模拟体系的Rg 结果分析Fig.14 Analysis of Rg results for each simulation system

3 结论

本研究基于PatC 与E-ascladiol 结合的前提下，运用生物信息学、分子对接和分子动力学模拟方法预测蛋白与分子之间的重要作用位点，以研究PatC 转运机制。生物信息学结果表明，该蛋白的二级结构以为α-螺旋和无规则卷曲为主，三级空间构象含有14 个跨膜螺旋，含有2 个结构域：N 端结构域和C 端结构域，形成独特的MFS 折叠。根据以上结果说明PatC 作为毒素输出泵转运PAT，但其理化性质也需进一步的验证。分子对接结果表明PatC与E-ascladiol 之间存在多个结合位点与多种相互作用力，与RMSF 数据结合分析得出P188 作为重要结合位点并对其进行模拟突变。此外，本研究通过分子动力学方法对两种复合体系（野生型Wild 和突变体P188A）进行各参数分析，结果说明P188 氨基酸位点可能是PatC 的关键作用位点，可为后续的微生物分子实验提供基础论证。综上所述，本文为棒曲霉素合成途径中的MFS 蛋白的转运机制（毒素外排作用）提供一定的基础，并进一步为毒素输出泵的生物学功能与转运机制研究提供参考。