楊菊,毛銀,黃曉強,周勝虎,鄧禹*

1(江南大學 生物工程學院,江蘇 無錫,214122)2(密歇根大學 計算醫(yī)學與生物信息系,密歇根 安娜堡,48109)

己二酸是一種重要的二元羧酸,廣泛應用于醫(yī)療衛(wèi)生、食品、化工等行業(yè)[1],據(jù)統(tǒng)計,超過60%的己二酸被用于生產(chǎn)尼龍類纖維(如尼龍6-6)[2-3]。目前,己二酸的大規(guī)模生產(chǎn)仍依賴化學合成,主要通過硝酸對環(huán)己醇-環(huán)己酮的混合物(醇酮油)進行氧化反應制取[4-5]?；瘜W合成法存在著工藝流程長、副產(chǎn)物(尤其是氮氧化物)排放多、產(chǎn)品收率低等問題。因此,人們迫切希望尋找可替代化學合成己二酸的新方法。目前,全生物法合成己二酸已成為可能[6]。

本課題組于2015年發(fā)現(xiàn)天然菌株ThermobifidafuscaB6中存在一條全生物合成己二酸的代謝途徑(圖1),并且將該途徑可移植到大腸桿菌中[7]。該途徑涉及5個酶：β-酮硫解酶(β-ketothiolase)、3-羥?；?輔酶A脫氫酶(3-hydroxyacyl-CoA dehydrogenase)、3-羥基己二酰-輔酶A脫氫酶(3-hydroxyadipyl-CoA dehydrogenase)、5-羰基-2-戊烯酰基-輔酶A還原酶(5-carboxy-2-pentenoyl-CoA reductase)和己二酰輔酶A合成酶(adipyl-CoA synthetase)[7]。總體來講,該代謝過程能通過消耗D-葡萄糖來合成己二酸,但產(chǎn)量較低(0.2 g/L)[8]。其中,5-羰基-2-戊烯?；?輔酶A還原酶是該途徑的限速酶,其催化5-羰基-2-戊烯?；?輔酶A(5-carboxy-2-pentenoyl-CoA)反應生成己二酰輔酶A(adipyl-CoA)的活性很低[7]。目前,關于改造5-羰基-2-戊烯?；?輔酶A還原酶以促進己二酸合成的研究尚未見報道。

近年來,基于結(jié)構(gòu)的計算酶設計方法被廣泛應用于改造天然酶以適應工業(yè)需要[9-14],甚至是從頭設計具有新功能的酶[15-17]。T.fusca5-羧基-2-戊烯酰-輔酶A還原酶的結(jié)構(gòu)尚未通過實驗解析出來。本研究首先通過同源建模構(gòu)建5-羧基-2-戊烯酰-輔酶A還原酶的結(jié)構(gòu)模型并對其進行分析,進而基于酶-配體復合物的結(jié)構(gòu),通過計算酶設計方法來改造5-羧基-2-戊烯酰-輔酶A還原酶的結(jié)合殘基,嘗試引入氫鍵網(wǎng)絡來增強突變酶與目的底物5-羧基-2-戊烯酰-輔酶A的結(jié)合能力以提高酶催化效率,并通過分子動力學模擬來檢驗設計的合理性。

圖1 己二酸的生物合成途徑Fig.1 Biosynthetic pathway for producing adipic acid

1 材料和方法

1.1 蛋白質(zhì)序列

5-羧基-2-戊烯酰-輔酶A還原酶包含385個氨基酸,其蛋白質(zhì)序列為：

MSDFDLYRPTEEHEALREAIRSVAEDKIAPHAADVDEQ-SRFPQEAYEALRASDFHAPHVAEEYGGVGADALATCIVIEEIA-RVCASSSLIPAVNKLGSMPLILSGSDEVKQRYLPELASGEAMFS-YGLSEREAGSDTASMRTRAVRDGDDWILNGQKSWITNAGISK-YYTVMAVTDPDGPRGRNISAFVVHIDDPGFSFGEPERKLGIKG-SPTRELIFDNVRIPGDRLVGKVGEGLRTALRTLDHTRVTIGAQ-AVGIAQGALDYALGYVKERKQFGKAIADFQGIQFMLADMAM-KLEAARQMVYVAAAKSERDDADLSFYGAAAKCFASDVAMEI-TTDAVQLLGGYGYTRDYPVERMMRDAKITQIYEGTNQIQRVV-MARQLLKK

1.2 同源建模

蛋白質(zhì)數(shù)據(jù)庫(protein data bank,PDB)[18]中沒有5-羧基-2-戊烯酰-輔酶A還原酶的實驗結(jié)構(gòu),因此本文通過同源建模來構(gòu)建其結(jié)構(gòu)模型。對上述序列進行BLASTp[19]搜索(采用默認參數(shù)),發(fā)現(xiàn)PDB中存在多個同源結(jié)構(gòu)。因為計算酶設計過程中需要同時對配體分子(FAD和5-羧基-2-戊烯酰-輔酶A)進行建模,所以本文選擇包含了FAD和類似于5-羧基-2-戊烯酰-輔酶A配體的實驗結(jié)構(gòu)作為模板。將上述蛋白質(zhì)序列提交到SWISS-MODEL服務器[20],搜索模板,結(jié)果顯示有超過50個結(jié)構(gòu)模板可供使用。本文選擇分辨率1.8?的晶體結(jié)構(gòu)4L1F_A作為模板建模[21],因為該結(jié)構(gòu)中包含了配體分子FAD和過硫化輔酶A(CoA persulfide,PDB文件中對應名字為COS,與目標底物5-羧基-2-戊烯酰-輔酶A結(jié)構(gòu)相似),并且與5-羧基-2-戊烯酰-輔酶A還原酶具有46.28%的序列相同性。

1.3 分子動力學模擬

本文使用GROMACS 5.4.0軟件進行分子動力學模擬研究[22]。首先使用其中的pdb2 gmx模塊構(gòu)建蛋白模型的拓撲文件,蛋白的力場選擇GROMOS96 43A1聯(lián)合原子力場。然后將蛋白-配體復合物模型置于一個充滿TIP3P水模型的十二面體盒子中心,盒子邊緣距離蛋白分子邊緣為1 nm。向盒子中添加適量的Na+以平衡體系電荷。然后利用grompp模塊完成體系的能量最小化。先將體系至于NVT(T=300 K)條件下平衡100 ps,然后置于NPT條件下平衡100 ps。在NVT和NPT過程中,對蛋白和配體分子上的重原子加入位置約束(1 000 kJ/mol)[23]。最后,參考文獻[24-25],執(zhí)行無約束分子動力學模擬,時長16 ns。

2 結(jié)果與分析

2.1 構(gòu)建蛋白-配體復合物結(jié)構(gòu)模型

由圖2可知,本研究構(gòu)建的5-羧基-2-戊烯酰-輔酶A還原酶模型(藍綠色)具有典型的α/β折疊結(jié)構(gòu)。因為該結(jié)構(gòu)模型與模板4L1F_A處于相同的坐標系,將4L1F_A中的FAD和COS配體分子復制到5-羧基-2-戊烯酰-輔酶A還原酶結(jié)構(gòu)模型中。需要指出的是,COS并不能參與反應,它不是一個可以參與反應的底物。如圖3所示,COS與目的底物5-羧基-2-戊烯酰-輔酶A具有相同的母核結(jié)構(gòu),但其側(cè)鏈不同。COS的側(cè)鏈為巰基,而5-羧基-2-戊烯酰-輔酶A的側(cè)鏈為5-羧基-2-戊烯酰基。本文使用蛋白質(zhì)設計工具EvoEF2[26]的配體生成模塊將配體COS修改為底物5-carboxy-2-pentenoyl-CoA(圖2)。考慮到5-羧基-2-戊烯?；鶄?cè)鏈具有較大的柔性,配體側(cè)鏈生成過程中對可自由旋轉(zhuǎn)的二面角進行離散(-180°～180°,間隔30°離散),一共生成20 726個底物分子構(gòu)象。放置過程中，同時計算出每個底物構(gòu)象的內(nèi)能及其與蛋白質(zhì)骨架的范德華排斥能(兩者均為正值,越小越好)。接下來對所有構(gòu)象分別按照2種能量由低到高排序,共有1 728個構(gòu)象，在2種能量排序中都處于前25%。最后,按照適當?shù)木礁嚯x偏差(root mean squared error,RMSD,單位為?)對1 728個構(gòu)象進行篩選去重。本文采用的RMSD閾值為0.2 ?,篩選后得到48個底物分子構(gòu)象；其中與酶骨架排斥能最小的分子如圖2所示。

圖2 5-羧基-2-戊烯酰-輔酶A還原酶的模型Fig.2 Structure model of 5-carboxy-2-pentenoyl-CoA reductase in complex with FAD and 5-carboxy-2-pentenoyl-CoA注：藍綠色為5-羧基-2-戊烯酰-輔酶A還原酶,粉紅色球棍模型為FAD配體,淺藍色球棍模型為5-羧基-2-戊烯酰-輔酶A,綠色棍狀模型為結(jié)合殘基

圖3 三種輔酶A配體的比較Fig.3 Comparison of three CoA ligands注:過硫化輔酶A(CoA persulfide,COS)為PDB結(jié)構(gòu)4L1F中的配體；乙酰輔酶A(acetyl-CoA)為推測的天然底物；5-羧基-2-戊烯酰-輔酶A(5-carboxy-2-pentenoyl-CoA)為目的底物

同源建模過程中搜索到的結(jié)構(gòu)模板絕大部分被注釋為乙酰輔酶A脫氫酶(acyl-CoA dehydrogenase)。其中與5-羧基-2-戊烯酰-輔酶A還原酶序列最相近為Mycobacteriumthermoresistibile乙酰輔酶A脫氫酶(PDB ID:3PFD)[27]；兩者的序列相同性達到73.35%。據(jù)此推斷,T.fusca5-羧基-2-戊烯酰-輔酶A還原酶的天然底物也可能為乙酰輔酶A或其類似物。從圖3可以看出,目的底物5-羧基-2-戊烯酰-輔酶A與乙酰輔酶A的側(cè)鏈有較大差異,5-羧基-2-戊烯?；^長,有更強的柔性,且包含強極性的羧基。從圖2可知,側(cè)鏈附近的結(jié)合殘基大多數(shù)為疏水殘基。因此,可推測這些疏水基團不能與5-羧基-2-戊烯酰基側(cè)鏈很好地結(jié)合,從而導致酶催化5-羧基-2-戊烯酰-輔酶A反應的活性較低。

2.2 計算設計改造5-羧基-2-戊烯酰-輔酶A還原酶

基于上述分析,本研究選擇與5-羧基-2-戊烯酰-輔酶A側(cè)鏈直接接觸(5?范圍內(nèi))的氨基酸殘基作為設計位點。這樣的殘基共有15個,分別為Thr 161、Ile 90、Leu 89、Ser 88、Pro 91、Ala 92、Val 93、Lys 95、Leu 96、Ile 250、Gln 253、Thr 249、Tyr 367、Thr 246和Glu 368(圖4)。其中Glu 368為催化殘基,設計過程中不改變其殘基類型。此外,本文選擇距離5-羧基-2-戊烯?；??范圍內(nèi)的其他氨基酸殘基作為側(cè)鏈重新安裝的位點(未顯示)；這些殘基的側(cè)鏈構(gòu)象可改變,但氨基酸類型不變。從圖4可以看出,14個設計位點中只有Ser 88、Lys 95、Thr 161、Thr 246、Thr 249 和 Gln 253為親水性殘基,其余9個位點為疏水性殘基,并且所有殘基均未能與5-羧基-2-戊烯酰基側(cè)鏈的羧基官能團形成氫鍵或鹽橋。

圖4 5-羧基-2-戊烯酰-輔酶A側(cè)鏈5?范圍內(nèi)的氨基酸殘基Fig.4 Amino acids within 5? to the side chain of 5-carboxy-2-pentenoyl-CoA

鑒于5-羧基-2-戊烯酰-輔酶A與乙酰輔酶A的差異(圖3),本研究嘗試通過計算設計引入氫鍵網(wǎng)絡來穩(wěn)定極性的5-羧基-2-戊烯?；ＴO計過程通過計算蛋白質(zhì)設計工具EvoEF2[26, 28]完成。每個蛋白質(zhì)設計位點均可在20種天然氨基酸類型中自由選擇,每個氨基酸的側(cè)鏈構(gòu)象取自SHAPOVALOV等[29]開發(fā)的蛋白質(zhì)骨架依賴旋轉(zhuǎn)異構(gòu)體庫(backbone-dependent rotamer library)；而5-羧基-2-戊烯酰-輔酶A配體的空間結(jié)構(gòu)從前面生成的48個構(gòu)象中自由選擇。計算酶設計的過程相當于求解組合優(yōu)化問題,即從每個位點(包括配體分子)選擇一種殘基類型(或構(gòu)象),使得酶設計體系的整體能量達到最低。由于需要搜索的構(gòu)象空間巨大,EvoEF2采用模擬退火優(yōu)化算法來搜索低能量的氨基酸序列,但并不一定保證能取得全局能量最低序列。由于設計過程存在一定隨機性,每次設計的低能量序列可能并不完全相同,因此可以通過多次獨立運行EvoEF2計算程序來產(chǎn)生許多備選序列進行分析。EvoEF2自動輸出每個設計的結(jié)構(gòu)和體系能量(單位為EvoEF2 energy unit,簡記為EEU)。

如表1所示,本研究生成10個不同的設計(記為Des0～Des9)作為對照,對天然序列(WT)的設計位點的側(cè)鏈構(gòu)象進行重新安裝,并計算其體系能量。不同設計的能量值存在一定差異,說明EvoEF2每次優(yōu)化求解的過程可能收斂于不同的解,但能量差異并不大,說明EvoEF2優(yōu)化算法收斂較好。10個設計的體系能量(-1462.58～-1452.97 EEU)明顯低于WT的能量(-1247.53 EEU)。這說明天然酶與底物5-羧基-2-戊烯酰-輔酶A的結(jié)合作用可能較弱,而通過設計產(chǎn)生突變則可能引入一些有利的結(jié)合作用力。從表1來看,除了92、95和246位點十分保守(始終選擇對應的WT氨基酸類型)之外,其余位點都是可變的(可能選擇非WT氨基酸類型)。

表1 基于計算改造5-羧基-2-戊烯酰-輔酶A還原酶產(chǎn)生的10個設計Table 1 Ten designs for computational engineering the 5-carboxy-2-pentenoyl-CoA reductase

仔細分析設計的結(jié)構(gòu)發(fā)現(xiàn)Des0、Des3、Des4和Des9這4個設計中底物分子的羧基可與氨基酸殘基形成氫鍵(或鹽橋)網(wǎng)絡。如圖5所示,Des0中形成了如下氫鍵：R253的胍基和Q250側(cè)鏈的胺基可與配體羧基的其中一個氧原子分別形成距離為3.0?和2.8?的鹽橋或氫鍵,S89的羥基可與R253的胍基形成距離為2.6?的氫鍵,并且R253胍基可與T364(非設計位點)的羥基形成一個較弱的氫鍵(3.7?)。Des9中引入的氫鍵網(wǎng)絡與Des0類似；形成氫鍵網(wǎng)絡的氨基酸殘基完全相同,配體5-carboxy-2-pentenoyl-CoA的羧基側(cè)鏈構(gòu)象稍有不同。Des3和Des4中,位點89選擇為蘇氨酸(T)；與Des0和Des9相比,Des3和Des4中的T89取代S89與R253的胍基形成氫鍵。

圖5 Des0設計中的氫鍵網(wǎng)絡Fig.5 Hydrogen-bonding network in Des0

2.3 基于分子動力學模擬檢驗設計的氫鍵網(wǎng)絡

計算酶設計過程中假定蛋白質(zhì)骨架固定不變,但在實際酶反應過程中,酶、配體及其相互作用都可動態(tài)變化。本研究通過分子動力學模擬來檢驗設計的氫鍵在動態(tài)過程中能否穩(wěn)定存在,以此來評估設計的優(yōu)劣。統(tǒng)計分析設計的氫鍵在整個動力學模擬過程中的變化,氫鍵具有方向性,統(tǒng)計氫鍵一般考慮3個參數(shù)：XHY鍵角α,HY距離d1和XY距離d2(圖6)。在一定范圍內(nèi),α越大、d1和d2越小氫鍵越強。通常認為α>130°,d1<2.2?,d2<3.2?可算作穩(wěn)定的氫鍵[24-25]。

圖6 氫鍵及相關幾何參數(shù)圖示Fig.6 Schematic of hydrogen bond and related geometric parameters

對分子動力學模擬過程進行統(tǒng)計分析表明,Des0設計中構(gòu)建的氫鍵網(wǎng)絡較穩(wěn)定。Des0中引入的4個氫鍵統(tǒng)計結(jié)果如圖7所示。

a～d-不同氫鍵的分析結(jié)果：第一列描述所要分析的氫鍵,第二列描述氫鍵距離隨模擬時間的變化情況,第三列描述氫鍵距離與角度的關系圖7 Des0設計中引入的氫鍵在分子動力學模擬過程中的統(tǒng)計結(jié)果Fig.7 Statistical result of hydrogen bonds introduced in Des0 during 16 ns MD simulation

在16 ns的模擬過程中,5-carboxy-2-pentenoyl-CoA羧基的其中一個氧原子可以與R253的胍基以及Q250的胺基形成穩(wěn)定的氫鍵(圖7-a和b,d1平均值約為2?,α平均值約為160°)。相較之下,S89的羥基與R253的胍基之間的氫鍵較弱,d1平均值約為3?,α平均值約為110°(圖7-c)。這可能是因為T364的羥基可競爭性地與R253的胍基形成較強的氫鍵(圖7-d,d1平均值約為2?,α平均值約為170°)。從圖5來看,Des0設計中R253可與S89形成的氫鍵較強,而與T364的氫鍵較弱。在動力學模擬過程中,這2個氫鍵的強弱雖然發(fā)生了轉(zhuǎn)換,但R253至少可與S89和T364其中之一形成較強的氫鍵。因此,分子動力學模擬結(jié)果表明部分設計(如Des0)中的氫鍵網(wǎng)絡可以起到穩(wěn)定5-羧基-2-戊烯酰-輔酶A柔性側(cè)鏈的作用。這些設計可用于進一步實驗驗證。

3 結(jié)論

菌株ThermobifidafuscaB6中的野生型5-羰基-2-戊烯?；?輔酶A還原酶(5-carboxy-2-pentenoyl-CoA reductase)為己二酸生物合成途徑中的決速酶。其主要原因是目標底物5-羰基-2-戊烯酰基-輔酶A與其天然底物乙酰輔酶A相比柔性更大、極性更強。本文嘗試采用基于結(jié)構(gòu)的計算酶設計方法對野生型5-羰基-2-戊烯?；?輔酶A還原酶的結(jié)合位點進行改造,試圖通過設計引入氫鍵網(wǎng)絡來更好地結(jié)合5-羰基-2-戊烯?；?輔酶A的極性側(cè)鏈,目的是降低酶促反應的活化能以提高催化效率。為檢驗設計的合理性,本研究通過分子動力學模擬來觀察10個設計中氫鍵網(wǎng)絡的變化情況,結(jié)果發(fā)現(xiàn)Des0設計中的氫鍵網(wǎng)絡很穩(wěn)定,與野生型酶相比,可增強與5-羰基-2-戊烯?；?輔酶A的結(jié)合作用。因此,Des0突變可用于進一步實驗驗證。