2021, Vol. 37中国科学院微生物研究所、中国微生物学会主办
文章信息
- 陈权, 吕成, 许菲
- Chen Quan, Lü Cheng, Xu Fei
- 基于柔性区域的计算设计改造玉米赤霉烯酮水解酶热稳定性
- Computation-aided design of the flexible region of zearalenone hydrolase improves its thermal stability
- 生物工程学报, 2021, 37(12): 4415-4429
- Chinese Journal of Biotechnology, 2021, 37(12): 4415-4429
- 10.13345/j.cjb.210067
-
文章历史
- Received: January 21, 2021
- Accepted: April 8, 2021
引用本文
|
玉米赤霉烯酮(Zearalenone，ZEN) 是造成全球粮食污染的霉菌毒素之一，对饲料行业和养殖行业带来巨大损失[1]。此外，ZEN的生殖发育毒性[2]和致癌性[3]严重危害牲畜以及人类健康。2013年，程传民等[4]对我国玉米副产物中ZEN的含量进行了抽样调查，发现其平均含量为678 μg/kg，远高于饲料卫生国家标准(GB13078-2017) 中[5] ZEN的最高允许量500 μg/kg。同时，玉米副产物中ZEN的检出率高达95.6%，其中中度污染率(检出量为500–3 000 μg/kg) 达到38.1%。由于ZEN广泛存在于玉米、小麦等谷物的存贮过程中，对谷物农副产品及饲料进行ZEN脱毒是防治其危害的重要措施之一。然而，传统的基于物理吸附[6]和化学降解[7]的脱毒策略对农副产品及饲料中营养物质都有着不可避免的破坏。物理吸附剂会因为在脱毒的同时对营养物质的吸附导致营养物质的流失，而化学法脱毒会因为化学试剂的引入造成营养物质的二次污染。由于生物脱毒法是在玉米加工和生产的任一环节添加活菌或微生物制备的酶制剂，对毒素进行降解，反应条件相对温和[8]，所以基于活菌降解或酶降解的生物脱毒方式在防治ZEN污染的研究中逐渐得到重视[9]。目前，酵母菌[10]、芽孢杆菌[11]和乳杆菌[12]等都可以实现ZEN毒素的有效降解。此外，来源于Clonostachys rosea的玉米赤霉烯酮水解酶(Zearalenone hydrolase，ZHD) 可以破坏ZEN的内酯环，实现降解[13]。由于活菌降解存在菌株培养周期长、降解活性低等不利因素，因此酶降解法在ZEN生物脱毒方面具有更重要的工业应用价值。2002年，Takahashi-Ando等[14]首次将Clonostachys rosea中提取的碱性水解酶ZHD101测序并在裂殖酵母和大肠杆菌中实现异源表达。2014年，Peng等[15]解析了ZHD101水解酶(PDB ID: 3WZL) 及其与ZEN复合体(PDB ID: 3WZM) 的晶体结构，并确认了影响ZHD101水解ZEN能力的关键残基与底物结合位点，为ZHD101水解能力的进一步改造以及应用打下了基础。然而，ZHD101酶学功能的最适温度在37–45 ℃左右，当温度升高至50 ℃时，ZHD101就会丧失大部分活性[16]。较弱的热稳定性阻碍了ZHD101在工业中的应用。此外，由于ZHD101催化ZEN的反应没有吸光值的变化，难以进行高通量筛选，这也是进一步改造ZHD101获得优良突变体的主要障碍。
近年来，计算机运算能力的持续提升和先进算法的不断涌现推动着蛋白质工程的发展[17-18]，多种计算设计手段成功地实现了酶热稳定性改造[17]。其中，分子动力学(Molecular dynamics，MD) 模拟可以通过提供动态的分子相互作用信息指导有助于结构稳定的蛋白质设计[19]。结合分子动力学模拟，降低柔性位点的波动性是提高蛋白质热稳定性的常见策略[20]。例如，使用刚性较强的脯氨酸取代甘氨酸[21]、或引入二硫键降低柔性区域的波动性[22]以及优化蛋白表面静电相互作用，都可以有效实现蛋白质稳定性的提升[23]。但是在通常情况下，基于MD模拟的分子间作用力分析提供的突变设计数量较少。最近，通过优化能量函数进行序列自动设计逐渐成为蛋白质计算设计的主流策略[24-25]。该设计方法的覆盖范围广且会得到较多潜在的突变体；但是相应的计算量也较为庞大，而且存在“假阳性”设计[26]。将蛋白质计算设计和同源蛋白的多序列比对相结合是一种有效的改善方法。这种方法不仅可以有效地减少突变设计文库，同时还可以避免突变设计对酶的活性位点构型的影响。Meng等[27]针对二聚体亚基界面敏感区域，结合二硫键挖掘、FRESCO策略计算和保守序列分析等多种手段进行突变预测，突变体有效改善了耶式假单胞菌转氨酶的稳定性，同时提高了生产纯胺的产率，其中的正向突变设计比例高达56%。Xiong等[28]建立的ABACUS策略通过统计能量函数进行序列设计，得到的人工蛋白的热稳定性远超天然蛋白，Tm值可达123.3 ℃。Goldenzweig等[26]基于多序列比对的位置特异性信息与Rosetta设计模拟开发了PROSS自动计算策略，并成功设计了带有51个突变的乙酰胆碱酯酶突变体，与野生型相比，该突变体表现出20 ℃的热稳定性增强，突变的引入有效改善了蛋白的核心堆积、表面极性和骨架刚性等，同时提高了重组蛋白的表达水平，降低了酶分子的聚集倾向。
本研究借助分子动力学模拟从ZHD101的柔性区域中选取了32个突变位点，由于ZHD101催化ZEN反应的吸光值没有变化，难以进行高通量筛选，因此选取其作为模式蛋白，进行虚拟饱和突变(共计608个)，得到了12个突变体。其中3个突变体设计可以有效实现热稳定性的提高，并保留甚至提高了对ZEN的水解活性。以ZHD101为模式蛋白，可以探索基于柔性区域的虚拟饱和突变的可行性，并为ZHD101以后的工业应用奠定基础。
1 材料与方法 1.1 材料 1.1.1 菌株与质粒Escherichia coli BL21 (DE3) 菌株与表达载体pET-28a (+) 购自Novagen公司，玉米赤霉烯酮水解酶基因zhd101由苏州金唯智生物技术有限公司合成。
1.1.2 主要试剂PrimeSTAR HS (Premix) 高保真PCR酶，限制性核酸内切酶Hind Ⅲ，NcoⅠ，DL 10 000 DNA marker以及Premixed Protein marker购买自宝日医生物技术有限公司。质粒DNA提取试剂盒以及DNA胶回收试剂盒购自爱思进生物技术有限公司。BCA蛋白浓度测定试剂盒购自碧云天生物技术有限公司。玉米赤霉烯酮购自Sigma- Aldrich公司。甲醇与乙腈为市售高效液相色谱纯，其余试剂均为市售分析纯。
1.2 方法 1.2.1 分子动力学模拟以ZHD101的晶体结构(PDB ID: 3WZL) 作为分子动力学模拟的初始结构，突变体初始构象通过氨基酸替换获得。本研究所有的分子动力学模拟均使用AMBER99SB-ILDN力场在Gromacs 5.0.4中实现[29]。酶分子被置于立方体的水盒子中，到盒子边缘的距离至少为1.0 nm，水分子采用的是TIP3P模型[30]。静电相互作用和范德华相互作用的最小截止距离设置为1.2 nm。模拟的时间步长设置为2 fs。采用steepest descent方法[31]对系统进行能量最小化。随后通过velocity rescaling算法保持体系温度，并且使用Berendsen的弱耦合方法控制压力为0.1 MPa[32]，先后进行50 ns NPT的模型优化和200 ns NVT平衡模拟。我们对野生型(298 K) 和3个突变体进行了重复实验，为了增大采样的随机性，分别在NPT模拟不同时间段提取结构作为初始构象，并设置了不同的随机数种子。轨迹1和轨迹2的随机种子数分别为18 777和20 289 (野生型)，160 866和115 386 (N156F)，291 715和22 903 (S194T)，5 934和24 158 (T259F)。
根据Barlow等[33]的定义，盐桥相互作用由带负电氨基酸(Glu或Asp) 羧酸酯基团中的氧原子与带正电氨基酸(Arg或Lys) 酰胺基团中的氮原子之间的距离确定，当带电氨基酸之间原子间最小距离小于4 Å时带电残基之间形成盐桥相互作用。
在分子动力学模拟过程中，残基的二级结构类型(Definition of secondary structure of proteins，DSSP) 会随作用力以及位置发生变化。其中，较为稳定的残基的二级结构会基本保持该残基初始结构中的类型；而较为灵活或不够稳定的残基的二级结构类型会在模拟过程中不断改变，很难一直维持初始结构。本研究统计了模拟过程中残基维持初始结构中二级结构类型的步数占总步数的比例，即二级结构倾向性(PDSSP，propensity of DSSP)。该信息有助于选取出较为灵活的残基。
1.2.2 构象自由能计算使用Rosetta Cartesian DDG对结构进行优化和构象自由能的计算，以评估野生型与突变体的稳定性变化。目标位点的残基分别突变为其他19种氨基酸，突变体构象自由能变化值(ΔΔG) 越小表示设计的突变体越稳定[34]。使用FastRelax流程在笛卡尔空间中优化野生型蛋白质，Rosetta Cartesian DDG在允许小幅度的主链运动的前提下，确定野生型和每个突变的最佳旋转异构体侧链构象[35]。与传统的自由能计算方法相比，Rosetta Cartesian DDG的计算速度更快[34]。由于在本研究中经过多轮筛选后剩余的突变设计相对较少，同时为了富集正向突变，本研究选取ΔΔG < –0.45 kcal/mol作为筛选标准。
1.2.3 位置特异性评分矩阵利用NBCI数据库搜寻ZHD101的同源序列进行多序列比对(Multiple sequence alignment，MSA)，通过MSA产生位置特异性评分矩阵(Position specific scoring matrix，PSSM)[36]。该矩阵根据残基在MSA中出现的可能性对氨基酸同一性进行加权，对每个位置中20种氨基酸出现的对数概率进行统计，以表示每个位点中所有20种氨基酸各自的对数似然性[26]。出现概率高的氨基酸会有较高的PSSM得分，出现概率较低的氨基酸PSSM得分较低。PSSM比对扫描可以帮助消除有害突变以及自然界中罕见或从未出现过的突变。PSSM得分大于0为有利突变，而PSSM得分小于0则为不利突变。本研究选取PSSM得分大于0作为筛选标准，实现突变为最常见的氨基酸。
1.2.4 定点突变及表达纯化以质粒pET28a-zhd101为模板，根据突变位点设计引物并进行全质粒PCR，将PCR产物进行回收后转化至E. coli BL21 (DE3)感受态中。本研究所使用的引物如表 1所示。
| Primer names | Primer sequences (5′–3′) | Size (bp) |
| N152V-F | CCTGGCCGTTGTTATGCTGAAC | 22 |
| N152V-R | GCTGATTTCCTCGTCTTCCAGC | 22 |
| V153F-F | GCCAACTTTATGCTGAACGAC | 21 |
| V153F-R | CTTGCTGATTTCCTCGTCTTC | 21 |
| M154F-F | CAACGTTTTCCTGAACGACG | 20 |
| M154F-R | GATCTTGCTGATTTCCTCGTC | 21 |
| N156F-F | GTTATGCTGTTTGACGTGAGC | 21 |
| N156F-R | CAGGATCTTGCTGATTTCCTC | 21 |
| D157Q-F | GAACCAGGTGAGCGGTGGTAG | 21 |
| D157Q-R | CATAACGTTGGCCAGGATCTT | 21 |
| Q166F-F | GAAGCCTGGTTTGCAATGGGC | 21 |
| Q166F-R | GCTACCACCGCTCACGTCGTT | 21 |
| S194T-F | CATTCCTCCGACCGCACCGGT | 21 |
| S194T-R | CGGATAACCACGGGCCCACAC | 21 |
| T259F-F | CAAGTATGTGGTGGAATTTACC | 22 |
| T259F-R | CTTGGCGAACACGTCCGGATG | 21 |
| T260V-F | GGTGGAAACCGTTCAGAAAC | 20 |
| T260V-R | CTTGGCGAACACGTCCGGATG | 21 |
| Q261R-F | GTGGAAACCACCCGTAAACAT | 21 |
| Q261R-R | CACATACTTGGCGAACACGTC | 21 |
| K262R-F | GTGGAAACCACCCAGCGTCAT | 21 |
| K262R-R | CCACATACTTGGCGAACACGT | 21 |
| H263Y-F | GTGGAAACCACCCAGAAATAT | 21 |
| H263Y-R | CACCACATACTTGGCGAACAC | 21 |
将E. coli BL21 (DE3)/pET28a-zhd101及其突变菌株平板划线于LB固体培养基，37 ℃恒温培养8 h，挑取单菌落接种至20 mL LB培养基中，37 ℃、200 r/min摇床培养12 h后，按1%接种量转接至100 mL TB培养基，37 ℃、200 r/min摇床培养约2 h至OD600为0.6-0.8后加入诱导剂，25 ℃、200 r/min诱导发酵24 h。
离心收集发酵菌体后重悬菌体，超声破碎20 min后，离心20 min收集破壁上清，用0.22 μm滤膜过滤，进行Ni2+柱亲和层析纯化：先用A液(0.02 mol/L Tris，0.5 mol/L NaCl，0.02 mol/L咪唑)平衡Ni2+亲和层析柱15个柱体积，再将破壁上清以1 mL/min的流速上样，上样结束后，用10%的B液(0.02 mol/L Tris，0.5 mol/L NaCl，0.5 mol/L咪唑)除去杂蛋白，然后进行梯度洗脱(25%和100%的B液) 获得目的蛋白，再用Desalting凝胶柱对目的蛋白进行脱盐处理。最后将目的蛋白稀释至0.2 mg/mL，对分离纯化后产物，用SDS-PAGE电泳进行检验。
1.2.5 圆二色谱实验圆二色谱(Circular dichroism，CD) 实验可测定蛋白二级结构及热熔融温度(Tm)。使用Chirascan圆二色谱仪和光程为0.1 cm的石英比色皿进行CD实验。将纯化脱盐后的酶液浓度均调节至0.2 mg/mL，取200 μL酶液在4 ℃，190– 260 nm波长范围内进行全波长扫描，每步增量为0.5 nm，平均扫描时间为2 ns。取最低吸收值处波长220 nm处作为热变吸光值进行热变实验，热变温度范围4–80 ℃，升温梯度为1 ℃/min，测定220 nm处的椭圆率。热熔融温度Tm通过以下公式估算：
|
其中θ(T)是观测到的椭圆率，θF(T)和θU(T)是对折叠和展开的基线进行线性拟合得到的估计椭圆率。热熔融温度Tm为T，其中F(T)=0.5。
1.2.6 酶活性测定ZHD101活性测定方法由底物减少量来表征。在100 μL的Tris-HCl (0.05 mol/L，pH 7.5)缓冲液中加入20 μL ZEN标准品，混匀后置于37 ℃的反应温度预热10 min，加入40 μL酶液，反应3 min后，加入200 μL纯甲醇终止反应。反应结束后经0.22 μm滤膜过滤后，取200 μL样品用于高效液相色谱检测。样品被60%的乙腈水溶液以1 mL/min的流速洗脱下来，在254 nm处测定吸光值，剩余底物的量在HPLC曲线下根据峰面积计算得出。
此外，将野生型和突变体的酶液分别在50 ℃水浴加热10 min，然后在冰上冷却10 min，按照反应体系进行反应，测定残余活性，以表征高温处理之后酶活性的动力学稳定性。
2 结果与分析 2.1 柔性区域的选取柔性区域和热敏感区域容易触发蛋白质解折叠，可以作为改造热稳定性的靶标。ZHD101分子的柔性区域，主要体现在晶体结构中温度因子(B-factor) 数值较大和蛋白结构动态波动性RMSF (均方根涨落，root mean square fluctuation) 数值较大的区域。借助B-FITTER软件[37]，每个残基标准化后的B-factor可以从ZHD101晶体结构数据中提取并统计(图 1A)[15]。ZHD101结构中表现出较高模糊性的区域主要包括帽子区域(132-187，B-factor可达91.1 Å2)，C-末端区域(250-264，B-factor可达71.3 Å2) 以及连接帽子区域和核心区的α9-α10螺旋(191-203，B-factor可达68.1 Å2)。同时，我们对ZHD101晶体结构在298 K温度下进行了200 ns的分子动力学模拟。结构的波动特征通过统计模拟过程中残基的位置偏移，即均方根涨落来反映，RMSF越大表明该区域的构象波动性越大。在ZHD101结构中具有较高波动性的区域位于帽子区域和α9-α10螺旋。RMSF在预测ZHD101的柔性区域方面表现出与B-factor分析较高的一致性。
|
| 图 1 ZHD101中柔性和热敏感区域的选取 Fig. 1 Identification of flexible and heat sensitive regions in ZHD101. (A) B-factor and RMSF (298 K) of wild type. (B) Shift of RMSF and DSSP under 298 K and 323 K, ΔRMSF=RMSF323K-RMSF298K, ΔPDSSP=PDSSP 323K-PDSSP 298K. The shaded regions are α-helix structures. (C) Mutation sites are shown in the structure of the wild type. Red spheres indicate flexible sites meet the ΔRMSF and ΔPDSSP. Yellow spheres indicate mutational sites meet the filter, PSSM score and ΔΔG calculation. |
| |
高温下的分子动力学模拟可以帮助选取蛋白质解折叠的热敏感区域，对ZHD101晶体结构在323 K温度下进行了200 ns的分子动力学模拟。比对298 K和323 K温度下ZHD101构象动力学特征，帽子区域的RMSF以及二级结构倾向性(PDSSP) 发生显著的变化(图 1B)，ΔRMSF最高可达3.3 Å，ΔPDSSP最低可达–55.8%。虽然C-末端区域在298K的波动性较弱(RMSF298K最高为1.0 Å)，但是该区域对温度表现出较高的敏感性。当温度升高时，该区域残基的RMSF有明显上升(ΔRMSF最高为1.3 Å)，同时PDSSP下降明显(ΔPDSSP最低可达–20.5%)。
当温度升高时，ZHD101结构中热敏感区域与柔性区域(如图 1B中红色位点所示) 高度重叠。在柔性区域内，ΔRMSF > 0.5 Å和ΔPDSSP < –10%作为筛选标准筛选出32个同时表现出较高柔性和热敏感性的位点(图 1C)。其中有19个位点位于帽子区域(L141，D143，E144，A151，N152，V153，M154，L155，N156，D157，V158，S159，E163，A164，W165，Q166，A167，H177，K178，Y187)，5个位点分布于C末端的α-螺旋区域(T259，T260，Q261，K262，H263)，此外还有8个位点分布于帽子区域下方的loop2 (I13)、loop8 (T73，E74)、loop14 (R189，T190)、α9螺旋(S194) 以及α10螺旋(D199)。
2.3 突变体的计算筛选在选取的32个突变位点上，依据突变所引起的构象自由能变化和同源序列比对中获得的氨基酸保守性进行了虚拟饱和突变(共608个)。同时满足PSSM > 0和ΔΔG < 0 kcal/mol的突变设计共有22个，ΔΔG的平均值为(–1.63±1.20) kcal/mol。为了富集正向突变，本研究选取了ΔΔG < –0.45 kcal/mol的12个突变设计进行实验验证。
结合自由能计算(ΔΔG < –0.45 kcal/mol) 和序列保守性(PSSM > 0) 的筛选，在32个候选位点的608个突变中共包含12个潜在的热稳定性突变设计(图 2，表 2)，包括N152V，V153F，M154F，N156F，D157Q (帽子区域)，T259F，T260V，Q261R，K262R，H263Y (C-末端区域)，S194T (α9螺旋) (图 1C，橘色球状位点)。其中，帽子区域和C-末端均有自由能降低较多(ΔΔG < –5.5 kcal/mol) 和PSSM得分较高(PSSM≥4) 的突变。在帽子区域中，V153F和M154F突变自由能降低最为明显(ΔΔG < –5.5 kcal/mol)，V153F和D157Q有较高的PSSM得分，而N152V和Q166F两个突变的自由能仅略微降低(–1 < ΔΔG < 0 kcal/mol)。C-末端中，T259F和H263Y的自由能变化远高于其他位点(T260V，Q261R和K262R) 的自由能降低幅度(ΔΔG < –4.4 kcal/mol)，同时T259F有着较高的PSSM得分。在α9和α10螺旋区域，只有S194T符合筛选的条件(PSSM=1，ΔΔG=–1.781 kcal/mol)。
|
| 图 2 突变位点的计算设计 Fig. 2 Computational design of mutant sites. (A) Stability prediction of candidate sites (ΔΔG) by Rosetta Cartesian DDG. (B) PSSM scores of candidate sites. Black squares highlight the mutants meet both two filters. |
| |
| Mutants | ΔRMSF (Å) | ΔPDSSP (%) | PSSM | ΔΔG (kcal/mol) | Tm (℃) | ΔTm (℃) | Activity (%) | |
| Relative | Residual | |||||||
| WT | – | – | – | – | 43.0 | – | 100 | 37.1 |
| N152V | 0.892 | –30.8 | 1 | –0.568 | 49.8 | 6.8 | 79.1 | 56.9 |
| V153F | 1.014 | –26.5 | 6 | –5.616 | 48.8 | 5.8 | 72.5 | 70.0 |
| M154F | 1.276 | –11.3 | 2 | –5.836 | 47.5 | 4.5 | 0.3 | 0.1 |
| N156F | 1.599 | –20.9 | 1 | –2.441 | 47.8 | 4.8 | 95.8 | 87.8 |
| D157Q | 0.815 | –31.0 | 5 | –1.524 | 45.2 | 2.2 | 46.8 | 53.2 |
| Q166F | 1.753 | –30.1 | 1 | –0.800 | 44.6 | 1.6 | 150.3 | 9.2 |
| S194T | 1.122 | –13.4 | 1 | –1.781 | 47.2 | 4.2 | 131.6 | 48.3 |
| T259F | 0.697 | –14.8 | 6 | –4.458 | 47.6 | 4.6 | 169.0 | 81.4 |
| T260V | 0.795 | –14.1 | 1 | –1.362 | 47.5 | 4.5 | 97.6 | 9.8 |
| Q261R | 0.994 | –18.3 | 2 | –1.946 | 43.1 | 0.1 | 98.6 | 30.1 |
| K262R | 1.103 | –15.6 | 5 | –1.671 | 44.8 | 1.8 | 110.9 | 79.2 |
| H263Y | 0.999 | –20.5 | 4 | –5.557 | 46.0 | 3.0 | 108.1 | 24.6 |
突变体使用表 1中的引物进行全质粒PCR所构建，突变后的质粒转化至E. coli BL21 (DE3) 后进行全基因测序，结果正确后进行表达纯化。将野生型玉米赤霉烯酮水解酶ZHD101表达纯化后，用Tris-HCl缓冲液稀释至浓度为0.2 mg/mL，取20 μL酶液上样至SDS-PAGE进行分子量分析。SDS-PAGE中目的蛋白条带与29.0 kDa的Marker条带齐平(图 3，泳道1)，与ZHD101理论分子量(28.75 kDa) 相符合。所有的突变体均能在大肠杆菌中可溶性表达，且条带对应的分子量大小正确，纯度在90%以上(图 3，泳道2-17)。
|
| 图 3 野生型ZHD101及其突变体在大肠杆菌中的表达纯化 Fig. 3 Expression and purification of wild type ZHD and its mutants expressed in E. coli BL21 (DE3). M: premixed protein marker; 1: purified wild type ZHD; 2-17: purified protein of mutant N152V, V153F, M154F, N156F, D157Q, Q166F, S194T, T259F, T260V, Q261R, K262R, H263Y, N156F/S194T, N156F/T259F, S194T/T259F, N156F/S194T/T259F. |
| |
野生型和突变体在4 ℃进行190-260 nm的全波长扫描(图 4A)。野生型与突变体在220 nm处均存在最低吸收值，所有蛋白全部折叠形成了正确的二级结构。在4-80 ℃的范围内测定蛋白质热熔融过程中椭圆率的变化情况。蛋白的热熔融温度(Tm) 由220 ns波长下的热变曲线获得(图 4B)。
|
| 图 4 野生型及突变体的热力学稳定性表征 Fig. 4 Far-UV CD spectra and thermal denaturation of wild type and mutants. (A) CD spectra at 4 ℃ in 0.02 mol/L Tris-HCl buffer (pH 7). (B) Thermal melting curves monitored at 220 nm in Tris-HCl buffer. The samples were shown in different colors as labeled in the plots. |
| |
同时，本研究测定了野生型与突变体的相对酶活性和在50 ℃热处理10 min后的残余酶活性(表 2) 以表征突变设计对酶活性的影响。
位于帽子区域的6个突变体有4个热稳定性有一定程度的提升，包括N152V (ΔTm=6.8 ℃)、V153F (ΔTm=5.8 ℃)、M154F (ΔTm=4.5 ℃)和N156F (ΔTm=4.8 ℃)，其中只有N156F突变保持了酶活性(相对酶活性为95.8%，以野生型未热处理的酶活性为相对标准)，并且在加热处理后残余酶活性(87.8%) 也得到了明显改善；N152V和V153F的相对酶活性(< 80%) 受到一定程度的损失；而M154F突变导致酶完全失活(表 2)。这是由于帽子区域与底物结合有关，突变可能对酶的催化过程有一定影响，导致了某些突变体相对酶活性的损失。
位于C-末端的5个突变体中，T259F和T260V热稳定性的提高相对较多(ΔTm > 4.0 ℃)。由于位于C-末端，远离催化中心，对酶活性基本没有影响，并且T259F的相对酶活性(169.0%) 有一定的提高。
此外，S194T作为α9螺旋中唯一存在的位点，热稳定性有所提升(ΔTm=4.2 ℃)，且表现出相对酶活性与残余酶活性均高于野生型(131.6%和48.3%)。
2.6 突变体的分子动力学模拟对于热稳定性有所提升(ΔTm > 4.0 ℃)，并且酶活性与野生型相似甚至高于野生型(相对酶活性 > 95%，残余酶活性 > 40%) 的3个突变体N156F、S194T和T259F进行了200 ns的动力学模拟。从主链原子的均方根偏差(Root mean square deviation，RMSD) 的变化趋势可以看出，野生型及突变体的构象在60 ns后达到平衡状态(图 5A)。在轨迹1中，达到平衡后3个突变体的RMSD值(相对结构为NVT模拟的初始结构) 均低于野生型，与野生型相比(2.28±0.12 Å)，N156F、S194T和T259F在后100 ns轨迹内的平均RMSD分别下降至(1.85±0.14) Å、(1.82±0.10) Å和(1.69±0.11) Å，轨迹2的数值与轨迹1相似，相较于野生型，突变引入减弱了构象波动性，更有利于ZHD101结构稳定性的提高。此外，对野生型及3个突变体后100 ns的轨迹进行RMSF分析可以发现，相对于野生型帽子区域较高的波动性(图 5B)，3个突变体位于帽子区域残基的RMSF都有明显降低。
|
| 图 5 野生型ZHD101 (黑色) 及突变体N156F (红色)、S194T (蓝色) 和T259F (橙色) 的分子动力学模拟分析(A：均方根偏差；B：均方根波动. 黑框标出的为帽子区域) Fig. 5 Molecular dynamics simulation analysis for wild type ZHD101 (black) and mutant N156F (red), S194T (blue) and T259F (orange). (A) The root mean square deviation (RMSD). (B) The root mean square fluctuation (RMSF). Black squares highlight the cap domain. |
| |
S194T突变的引入并没有实现局部柔性区域的刚性化。相反，该突变导致附近残基的二级结构倾向性均有一定程度的降低，在构象上可能表现出更大的自由度。这种区域的结构特征变化导致位于α9螺旋的R189与位于帽子区域的E142之间的盐桥相互作用被破坏(图 6A)。该区域强相互作用力的丢失有一定概率引起突变位点附近区域以及帽子区域的构象以及盐桥作用力的调整(图 6B)，其中R189与E74以及E142与R185之间的盐桥成键概率分别增加了13.1% (轨迹2：24.2%) 和33.0% (轨迹2：52.9%)。由S194T突变引起的区域盐桥相互作用力重排，对构象稳定性的提高可能有所帮助。此外，在野生型ZHD101中，帽子区域的残基与其他区域的多个残基之间存在显著的动力学相关性。当S194T突变后，由于帽子区域的E142残基与本体R189残基之间的相互作用的丢失切断了α9螺旋与帽子区域的联系，我们推测区域内的构象变化与作用力的重排进一步减弱了区域间残基之间的动力学相关性。
|
| 图 6 N156F、S194T和T259F热稳定性提高的机理分析 Fig. 6 Mechanism analysis of improved thermal stability of N156F, S194T and T259F. (A) Local salt bridge changes caused by S194T mutation. (B) Comparison of the distance between the oxygen atom in the side chain of E142 and the nitrogen atom in the side chain of R185 in wild type and S194T. (C) Atypical non-covalent effects introduced by the N165F mutation. (D) The comparison of the distance between the oxygen atom of N156 side chain and the centroid of benzene ring of F156 side chain and the nitrogen atom of N152 side chain. (E) The T259F mutation fills the cavities on the surface of the protein, and the black and red squares highlight the cavities. (F) The comparison of the residue solvent accessible surface area (SASA) between WT and the T259F mutant. |
| |
值得注意的是，N156F突变设计位于结构的表面，在通常情况下并不利于蛋白质的折叠。由此推断，N156F可能与邻近的残基形成非典型性的非共价相互作用，比如与L155形成CH-π相互作用[38]，或者与N152形成NH-π相互作用(图 6C)[39]。这些非共价相互作用在蛋白结构中常有发生，并对蛋白的热稳定性具有一定的贡献[40]。从图 6D中可以看出，N156与N152侧链之间的平均距离为9.1 ű2.7 Å (轨迹2：8.4 ű2.1 Å)。而N156突变为苯丙氨酸后，其侧链苯环质心与N152侧链氮原子的平均距离仅为3.4 ű0.7Å (轨迹2：3.7 ű0.3 Å)，可形成NH-π相互作用。
T259位于蛋白质表面的空穴中，突变位点由苏氨酸突变成苯丙氨酸后，残基侧链长度和疏水性的增加可以将侧链埋藏于蛋白内部，并填充表面空穴，从而提高蛋白质的热稳定性(图 6E)[41]。将位于表面的疏水残基的侧链隐藏在蛋白内部已经被证明可以提高天然蛋白的稳定性[42]。从图 6F中可以看出，与野生型相比(11, 339.9±216.8 Å2)，T259F的溶剂可及表面积(SASA，solvent accessible surface area) 有所减小(轨迹1：11 068.6 ű172.9 Å2，轨迹2：11 104.6 ű218.5 Å2)，这对热稳定性的提高可能具有一定的贡献。
2.8 组合突变体的性质表征为了获得性质较好的突变体，并考察突变位点的上位效应，对N156F、S194T和T259F进行迭代组合突变，分别为N156F/S194T、N156F/T259F、S194T/T259F、N156F/S194T/T259F。全波长扫描结果显示，这4个组合突变体折叠形成了正确的二级结构，在220 nm处存在最低吸收峰(图 7A)。3个组合突变体热的熔融温度(Tm) 由220 ns波长下的热变曲线获得(图 7B)。其中热稳定性提高较多的为N156F/S194T (ΔTm=6.7 ℃)，所有组合突变体的Tm提高效果都弱于单点突变的Tm叠加效果，这说明了这3个突变之间可能存在负上位效应。与野生型相比，4个组合突变的相对酶活性(> 100%) 和残余酶活性均有所提升(> 50%) (表 3)。
|
| 图 7 野生型及组合突变体的热力学稳定性表征 Fig. 7 Thermal denaturation of wild type and combined mutants. (A) CD spectra at 4 ℃ in 0.02 mol/L Tris-HCl buffer (pH 7). (B) Thermal melting curves monitored at 220 nm in Tris-HCl buffer. The samples were shown in different colors as labeled in the plots. |
| |
| Mutants | Tm (℃) | ΔTm (℃) | Activity (%) | |
| Relative | Residual | |||
| WT | 43.0 | – | 100 | 37.1 |
| N156F/S194T | 49.7 | 6.7 | 103.6 | 89.7 |
| N156F/T259F | 46.7 | 3.7 | 150.1 | 60.3 |
| S194T/T259F | 46.3 | 3.3 | 171.4 | 54.6 |
| N156F/S194T/T259F | 49.1 | 6.1 | 152.9 | 80.2 |
这3个位点相互之间距离较远，稳定蛋白质的机理各不相同，且组合后表现出不同效果。N156F和S194T热稳定性的提高是可叠加的(ΔTm≈7 ℃)，T259F与N156F，S194T形成的双点突变的Tm与单点突变类似(ΔTm≈4 ℃)，T259F/N156F/S194T的三点突变与双点突变类似(ΔTm≈7 ℃)，似乎T259F对稳定性的增强作用被抑制了。N156F和S194T是通过作用力改变提升热稳定性，T259F通过表面填充改变蛋白质疏水特性提升热稳定性。T259F位于C-末端区域，可能该区域构象多变并且对蛋白质整体构象依赖性高，这可能是在组合突变中T259F稳定性被抑制的一个原因。
3 结论本研究通过柔性和热敏感区域的分析，确定了3个区域的32个突变位点。因此，对于每个位点进行了虚拟饱和突变，共计608个突变体。依据构象自由能变化和多序列比对保守性筛选出了12个合理的突变体。在12个突变体中，其中7个热稳定性有一定程度的提升(ΔTm > 4 ℃)。在热稳定性提高显著的突变体中，N156F、S194T和T259F这3个突变体保留了野生型的酶活性，甚至有一定幅度的提高(相对酶活性分别为95.8%、131.6%和169.0%)。
S194T、N156F和T259F可能分别通过盐桥重排，形成NH-π和填补蛋白表面空腔提高热稳定性。通常情况下，热稳定性的提高以酶活性作为补偿。但是这3个突变体保留甚至提高了酶活性，它们对催化作用的影响机理可能比较复杂。以S194T为例，S194T突变与催化中心距离较远，但是该突变可能会导致附近残基在构象上发生变化，区域盐桥相互作用力重排，可能影响了催化中心与ZEN的结合，对酶活性产生了影响，这种别构效应在功能蛋白中普遍存在，比如在酪氨酸磷脂酶蛋白中被观察到[43]。因此，选取ZHD101作为模式蛋白，通过分子动力学模拟结合虚拟饱和突变的策略，可有效提高ZHD101的热稳定性及其在工业应用中的潜力。
| [1] |
Zinedine A, Soriano JM, Moltó JC, et al. Review on the toxicity, occurrence, metabolism, detoxification, regulations and intake of Zearalenone: an oestrogenic mycotoxin. Food Chem Toxicol, 2007, 45(1): 1-18. DOI:10.1016/j.fct.2006.07.030
|
| [2] |
D'Mello JPF, Placinta CM, MacDonald AMC. Fusarium mycotoxins: a review of global implications for animal health, welfare and productivity. Animal Feed Sci Technol, 1999, 80(3/4): 183-205.
|
| [3] |
Berek L, Petri IB, Mesterházy A, et al. Effects of mycotoxins on human immune functions in vitro. Toxicol Vitro, 2001, 15(1): 25-30. DOI:10.1016/S0887-2333(00)00055-2
|
| [4] |
程传民, 柏凡, 李云, 等. 2013年玉米赤霉烯酮在饲料原料中的污染分布规律. 中国畜牧杂志, 2014, 50(16): 68-72, 77. Cheng CM, Bai F, Li Y, et al. Distribution of zearalenone contamination of feed raw materials in 2013. Chin J Animal Sci, 2014, 50(16): 68-72, 77 (in Chinese). DOI:10.3969/j.issn.0258-7033.2014.16.014 |
| [5] |
国家标准化管理委员会, 国家质量监督检验检疫总局. 饲料卫生标准: GB 13078-2017. Standardization Administration, State Administration for Market Regulation. Hygienical standard for feeds: GB 13078-2017 (in Chinese). |
| [6] |
Bueno DJ, Di Marco L, Oliver G, et al. In vitro binding of zearalenone to different adsorbents. J Food Prot, 2005, 68(3): 613-615. DOI:10.4315/0362-028X-68.3.613
|
| [7] |
Abd Alla ES. Zearalenone: incidence, toxigenic fungi and chemical decontamination in Egyptian cereals. Nahrung, 1997, 41(6): 362-365. DOI:10.1002/food.19970410610
|
| [8] |
俞建良, 苏会波, 李凡, 等. 玉米干法酒精生产过程中真菌毒素降解技术的探讨. 酿酒科技, 2016(9): 59-64. Yu JL, Su HB, Li F, et al. Mycotoxin degradation technology in ethanol production by dry milling corn. Liquor-Mak Sci Technol, 2016(9): 59-64 (in Chinese). |
| [9] |
赵仁勇, 宋斌. 玉米赤霉烯酮生物脱毒研究进展. 河南工业大学学报(自然科学版), 2018, 39(2): 113-121. Zhao RY, Song B. Advances in biological detoxification of zearalenone. J Henan Univ Technol (Nat Sci Ed), 2018, 39(2): 113-121 (in Chinese). DOI:10.3969/j.issn.1673-2383.2018.02.020 |
| [10] |
Vekiru E, Hametner C, Mitterbauer R, et al. Cleavage of zearalenone by Trichosporon mycotoxinivorans to a novel nonestrogenic metabolite. Appl Environ Microbiol, 2010, 76(7): 2353-2359. DOI:10.1128/AEM.01438-09
|
| [11] |
Lei YP, Zhao LH, Ma QG, et al. Degradation of zearalenone in swine feed and feed ingredients by Bacillus subtilis ANSB01G. World Mycotoxin J, 2014, 7(2): 143-151. DOI:10.3920/WMJ2013.1623
|
| [12] |
El-Nezami H, Polychronaki N, Lee YK, et al. Chemical moieties and interactions involved in the binding of zearalenone to the surface of Lactobacillus rhamnosus strains GG. J Agric Food Chem, 2004, 52(14): 4577-4581. DOI:10.1021/jf049924m
|
| [13] |
Kakeya H, Takahashi-Ando N, Kimura M, et al. Biotransformation of the mycotoxin, zearalenone, to a non-estrogenic compound by a fungal strain of Clonostachys sp.. Biosci Biotechnol Biochem, 2002, 66(12): 2723-2726. DOI:10.1271/bbb.66.2723
|
| [14] |
Takahashi-Ando N, Kimura M, Kakeya H, et al. A novel lactonohydrolase responsible for the detoxification of zearalenone: enzyme purification and gene cloning. Biochem J, 2002, 365(pt 1): 1-6.
|
| [15] |
Peng W, Ko TP, Yang YY, et al. Crystal structure and substrate-binding mode of the mycoestrogen-detoxifying lactonase ZHD from Clonostachys rosea. RSC Adv, 2014, 4(107): 62321-62325. DOI:10.1039/C4RA12111B
|
| [16] |
Takahashi-Ando N, Ohsato S, Shibata T, et al. Metabolism of zearalenone by genetically modified organisms expressing the detoxification gene from Clonostachys rosea. Appl Environ Microbiol, 2004, 70(6): 3239-3245. DOI:10.1128/AEM.70.6.3239-3245.2004
|
| [17] |
曲戈, 朱彤, 蒋迎迎, 等. 蛋白质工程: 从定向进化到计算设计. 生物工程学报, 2019, 35(10): 1843-1856. Qu G, Zhu T, Jiang YY, et al. Protein engineering: from directed evolution to computational design. Chin J Biotech, 2019, 35(10): 1843-1856 (in Chinese). |
| [18] |
田健, 王平, 伍宁丰, 等. 理性设计提高蛋白质热稳定性的研究进展. 生物技术进展, 2012, 2(4): 233-239. Tian J, Wang P, Wu NF, et al. Recent advances in the rational design to improve the protein thermostability. Curr Biotechnol, 2012, 2(4): 233-239 (in Chinese). DOI:10.3969/j.issn.2095-2341.2012.04.01 |
| [19] |
Childers MC, Daggett V. Insights from molecular dynamics simulations for computational protein design. Mol Syst Des Eng, 2017, 2(1): 9-33. DOI:10.1039/C6ME00083E
|
| [20] |
Yu H, Huang H. Engineering proteins for thermostability through rigidifying flexible sites. Biotechnol Adv, 2014, 32(2): 308-315. DOI:10.1016/j.biotechadv.2013.10.012
|
| [21] |
Tian J, Wang P, Gao S, et al. Enhanced thermostability of methyl parathion hydrolase from Ochrobactrum sp. M231 by rational engineering of a glycine to proline mutation. FEBS J, 2010, 277(23): 4901-4908. DOI:10.1111/j.1742-4658.2010.07895.x
|
| [22] |
Pikkemaat MG, Linssen ABM, Berendsen HJC, et al. Molecular dynamics simulations as a tool for improving protein stability. Protein Eng Des Sel, 2002, 15(3): 185-192. DOI:10.1093/protein/15.3.185
|
| [23] |
Strickler SS, Gribenko AV, Gribenko AV, et al. Protein stability and surface electrostatics: a charged relationship. Biochemistry, 2006, 45(9): 2761-2766. DOI:10.1021/bi0600143
|
| [24] |
Rohl CA, Strauss CE, Misura KM, et al. Protein structure prediction using Rosetta. Methods Enzymol, 2004, 383: 66-93.
|
| [25] |
Schymkowitz J, Borg J, Stricher F, et al. The FoldX web server: an online force field. Nucleic Acids Res, 2005, 33(web server issue): W382-W388.
|
| [26] |
Goldenzweig A, Goldsmith M, Hill SE, et al. Automated structure-and sequence-based design of proteins for high bacterial expression and stability. Mol Cell, 2016, 63(2): 337-346. DOI:10.1016/j.molcel.2016.06.012
|
| [27] |
Meng QL, Capra N, Palacio CM, et al. Robust ω-transaminases by computational stabilization of the subunit interface. ACS Catal, 2020, 10(5): 2915-2928. DOI:10.1021/acscatal.9b05223
|
| [28] |
Xiong P, Wang M, Zhou X, et al. Protein design with a comprehensive statistical energy function and boosted by experimental selection for foldability. Nat Commun, 2014, 5: 5330. DOI:10.1038/ncomms6330
|
| [29] |
Lindorff-Larsen K, Piana S, Palmo K, et al. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins, 2010, 78(8): 1950-1958. DOI:10.1002/prot.22711
|
| [30] |
Jorgensen WL, Chandrasekhar J, Madura JD, et al. Comparison of simple potential functions for simulating liquid water. J Chem Phys, 1983, 79(2): 926-935. DOI:10.1063/1.445869
|
| [31] |
Bussi G, Donadio D, Parrinello M. Canonical sampling through velocity rescaling. J Chem Phys, 2007, 126(1): 014101. DOI:10.1063/1.2408420
|
| [32] |
Berendsen HJC, Postma JPM, Van Gunsteren WF, et al. Molecular dynamics with coupling to an external bath. J Chem Phys, 1984, 81(8): 3684-3690. DOI:10.1063/1.448118
|
| [33] |
Barlow DJ, Thornton JM. Ion-pairs in proteins. J Mol Biol, 1983, 168(4): 867-885. DOI:10.1016/S0022-2836(83)80079-5
|
| [34] |
Park H, Bradley P, Greisen P, et al. Simultaneous optimization of biomolecular energy functions on features from small molecules and macromolecules. J Chem Theory Comput, 2016, 12(12): 6201-6212. DOI:10.1021/acs.jctc.6b00819
|
| [35] |
Kellogg EH, Leaver-Fay A, Baker D. Role of conformational sampling in computing mutation-induced changes in protein structure and stability. Proteins, 2011, 79(3): 830-838. DOI:10.1002/prot.22921
|
| [36] |
Altschul SF, Gertz EM, Agarwala R, et al. PSI-BLAST pseudocounts and the minimum description length principle. Nucleic Acids Res, 2009, 37(3): 815-824. DOI:10.1093/nar/gkn981
|
| [37] |
Reetz MT, Carballeira JD, Vogel A. Iterative saturation mutagenesis on the basis of B factors as a strategy for increasing protein thermostability. Angew Chem Int Ed Engl, 2006, 45(46): 7745-7751. DOI:10.1002/anie.200602795
|
| [38] |
Scheiner S, Kar T, Gu Y. Strength of the Calpha H..O hydrogen bond of amino acid residues. J Biol Chem, 2001, 276(13): 9832-9837. DOI:10.1074/jbc.M010770200
|
| [39] |
Steiner T, Koellner G. Hydrogen bonds with pi-acceptors in proteins: frequencies and role in stabilizing local 3D structures. J Mol Biol, 2001, 305(3): 535-557. DOI:10.1006/jmbi.2000.4301
|
| [40] |
Hughes RM, Waters ML. Influence of N-methylation on a cation-π interaction produces a remarkably stable β-hairpin peptide. J Am Chem Soc, 2005, 127(18): 6518-6519. DOI:10.1021/ja0507259
|
| [41] |
Borgo B, Havranek JJ. Automated selection of stabilizing mutations in designed and natural proteins. Proc Natl Acad Sci USA, 2012, 109(5): 1494-1499. DOI:10.1073/pnas.1115172109
|
| [42] |
Munson M, Balasubramanian S, Fleming KG, et al. What makes a protein a protein? Hydrophobic core designs that specify stability and structural properties. Protein Sci, 1996, 5(8): 1584-1593. DOI:10.1002/pro.5560050813
|
| [43] |
Lu C, Knecht V, Stock G. Long-range conformational response of a PDZ domain to ligand binding and release: a molecular dynamics study. J Chem Theory Comput, 2016, 12(2): 870-878. DOI:10.1021/acs.jctc.5b01009
|