靶點?
RNA 聚合酶II
計算方法
QM/MM MD
計算軟件
QChem/Amber
作者與單位
香港理工大學黃旭輝教授課題組（http://compbio.ust.hk）
發表雜志
Nature?Catalysis
前言
酶催化機制的理論研究一直以來是計算化學、計算生物學的重要方向之一，其中QM/MM MD是當前最為廣泛應用且較為認可的計算模擬方法。而計算與實驗相結合，往往能更加深入探究復雜的酶催化機制。在這分享香港理工大學黃旭輝教授在2019 Nature Catalysis上一篇關于RNA聚合酶II 催化內裂過程的理論模擬與實驗驗證相結合的文章。[1]
背景
在基因表達中，為了保證較高的轉錄保真度，RNA聚合酶(RANPs)能夠通過其內在的核酸酶活性來校對新生RNA的轉錄，這是一種水解裂解反應，可以進一步受到其他轉錄因子的刺激。[2]而在一些物種中，轉錄因子的刺激并非必需，如釀酒酵母RNA 聚合酶II (Pol II)?（Fig.1a）的內在裂解活性已經足以維持細胞的活力。
RNA聚合酶利用相同的活性位點進行核苷酸聚合和內在的裂解反應(下簡稱內裂)。在單個核苷酸發生錯配時，RNA聚合酶采取單-核苷酸-回溯態(one-nucleotide-backtracked state)，在無轉錄因子參與下通過水解磷酸二酯鍵發生內裂并移除3’-端的核苷酸（Fig.1b）。研究結果表明，RNA聚合酶內裂采取了與其他酶中（如多種DNA聚合酶）類似的雙金屬機制。然而對于RNA聚合酶來說，核心問題仍然是確定這個裂解反應所涉及的廣義酸堿(general acid/base)問題。
Fig. 1 ?a, QM/MM-MD模型(PDB:3GTG)及QM區域劃分。b, 反應機制示意圖, 步驟1（左）:核苷酸內切斷裂，以水為親核試劑，端位磷酸O作為廣義堿；步驟2（右）：水介導的質子轉移過程。
圖片來源Nature Catalysis
?
雖然時間分辨x射線實驗研究已經成為實時監測酶反應過程的一項強有力的技術，[3]但由于晶體結構分辨率的限制，氫原子無法準確定位，仍然無法完整揭示RNA 聚合酶II (Pol II)中的反應機制。尤其是對于參與激活親核進攻的廣義堿，以及參與質子化O離子離去基團的廣義酸的分析研究。之前的研究推測末端核苷酸的堿基或參與進來的氨基酸殘基（如His等）可能作為廣義堿，[4]而對于廣義酸則仍不清楚。考慮到實驗研究的局限，計算研究可以作為強大的與實驗互補的方法，來揭示反應機制中的細節。[5]
?于是研究人員采用從頭算量子力學/分子力學的分子動力學模擬（ab initio quantum mechanics/molecular mechanics– molecular dynamics, aiQM/MM–MD）方法，解析Pol II 回溯過程中的內裂機制。研究結果表明，端部錯配核苷酸上的磷酸O原子出人意料的充當廣義堿，并協助進攻水（attacking water）的去質子化。而暴露的3 ‘ –端氧陰離子在無廣義酸的情況下，通過水介導的質子轉移進行質子化。上述理論計算預測得到了生化實驗的驗證。
?
研究內容
RNA 聚合酶II (Pol II)內裂的協同反應機制
研究人員基于回溯態的Pol II體系進行了aiQM/ MM-MD模擬，體系包含103 QM區原子和372091 MM區原子（Fig. 1），使用方法基組為B3LYP/6-31G*。在典型的aiQM/ MM-MD模擬中，一般采用球邊界方法將中心區域封閉在截斷的球體中(如Qchem-Tinker)。而該研究采用了新近開發的Qchem-Amber考慮了周期性邊界條件，可以更準確地描述遠距離靜電相互作用。[6]由于體系中含有兩個金屬Mg2+離子，在QM/ MM-MD模擬前首先仔細驗證了初始結構，確認催化所需的重要結構保持得很好，包括兩個Mg2+離子之間的距離以及它們的配位。
?
通過MD模擬得到合理的復合物初始模型后，沿反應坐標進行QM/MM計算，先后獲得一維(1D)和二維(2D)勢能分布圖，然后以傘形采樣(Umbrella Sampling)的方式基于QM/ MM MD模擬進行充分的采樣，沿二維勢能面上最小勢能路徑進行采樣，共收集了115個窗口（每個窗口收集7ps軌跡），利用加權直方圖分析方法(weighted histogram analysis method，WHAM)計算平均力勢(PMF)即獲得反應自由能曲線。同時通過將每個窗口的QM/ MM MD模擬時長加倍來驗證QM/MM MD自由能模擬是否收斂。
結果表明，回溯態Pol II中錯配核苷酸的磷酸二酯鍵斷裂過程采取協同機制。研究人員通過系統的調查反應中心附近可作為潛在廣義堿來激活親核試劑的基團/原子，發現并不是Pol II中的氨基酸殘基而是下游的磷酸O (G12 pro-Rp)作為廣義堿來激活MgB-結合的親核進攻試劑(活化能壘~29 kcal mol?1, 見Fig.2)。反應過程所歷經的態分別為反應態RS(自由能0 kcal mol?1)，過渡態TS1(29.18 ± 0.58 kcal mol?1)，中間體態IS(20.43 ± 0.99 kcal mol?1)，過渡態TS2(28.66 ± 0.31 kcal mol?1)和產物態PS(0.77 ± 0.91 kcal mol?1)。
?
Fig. 2 ?a,一維自由能曲線。b, 步驟1（左）和步驟2（右）的二維反應坐標及自由能面。
圖片來源Nature Catalysis
?
QM/ MM-MD模擬表明，內部裂解的協同機制與過渡態(TS1)有關，該過渡態包含與親核進攻試劑和離去基的部分鍵合(Fig. 3)。這會導致亞穩態IS出現，其中包含一個暴露的A10 (左側核苷酸)3 ‘-端，并被MgA的配位水進一步質子化。
?
Fig. 3 反應途徑中各個態的結構及反應歷程示意圖。
圖片來源Nature Catalysis
?
反應的最后一步是完成質子轉移過程，質子轉移將使活性部位的水分子恢復中性。為了驗證TS的結構，將TS1和TS2進行弛豫(注：TS結構很不穩定，傘型采樣中通過外加限制力才可對不穩定的TS進行采樣，此處弛豫是將TS的限制放開，進行QM/MM MD采樣。理論上采樣可以夠捕捉到由TS向兩側局部能量最小態的運動，從而可輔助驗證TS和反應途徑的合理性)，弛豫途徑與所確定的機制完全一致:即進攻水被末端磷酸O去質子化。同時還檢測了IS的穩定性，其結構包含一個質子化的末端磷酸和一個裂解的磷酸二酯鍵，這進一步支持了3’-端核苷酸磷酸O接受質子并作為廣義堿的能力(Fig. 3)。研究所提出的下游磷酸(G12 pro-Rp) O原子作為內裂過程的廣義堿?，得到了釀酒酵母Pol II和嗜熱菌RNAP的單-核苷酸-回溯態晶體結構的支持。這二者結構中都包含扭曲骨架結構（此處指核酸骨架或核酸主鏈），位于回溯核苷酸和上游核苷酸之間，使前者的pro-Rp O原子能夠靠近活性位點(Fig. 4)。
?
Fig. 4 ?釀酒酵母Pol II和嗜熱菌RNAP的單-核苷酸-回溯態晶體結構。二者都包含扭曲的骨架結構使pro-Rp 原子靠近活性位點。
圖片來源Nature Catalysis
?
3′-RNA核苷酸磷酸O原子作為廣義堿(實驗驗證部分)
上文提到研究所提出的機制中，G12的pro-Rp O原子作為激發進攻水裂解的廣義堿。因此G12 pro-Rp O原子與進攻水之間的距離是激活的重要因素。為了直接驗證這一機制，后續采用了兩種不同的方案來研究pro-Rp O在內裂中的作用。
?在第一種方案中，通過將末端G12和U11之間的3 ‘ -5 ‘磷酸二酯鍵替換為2 ‘ -5 ‘磷酸二酯鍵，增加了G12的pro-Rp O和攻擊水的距離。其結果使得G12非橋接磷酸O與進攻水之間的距離大大增加，從而使其不能再作為廣義堿。為了證明這一點，對2 ‘ -5 ‘RNA進行MD模擬，結果顯示G12 pro-Rp O與進攻水上氫平均距離為5.47±0.29 A，大于3’ -5的3.44±0.29 A，嚴重影響催化活性。
為進一步驗證上述理論觀察，研究人員使用Fig. 5a中描述的骨架進行了內裂實驗，RNA的G12呈3 ‘ -5 ‘或2 ‘ -5 ‘連接方向。通過對Pol II的內裂實驗分析，可以清楚地看出U11和G12核苷酸兩者之間的鍵合發生了變化，并大大降低了裂解率。在量化裂解速率時，研究人員發現雙相離解方程(double-phase dissociation equations)比單相(single-phase)更適合擬合實驗數據(Fig. 6)。值得注意的是，2 ‘ -5 ‘連接的裂解速率(kfast)比3 ‘ -5 ‘連接裂解速率(慢14倍以上(0.0091±0.0022 min-1 ?vs ?0.14±0.02 min-1)。此外, 3 ‘ -5 ‘連接骨架的快速相（fast phase, 注：雙相法分為快速和慢速相，比例根據實驗擬合結果得出）比例（69.5±6.0%）也高于2 ‘ -5 ‘連接支架（32.1 ± 8.6%）。
?
Fig. 5 ?改變RNA 3 ‘-端的連接降低內裂比率。a,內裂實驗所用核酸序列。b,c, 與3’-5’ 連接相比，2’-5’連接將會大大降低內裂速率。d, pH對反應的影響。e, G12 pro-Rp O原子與進攻水之間的平均距離。f, 三個2 ‘ -5 ‘連接骨架的代表性類反應物結構，標注為G12 pro-Rp O原子與進攻水距離。
圖片來源Nature Catalysis
?
Fig. 6 ?用兩種方法擬合內裂反應產物曲線（單相(a) vs雙相(c)）。右圖為早期時間點的對應放大圖。
圖片來源Nature Catalysis
這些實驗揭示了兩個重要的觀察結果，與G12在回溯內裂過程中的作用有關。首先，雙相擬合結果表明，回溯過程可能同時存在兩種不同的內裂機制。
其次， kfast受回溯核苷酸磷酸二酯鍵連接方式的高度(3’-5’ vs 2’-5’連接)。這一結果明確表明，回溯核苷酸的連接改變會強烈影響其內在的裂解活性。這種改變直接影響了磷酸基的位置，(Fig. 5f)，因此，裂解活性降低了14倍，并且快速相比例也降低。這些結果支持理論計算中提出的磷酸鹽基團在內裂過程中的促進作用。
為了檢驗2 ‘ -5 ‘和3 ‘ -5 ‘連接骨架的內裂機制是否相同，研究人員比較了2 ‘ -5 ‘和3 ‘ -5 ‘連接的pH譜。高pH值能通過去質子化進攻水來刺激磷酸二酯水解，因此pH-依賴的裂解曲線可以為反應機制提供重要信息。實驗獲得的pH譜非常相似(Fig. 5d)，這表明是相同的一組原子參與了裂解反應，這個結果與前面獲得的結論一致。
為了進一步證實G12 pro-Rp O原子的參與，再采用第二種方案進行實驗，在不顯著改變活性位點原子結構的情況下，使用S原子來取代這個關鍵O原子(Fig. 7a)。含有硫代取代物的寡核苷酸會產生兩個立體異構體，即(pro-)Rp-S和Sp-S(Fig. 7a)。采用高效液相色譜法進行分離，從混合物中提取Rp-S和Sp-S寡核苷酸，并測定分離組分的裂解速率。在得到的兩個峰(peak1和peak2)中，peak2的裂解速率顯著降低，而peak1的裂解變化不大(Fig. 5b, c)。這一結果與提出的反應機理是一致的。一方面，立體定向的Rp-S取代了末端磷酸O (Fig. 7a左)會降低去質子化水的能力，造成裂解速率降低，因此實驗Rp-S取代對應peak2。另一方面，在Sp-S取代過程中，O仍然是質子的受體，裂解速率受到的影響很小(Fig. 7a右)，因此認為Sp-S取代對應于peak1。
?
Fig. 7 ?S取代磷酸O對內裂的影響。a,兩種不同的硫取代異構體結構。b,c, 硫代異構體對內裂速率的影響。
圖片來源Nature Catalysis
?
3′-RNA 核苷酸堿基作為可選廣義堿的有限作用
之前的研究表明，雖然細菌和真核生物之間存在差異，但回溯核苷酸的堿基可能直接作為廣義酸/堿參與內裂過程。為了研究鳥嘌呤(G)堿基是否能作為廣義酸/堿參與內裂，實驗將鳥嘌呤與取代N7原子鳥嘌呤的酶催化內裂速率進行比較。結果表明二者對內裂速率的影響都非常小。與這些實驗相一致，以G12的N7原子為廣義堿，進行QM/MM計算的結果顯示活化能勢壘增加了20 kcal mol-1。對G12的O6原子作為廣義堿的可能性同樣進行了QM/MM計算，結果顯示活化能勢壘同樣增加了~20 kcal mol-1(Fig. 8)。綜上所述，G12的堿基對Pol II內裂的貢獻非常小。
?
Fig. 8 ?QM/MM計算排除三個可能作為廣義堿的堿基原子。
圖片來源Nature Catalysis
?
Pol II?殘基不太可能作為廣義酸
考慮到天冬氨酸D485是距離新出現的A10 3’-端O原子最近的蛋白殘基，研究者采用質子化的D485進行MD模擬檢測D485是否可以作為廣義酸，以及質子化的D485在回溯態Pol II的活性位點是否穩定。在所獲得的三段MD軌跡中，其中一個呈現了一個扭曲的活性位構型，U11不與MgA配位，而另外兩軌跡中質子化的D485保持穩定并指向口袋外側(Fig. 9a)。為了進一步研究在去除最近的水分子后，D485的質子是否能夠被成功的容納，計算使用五配位MgA重復了MD。結果表明，在位置約束平衡過程中，質子在離活性位點較遠的D485的穩定性并沒有增加(Fig. 9b)。為了消除源自力場的不穩定因素，在對D485質子進行限制的情況下重復MD模擬。在所有結果中，MgA-MgB的距離在平衡過程中增加到5A附近(Fig. 9c)。因此可以得出結論，對于回溯態的活性位點，質子化的D485不是一個穩定存在的狀態，并且不能作為廣義酸。所以，結果表明Pol II中的蛋白殘基在裂解反應中不直接作為廣義酸參與質子轉移。
?
Fig. 9 質子化D485的三組模型MD模擬。a, 6配位MgA; ?b, 5配位MgA; ?c, 對5 配位MgA及D485質子進行限制。
圖片來源Nature Catalysis
結語
這篇催化機制研究案例文章是非常典型的計算結合實驗的工作。RNA聚合酶II內裂機制中，由于不確定參與反應的廣義酸/堿，一直存在爭議。研究工作在獲得可靠的模型，確定QM區原子及選擇反應坐標后，通過QM/MM-MD模擬找出了可行的RNA聚合酶II內裂機制，分析了磷酸O原子作為廣義堿參與反應的過程，同時將其他的可能性進行逐一排除。而后根據計算結果，巧妙的設計了兩個實驗方案（改變核酸骨架影響磷酸位置和S代磷酸O原子），給出了有力的實驗證據，對理論計算提出的反應機制進行了驗證。從中我們可以學習到在對酶反應機制進行研究時，如何瞄準關鍵科學問題，在一開始構思合理的計算模型，使得計算結果能夠闡明催化機制的同時，更好的配合實驗或在一定程度上指導實驗方案設計。但文中提出的反應機制能壘偏高(29kcal/mol)，作者認為這與RNA聚合酶II內裂速率較慢相符(比受轉錄因子刺激時慢~400倍)，并通過進一步計算，利用類似排除法的方式將其他可能反應途徑進行排除。不過值得我們注意的是在反應途徑能壘偏高時，同樣要考慮到初始模型、反應坐標、所用QM方法和QM區選擇等影響因素，確認計算的合理性，再進行實驗驗證。
?
參考文獻
(1) Tse, C. K. M.; Xu, J.; Xu, L.; Sheong, F. K.; Wang, S.; Chow, H. Y.; Gao, X.; Li, X.; Cheung, P. P.-H.; Wang, D.; Zhang, Y.; Huang, X. Intrinsic cleavage of RNA polymerase II adopts a nucleobase-independent mechanism assisted by transcript phosphate. Nature Catalysis 2019,?2, 228-235.
(2) Udo, L.; Winfried, H. Transcriptional fidelity and proofreading in Archaea and implications for the mechanism of TFS-induced RNA cleavage. Mol. Microbiol. 2010,?52, 1133-1143.
(3) Samara, N. L.; Yang, W. Cation trafficking propels RNA hydrolysis. Nature Structural & Molecular Biology 2018,?25, 715-721.
(4) Mishanina, T. V.; Palo, M. Z.; Nayak, D.; Mooney, R. A.; Landick, R. Trigger loop of RNA polymerase is a positional, not acid-base, catalyst for both transcription and proofreading. Proceedings of the National Academy of Sciences of the United States of America 2017,?114, E5103.
(5) Lu, X.; Fang, D.; Ito, S.; Okamoto, Y.; Ovchinnikov, V.; Cui, Q. QM/MM free energy simulations: recent progress and challenges. Mol Simul 2016,?42, 1056-1078.
(6) Zhou, Y.; Wang, S.; Li, Y.; Zhang, Y. Born-Oppenheimer Ab Initio QM/MM Molecular Dynamics Simulations of Enzyme Reactions. Methods Enzymol. 2016,?577, 105-118.
內容來源自QMCLab，唯信計算整理發布