計算設計改造Thermobifida fusca 5-羧基-2-戊烯酰-輔酶A還原酶促進己二酸生產

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

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

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

本課題組于2015年發現天然菌株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-葡萄糖來合成己二酸,但產量較低(0.2 g/L)[8]。其中,5-羰基-2-戊烯酰基-輔酶A還原酶是該途徑的限速酶,其催化5-羰基-2-戊烯酰基-輔酶A(5-carboxy-2-pentenoyl-CoA)反應生成己二酰輔酶A(adipyl-CoA)的活性很低[7]。目前,關于改造5-羰基-2-戊烯酰基-輔酶A還原酶以促進己二酸合成的研究尚未見報道。

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

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

1 材料和方法

1.1 蛋白質序列

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

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

1.2 同源建模

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

1.3 分子動力學模擬

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

2 結果與分析

2.1 構建蛋白-配體復合物結構模型

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

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

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

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

基于上述分析,本研究選擇與5-羧基-2-戊烯酰-輔酶A側鏈直接接觸(5?范圍內)的氨基酸殘基作為設計位點。這樣的殘基共有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-戊烯酰基7?范圍內的其他氨基酸殘基作為側鏈重新安裝的位點(未顯示)；這些殘基的側鏈構象可改變,但氨基酸類型不變。從圖4可以看出,14個設計位點中只有Ser 88、Lys 95、Thr 161、Thr 246、Thr 249 和 Gln 253為親水性殘基,其余9個位點為疏水性殘基,并且所有殘基均未能與5-羧基-2-戊烯酰基側鏈的羧基官能團形成氫鍵或鹽橋。

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

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

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

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

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

圖5 Des0設計中的氫鍵網絡Fig.5 Hydrogen-bonding network in Des0

2.3 基于分子動力學模擬檢驗設計的氫鍵網絡

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

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

對分子動力學模擬過程進行統計分析表明,Des0設計中構建的氫鍵網絡較穩定。Des0中引入的4個氫鍵統計結果如圖7所示。

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

在16 ns的模擬過程中,5-carboxy-2-pentenoyl-CoA羧基的其中一個氧原子可以與R253的胍基以及Q250的胺基形成穩定的氫鍵(圖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個氫鍵的強弱雖然發生了轉換,但R253至少可與S89和T364其中之一形成較強的氫鍵。因此,分子動力學模擬結果表明部分設計(如Des0)中的氫鍵網絡可以起到穩定5-羧基-2-戊烯酰-輔酶A柔性側鏈的作用。這些設計可用于進一步實驗驗證。

3 結論

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