前言

        小编要分享的这篇文章是2020年10月19发表在Nature Communications杂志上，很荣幸征得作者本人（北京林业大学林木分子设计高精尖创新中心博士生 聂帅）的同意写一篇解读博客，在此表示感谢，小编水平有限，要是有解读不准确的地方欢迎广大朋友指出。文章研究的是杜鹃花，主要做了三部分内容：（1）基因组组装注释 （2）基因组进化分析 （3）转录组数据表达分析。文章大体的研究思路和常见的基因组文章差不多，甚至没有做常见的群体进化分析，但是它尽然能发表在NC上，这和该文章的分析方法创新有关，文章中用到了一些比较新的方法，并且得到了不错的结果，小编断定该文章会出现很多的模仿者

文章链接：Chromosome-level genome assembly of a parent species of widely cultivated azaleas | Nature Communications

摘要

        杜鹃花科（Ericaceae）的杜鹃（Azaleas）构成了最多样化的园艺植物之一，因其文化和经济重要性而闻名于世。本文中，作者报道了杜鹃（Rhododendron simsii）的染色体级别基因组组装，而R. Simsii是目前栽培种杜鹃的主要祖先。基因组分析显示在杜鹃花科植物的辐射扩散之前发生过一次全基因组复制事件，并且该WGD事件可能作用于杜鹃开花时间的基因组结构。小范围的基因复制事件作用于参与杜鹃花色素合成的基因家族的扩张。作者通过详细的时间序列基因共表达网络分析，重构了杜鹃花青素和类胡萝卜素的代谢通路，并揭示了其潜在的调控网络。MYB、bHLH和WD40转录因子可能共同作用于杜鹃中花青素积累的调控，尤其是花着色的起始阶段；另外，WRKY转录因子在后期控制花的渐进着色。本文的工作为理解开花时间和花着色调控的遗传机制提供了基石，可能有助于未来加速杜鹃花的选择育种

引自：Nature Communications：杜鹃花基因组

背景

        杜鹃花（Rhododendron simsii），又名映山红，如下图，双子叶植物，我国十大名花之一，东亚特有，世界范围上种植最广泛的杜鹃花种类

        尽管之前研究发表过两版杜鹃花的基因组，但是都只做了二代测序，做到contig版本水平，该文章用了pacbio数据，并且挂载了Hi-C，最后得到13条染色体，为杜鹃花关键性状遗传机制研究和基因组编辑等高效分子设计育种的建立提供了重要基础

组装注释

如下图

（1）与另外两版发表的基因组比，该基因组用了pacbio数据和Hi-C数据，最后挂载到13条染色体上，全基因组大小为528.6M，共552条scaffold

（2）contig N50为2.23M，BUSCO评估93.68%，LAI评估18.10，这几个常见的基因组评估指标都挺高的

（3）注释得到34170个编码基因，47.48%重复序列，rcRNA的注释（482 tRNAs,64 rRNAs,211 miRNAs, 16 tRNAs and 158 snoRNAs），并给出了一些常用数据路能注释到的基因的比例（NR，SwissProt,TrEMBL,Pfam,GO, KO）

注：基因组注释主要包括重复序列注释，基因结构注释，rcRNA注释，重复序列注释的时候作者详细鉴定了LTR-TRs和solo-LTR，Gypsy 占基因组11.90% ，Copia 占基因组4.00%，整个注释工作比较复杂，用到的软件非常的多，小编这里不进行讲述

基因组进化分析

（1）比较基因组分析

图1a：文章用了15个物种加2个外群，总共聚类得到806个orthogroups单拷贝，然后构建了ML树，并计算了一下分歧时间和扩张收缩，将大分支Ericales上发生的WGD以及γ-WGT事件用不同的图标进行了标注，可以看到，文章研究的杜鹃花和另外两种杜鹃花在14.54个百万年发生分化，在大约78个百万年的时候有一次WGD事件（使用软件：OrthoFinder聚类，MUSCLE比对，IQ-TREE建树，MCMCTree算分歧时间，CAFÉ计算扩张收缩）

图1b：计算Ks值推测WGD事件，下边的图计算了五个物种自己的Ks值，可以看到，在0.65的时候三种杜鹃花有一个峰值，在此发生的WGD事件，右上角的图是葡萄与五个物种blast比对得到一对一的orthogroups，然后计算得到，计算方法见原文，推算成百万年后结果和下边的图大概是一致的

图1c：13条染色体circos图，由外之内分别是Ogre（a family of Gypsy LTRRT）、Tekay（a family of Gypsy LTR-RT）、Gypsy LTR-RT、Copia LTR-RT、基因密度、GC含量、共线性区块，该图特点是将每条染色体的着丝粒位置标出来了，并且将共线性区块分成了WGD（红色）和WGT（蓝色）（咨询了作者，杜鹃花的共线性区块基本上来自WGD和WGT，每个区块到底是WGD还是WGT没有一个标准，存在主观性，个人以为很多物种的共线性区块不仅仅来自加倍化事件，作者在这拓展了很好的分析思路）

（2）基因复制分析

    文章这里用了DupGen_finder软件【DupGen_finder】，该软件能将所有的基因分为五类复制基因：

WGD：全基因组复制

TD：串联重复（相邻的两个重复基因）

PD：近端重复（相隔10个以内基因的重复基因）

TRD：转置重复（祖先和新基因座组成的重复基因）

DSD：分散重复（不相邻也不共线性的重复基因）

SL：单拷贝

对这几类概念感兴趣的可以参考：Gene duplication and evolution in recurring polyploidization–diploidization cycles in plants | Genome Biology | Full Text

    可以看到，全基因组所有的基因有89%可以分为这五类，不足的是这个软件分类得到的结果假阳性也挺高的，该文章验证了假阳性结果是可以接受的，小编认为这篇文章最大的亮点也是最大的创新在于将基因分成了五类复制然后进行接下里的分析

图2a：五类复制的ka/ks值箱线图，可以发现PD和TD的值较其他三类明显要高，说明这两类复制基因发生了更多的正选择，这两类基因在植物应对逆境环境中有重要的作用，尤其是TD，所以文章也重点关注了这两类复制基因

图2b：扩张基因（6,754个）、复制基因、正选择基因（1657个）的venn图，先看最外去个大圈和最中间的黄色圈，这是正选择（EGFs）和五类复制的venn图，能看到每一类复制并且扩张基因的个数，然后看大圈里边的小圈，黄线以外是正选择和复制基因的交集，黄线以内是正选择和扩张基因的交集，最后发现所有发生正选择和扩张的基因中有57.6%由TD和PD引起

图2c:对发生扩张基因的五类复制进行了GO富集，可以看到，尤其是TD，富集得到一些重要的途径

图2d:通过FLOR-ID数据库得到与花周期相关的基因家族序列，然后鉴定得到杜鹃花有424个基因与花周期相关，图片展示了五类复制中与花周期有关的基因个数，最后发现了一个OG0000614基因家族（高亲和性蔗糖转运蛋白，转录表达证明在花发育过程中持续高表达）

转录组数据表达分析

（1）类胡萝卜素和花青素/黄酮醇生物合成基因家族转录表达分析

    文章用了五个不同时期的花组织测了表达量（图4a），大多数基因组文章测的是不同组织器官上的表达，文章后边用到了一个时间序列分析软件来分析这有连续时间差异的五个时期表达量，这里也是文章的一个亮点和创新点

图3：文章鉴定了和类胡萝卜素、花青素/黄酮醇相关通路的基因家族，上图是花青素/黄酮醇生物合成的三个重要途径（Anthocyanin 花青素、Flavonoids类黄酮素、Flavonol黄酮醇）的基因家族表达分析，蓝色ID标签是TD，棕色标签是PD复制，在所有这些和花青素/黄酮醇生物合成相关的基因中TD/PD占了42%，可见这两类基因在花色调节过程中有很重要的作用。但是有意思的是，TD/PD在所有花周期基因家族中的占比很小（图2d），总共只有13.66%，花周期大多数由WGD引起。

 （2）时间序列共表达网络分析

    文章这里用了TO-GCNs（time-ordered genecoexpression networks）软件，，时间序列基因共表达分析软件，2019年发表在PNAS上，文章链接【TO-GCNs】，小编没有用过，听作者说挺好用的，输入文件为全部的表达基因和关注的表达基因，就能得到共表达网络。文章发现了8067个高表达基因（618个转录因子TFs和7,449个结构基因），用来构建共表达网络

图4a：五个不同时期的花组织，又分成了三个大的时间段initial：T1 、transitional：T2,T3 、terminal：T4,T5

图4b：TO-GCNs软件做出来的时间序列表达网络，文章做的主要是618个TFs的共表达网络，分成了8个系列（结合表达热图定义成8个系列），粉色圈内的点是五类转录因子家族（紫色点：MYB，橄榄色点：bHLH，绿宝石色：WD40，蓝色：WRKY ，绿色：ERF）的共表达情况，可以看到尤其在initial阶段共表达特别明显

图c：五个不同时期8个系列的表达热图，可以看到分成八个系列能够很好的展示出时间序列共表达情况

（3）initial阶段花色有关基因共表达网络分析

    上边的时间序列共表达看到，在initial阶段共表达特别明显，所以文章对这个阶段进行了与花色有关基因的共表达分析，最后推测MYB、bHLH、WD40三个转录因子家族可能构成复合体从而共同决定着花色形成

图5a：类胡萝卜素生物合成共表达网络

图5b：花青素/黄酮醇生物合成共表达网络

图5c：花青素/黄酮醇相关基因和MYB、bHLH、WD40基因家族的共表达

图5d：F3H、MYB、bHLH、WD40在五个时期的表达热图

（4）terminal阶段花色有关基因共表达网络分析

    除了initial阶段，terminal阶段也有很强的共表达，所以最后还做了这个阶段与花色有关的共表达分析

图6a：花青素/黄酮醇生物合成共表达网络

图6b：F3H、MYB、C2H2、C3H、GRAS五类转录因子基因上游2kb DNA结合位点情况

图6c：F3H基因的分级调控

注：对于各个表达网络更深层次的研究意义小编也是看的很糊涂，大家感兴趣可以仔细研究原文，能力有限，不再扩展。

总结

（1）该研究结合 Pacbio 长读长测序与Hi-C技术，获得了映山红杜鹃高质量的基因组组装和功能注释

（2）分析了基因组组成特点，识别了开花调控基因和花色相关代谢通路

（3）基于时间序列的基因共表达分析揭示出与花色决定相关等级基因调控网络，识别到 MYB, bHLH, 和 WD40 三个转录因子家族成员可能构成复合体从而共同决定着花色形成