基于分子動力學模擬的擴展青霉棒曲霉素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 蛋白的轉運機制（毒素外排作用）提供一定的基礎，并進一步為毒素輸出泵的生物學功能與轉運機制研究提供參考。