基于分子動力學模擬的擴展青霉棒曲霉素MFS 蛋白轉(zhuǎn)運機制研究

楊 琪，王艷玲

（蘭州理工大學生命科學與工程學院，甘肅蘭州 730050）

青霉病是全球最常見且經(jīng)濟上最重要的采后水果腐爛病之一[1]，其中擴展青霉（Penicillium expansum，P.expansum） 是引起果實青霉病或軟腐病的主要病原菌[2]。P. expansum可在水果作物（蘋果、梨、柑橘和櫻桃等）及其衍生品、加工產(chǎn)品的采后、貯藏、運輸和加工期間侵襲[3-7]。同時，P. expansum還可產(chǎn)生具有危害的次級代謝產(chǎn)物——棒曲霉素（Patulin，PAT）[8-10]，通過誘導DNA 鏈斷裂、氧化損傷以及DNA 鏈之間相互交聯(lián)、產(chǎn)生核質(zhì)橋等機制影響微核形成，從而引起DNA 結(jié)構(gòu)性損傷[11]，產(chǎn)生遺傳毒性、細胞毒性和免疫毒性作用等[12]。

20 世紀80 年代，Puel 等提出PAT 的生物合成途徑是一個酶級聯(lián)反應(yīng)[13-14]。近幾年，Li 等[15]鑒定得到參與PAT 生物合成的基因簇，包括15 個基因（PatA～PatO），具體合成途徑見圖1。其中PatC基因編碼MFS 超家族轉(zhuǎn)運蛋白（Major Facilitator Superfamily，MFS）。MFS 的主要功能是將PAT 的前體E-ascladiol 轉(zhuǎn)運至胞外，最終合成PAT[16]，但具體的轉(zhuǎn)運機制尚不清楚。研究表明，E-ascladiol 對來自人類肝、腎、腸和免疫系統(tǒng)的細胞無毒性，因此，增加E-ascladiol 累積以減少棒曲霉素含量是控制棒曲霉素風險的有效手段[17]。目前，MFS 轉(zhuǎn)運蛋白包括105 個家族，是最大的次級跨膜轉(zhuǎn)運蛋白超家族[18]。MFS 轉(zhuǎn)運蛋白大多數(shù)都含有400～600 個氨基酸殘基[19]，具有12、14 或24 個不規(guī)則跨膜α-螺旋（TMSs）包圍形成深中央親水腔[20-21]，催化底物在兩個方向上進行跨膜運輸，可以是運輸單個底物的單轉(zhuǎn)運子，可以是將底物和偶聯(lián)離子運輸在一起的同轉(zhuǎn)運子，也可以是按相反方向依次運輸兩個不同分子的反轉(zhuǎn)運子（圖2a）。根據(jù)交替開放模型，MFS 能夠在向胞內(nèi)開放，閉塞和向外開放構(gòu)象之間變化（圖2b），轉(zhuǎn)運蛋白上的底物結(jié)合位點暴露在細胞膜的不同側(cè)，底物從細胞膜的一側(cè)結(jié)合到蛋白上，而從另一側(cè)被釋放（圖2c）。此外，MFS 轉(zhuǎn)運蛋白還可以間接調(diào)節(jié)真菌的內(nèi)部pH和應(yīng)激反應(yīng)機制[22]，參與多種真菌毒素的輸出。如擬枝孢鐮刀菌（Fusarium sporotrichioides）的TRI12和大豆紫斑菌（Cercospora kikuchii）的CFP[23-24]，兩者均編碼MFS 轉(zhuǎn)運蛋白，分別參與真菌毒素-單端孢霉烯（Trichothecenes）和尾孢菌素（Cercosporin）的分泌。植物病原體，玉米圓斑真菌（Cochliobolus carbonum）的TOXA也與HC 毒素的宿主特異性毒素的分泌有關(guān)[25-26]。TOXA、Tri12及CFP缺失突變體都顯示出毒素產(chǎn)生的急劇減少，Tri12缺失突變體也觀察到生長遲緩，但CFP缺失突變體與野生型菌株一樣正常生長，只是對大豆的毒力降低[26-28]。此外，在攜帶毒素合成基因的基因簇中也發(fā)現(xiàn)了一些編碼MFS 的基因[29]。然而，它們在運輸特定代謝產(chǎn)物中的作用仍待闡明。可惜的是，對MFS 轉(zhuǎn)運機制的研究甚少，目前只解析得到幾個細菌蛋白的完整拓撲結(jié)構(gòu)。相關(guān)PAT 毒素輸出的研究更是沒有。因此致力于探索PatC 的轉(zhuǎn)運機制，有利于蘋果青霉病病原菌的防治的研究。

圖1 PAT 的生物合成途徑[30]Fig.1 PAT biosynthetic pathway[30]

圖2 MFS 蛋白的運輸機制[31]Fig.2 Mechanism of transport of MFS[31]

隨著計算生物學的發(fā)展，越來越多的研究者采用分子動力學模擬生物大分子在三維空間的運動狀態(tài)，從而分析其發(fā)揮生物功能的作用機制。本研究采用生物信息學及分子動力學模擬方法探索擴展青霉棒曲霉素PatC 的轉(zhuǎn)運機制，研究其與底物分子Eascladiol 的結(jié)合位點，并針對PRO188 定點突變，分析復(fù)合體及突變體的RMSD、Rg、RMSF 等變化，以期為研究該蛋白的結(jié)構(gòu)和轉(zhuǎn)運機制奠定基礎(chǔ)。

1 材料與方法

1.1 生物信息學分析

通過美國國家生物技術(shù)信息中心（National Center for Biotechnology Information，NCBI）搜索Penicillium expansumT01（登錄號：AYHP00000000.1），下載Penicillium expansumT01 中全部527 段鳥槍法測序結(jié)果。利用DNAMAN 生物信息軟件BLAST其他亞種Penicillium expansum strainNRRL 35695（登錄號：KF899892.1）中PatC基因及其編碼的氨基酸信息，結(jié)果與Penicillium expansumT01 中PatC基因序列完全一致。采用Mega 5.0 鄰近法（Neighbor Joining，NJ）進行氨基酸序列同源性分析及構(gòu)建系統(tǒng)進化樹；通過在線網(wǎng)站分析PatC 的理化性質(zhì)、跨膜結(jié)構(gòu)域、蛋白質(zhì)功能結(jié)構(gòu)域、二級結(jié)構(gòu)及三級結(jié)構(gòu)等，具體所用網(wǎng)站見表1。

表1 生物信息學相關(guān)分析工具Table 1 Bioinformatics correlation analysis tools

1.2 分子對接

配體分子和蛋白三維結(jié)構(gòu)的優(yōu)化：分子對接用于驗證活性成分與關(guān)鍵靶點之間的結(jié)合活性。通過AlphaFold 蛋白質(zhì)三維結(jié)構(gòu)數(shù)據(jù)庫在線搜索PatC 的三維結(jié)構(gòu)，再對其三維模型導入AutoDocktools（v1.5.6）進行加氫、計算電荷、分配電荷、指定原子類型等優(yōu)化處理；從Pubchem 數(shù)據(jù)庫下載小分子或蛋白的3D 結(jié)構(gòu)，導入ChemBio3D Ultra 14.0 進行能量最小化。將優(yōu)化好的小分子導入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 對接的最佳結(jié)合構(gòu)象，作為初始復(fù)合物結(jié)構(gòu)進行MD 模擬。隨后使用蛋白定點突變軟件，構(gòu)建新型PatC 突變體結(jié)構(gòu)（P188A）觀察突變后的蛋白與底物對接后復(fù)合物的活性。

具體操作如下：使用分子動力學軟件Gromacs2020.6 進行模擬。復(fù)合物小分子在GAFF 力場下生成拓撲文件，蛋白在AMBER99SB-ILDN 力場生成相應(yīng)參數(shù)化文件。以復(fù)合物為中心，添加SPC/E 水模型分子，生成立方體水盒子，向復(fù)合物水盒子里添加4 個CL-離子保持模擬體系的電中性。在體系完成后，使用最陡能量下降法進行50000 步能量最小化，使體系充分穩(wěn)定。采用正則系綜（NVT），Berendsen 熱浴法控制溫度為298.15 K，模擬50000 步，步長2 fs。采用等溫等壓系綜（NPT），Parrinello-Rahman 控壓，1 個大氣壓，模擬50000 步，步長2 fs。進行mdrun，步長2 fs，模擬10000000，總計20 ns 的模擬，對模擬后的數(shù)據(jù)計算MMPBSA。

1.4 數(shù)據(jù)處理

利用Discovery studio 2016 對分子對接結(jié)果進行可視化，且利用PyMOL（4.3.0 版）軟件繪制蛋白與分子之間相互作用的氨基酸殘基。分子動力學模擬使用軟件gromacs 分析各組參數(shù)變化，最終通過程序統(tǒng)計結(jié)果文件計算出體系各種宏觀性質(zhì)和微觀性質(zhì)，并使用xmgrace 繪制相應(yīng)的數(shù)據(jù)圖，如RMSD、RMSF、SASA 及氫鍵變化等。PAT 生物合成圖（圖1）使用在線網(wǎng)站FigDraw（https://www.figdraw.com）繪制，圖2 使用Adobe Illustrator 軟件繪制。

2 結(jié)果與分析

2.1 構(gòu)建系統(tǒng)進化樹

選取與PatC 同源性較高的多條已知菌株的氨基酸序列（表2），利用MEGA5.0 軟件構(gòu)建NJ 系統(tǒng)進化樹（圖4）。所選的11 個物種序列與PatC 的氨基酸序列一致性較高。系統(tǒng)進化樹結(jié)果表明，Penicillium expansum（擴展青霉）與Penicillium griseofulvum（灰黃青霉）處于同一分支，親緣關(guān)系較近。此現(xiàn)象可能是因為兩者的PAT 合成基因簇均包含相同的15 個基因，且基因排布順序相同。此外，與Penicillium digitatum（指狀青霉）和Penicillium italicum（意大利青霉）關(guān)系較遠，推測是因為兩者均只含有棒曲霉素生物合成基因簇的部分基因且不產(chǎn)生PAT[33]。

表2 11 個MFS 轉(zhuǎn)運蛋白登錄號Table 2 11 MFS transporters accession numbers

圖4 PatC 的系統(tǒng)發(fā)育樹Fig.4 Phylogenetic tree of PatC

2.2 蛋白結(jié)構(gòu)功能域預(yù)測

利用NCBI Domains & Structures 對PatC 進行功能結(jié)構(gòu)域預(yù)測分析。如圖5 結(jié)果表明，該蛋白氨基酸序列53aa～480aa 為MFS 結(jié)構(gòu)域，驗證PatC 屬于MFS 轉(zhuǎn)運蛋白超家族。Paulsen 等[34]研究MFS家族時發(fā)現(xiàn)，大多數(shù)MFS 蛋白的保守特征序列是基序A（GrLaDrfGrRRv）、基序B（LvaaRvlQGxGA）和基序C（gxxxGPxxxGGxl），其中位于TMS4 內(nèi)的基序B 參與大多數(shù)MFS 成員的質(zhì)子轉(zhuǎn)移。Jiang 等[35]提出基序A 是MFS 成員中高度保守的基序，位于TMS 2 和TMS 3 之間的環(huán)中，在構(gòu)象調(diào)節(jié)中起著至關(guān)重要的作用，特別是對外構(gòu)象的穩(wěn)定性。此外，Motif C 被稱為反轉(zhuǎn)運基序，在TMS5 內(nèi)形成紐扣，可能與底物結(jié)合位點的取向有關(guān)[36]。

圖5 PatC 的功能結(jié)構(gòu)域預(yù)測Fig.5 Prediction of functional domains of PatC

2.3 理化性質(zhì)分析

利用在線軟件ProtParam 對PatC 的理化性質(zhì)進行預(yù)測，結(jié)果顯示（表3），PatC 的氨基酸數(shù)量是546 個，分子量為58733.15；理論等電點PI：5.57，為酸性蛋白；該蛋白氨基酸序列中亮氨酸（Leu）、甘氨酸（Gly）及丙氨酸（Ala）三種氨基酸含量較高，所占比值分別為12.1%、10.1%和9.7%；蛋白不穩(wěn)定指數(shù)：40.3，為不穩(wěn)定蛋白；從脂肪系數(shù)數(shù)據(jù)分析表明該蛋白脂肪系數(shù)：113.77，為脂溶性蛋白。親水性平均系數(shù)（GRAVY）值為0.712，表明該蛋白為疏水性蛋白。

表3 PatC 的理化性質(zhì)分析Table 3 Physicochemical property analysis of PatC

2.4 跨膜結(jié)構(gòu)預(yù)測

采用在線軟件TMHMM 對PatC 的跨膜結(jié)構(gòu)進行預(yù)測，由圖6 和上述理化性質(zhì)結(jié)果表明該蛋白具有強疏水性氨基酸，具有能夠行使跨膜轉(zhuǎn)運的功能，符合轉(zhuǎn)運蛋白的特點。

圖6 PatC 跨膜區(qū)預(yù)測Fig.6 Transmembrane structure prediction of PatC

2.5 蛋白質(zhì)二級結(jié)構(gòu)預(yù)測

利用在線SOPMA 軟件對PatC 二級結(jié)構(gòu)進行預(yù)測和分析（圖7），PatC 中二級結(jié)構(gòu)主要存在4 種：α-螺旋、β-轉(zhuǎn)角、β-折疊和無規(guī)則卷曲。二級結(jié)構(gòu)元件中α-螺旋所占比例最高，此外無規(guī)則卷曲、β-折疊和β-轉(zhuǎn)角所占比例依次減少分別為30.95%、24.73%和6.23%。

圖7 PatC 的二級結(jié)構(gòu)預(yù)測Fig.7 Secondary structure prediction of PatC

2.6 蛋白質(zhì)三級結(jié)構(gòu)預(yù)測

PatC 同源性蛋白較少，無法利用SWISS MODEL同源建模方法得到其三級結(jié)構(gòu)。因此通過AlphaFold預(yù)測PatC 的三級結(jié)構(gòu)。由圖8 說明該蛋白是14 次α-螺旋結(jié)構(gòu)（與跨膜結(jié)構(gòu)預(yù)測相吻合），可能是進化過程中以12 次跨膜α-螺旋為基礎(chǔ)產(chǎn)生的[37]，分為2 個結(jié)構(gòu)域：N 端結(jié)構(gòu)域和C 端結(jié)構(gòu)域。每個結(jié)構(gòu)域都由6 個α-螺旋組成，2 個結(jié)構(gòu)域呈現(xiàn)二次贗對稱（twofold psudosymmetry）[20]。與大腸桿菌中的二肽轉(zhuǎn)運蛋白（Escherichia coli YbgH，EcYbgH）結(jié)構(gòu)類似，為經(jīng)典的12 次跨膜α-螺旋構(gòu)象加兩個TM 螺旋的插入（HA 和HB）[38]。

圖8 PatC 的三級結(jié)構(gòu)預(yù)測Fig.8 Tertiary structure prediction of PatC

2.7 分子對接

采用AutoDock Vina1.1.2 將PatC 與E-ascladiol進行半柔性對接，分子對接結(jié)果如圖9，其之間的平均結(jié)合能為-5.0 kcal/mol（＜0），說明配體和受體蛋白可以自發(fā)結(jié)合。PatC 與E-ascladiol 之間主要以傳統(tǒng)氫鍵（Conventional hydrogen bond）、堆積鍵（Pi-Pi Tshaped）及疏水作用力（Pi-Alkyl）相互作用。具體相互結(jié)合模式：E-ascladiol 與PatC 結(jié)合至蛋白空腔，且與該蛋白受體的氨基酸殘基SER353、TYR336、PRO339、及PRO188 分別形成2.74?氫鍵，2.96?Pi-Pi 堆積作用力，5.21?和529?疏水作用力，這些相互作用使得蛋白質(zhì)與小分子緊密結(jié)合。

圖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）是用來衡量一段時間內(nèi)分子結(jié)構(gòu)中原子位置的平均移動，是目的結(jié)構(gòu)特定原子相對于參照結(jié)構(gòu)在整個時間范圍內(nèi)的距離均方根，是判斷蛋白質(zhì)結(jié)構(gòu)結(jié)構(gòu)變化的重要參數(shù)之一[39]。

通過Protein（純蛋白）與Wild（野生型復(fù)合體）的RMSF 數(shù)據(jù)分析，如圖10 所示，PatC 在多個氨基酸位點發(fā)生柔性變化。與野生型復(fù)合物對比，純PatC 分別在第127、181、188～197、231～241、273、296、333～341、457～463 位殘基區(qū)域即（Ala127、Gly181、Pro188～Ser197、Gly231～Val241、Asp273、Phe296、Val333～Try341、Gln457～Ala463）發(fā)生了形變拉伸或縮短。其中位于Pro188～Ser197aa、Gly231～Val241aa 氨基酸區(qū)域內(nèi)，結(jié)合底物分子后的PatC 發(fā)生柔性變化，在一定程度上增加其靈活度，改變其空間構(gòu)象。因此可推測出上述氨基酸區(qū)域內(nèi)可能存在PatC 與E-ascladiol 之間的結(jié)合位點。

圖10 各模擬體系的RMSF 結(jié)果分析Fig.10 Analysis of RMSF results for each simulation system

2.8.2 點突變分析 通過上述分子對接和復(fù)合物分子動力學模擬兩種方法得出不同的結(jié)合位點。由圖11可知，MD 計算得到6 個結(jié)合位點，其中第188 位氨基酸與分子對接結(jié)果重疊，且與RMSF 數(shù)據(jù)相吻合，因此使用Pymol 軟件對PatC 中第188 位氨基酸殘基位點突變，將原有的脯氨酸（P）突變?yōu)楸彼幔ˋ），即產(chǎn)生新的P188A 的PatC 突變體。通過將PatC 與P188A 的突變位點進行構(gòu)象對比，如圖12，發(fā)現(xiàn)在第188 位氨基酸位置上的點突變并未使該位點的二級結(jié)構(gòu)出現(xiàn)明顯差別，且不影響整體蛋白質(zhì)構(gòu)型，因此繼續(xù)使用此結(jié)果進行后續(xù)實驗研究。

圖11 各個結(jié)合位點的總結(jié)合能Fig.11 Total binding energy of each binding site

圖12 PatC 及其突變體結(jié)構(gòu)Fig.12 Structure of PatC and its mutant

2.8.3 RMSD 分析 RMSD（Root Mean Square Distance/Deviation，RMSD）是用來衡量同一分子結(jié)構(gòu)不同構(gòu)象之間的差異，常用來分析體系分子整體結(jié)構(gòu)的穩(wěn)定性，RMSD 變小表明結(jié)構(gòu)趨于穩(wěn)定，反之趨于不穩(wěn)定[39]。

將純Protein（純蛋白）、Wild（野生型復(fù)合體）、P188A（突變體復(fù)合體）三個復(fù)合體系體系分別進行了20 ns 的分子動力學模擬實驗。在對其數(shù)據(jù)進行采集和分析中，三個模擬體系結(jié)構(gòu)中的原子骨架的根均方偏差RMSD 值如圖13 所示。當反應(yīng)進行至15 ns時，三個反應(yīng)體系曲線趨于平穩(wěn)，雖有一定的波動，但范圍在0.05 nm 之內(nèi)。說明三個體系均達到了平衡狀態(tài)，且三個反應(yīng)體系的RMSD 終值都停留在0.3 nm以下，其中Protein、Wild 復(fù)合體系和Mutant 復(fù)合體系的RMSD 平均值分別為0.17、0.23 和0.18 nm 左右，可看出三者在模擬過程中結(jié)構(gòu)均保持穩(wěn)定。且相較于Wild，Mutant 穩(wěn)定性更佳。

圖13 各模擬體系的RMSD 結(jié)果分析Fig.13 Analysis of RMSD results for each simulation system

2.8.4 Rg 分析 Rg（Radius of Gyratio，Rg）回旋半徑用來衡量分子結(jié)構(gòu)的緊密度或折疊程度，Rg 考察主要選擇平衡時狀態(tài)進行[39]，因此主要選擇后10 ns進行分析，由圖14 可得，后10 ns 兩個反應(yīng)體系中Rg 變化均很小，表明這三個反應(yīng)體系穩(wěn)定且在反應(yīng)過程中均為平衡狀態(tài)且復(fù)合物結(jié)構(gòu)致密，且三個反應(yīng)體系的平均值均在2.5 nm 以下，并且空白蛋白的Rg值最高，兩種結(jié)合底物分子的復(fù)合物體系中Rg 值均有明顯減少，由此可見在PatC 結(jié)合了底物分子后，整體構(gòu)象進行了收縮。

圖14 各模擬體系的Rg 結(jié)果分析Fig.14 Analysis of Rg results for each simulation system

3 結(jié)論

本研究基于PatC 與E-ascladiol 結(jié)合的前提下，運用生物信息學、分子對接和分子動力學模擬方法預(yù)測蛋白與分子之間的重要作用位點，以研究PatC 轉(zhuǎn)運機制。生物信息學結(jié)果表明，該蛋白的二級結(jié)構(gòu)以為α-螺旋和無規(guī)則卷曲為主，三級空間構(gòu)象含有14 個跨膜螺旋，含有2 個結(jié)構(gòu)域：N 端結(jié)構(gòu)域和C 端結(jié)構(gòu)域，形成獨特的MFS 折疊。根據(jù)以上結(jié)果說明PatC 作為毒素輸出泵轉(zhuǎn)運PAT，但其理化性質(zhì)也需進一步的驗證。分子對接結(jié)果表明PatC與E-ascladiol 之間存在多個結(jié)合位點與多種相互作用力，與RMSF 數(shù)據(jù)結(jié)合分析得出P188 作為重要結(jié)合位點并對其進行模擬突變。此外，本研究通過分子動力學方法對兩種復(fù)合體系（野生型Wild 和突變體P188A）進行各參數(shù)分析，結(jié)果說明P188 氨基酸位點可能是PatC 的關(guān)鍵作用位點，可為后續(xù)的微生物分子實驗提供基礎(chǔ)論證。綜上所述，本文為棒曲霉素合成途徑中的MFS 蛋白的轉(zhuǎn)運機制（毒素外排作用）提供一定的基礎(chǔ)，并進一步為毒素輸出泵的生物學功能與轉(zhuǎn)運機制研究提供參考。