客户文章 | Plant Physiology玉米m6A进化由基因组重复事件介导
论文标题：Evolution of the RNA N6-methyladenosine methylome mediated by genomic duplication
刊登日期：2019年8月
发表杂志：Plant Physiology
影响因子：7.394
研究机构：西北农林科技大学马闯课题组
RNA上的m6A修饰是真核生物中最丰富的RNA表观遗传修饰形式。鉴于m6A进化与哺乳动物基因组中核苷酸序列的选择性限制有关，作者假设至少部分程度上m6A进化与复杂多倍体植物基因组中的基因组重复事件有关。
为了验证这一假设，作者对两个基因组亚型的玉米叶片与幼苗（干旱胁迫处理）进行了m6A测序及转录组测序，并从两个亚基因组（分别为maize1和maize2）中分别鉴定了5893和3811个基因，这些基因总共鉴定了11968个unique peaks。
这些亚基因组中的每一个样本都覆盖了多达2200个重复基因。在这些重复基因中，那些携带m6A Peak的基因在保留率(retention rate)上有显著差异。
m6A甲基化基因对出现退化并在这种在亚基因组上出现的功能分离现象与多重序列特征相关，并受不对称进化速率的影响。
此外作者还描述了由全基因组重复及串联复制事件（whole genome duplication and tandem duplication）介导m6A甲基化基因和转座因子的共同进化模式，揭示了m6A重复功能因子的进化保守性和多样性，以及玉米RNA上的m6A修饰在干旱胁迫中的潜在作用。
这篇论文强调了m6A修饰与基因组复制事件之间的复杂互作关系，为理解由基因复制事件介导的多倍体植物m6A修饰进化机制提供了参考。
在植物中，基因组重复事件是植物基因组进化的主导力量。基因组复制事件包括全基因组复制、片段复制、串联复制和易位复制，来产生特定基因组背景（即复制区域）的来源。
在基因组重复事件之后，重复的基因受到不同程度的纯化选择，其中一部分在一个被称为分离（fractionation）的过程中丢失。
基因组中还有许多重复基因作为paralog pairs被保留了下来，其中单个基因可能被亚功能化（分割和共享原始基因功能）或获得新的功能。
据报道，分离偏差和复制基因的功能差异与DNA甲基化、转座子的移动率、基因表达和转录后调控的差异有关。
因此，基因组重复事件提供了一个特定的基因组背景来源，可能对转录调控和转录后调控产生深远的影响。
作者采用Well Water即WW以及干旱胁迫DS一共两种处理方式，使用玉米品种为B73和Han21，每组三个生物学重复总计12个IP文库12个Input文库。
其中，Han21转录组及m6A-seq数据比对的是Han21 Pseudogenome假基因组。即通过STAR比对到B73基因组后，用GATK进行SNP Calling获得总计94761个exon区域SNP后对原有的B73基因组进行替换，获得Han21假基因组。
我们可以看到12个样本peak归一化后整体趋势图，基本上都集中在3’ UTR和终止密码子附近。
而这12个样本获得总计11219个unique genes中的11968个unique peaks。
无论是B73还是Han21，超过70%的Peaks都落在了3’ UTR区。
另外为了保证准确性，作者使用了MEME和HOMER两款软件进行motif分析。最终结果表明有大约超过90%的peaks含有m6A的经典motif RRACH（R=A/G; H=A/C/U）。此外另一个属于植物m6A特有的motif URUAY（R=A/G; Y=C/U）也被鉴定了出来。
值得注意的是，URUAY并不是在本文中首次被报道。早在先前拟南芥m6A的文章中，一种甲基化阅读蛋白ECT2中就有很典型的URUAY motif（案例解析：拟南芥去甲基化酶ALKBH10调控成花转变 | m6A专题；案例解析：m6A-YTH组件调控拟南芥叶片发育时间和形态发生 | m6A专题）。
Dot blot实验表明，m6A抗体能够与带有m6A修饰的URUAY motif结构的oligo特异性结合。
具体来说就是人工体外合成一段带有m6A修饰的oligo且具有URUAY motif结构，而在碱基序列不变的情况下没有m6A修饰的oligo作为对照。下一步就是与m6A抗体进行binding实验。
作者把一个含有2个m6A peak以上的基因称之为高甲基化基因（总计660个unique），1个m6A peak为低甲基化基因（10559个），没有m6A peak则是非甲基化基因（总计28105个）。
另外比较重要的是，m6A基因在SNP数量显著低于非m6A基因，这从另一个方面证实了m6A进化上的保守性。
进化过程中，玉米经过基因组重复事件后从高粱开始出现了分化。当一对复制基因其中一个拷贝丢失了之后，另一个基因会进化成singleton。
由于基因biased fractionation存在，两个玉米基因型之间的重复基因丢失并不同步或均匀，这使得其中一个玉米基因组中比另一个玉米基因组中保留了更多的基因。
显然m6A基因在Maize2 singleton中频率显著高于Maize1。
m6A基因在两个亚基因组之间的这些差异很可能是亚基因组biased fractionation的进化结果。
Maize1中m6A singletons比在Maize2显著提高。这表明m6A基因的biased fraction可能与基因长度有关。
下一步作者比较了Maize1和Maize2两个亚基因组的进化速率ω。在高粱和玉米两个物种之间针对基因同源序列进行Ka/Ks分析。m6A基因的进化速率显著大于非m6A基因。m6A基因的Ks值显著小于非m6A基因。这表明m6A基因的进化速度快，进化时间短。
作者还发现，Maize1的m6A Singleton的进化率明显低于Maize2的m6A Singleton，但Maize1和Maize2的非m6A Singleton的进化率没有显著差异。
结果表明，Maize1的m6A Singleton比Maize2的m6A Singleton有更高的纯化选择度。这种不对称的纯化选择可能会对m6A Singleton的biased fractionation产生影响。
相比之下，Maize1和Maize2中的m6A duplicates在相同的纯化选择水平下进化，但Maize1中非m6A duplicates的进化速率明显低于Maize2中非m6A duplicates的进化速率。
这表明m6A的修饰可能与两个亚基因组中重复基因间进化速率的差异有关。
对m6A duplicates的进一步分析表明，进化时间（Ks）与具有不同m6A模式的重复基因对有显著差异。
基因转座转位（Gene transposition）可能导致syntenic duplicates分离为两个Singleton。转座基因对m6A的分化比例明显高于未转座的同基因重复基因对。
这些结果表明，新的duplicates对m6A的修饰差异小于旧的duplicates对m6A的修饰差异，基因转位转座可以增强duplicates pairs之间m6A的分化程度。
对转录组和m6A-seq进行联合分析表明，m6A基因在转录组里表达量显著高于非m6A基因，m6A基因的Singleton duplicates率显著高于非m6A基因。另外Singleton的基因表达量显著高于duplicates的基因。
作者定义IM为为identical-m6A（两个duplicate gene pair或partners都有m6A peak），DM为diverged-m6A（其中一个gene duplicate pair或partner带有m6A peak），而Non-MP则是都不带有m6A peak。
总的来说，基因表达分化程度DM显著大于IM，同样基因表达丰度也是DM显著大于Non-MP。这些结果表明m6A修饰更可能发生在主动转录的基因上，并且可能与重复基因对的保留率和表达差异有关。
另一个有意思的点是，m6A基因更加靠近转座子TE且与TE相关的基因m6A水平更高。此外m6A基因进化速度更快，TE相关m6A基因进化速率显著高于非TE相关m6A基因。
这些证据表明经过基因组重复事件后m6A基因出现了较为宽松的选择，开始更加靠近转座子TE区域并聚集。这一趋势预示着m6A基因可能与TE协同进化。
尽管都参与到串联重复事件的m6A基因与非m6A基因都有较高的替代率（substitution rate），但是m6A基因有更强的选择性。
之后作者发现串联重复TD基因中m6A的频率明显低于非TD基因。TD基因簇中DM模式与IM模式的比率明显高于基因组重复事件duplicates中观察到的比率。
总之，m6A相关TD基因与转座子的距离明显小于非TD基因；m6A基因是转座子TE相关基因。这一比率在非TD基因中显著降低。
这些结果表明，m6A TD基因的进化过程伴随着TD事件过程中偏向转座子的富集积累。
作者对将拟南芥、水稻、人、酵母中的writers、erasers、readers与玉米中的基因进行了进化分析后发现，玉米中MTA、MTB、HAKAI家族、FIP37家族、ALKBH9以及ECT家族都与拟南芥水稻等物种在进化上高度保守。
同时作者鉴定了基因组重复事件的五个同源家族基因对：ZmFIP37a/ZmFIP37b、ZmHAKAIa/ZmHAKAIb、ZmECTa/ZmeECTb、ZmECTd/ZmECTe以及ZmALKBH10a/ZmALKBH10b。
对于ZmHAKAIa/ZmHAKAIb和ZmECTd/ZmECTe来说，这些基因彼此之间进化速率相似。相反的是，另外三对基因进化速率相差很大。
如此反差大的结果表明，m6A功能基因在进化上保守性和多样性都有可能被基因组重复事件抑制。
干旱胁迫处理后，去甲基化酶ALKBH10家族以及ECT家族的mRNA表达量都出现了显著上调，m6A整体水平下降。这表明干旱胁迫诱导去甲基化酶高表达，从而引起了RNA整体m6A水平出现下降。
无论是Han21还是B73，差异甲基化Peak（DMP）绝大部分都富集在3’ UTR。
B73 DMP GO富集分析表明，3’ UTR DMP大量富集在蛋白磷酸化（protein phosphorylation）和细胞凋亡调控（regulation of apoptotic process）上，终止密码子DMP大量富集在组蛋白替换（histone exchange）和细胞质翻译（cytoplasmic translation）上。
Han21 DMP GO富集分析表明，3’ UTR DMP大量富集在转录调控（regulation of transcription）和蛋白转运（protein transport）上，DMP附近的基因则大量富集在光修复/光复活（photoreactive repair）以及siderophore biosynthetic process。
以上结果表明B73和Han21的m6A基因差异较大，调控上有各自独有的调控模式。
Han21对干旱胁迫耐受性显著高于B73，表型差异极其明显，具体包括株高增加、含水量失水率变化以及根系发育等。
下一步，作者对低甲基化基因进行GO富集分析。在B73中，这些基因主要集中在表皮发育、凋亡信号负调控通路、脱落酸ABA等。
而在Han21中，这些基因大部分富集在非生物应激功能如氧化应激反应、反渗透应激反应、乙酰辅酶A代谢、乙烯介导信号通路传导等，其他还包括胚胎发育、葡萄糖与淀粉代谢、种子休眠等。这些基因功能与玉米自交系的表型反应有明显的对应关系。
VACUOLAR INVERTASE 2（VI2）是拟南芥根系伸长的一个调节因子。在干旱胁迫期间，Han21中VI2上m6A peak甲基化水平降低了两倍以上。
参与拟南芥根系生长、细胞分裂和根系结构的肌动蛋白基因家族蛋白Actin-7（ACT7）m6A peak在B73干旱胁迫样品中也呈现甲基化。
此外，表皮蜡质层积累越多越能限制水分流失，增加缺水耐受性。B73和Han21之间与蜡质层富集相关的基因CER4和CER10（两个编码醇形成脂肪酰辅酶A还原酶的基因）的m6A peak甲基化水平的显著差异。
此外，在编码拟南芥茎蜡油酯生物合成所需的二酰甘油酰基转移酶基因WSD1也受到旱灾胁迫的严重抑制。
最后，作者单独进行了m6A-IP-qPCR验证后发现，与测序结果吻合一致。
知道这些数据库，轻松玩转microRNA靶基因预测 | miRNA专题
结核标记物miR-423-5p升高通过抑制自噬体-溶酶体融合在开放性肺结核发展中起重要作用 | 用户文章
Nat Commun：利用小RNA测序揭示植物抗病基因与miRNA以及phaisRNA之间的调控关系 | miRNA专题
一篇miRNA跨界调控研究思路的解析——植物miR159抑制乳腺癌 | miRNA专栏