Nat Commun | 含m6A的RNA二级结构预测
本次给大家分享是2022年3月发表于Nature Communication的一项研究，论文标题为Secondary Structure prediction of RNA sequences including N6-methladenosine。
课题组介绍
通讯作者是美国罗契斯特大学的David H. Mathew教授，主要研究方向为结合计算机和生物信息学，开发靶向RNA的药物和基于RNA的分子工具，目前已与其他课题组合作开发了预测RNA二级结构的RNAstructure软件。通过该软件，73%的碱基对在一组具有已知结构的不同序列被正确预测。
David H. Mathew 教授
研究背景
RNA分子二级结构预测已用于：①确定miRNA结合位点；②设计siRNA；③确定RNA与蛋白质的结合位点；④发现功能性RNA结构。但RNA具有包括肌苷的脱氨、假尿苷化、胞嘧啶的5-甲基化和腺苷的N6-甲基化等多种修饰，其中腺苷的N6-甲基化（m6A）在mRNA和lncRNA中十分普遍，它具有：①引起RNA结构转换；②暴露蛋白质结合位点；③调节剪接等功能。目前，由于缺乏相关热力学参数及配套软件，RNA二级结构预测研究中无法兼顾m6A修饰的影响。
工作概述
本文的研究思路为：①利用对RNA测得的热力学参数，为m6A、A、C、G和U的折叠字母制定了一整套近邻参数（近邻模型：假设一个碱基对对热力学性质的贡献只取决于两个相邻的碱基对），本文着重m6A在螺旋的中间、悬垂端和末端错配的近邻参数上；②利用RNAstructure软件预测含m6A的RNA对应的二级结构（图1）。作者还将修改后的该软件运用于已知的lncRNA转移相关的肺腺癌转录物（MALAT1），为其N6-甲基腺苷引起的蛋白结合亲和力的变化提供计算模型。
图1 研究方法概览
结果与讨论
围绕Watson-Crick碱基对、G-U和m6A-U，进行近邻参数拟合。通过获取29条含m6A的RNA双链的热力学参数，与非甲基化的相同序列相对比，考察甲基化对折叠稳定性的影响。图2显示甲基化的ΔΔG°37高度依赖于序列(ΔΔG°=甲基化序列ΔG°-未甲基化序列ΔG°)，可以看出在+2.1 - -0.1 kcal/mol不等，尤其是当m6A在螺旋的中间时，ΔΔG°值都是正的，说明m6A对折叠是不稳定的，而在末端时则基本没有区别甚至有助于折叠的稳定。
图2 甲基化ΔΔG°37高度依赖于序列
通过线性回归拟合，图3A显示有无甲基化A-U对的相同堆积相比的增量，根据图1，作者预期m6A-U对的近邻堆积比A-U对的堆积要不稳定，最不稳定的是蓝圈标注的2个连续的m6A-U。从整体看m6A-U即蓝色柱状图，不难发现m6A-U与GC相邻时比跟其他碱基对相邻更加稳定。图3B 显示了4条RNA duplex的序列（互补链未标出），它们都具有相同的碱基对堆积，下面2条含m6A的相同序列自由能差值仅为0.18 kcal/mol，仅为总自由能变化的4%，表明位于末端的m6A-U并不破坏稳定。
图3 m6A-U与A-U碱基配对稳定性的比较
对于二级结构的预测，如图4所示，m6A作为悬垂端或作为末端错配中的一个成分，比A能更大程度地稳定二级结构的形成。在3′和5′悬垂端，m6A比A要稳定-0.43 ± 0.15千卡/摩尔。在末端错配中，m6A比A平均要稳定-0.28 kcal/mol，这种稳定性取决于序列，ΔΔG°37的范围从-0.74 kcal/mol到+0.02 kcal/mol。
图4 m6A比A堆积在螺旋上更能稳定RNA二级结构
为了测试获得的近邻参数，作者开展了光学熔融实验。从测试中可以看出近邻参数是准确的，偏差都在实验和近邻参数的不确定性范围内，可以用于RNA二级结构预测。作者指出近邻参数存在一定局限性，即不考虑单个凸起核苷酸的序列特征，因此对其进行建模是在当前的功能形式之外。
作者将一整套近邻参数运用于RNAstructure软件，测定了lncRNA转移相关的MALAT1的N6-甲基腺苷引起的蛋白结合亲和力的变化。MALAT1是核不均一核糖核蛋白C (hnRNP C)的结合靶点，已经证实未甲基化的该RNA且形成单一的茎环结构，m6A可以使该RNA与hnRNP C更容易结合。
图5 m6A最近邻参数测试及RNAstructure软件测试
作者预测32个核苷酸RNA的最低自由能结构，为一个发夹结构，与此前研究的未甲基化的序列一致。同时，软件将甲基化位点预测在A22位点上，且有无甲基化都没有改变RNA的二级结构。对于hnRNP C的结合，由第10位开始的五个Us组成的位点必须是不配对的。因此，作者用邻近发夹环的三个碱基对断裂的ΔG°37来估计可及性（图5A，开放结构）。ΔΔG°37是指在第22位有无N6-甲基化的蛋白质可及性差异为0.6±0.4千卡/摩尔，作者提出的模型是MALAT1发夹的A22 N6-甲基化并没有打开结合位点，而是通过去稳作用促进了hnRNP C蛋白介导的结合位点碱基对的断裂。为了验证这点作者还用chemical mapping和核磁共振谱图并结合其他研究的enzymatic mapping等手段确定该RNA有无发生甲基化都是形成相同的发夹结构（图6）。
图6 MALAT1发夹的酶学图谱、化学图谱和核磁共振谱图
为了进一步测试m6A近邻参数和软件，作者预测了18,026个mRNA的结构，这些mRNA被转录组测序确定含有m6A，并且有PARS结构图的数据。作者使用近邻参数和RNAstructure来估计甲基化位点被埋在螺旋中的概率，即在两个其他碱基对之间堆积的碱基对中，使用局部序列的800个核苷酸片段来估计配对概率，因为作者以前发现800个核苷酸片段的配对概率估计与全局二级结构预测的概率密切相关。
作者发现甲基化位点的未甲基化的A比相邻的核苷酸更不可能被埋在一个螺旋中，相对于A，m6A被埋在螺旋中的概率有很大的变化（A为21%，m6A为13%）（图5B，位置0为甲基化位点）。这表明可能有广泛的结构转换受到了N6-甲基腺苷的影响。我们还可以将我们的结果与同一序列的PARS数据进行比较（图5C）。PARS分数量化了局部配对的酶切估计。较低的PARS分数表明核酸酶S1相对于核酸酶V1有更大的裂解作用，因此有更大程度的不配对，因为核酸酶S1对环有特异性，而核酸酶V1对螺旋有特异性。对这种差异的一个可能性是：PARS将S1裂解归于裂解位点的5′的碱基，假设裂解的5′的碱基是未配对的。当裂解位点的碱基3′未配对时也会发生裂解，因此PARS对甲基化位点5′的评分可能高估了未配对的倾向，因为有些未配对的倾向应该归于甲基化位点。例如，先前的5S rRNA结构的S1图谱与未配对核苷酸的5′和3′的裂解一致。之前对甲基化位点的PARS评分分析也得出结论，数据表明m6A位于碱基配对区和环区之间的过渡结构中，与结构预测估计一致。
总结
研究团队证明m6A在结构中的位置决定了相对于有A的相同结构，折叠稳定性的变化，现在获得堆积参数可以量化这种序列依赖性的变化。此外，通过RNAstructure（https://rna.urmc.rochester.edu/RNAstructure.html）可以对具有m6A的序列进行结构和折叠稳定性的定量预测，作者未来还将扩展到其他修饰碱基上。
