利用同源建模、分子模拟和分子对接技术分析家蚕中碳酸酐酶的结构和催化机制

陈胤熹，黎菁菁，罗佳伟，郑少鹏，余洁婷，黄嘉惠，李鑫尧，余铭怡，郝锦亨，黎佩瑜，古伟明，吴易达，曹诗林,，赖林浩

（1.佛山科学技术学院食品科学与工程学院，广东佛山 528000）（2.华南农业大学食品学院，广东广州 510641）（3.广东工业大学轻工化工学院，广东广州 510006）（4.暨南大学理工学院，广东广州 510640）（5.华南理工大学食品科学与工程学院，广东广州 510640）

碳酸酐酶（Carbonic Anhydrase，碳酸酐酶）是一种重要的含锌金属酶，广泛存在于动、植物及微生物体中。碳酸酐酶不仅能够快速、高效地催化二氧化碳和碳酸氢盐的可逆转化而且还具有维持生物体内环境稳态、参与光合作用、钙化作用、CO2及碳酸盐的转运等生命功能[1]。由于碳酸酐酶的这些特性，它既可被应用于生物检测，CO2捕集和生理诊断等领域，也被广泛应用于食品和药品工业中，因此具有良好的经济价值和应用前景[2-3]。

根据碳酸酐酶的氨基酸序列不同，一般可分为α、β、γ、δ、ε、η、θ和ι等八种类型[4]，游离的碳酸酐酶普遍存在着容易失活、热稳定性低、可重复利用率低等问题，难以得到大规模应用[5]，因此，开拓更多来源的碳酸酐酶酶，并对其进行分子改造和设计，是解决该问题的潜在途径。获得准确的碳酸酐酶的三维结构是理性改造的基础。

目前已知的获得蛋白质结构的方法有很多，例如：核磁共振，X射线单晶衍射、冷冻电镜等实验方法，但是这些方法成本高，耗时长，难以得到更广泛的应用。另一方面，可以使用蛋白质建模技术预测蛋白质的结构信息。蛋白质建模的方法包括同源建模和非同源建模。若在现有蛋白质结构数据库可检索出与目的蛋白具有序列同源性的模板，则可进行同源建模。若无法得到同源模板，或同源模板无法覆盖蛋白全序列，则需要通过穿针引线、从头算或分段建模和结构域组装等方法进行建模。通过从NCBI与SilkDB3.0中获得大量的数据[6,7]，使用同源建模将已知结构的蛋白质的序列与同源蛋白质序列进行比对，通过分析蛋白质序列中的保守区域来识别蛋白质活性区域，可望快速预测出蛋白质的结构信息[8]。

家蚕（Bombyxmori，Bm）的第五龄幼虫的丝腺中存在着四种不同的上皮细胞，其中一种产生活性碳酸酐酶，用于调节家蚕内环境的 pH，维持管腔内的pH梯度[9]。本文以家蚕碳酸酐酶为研究的对象，通过蛋白建模、结构分析、分子动力学模拟、分子对接等手段，建立和分析了家蚕碳酸酐酶的蛋白质结构模型、活性区域和催化作用机制，并对所建立的模型进行了一系列的计算评分，最后采用分子动力学模拟对结果进行了验证。本文研究为家蚕碳酸酐酶后续的分子改造提供了理论基础。

1 材料与方法

1.1 数据库与软件

家蚕碳酸酐酶的基因和蛋白质序列信息来源于NCBI与SilkDB3.0数据库[10]。NCBI是美国国家分子生物学信息资源中心，具有广大的公共数据库，可被用于检索和分析基因组数据等。而SilkDB3.0是由西南大学蚕丝科学与技术研究团队维护的家蚕的多组学数据库，可进行可视化的多组学分析。蛋白质立体结构信息从蛋白质结构数据库RCSB PDB获得，该数据库保存了通过X射线衍射、核磁共振、冷冻电镜等实验手段获得蛋白质等生物大分子的三维结构数据。利用SWISS-MODEL计算平台和Modeller软件进行蛋白质建模。其中SWISS-MODEL是一个自动化的蛋白质比较建模服务器。Modeller是一款同源建模软件，主要用于家蚕碳酸酐酶的优化与评价[11]。利用PROCHECK评价蛋白质结构合理性[12]。利用分子三维结构显示软件PyMol对蛋白质三维结构进行可视化分析[13]。通过Gromacs软件[14]考察酶在催化过程中结构变化及反应体系的状态变化。AutoDock-Vina是由Scripps研究所的Olson实验室开发与维护的分子对接的开源软件[15]，用于对乙酸对硝基苯酯底物与家蚕碳酸酐酶分子对接计算，获得最优的对接结果，获得对接模型后，结合gromacs的MM/PBSA分析插件，共同分析家蚕碳酸酐酶与底物间相互作用的具体数值和具体部位。MM/PBSA[16]是利用分子力学和连续介质模型来计算蛋白和配体之间存在的结合自由能，其中配体可以是蛋白、小分子或者多肽，在该方法中，计算分子力学（MM）、极性溶剂化能（PB）、非极性溶剂化能（SA）的相关能量得以被分解计算及分析。其中PROCHECK、AutoDock-Vina以及Modeller已经被本团队进行二次开发并部署于互联网在线平台中（https://atomevo.com/），开放给研究人员免费使用。

1.2 方法

1.2.1 碳酸酐酶同源建模

在NCBI数据库和SilkDB3.0中获得家蚕碳酸酐酶的信息（XP_004922882.1）序列为：MDNRTVKID PKFLSAQPKKTSSDAEHISRLRPSQSPIAISLSRCPT WSSLDPLKFKGYWDSNANAILLNNGSTAYFTFND ASVRPTLSGGPLIGEYIFEQMHFHWSVDDFTGCEH VLDGHGYAAECHFVHYNSKYESLETAVGHPDGLA VVGFLLETVDAPNPRFDRLVQGLEGIQKRESVMN VTSESLLWMDREDLQIGNYVTYKGSLTTPPYTEC VTWIIYEKPVQIGSEQLGLLRQLEGPDSQPIERNVR PTQRHPPGHSVIYVKQVRSKL。

首先使用 SWISS-MODEL[17]对所得到的家蚕碳酸酐酶蛋白质序列进行模板搜索，筛选出酶种类接近、同源性高、序列覆盖度高的模板，并进行同源建模，其中SWISS-MODEL服务器模板数据库ExPDB是从PDB中提取的[18]。

如出现单个模板无法覆盖整个目标序列的情况，则先进行分段建模，再进行穿针引线拼接。具体方法是：（1）若结构域存在同源性高的模板，则利用SWISS-MODEL进行同源建模，若部分结构域无同源模板，则使用trrosetta进行建模。（2）将得到的各结构域模型作为模板，使用本团队自行开发的Modeller软件的穿针引线拼接程序进行拼接，从而得到目标蛋白的结构。

1.2.2 碳酸酐酶潜在活性区域

综合分析RCSB PDB数据库中的模板，以及对应结构的参考文献中的活性区域信息进行分析。利用活性区域的信息包括：易形成疏水口袋、组氨酸聚集、活性中心金属离子临近区域。再通过比对 SWISSMODEL中最优结果酶蛋白氨基酸序列与RCSB PDB得到同样种属的酶蛋白的氨基酸序列，得出酶潜在的活性区域。

1.2.3 碳酸酐酶活性区域与底物的分子对接

以乙酸对硝基苯酯作为模型底物，以碳酸酐酶模型作为受体结构，利用Auto Dock-Vina模块中将受体与配体进行对接。根据1.2.2得到碳酸酐酶活性区域的空间位置确立一个长方形盒子限制计算范围，计算得到盒子参数后与受体的PDB文件、配体的PDB文件共同上传至 Auto Dock-Vina模块进行计算，并利用Vina-XScore连用模块，对分子对接结果进行评分，找出受体与配体结合模式[15]。

1.2.4 家蚕碳酸酐酶的分子动力学模拟

本课题利用 Gromacs软件进行分子动力学模拟[19]，以Amber03力场。模拟步骤如下：首先，输入在1.2.3中对接后所输出的蛋白结构和力场参数的相关文件（电荷，键合参数和非键参数等势能函数）。随后，确定其体系的大小并且对体系进行能量最小化，再进行NVT与NPT平衡。随后开始100 ns分子模拟，并且分析模拟输出关于溶剂可及面积、氢键数目和径向分布函数的数据，并进行MM/PBSA分析。

2 结果与分析

2.1 家蚕碳酸酐酶同源建模及结构评价

首先通过SWISS-MODEL进行同源建模，结果显示家蚕碳酸酐酶的主要区域（25E-268K）可以利用Protein Data Bank中4ZX1模型为模板进行建模，但是1M-24A以及269Q-274L这两个区域缺少模板。因此先通过SWISS-MODEL同源建模得到主体结构区域。1M-24A以及269Q-274L区域通过trRosetta进行非同源建模。随后利用本课题组自行开发的基于Modeller的穿针引线建模方法，将三个模板进行穿针引线拼接，得到家蚕碳酸酐酶模型。

Ramachandran plot用于描述蛋白质结构中氨基酸残基二面角φ和Ψ是否处于合理区域，可用于表征同源建模结果的合理性[20]。结果表明（图 1）其完全允许区域（89.3%）与允许区域（10.3%），总和超过了99%[21]，表明本研究所得的家蚕碳酸酐酶结构合理。

2.2 家蚕碳酸酐酶潜在活性部位分析

通过前人对与家蚕碳酸酐酶同种属的碳酸酐酶 IX（PDB代号：4XZ1）的研究表明[22]，4ZX1的活性区域为 62N～65S，89L～91Q，120L～122H，131V～143V，198L～206C。利用分子叠合技术对家蚕碳酸酐酶与4ZX1的蛋白结构进行叠合对比（图 2a），确定家蚕碳酸酐酶活性口袋的空间位置。通过蛋白序列进行比对分析（图2b），得到家蚕碳酸酐酶的潜在活性口袋为：70N～72S（区域A），98F～100Q（区域B），128F～130H（区域C），138L～150V（区域D），209L～217C（区域E）。

以家蚕碳酸酐酶蛋白的活性口袋周边（空间中心坐标为 X=0.1、Y=1.4、Z=88）为对接位点，通过Autodock-Vina进行蛋白与对乙酸对硝基苯酯之间的分子对接。对接后得到最优对接结果表明，家蚕碳酸酐酶与底物结合能为-6.1 Kcal/mol。

2.3 家蚕碳酸酐酶的分子动力学模拟数据分析

均方根偏差（RMSD）分析可以揭示家蚕碳酸酐酶整体构象的稳定性。该结果显示（图 3a），在模拟的前10 ns，家蚕碳酸酐酶的RMSD迅速从0 nm上升至0.3 nm，随后在0.3 nm～0.35 nm左右波动。家蚕碳酸酐酶最后50 ns的RMSD值约为0.35 nm，表明该酶的结构稳定。家蚕碳酸酐酶最后50 ns的平均可及面积（图3b）为142 nm2。家蚕碳酸酐酶蛋白分子内平均氢键数（HBN）为181，家蚕碳酸酐酶与水溶液间的平均氢键数为587。

径向分布函数可以揭示酶蛋白与水的分布距离以及对应的概率密度，对模拟输出的数据作图分析，研究表明水与家蚕碳酸酐酶在距离＜0.35 nm区间形成一个聚集峰，该峰主要归因于水和蛋白之间的氢键作用（图4）。

2.4 家蚕碳酸酐酶MM/PBSA分析

在动力学分子模拟后，通过MM/PBSA法进一步分析对接模型的结合自由能，得到总结合能为-40.72 kJ/mol。其中，对自由能的贡献主要是范德华力的相互作用-57.86 kJ/mol。而极性溶剂化21.39 kJ/mol，对结合有显著的反作用。通过分析蛋白质氨基酸残基结合能（图4），得知家蚕碳酸酐酶对接后的模型中，分子间相互作用力主要存在于家蚕碳酸酐酶活性口袋中的D、E区域。

表1 家蚕碳酸酐酶的主要氨基酸残基与底物乙酸对硝基苯酯的g_mmpbsa分析结果Table 1 The g_mmpbsa analysis of 4-nitrophenyl acetate and the main residues of BmCA

3 结论

3.1 本文通过同源建模建立结构合理的家蚕碳酸酐酶模型，结果显示其氨基酸残基处于完全允许区域（89.3%）与允许区域（10.3%）的总和超过了99%。2010年邓秋红等预测甘蓝型油菜碳酸酐酶结构，其模型的完全允许区域（87.0%）与允许区域（12.4%）总和超过99%[23]，通过Ramachandran plot分析的结果与前人的研究数据对比，证明该酶结构模型可靠性良好。2015年Mahon Brian P等研究碳酸酐酶IX（PDB代号：4ZX1），获得了潜在活性区域为：62N～65S、89L～92Q、120L～122H、131V～143V、198L～206C。利用分子叠合技术对家蚕碳酸酐酶与 4ZX1的蛋白结构叠合对比，获得了家蚕碳酸酐酶潜在的活性口袋：70N～72S、98F～100Q、128F～130H、138L～150V、209L～217C。结果表明同源建模所得的家蚕碳酸酐酶潜在活性区域与4ZX1号蛋白活性区域是相吻合的。

3.2 随后，利用Autodock-Vina将家蚕碳酸酐酶与乙酸对硝基苯酯进行对接，发现最优对接结果的家蚕碳酸酐酶与底物的结合能为-6.1 kcal/mol，2007年向福利用分子对接技术研究与本文同种属的碳酸酐酶（PDB：2CAB）与乙酰唑胺进行分子对接，得到对接结合能为-2.1 Kcal/mol[24]，经对比表明家蚕碳酸酐酶与底物乙酸对硝基苯酯经对接得到的结果较为稳定。通过MM/PBSA探究家蚕碳酸酐酶与乙酸对硝基苯酯的相互作用，其结果表明：（1）范德华力在结合中占主导地位，而极性溶剂化对结合有显著的反作用；（2）家蚕碳酸酐酶和底物相互作用的区域为138L～150V（区域D）和209L～217C（区域E）。2013年孙维琦等运用分子动力学模拟和自由能计算方法，研究了苯磺酰胺分子从碳酸酐酶II的活性位点的相互作用，显示关键残基198L、199T和200T（该区域相当于本研究中的区域E）与苯磺酰胺的氢键作用阻碍了底物从酶中的脱离[25]。

3.3 综上所述，本文通过同源建模获得了一个良好的蛋白质结构模型，并通过分子模拟与分子对接技术探究其潜在活性位点与催化机制，为家蚕碳酸酐酶进一步改造探究以及实际生产应用提供理论依据。