黑曲霉α-L-鼠李糖苷酶與柚皮苷分子動力模擬研究

鞏建業(yè) ， 吳喆瑜 ， 劉嘉男 ， 廖 輝 ， 李利君 *，2， 倪 輝 ，2

（1.集美大學 食品與生物工程學院，福建 廈門361021；2.福建省食品微生物與酶工程重點實驗室，福建 廈門361021）

柑橘是全球第一大種植和消費水果，是僅次于小麥和玉米的第三大國際貿易農產品[1]。目前對于柑橘類水果的主要加工方式為外果皮生產蜜柚香精、皮囊渣等提取果膠和鮮切柚果等，但對比橙子、檸檬等柑橘類水果，產量最大、附加值最高的蜜柚果汁類產品始終無法深入開發(fā)，主要原因是蜜柚壓榨果汁具有明顯的苦味，這大大制約了蜜柚的高附加值利用[2]。柚皮苷是柑橘類果汁加工中后呈苦味的重要物質之一，而α-L-鼠李糖苷酶能高效降解柚皮苷，生成有醫(yī)用價值的普魯寧和L-鼠李糖[3-4]。但工業(yè)生產中往往需要性質更優(yōu)良的α-L-鼠李糖苷酶，這就需要進行蛋白質工程改造。研究酶與底物相互作用方式，可為工程改造奠定基礎[5]。

隨著核磁共振、X射線衍射、冷凍電鏡及同源建模、從頭計算等解析蛋白質結構技術的發(fā)展，越來越多的蛋白質三維結構被測出來。這為研究蛋白質與底物小分子奠定了基礎。分子對接是一種靜態(tài)的研究受體大分子與配體小分子結合作用的分子模擬方法，并通過對接相互作用能、關鍵殘基等方面評價受體與配體的相互作用[6]。分子動力模擬按照給定溫度下的麥克斯韋速率分布給原子賦予一定的初始速度，在這一定的初始條件下求解牛頓第二定律的微分方程，得到蛋白質中所有原子的運動軌跡[7]，MM-PBSA計算結合自由能可以在一定程度上反映受體與配體之間結合的主要推動力[8]。因此，在基于黑曲霉α-L-鼠李糖苷酶和柚皮苷結構的基礎上進行分子對接和分子動力模擬，進而利用MMPBSA法計算結合自由能，分析黑曲霉α-L-鼠李糖苷酶和柚皮苷結合主要推動力，并揭示結合過程中起關鍵作用的氨基酸，探討α-L-鼠李糖苷酶和柚皮苷的結合機理，可為開發(fā)更高效降解柚皮苷的α-L-鼠李糖苷酶提供理論依據(jù)。

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]已經(jīng)證實（α/α）6桶裝結構域是α-L-鼠李糖苷酶的催化結構域。在100次對接結果中，柚皮苷與α-L-鼠李糖苷酶的對接位置均在（α/α）6的桶裝結構域底部的loop區(qū)，這與Grandits等人[19]的分子對接結果完全吻合。圖2（a）顯示了柚皮苷與α-L-鼠李糖苷酶最佳對接結果得相對位置，圖2（b）顯示了柚皮苷范圍內可能有結合作用的的氨基酸。Autodoc4.2依據(jù)半經(jīng)驗的自由能計算方法來評價受體和配體之間的能量匹配，其值越小表明配體與受體之間的親和力越好[20]。在100次對接實驗結果中結合自由能均為負值，最佳對接結果的結合作用能為-6.81 kcal/mol，這表明柚皮苷與黑曲霉α-L-鼠李糖苷酶有良好的結合作用。

2.2 柚皮苷與α-L-鼠李糖苷酶復合物分子動力模擬

我們對柚皮苷與α-L-鼠李糖苷酶的模型復合物進行了20 ns的分子動力模擬，均方根偏差（RMSD）反應出2個分子結構相同原子間的距離，可以用來評價體系穩(wěn)定性，若均方根偏差隨時間的變化較小，則表明體系達到平衡狀態(tài)[21]。圖3顯示了在20 ns模擬過程中RMSD隨時間變化的結果，α-L-鼠李糖苷酶在8.0 ns左右達到了平衡狀態(tài)，柚皮苷在2.5 ns達到了平衡狀態(tài)。分子動力模擬平衡后的軌跡用于進一步分析。

圖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部分）三部分組成，計算所得的結合自由能可用來反映配體與受體結合的穩(wěn)定性[22]。表1列出了柚皮苷與α-L-鼠李糖苷酶對結合自由能貢獻的各項能量值。柚皮苷與α-L-鼠李糖苷酶之間的范德華力、庫倫靜電相互作用，極性溶劑化能和非極性溶劑化能分別為-285.931，-24.353，119.325 kJ/mol和-27.068 kJ/mol，這表明柚皮苷與α-L-鼠李糖苷酶相互作用以范德華作用力為主，靜電相互作用及非極性溶劑化能力對結合自由能貢獻較小，極性溶劑化能阻礙柚皮苷與α-L-鼠李糖苷酶的結合。此外，在開發(fā)酶抑制劑時，往往也需要考慮抑制劑與酶分子的相互作用力，計算柚皮苷與α-L-鼠李糖苷酶相互作用力為開發(fā)酶抑制劑提供了思路。

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

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

為了揭示黑曲霉α-L-鼠李糖苷酶中與柚皮苷結合的關鍵氨基酸，將MM-PBSA總結合自由能分解到各個氨基酸殘基上（見圖4），發(fā)現(xiàn)Trp236、Asp239、Pro276、Ser339、Ala340、Asp341、Tyr359、Phe461、Phe465、Phe501、Glu503、Tyr516、Pro520、Val522 這 14 個 殘基對總結合自由能貢獻值較大。表2列出了這些關鍵氨基酸的MM-PBSA結合能分解，除Glu503與柚皮苷的主要作用方式為極性溶劑化能外，其它殘基與柚皮苷的主要作用方式均為范德華力和庫倫電荷作用，另外Asp239也發(fā)揮了有利于柚皮苷結合的極性溶劑化能。這為蛋白質工程設計更高效α-L-鼠李糖苷酶分子提供了理論數(shù)據(jù)。

圖4 柚皮苷與α-L-鼠李糖苷酶復合體系結合自由能分解圖Fig.4 Binding free energy of complex of naringin and Aspergillus niger α-L-rhamnosidase

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

非極性相互作用能對結合自由能有促進作用，柚皮苷可與黑曲霉α-L-鼠李糖苷酶形成較多的疏水相互作用。圖5為通過Ligplot軟件統(tǒng)計顯示的柚皮苷與α-L-鼠李糖苷酶結合狀態(tài)下的疏水氨基酸殘基。

由圖5（a）可知在柚皮苷與α-L-鼠李糖苷酶結合 過 程 中 ，Trp236、Ala340、Ile462、Phe461、Tyr516、Val522、Trp528等殘基均與柚皮苷形成了疏水相互作用。將這些殘基對應于三維空間模型，可以發(fā)現(xiàn)這些氨基酸殘基在α-L-鼠李糖苷酶（α/α）6桶裝催化結構域底部形成了一個疏水性口袋，而底物柚皮苷被包裹在這個疏水口袋里（圖5（b））。由此可見疏水作用可能是α-L-鼠李糖苷酶（α/α）6桶能夠結合柚皮苷的重要原因。

2.6 柚皮苷與黑曲霉α-L-鼠李糖苷酶氫鍵作用分析

通過對柚皮苷與α-L-鼠李糖苷酶結合自由能和疏水作用的分析，發(fā)現(xiàn)極性相互作用雖然很大程度上表現(xiàn)為阻礙柚皮苷與α-L-鼠李糖苷酶，但并不是所有形成極性相互作用的氨基酸都不利于結合的形成，例如Ser286在結合過程中也表現(xiàn)出了有利于結合自由能的極性作用，因為在結合過程中它與柚皮苷形成了氫鍵。利用VMD軟件我們分析了分子動力模擬中平衡后的氫鍵，將體系中穩(wěn)定的氫鍵統(tǒng)計在表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 結 語

越來越多的人采用分子動力學研究生物分子相互作用時的動態(tài)變化過程。作為一個被廣泛應用的重要手段，分子動力學模擬在揭示配體-受體相互作用關系方面變得日益重要[23-24]。本研究表明范德華力是柚皮苷與黑曲霉α-L-鼠李糖苷酶結合的主要驅動力，庫倫靜電力和非極性溶劑化能對于結合作用較小。疏水作用分析發(fā)現(xiàn)柚皮苷結合在Trp236、Ala340、Ile462、Phe461、Tyr516、Val522、Trp528等殘基形成在一個疏水性口袋里面，而Ser286、Phe465、Pro520與柚皮苷形成了穩(wěn)定的3條氫鍵。這些位點都是黑曲霉α-L-鼠李糖苷酶與柚皮苷相互作用的關鍵位點，為深入研究α-L-鼠李糖苷酶與柚皮苷的構效關系確定了實驗靶點。因此本研究為了解α-L-鼠李糖苷酶與柚皮苷之間的相互作用提供了重要信息，同時也為α-L-鼠李糖苷酶的抑制劑的開發(fā)和α-L-鼠李糖苷酶的定向改造提供了理論支持。