PacBio技术破译首个木兰类植物——牛樟基因组

牛樟预估基因组大小为823.7± 58.2 Mb/1 C， 2n=24。本研究先用PacBio测序覆盖85×深度，reads N50为11.1kb，contig N50为0.9Mb，组装基因组728.3Mb。然后使用141×Illumina reads进行校正，接着用207בChicago’重组染色质和204×Hi-C双端reads进行染色体挂载。最终，整合基因组大小730.7Mb，含2153条scaffolds，占流式细胞仪估算基因组大小的91.3%，scaffold N50达到50.4Mb，占12条对应于牛樟染色体超过90%的pseudomolecules，且注释蛋白质BUSCO评估不低于89%。

杂合位点的空间分布高度可变，有23.9%的基因组每1kb仅有不到1个SNP位点，而有10%的基因组1kb有至少12.6个。杂合区域呈现随机分布，在第11号scaffold上达到最多20.2Mb（参见图1a），但在序列覆盖上较其他杂合区域更均衡，推断与选择性清除、同系繁殖或最近的种群瓶颈效应相关。同时，位于这些杂合区域上的基因在木质素生物合成过程与半乳糖代谢中尤为富集，可能在木质素-碳水化合物复合物的合成中起到作用(见图1b)。

图1 牛樟基因组杂合度(a.每100-kb非重叠窗口的杂合双等位基因的SNPs数目对最大的12条scaffolds作图，Indels被排除；b.使用PSMC方法推断的有效种群大小的历史；c.对每个100-kb非重叠窗口，分布从顶端到末端为基因密度、转录组和三类重复序列。红色T字母表示scaffold末端存在端粒重复簇，LINE表示长散在重复序列 )

转座元件和散在重复序列占到基因组的48%，最大的类别为长末端重复转座子（LTR RT）（25.53%），其次是DNA转座子占12.67%。在LTR RT中，又分别有40.75%为Ty3/Gypsy、23.88%为Ty1/Copia(图1c)。另外，LTR富集的区域较基因组其他区域平均多覆盖了35%，正是其在组装中被折叠导致了流式细胞仪和k-mer分析对基因组大小评估的差异。

再者，染色体水平的scaffolds在染色体中心区域呈现蛋白编码基因密度低、转座元件分布密度高的特点(图1c)。同时研究还发现了687kb核质体类DNA序列（NUPTs），其中96%是小于500bp的短片段。

研究从上述13个物种中搜寻到了211个单拷贝直系同源数据集，整合为超级矩阵后使用最大似然法进行系统发育树的构建，结果显示木兰亚纲和双子叶植物支系形成姐妹类群关系（见图2）。该系统树在加入额外22个木兰类物种的转录组数据后仍然保持原来的拓扑结构，只是自展支持率较之前稍低。

最后，通过化石年代的标定，研究估算出木兰亚纲和真双子叶谱系于136.0-209.4Ma分化出来，这一结果也得到了最近其他研究结果的支持（图2）。

图2 基于13个物种中211个单拷贝直系同源基因构建的物种系统发育树（+和-旁的数字分别表示基因家族的扩张和收缩；括弧里的绿色数字表示分歧时间的估算；如无额外说明，所有自展支持率为100%）

研究在72.7%的基因组中鉴定到16498个基因对分布于992个共线性区上，这些共线性区又有72.3%与基因组上超过一个位点呈共线性，也就是说：牛樟祖先发生了不止一次WGD事件（参见图3a）。染色体区域成对的广泛共线性以及每个区域与另外两个基因组片段显著、但非共线性的配对和两轮古代WGD事件息息相关(图3a)。

无油樟（Amborella trichopoda ）是其它现生所有被子植物姐妹类群的唯一物种，自现生开花植物谱系最近的共同祖先分化以来没有发生过WGD事件的证据。因此，研究做了一个假设：两轮WGD事件是牛樟祖先在牛樟和无油樟分化之后发生的，两者基因组的共线性分析发现牛樟1-4个片段匹配上了无油樟基因组的单个区域，证实了上述假设（图3b）。

那么两轮WGD事件什么演化节点上发生的呢？研究者又估算了基因组内和种间的同源染色体Ks(同义替换)分布。牛樟基因组内重复片段峰值在0.46和0.76附近(图4a)，正好和两轮WGD事件对应，研究据此来推断核型进化事件。洛矶山耧斗菜（Aquilegia coerulea）隶属于其它所有现生真双子叶植物姐妹类群的毛茛目，其与牛樟直系同源的分析显示主峰在Ks=1.41左右（图4a），而耧斗菜基因内重复片段Ks=1，由此推断在牛樟和耧斗菜谱系分化后发生了独立的WGD事件。通过挖掘17个樟目和木兰目转录组数据的信息，研究发现樟科Ks分布的两个峰与牛樟Ks分布一致，对应上了牛樟祖先WGD事件的两个基于共线性的推断(图4b和图3)。而其他樟目和木兰目物种仅发现一个Ks峰，推断WGD事件在其演化过程中仅发生一次。在耧斗菜数据中观测到的Ks峰可能对应的是真双子叶植物和木兰亚纲分化后毛茛目内部的WGD事件（图4a）。

图3 牛樟基因组的进化分析（a.牛樟基因组中的637个共线性区基因组内关系示意图；共线性区由五色线条表示，代表祖先染色体组型；紫色区块表示分配到第一种颜色区域上的共线性区；b.牛樟基因组内第一种颜色区域和无油樟对应的共线性关系）

图4 牛樟和其他物种之间同义替换的密度图（a.牛樟与洛杉矶耧斗菜及两者间共线性区鉴定到的成对直系同源重复序列；b.樟科和木兰目基因组内成对重复序列Ks，虚线为在牛樟上观察到的Ks峰，棕色和灰色的线条分别为牛樟和其他樟科Ks分布）

研究试图通过对蛋白家族（pfam）域进行注释、评估其在13个用于系统发育分析的种子植物基因组中的分布来鉴定牛樟特有的基因和蛋白域。分析发现，牛樟、真双子叶植物和单子叶植物之间有相当大的重合，说明三个谱系在分化后功能发生了显著的多样化。包括萜烯合成酶（tps）羧基末端结构域在内的蛋白质结构域的获得涉及植物蒸腾效率中的防御反应和富含亮氨酸的重复序列。有意思的是，研究者发现牛樟拥有21个EIL转录因子的拷贝，比之前报道的拥有最多拷贝数的香蕉基因组还多4个；同时EIL是通过激活乙烯应答因子（ERF）来启动乙烯信号应答的，牛樟中的ERF同样高度扩张。通过EILs的扩张来刺激ERF从而实现对下游效应的调控形成了牛樟特有的属性。

接下来，研究还评估了种子植物系统发育中直系同源类群的扩张和收缩。基因家族大小的演化在系统发育过程中是动态的，进化为牛樟的支系没有呈现出显著的扩张或收缩变化。GO富集分析揭示出牛樟不同的基因家族享有共同的功能或者单个基因家族经历了大规模的扩张。

牛樟基因组注释包含了387个R基因模型，其中82%都属于核苷酸结合位点上亮氨酸富集的重复序列（NBS-LRR）或者卷曲螺旋结构的NBS-LRR类型。在13个研究基因组中，牛樟在非栽培植物中含有最多数目的R基因。2465个NBS结构域的系统发育分析也显示基因家族内的分支在真双子叶植物、单子叶植物和木兰亚纲中是独立分化的。引人注意的一点是，牛樟NBS基因最分化的一支与双子叶植物NBS基因最保守的一支形成了姐妹类群。

图5 101个CkTPS基因的系统发育位置

对于牛樟基因组，最显著的特征是其庞大数量的TPS基因（CkTPS）。本研究在牛樟基因组中预测、注释了101个CkTPS，是迄今基因组中发现数目最多的。在加了两个木兰亚纲物种转录组数据并进行系统发育分析后在之前已描述过的7个种子植物TPS基因亚族中明确了6个亚族CkTPS的系统位置(图5和表1)。101个CkTPS基因有7个与催化形成20-碳异戊二烯类化合物的关键酶有关，另外94个很可能编码10-碳单萜合成酶、15-碳倍半萜烯合成酶以及其他20-碳双萜合成酶（表1）。

值得引起注意的是CkTPS基因系统树还解决了TPS-a、TPS-b、 TPS-f和TPS-g基因亚族内部樟科特异性的TPS基因分支。总之，分析表明不论是樟科起源之前还是之后，木兰亚纲TPS基因都处于不断分化之中。

最后研究还分析了不同染色体上TPS基因的分布密度。研究提到，TPS基因在染色体上是不均匀分布的，并且独立亚族中的成员聚类以串联重复的形式存在。在牛樟基因组最大的12条scaffolds上观察到了76个TPS基因，而7号scaffold包含了隶属于多个亚族的29个TPS基因，相比之下，1号scaffold仅含CKTPS-c的两个成员，再有24个CKTPS基因定位在其他小一些的scaffolds上等等。