Hsp70核苷酸結(jié)合域突變體影響酶活的分子動力學(xué)模擬研究

周小紅，徐利楠，薛友林，宋有濤，2*

(1.遼寧大學(xué)生命科學(xué)院，遼寧沈陽110036;2.遼寧大學(xué)環(huán)境學(xué)院，遼寧沈陽110036;3.遼寧大學(xué)輕型產(chǎn)業(yè)學(xué)院，遼寧沈陽110036)

Hsp70家族是細(xì)胞內(nèi)一類高度保守的分子伴侶蛋白，可以與其它分子伴侶共同參與多種重要的生命活動，如非天然多肽的折疊、多聚體的組裝、細(xì)胞器蛋白的轉(zhuǎn)運及錯誤折疊蛋白的降解等［1］。所有的Hsp70家族都由一個氮端44 kDa的核苷酸結(jié)合域(NBD，包括 NBD-I，NBD-II)和一個碳端25 kDa的底物結(jié)合域(SBD，Substrate-binding domain，包括SBD-α，SBD-β)組成［2］。

此前，McKay等人的研究發(fā)現(xiàn)牛Hsc70(Heat shock cognate 70 protein)的ATP酶活性與金屬離子及活性位點 T13、T204的自身磷酸化有重要關(guān)系［3］，其中，T13側(cè)鏈的羥基對 ATP酶活性具有重要作用［4］。此外，Chiappori等人研究指出，loop1的柔性對ATP酶活性具有重要作用，柔性越強ATP酶活性越高［5］。這些結(jié)果暗示著loop1的構(gòu)象改變可能引起Hsp70氮端ATP酶活的變化。隨后，Jones等人的研究指出酵母Hsp70氮端NBD存在大量影響ATP酶活性的位點［6-7］，Chang等人在 DnaK(大腸桿菌Hsp70)中得出相同結(jié)果［8］。但是，Jones和Chang等人都沒有從分子水平去探討loop1的構(gòu)象改變和ATP酶活的關(guān)系。

為了從分子水平上探討loop1構(gòu)象的改變與ATP酶活的關(guān)系，在已有研究表明同源的突變體功能相似的基礎(chǔ)上［9］，我們選用同源性很高的DnaK進(jìn)行動力學(xué)模擬(酵母Hsp70晶體結(jié)構(gòu)未知)，研究了點突變可能通過loop1的構(gòu)象改變而引起ATP酶活變化——酵母 Hsp70氮端 NBD中 A17V(DnaKA17V)、R23H(DnaKT23H)、G32S(DnaKG32S)增強ATP酶活、R167H(DnaKR167H)對 ATP酶活無影響［6］。該研究不僅是對Jones等人生化試驗數(shù)據(jù)的分子解釋，而且對后續(xù)Hsp70氮端NBD的點突變研究具有重要的借鑒意義。

1 實驗過程

1.1 蛋白質(zhì)的來源及模型構(gòu)建

酵母Hsp70氮端NBD初始三維結(jié)構(gòu)模型來自同源蛋白 E.coli DnaK(RCSB數(shù)據(jù)庫，PDB編號1DKG:D)［10］。將已知的野生型DnaK NBD晶體結(jié)構(gòu)模型與對應(yīng)的17位、23位、32位、167位突變體一級序列提交至 Swiss Model服務(wù)器上(http://swissmodel.Expasy.org)［11］，根據(jù)同源建模算法得到所有突變型蛋白的三維結(jié)構(gòu)，將該模型作為后續(xù)動力學(xué)模擬突變型A17V、R23H、G32S、R167H NBD的初始模型。

1.2 動力學(xué)模擬過程

本模擬利用由動力學(xué)模擬軟件Gromacs 4.5.5［12］和三維結(jié)構(gòu)分析軟件 PyMOL 進(jìn)行。各模型的模擬步驟及參數(shù)設(shè)定均相同。具體步驟及參數(shù)如下:

(1)力場選擇:根據(jù)不同力場偏向于不同的結(jié)構(gòu)狀態(tài)模擬，該模擬各體系均采用NPT系綜，計算力場選用GROMOS96 43a1。

(2)溶解蛋白:將結(jié)構(gòu)模型放進(jìn)一個立方體的周期性邊界盒子中并置于該盒子的中心，設(shè)定模型中任一原子距盒子邊緣的距離大于1.0 nm，以保證其不會和另一個盒子中的鏡像相互作用。水溶液環(huán)境采用SPC模型水分子(約含25 169個水分子)。

(3)平衡電荷:正常的生理條件下，體系的pH值應(yīng)接近中性，總電荷為0。向各溶液體系中加入不同數(shù)量的Na+或Cl-，保持體系 pH=7，總電荷為0。

(4)能量最小化:為了消除可能的原子碰撞，根據(jù)蛋白質(zhì)折疊的熱力學(xué)，先將各模型在真空中以最陡下降法(steep)優(yōu)化2 000步的能量使其達(dá)到結(jié)構(gòu)穩(wěn)定狀態(tài)。

(5)限制性模擬:根據(jù)固定主鏈，優(yōu)化側(cè)鏈的原則，為了使側(cè)鏈找到一個較穩(wěn)定的狀態(tài)避免出現(xiàn)溶劑分子分布不均的問題，進(jìn)行了80 ps的限制性模擬。

(6)最終的動力學(xué)模擬:達(dá)到上述平衡之后，各體系在恒定的溫度和壓力下進(jìn)行10 ns的分子動力學(xué)模擬。分子的運動軌跡每4 ps收集一次，以便用于后續(xù)的數(shù)據(jù)分析。

此外，為了保持整個模擬過程中體系溫度穩(wěn)定在300 K、壓力穩(wěn)定在 1.0×105Pa，分別采用“V-rescale”和“Parrinello-Rahman”方法對體系溫度和壓力進(jìn)行控制，耦合常數(shù)分別為0.1 ps和0.5 ps。模擬的時間積分步長為2fs，整個體系使用周期性邊界條件，采用LINCS(linear constraint solver)算法對體系中所有鍵的鍵長進(jìn)行約束，靜電相互作用采用PME(particle mesh Ewald)算法進(jìn)行估算，截斷半徑為1.0 nm。范德華力通過“cut-off”方法進(jìn)行估算，截斷半徑為1.0 nm。

2 結(jié)果與討論

2.1 結(jié)構(gòu)穩(wěn)定性

RMSD(Root mean square deviation，均方根偏差)可用來測量進(jìn)行比對的兩個蛋白質(zhì)相同位點的氨基酸殘基碳原子之間的距離，是衡量模擬過程中蛋白質(zhì)穩(wěn)定程度的重要參數(shù)，間接反應(yīng)蛋白質(zhì)靈活性。通過比較整體RMSD曲線的變化趨勢得出在7 ns時各模型均達(dá)到平衡狀態(tài)(見圖1a)，結(jié)構(gòu)穩(wěn)定。對loop1穩(wěn)定性分析發(fā)現(xiàn)，各模型在7 ns時達(dá)到平衡(見圖1b)，平均值分別為:WT(0.51 nm)、R167H(0.69 nm)、G32S(0.81 nm)、A17V(0.84 nm)、T23H(0.75 nm)。從這些數(shù)據(jù)可以明顯看出突變模型A17V、T23H、G32S中l(wèi)oop1的RMSD平均值均比WT、R167H高。表明突變模型 A17V、T23H、G32S的柔性比 WT、R167H 強，暗示 A17、T23、G32三個位點可能通過loop1構(gòu)象變化，引起酵母Hsp70氮端NBD的ATP酶活性變化。

圖1 (a)蛋白RMSD隨時間變化(b)loop1RMSD隨時間變化Fig.1 (a)The protein RMSD as a function of simulation time(b)The loop1 RMSD as a function of simulation time

2.2 Hsp70氮端ATP酶活性區(qū)域結(jié)構(gòu)對比

通過三維可視軟件PyMOL，標(biāo)示出loop1和重要殘基T11在E.coli DnaK氮端NBD的位置(見圖2)。

圖2 DnaK NBD三維結(jié)構(gòu)Fig.2 DnaK NBD three-dimensional structural

為進(jìn)一步分析loop1怎樣的構(gòu)象改變而引起Hsp70氮端NBD的ATP酶活性變化，我們對A17V、T23H、G32S、R167H、WT 進(jìn)行了 10 ns時整體蛋白三維結(jié)構(gòu)對比分析。結(jié)果顯示，突變模型R167H與WT在loop1處具有相似的構(gòu)象，重要殘基T11側(cè)鏈未向內(nèi)發(fā)生移動(見圖 3a)。突變模型 A17V、T23H、G32S與WT結(jié)構(gòu)對比發(fā)現(xiàn)，T11的側(cè)鏈均向內(nèi)發(fā)生不同程度移動(見圖3b，圖c，圖d)，預(yù)示著在酵母Hsp70 ATP水解過程中，A17V、T23H、G32S三個突變模型與WT和R167H相比能夠更快速的通過T11側(cè)鏈的羥基與ATP進(jìn)行相互作用，攻擊ATPγ-磷酸基團(tuán)，引起 ATP快速水解，從而增強ATP酶活性。

圖3 模擬后loop1結(jié)構(gòu)對比:WT(T11:Sticks)(a)R167H(T11:lines)(b)A17V(T11:lines)(c)T23H(T11:lines)(d)G32S(T11:lines)Fig.3 loop1 structural alignment after MD simulation:WT(T11:Sticks)(a)R167H(T11:lines)(b)A17V(T11:lines)(c)T23H(T11:lines)(d)G32S(T11:lines)

2.3 Hsp70 ATP酶活性區(qū)域距離分析

為了深入探討前面提到的構(gòu)象改變與功能變化關(guān)系，我們引入了一個ATP分子(RCSB數(shù)據(jù)庫，PDB編號:3LDL)進(jìn)行距離分析。McKay等人研究指出，Hsc70 T13(DnaKT11)殘基對ATP水解的重要作用，是通過T13側(cè)鏈的羥基對ATPγ-磷酸基團(tuán)作用，引起ATP酶水解速率變化［5］。數(shù)據(jù)顯示，10 ns時T11側(cè)鏈的羥基與γ-磷酸基團(tuán)距離(見圖4a，b，c，d，e)WT(5.5 ?)，R167H(6.0 ?)均大于突變模型 A17V(1.3 ?)、T23H(4.9 ?) 和 G32S(4.6 ?)。數(shù)據(jù)表明，T11側(cè)鏈的羥基與γ-磷酸基團(tuán)距離變近，其相互作用增強，可能引起ATP酶活性增強，與前面結(jié)論相一致。從分子水平初步解釋了生化試驗酵母Hsp70 A17V，R23H，G32S點突變增強ATP酶活結(jié)論。

2.4 loop1 氫鍵分析

為了探討引起loop1發(fā)生向內(nèi)移動的原因，我們對loop1與周圍氨基酸的氫鍵交互作用進(jìn)行了分析，通過計算得到各突變體模型10ns模擬過程中氫鍵數(shù)量的平均值分別為:WT(1.2個)、R167H(1.2個)、A17V(2.6 個)、T23H(2.4 個)、G32S(0.4個)。結(jié)果表明，A17V、T23H突變體模型在10ns模擬過程中的氫鍵數(shù)量顯著高于WT、R167H，然而有趣的是G32S突變體模型在10ns模擬過程中的氫鍵數(shù)量低于WT、R167H(見圖5a)。隨后氫鍵存活時間分析發(fā)現(xiàn)(見圖5b)，在10ns模擬過程中穩(wěn)定氫鍵數(shù)量分別為WT(1個)、R167H(1個)、A17V(3個)、T23H(3個)、G32S(1個)。以上數(shù)據(jù)表明A17V、T23H點突變后，引起loop1中與T12相互作用的氫鍵數(shù)量顯著增加，增強了loop1附近的氫鍵網(wǎng)絡(luò)，從而間接調(diào)控loop1中T11側(cè)鏈向內(nèi)移動，引起ATP酶活性的增強;值得注意的是，G32S點突變后盡管減少了loop1氫鍵數(shù)量，但是氫鍵存活時間表明，在沒有與T12相互作用的氫鍵存在時，每隔一段時間會重復(fù)出現(xiàn)T11-S38氫鍵，有規(guī)律的影響T11殘基的狀態(tài)，可能直接引起T11側(cè)鏈向內(nèi)移動，從而增強ATP酶活性。

圖4 模擬后T11殘基與γPi距離(單位:?)(a)WT(b)R167H(c)A17V(d)T23H(e)G32SFig.4 the distance of T11 and γPi after MD simulation(?)(a)WT(b)R167H(c)A17V(d)T23H(e)G32S

圖5 (a)模擬過程中l(wèi)oop1氫鍵數(shù)量(b)模擬過程中l(wèi)oop1氫鍵存活時間Fig.5 (a)the number of hydrogen bond during the MD simulation(b)Hydrogen bond existence map of loop1 during the MD simulation

3 結(jié)論

Hsp70分子伴侶SBD和NBD的功能既相互聯(lián)系又相互獨立，分析其中一個結(jié)構(gòu)域的構(gòu)象功能變化，能夠間接反應(yīng)出整個Hsp70分子伴侶的功能變化。在本研究中，A17V，T23H，G32S點突變引起Hsp70氮端NBD ATP酶活性區(qū)域的ATP結(jié)合口袋袋口Loop1柔性增強，重要殘基T11側(cè)鏈向內(nèi)移動，接近ATPγ-磷酸基團(tuán)，增強其與ATP的作用，提高ATP水解速率，從而增強ATP酶活性，引起Hsp70分子伴侶功能變化。

本研究為我們從分子水平上解釋Jones等人的生化試驗結(jié)果(A17V，R23H和G32S點增強ATP酶活，R167H與WT對 ATP酶活沒有影響［6］)提供了數(shù)據(jù)支持。同時也從分子水平對 Hsc70 T13(DnaKT11)側(cè)鏈的羥基對ATP酶具有重要作用［4］進(jìn)行了一定驗證，對后續(xù)Hsp70氮端NBD的突變研究具有重要的借鑒意義。

References)

［1］ Sharma D，Masison DC.Hsp70 Structure，F(xiàn)unction，Regulation and Influence on Yeast Prions［J］.Protein Peptide Letters，2009，16(6):571-581.

［3］ Johnson ER，McKay DB.Mapping the role of active site residues for transducing an ATP-induced conformational change in the bovine 70-kDa heat shock cognate protein［J］.Biochemistry，1999，38(33):10823-10830.

［4］ Sousa MC，McKay DB.The hydroxyl of threonine 13 of the bovine 70-kDa heat shock cognate protein is essential for transducing the ATP-induced conformational change［J］.Biochemistry，1998，37(44):15392-15399.

［5］ Chiappori F，Merelli I，Colombo G，et al.Molecular Mechanism of AllostericCommunication in Hsp70 Revealed by Molecular Dynamics Simulations［J］.PLoS computational biology，2012，8(12):e1002844.

［6］ Jones GW，Masison DC.Saccharomyces cerevisiae Hsp70 Mutations Affect［PSI+］Prion Propagation and Cell Growth Differently and Implicate Hsp40 and Tetratricopeptide Repeat Cochaperones in Impairment of［PSI+］［J］.Genetics，2003，163(2):495-506.

［7］ Loovers HM，Guinan E，Jones GW.Importance of the Hsp70 ATPase Domain in Yeast Prion Propagation［J］.Genetics，2007，175(2):621-630.

［8］ Chang L，Thompson AD，Ung P，et al.Mutagenesis reveals the complex relationships between ATPase rate and the chaperone activities of Escherichia coli heat shock protein 70(Hsp70/DnaK)［J］.Journal of Biological Chemistry，2010，285(28):21282-21291.

［9］ Xu L，Hasin N，Shen M，et al.Using Steered Molecular Dynamics to Predict and Assess Hsp70 Substrate-Binding Domain Mutants that Alter Prion Propagation［J］.PLoS computational biology，2013，9(1):e1002896.

［10］ Harrison CJ，Hayer-Hartl M，Di Liberto M，et al.Crystal structure of the nucleotide exchange factor GrpE bound to the ATPase domain of the molecular chaperone DnaK［J］.Science，1997，276(5311):431-435.

［11］ Kiefer F，Arnold K，Künzli M，et al.The SWISS-MODEL Repository and associated resources［J］.Nucleic acids research，2009，37(suppl 1):D387-D392.

［12］ Pronk S，Páll S，Schulz R，et al.GROMACS 4.5:a highthroughput and highly parallel open source molecular simulation toolkit［J］.Bioinformatics，2013，29(7):845-854.