日本朝颜（牵牛花）的基因组测序和分析
日本那部号称世界上最早长篇写实小说的《源氏物语》第四卷《夕颜》中有这样的一句话：名花褪色终难弃，爱煞朝颜欲折难！此处对“朝颜”有个注解：朝颜即牵牛花。呵呵，还挺诗意好啦，我们直接切入正题，看看Japanese都对牵牛花做了些啥子，竟然发了NC嘞~
篇幅有限，本文就不po原文啦，要看的自行百度唷
研究背景
知道么，牵牛花和番薯一样都属于番薯属喔，番薯属还是旋花科中最大的属！（厉害了wuli牵牛花）牵牛花（日本朝颜）一直以来都被作为模式生物来研究花卉性状的遗传基础，并且拥有超过1500个突变系。本文作者利用二代和三代测序数据，组装出牵牛花基因组，进而将scaffolds挂载到了15条染色体上。然后利用组装基因组检测Tpn1家族的转座子，这类转座子作为形成牵牛花的诱导因子，以及在分析矮杆基因CONTRACTED中的重要作用而广为人知（拓展：牵牛花属于旋花科，但矮牵牛就属于茄科啦）。最后对旋花科和茄科物种进行了比较基因组学分析。
材料
牵牛花TKS株系
测序方法
通过PacBio测序得到39.4 GB，共52.6×三代数据，read N50为10.3 kb，最长read长度为48.1 kb，平均read长度为6.8
kb。构建300 bp和500 bp的2个PE150小文库和插入片段分别为3 kb、5 kb、10 kb（*2）、15 kb及20 kb的6个大文库，利用Illumina HiSeq测序获得906×二代数据。
研究结果
1、基因组组装
作者先通过流式细胞术评估确定基因组大小在750 Mb左右，再利用PacBio数据初步组装到contig水平，得到基因组大小为736.4 Mb，contigN50为1.83 Mb。将二代小文库与基因组比对，共检测到1,532个SNP，20,479个deletion和6,549个insertion，进而得出组装的单碱基准确率达到99.99%，看来组装效果还是不错哒~接着利用二代小文库数据对组装基因组进行纠错，并去除了长达1.15 Mb的线粒体和叶绿体序列。然后利用二代大文库数据进行scaffold搭建，并利用PacBio数据进行补洞，最终组装基因组的scaffoldN50达到3.72 Mb。作者还利用Soapdenovo2进行了组装，得到的基因组大小为1.1 Gb，而当只考虑1 kb以上的scaffolds和contigs时，基因组大小则为768 Mb，scaffoldN50为3.5 Mb，contigN50为9.5 kb。
2、错拼检测和染色体构建
scaffold水平的基因组组装完成，接下来就该挂载染色体啦~
首先利用基于Illumina的RAD-seq（简化基因组测序）技术对两个亲本和207个子代样本进行测序，分别获得86.1 Mb的亲本reads，和562.2 Mb的子代reads。挂载前，先将reads与组装基因组进行BWA比对，过滤掉重复及没有限制性酶切位点的reads。为了使得最终80%以上的样品中含有marker，并且每个样品中至少含有80%的marker，作者又利用STACKS来检测SNP并过滤marker。接下来，就是利用Onemap构建遗传图谱了。
根据遗传图谱检测并纠正错拼造成的scaffold嵌合体。如果一条scaffold含有一串指向两个不同连锁群的连锁标记，并且两种标记之间以N（即gap）相连，就将其作为嵌合体在N处分开。对于contig水平的嵌合体，则将其分为三部分，头尾的部分分别属于两个不同的染色体，中间部分依旧作为嵌合体保留。在scaffold搭建完成后，进行了第一次scaffold嵌合体切割，断开了52条scaffolds，而补洞完成后，又断开了29条scaffolds。但是contigN50和scaffoldN50依旧能达到1.87 Mb和2.88 Mb。最后将scaffolds之间以gaps连接成染色体，然后根据marker的顺序确定scaffold方向，marker不足的忽略方向问题，但依旧会作为染色体的一部分。最终染色体覆盖率达到了91.42%，其中包括25.53%未定向的scaffolds。
3、组装评估
CEGMA和BUSCO都可以用来评估组装完整性，作者就利用这两种方法相互印证。CEGMA评估结果显示，组装完整性达到了94.35%，拥有99.60%的核心基因；BUSCO结果显示完整性也达到了95%。进一步评估表明，在93,691个ESTs中有99.11%的支持率，其中coverage在90%以上的EST占97.40%；而20,874个BAC对scaffolds和染色体的支持率则分别是94.92%和97.78%；将转录组数据与组装基因组比对，发现来自胚胎组织的转录组比对率为94.7%，而其余5个组织（根、茎、叶、花、胚胎、种子）的转录组比对率则为96%。这说明牵牛花中保留了大部分保守的核心基因，并且组装质量高。5条约100 kb的全长BAC序列都能完整覆盖到scaffolds上，其中一条BAC序列包含了12.6 kb的Tpn1家族的转座子——TpnA2，这也说明了长序列和高拷贝数的重复元件都已成功组装出来。
举个例子证明，作者在30条携带有端粒重复元件（AAACCCT）的scaffolds中检测到串联重复序列，其中13条因为完全由串联重复序列组成，而无法挂载到遗传图谱上。并且，作者在染色体2、6、8、14号的两端都检测到了串联重复序列，但在染色体3、4、5、9、10、12、13、15号中只在一端检测到。另外还发现18s、5.8s和25s都与串联重复区域中的NOR（nucleolar organizer regions）区域相关，但5s rDNA序列却离NOR区域很远。
4、重复序列分析及Tpn1转座子识别
完成了基因组组装，就可以进一步注释分析啦~
首先使用RepeatModeler预测重复序列，得出基因组重复序列含量为63.92%，其中LTR（长末端重复序列）占比最高。而LTR中，copia和gypsy元件分别占基因组的12.92%和14.46%。通过识别末端反向重复序列（terminal inverted repeat,TIR）和靶位点重复序列（target site duplication，TSD）检测到339个Tpn1家族的转座子。在这些转座子中，除了一个转座子含有的TSD序列是5 bp之外，其他全都只有3 bp。由BLAST比对结果显示，大部分Tpn1转座子在5’和3’末端区域都含有SRR序列。已知TIR和SRR是转座子转录的顺式调控元件（cis-requirements存在于DNA链上，对转录调控起作用；trans-不存在于DNA链，但与cis结合，对转录调控起作用的反式作用因子），那么可以推测Tpn1转座子很有可能是可以转座的。但同时发现32个Tpn1转座子的SRR区域中包含了大量重排，所以这些转座子也有可能已经失活了。有29个Tpn1转座子存在于基因的5’UTR和内含子区域，这有可能会干扰到基因功能。作者猜测Tpn1家族自发的转座子可能同时拥有TnpA和TnpD转座酶编码序列，如En/Spm。为了证实这个猜想，作者将玉米和金鱼草中的TnpA、TnpD序列与牵牛花中339个转座子序列BLAST比对，发现了两个与TnpD同源的转座子——TpnA3和TpnA4，但没有找到与TnpA的同源拷贝。在预测基因或转录本中也没有找到TnpA和TnpD对应的转录本序列，由此猜测转座酶在TKS株系中的转录是沉默的。在Q1072株系中Tpn1具有转录活性，并且存在与TnpA、TnpD同源的cDNA序列——TnpA1、TnpA2。因为TpnA1具有功能性的TIR和SRR序列，因此作者推测TpnA1是一个自发性转座元件，但并没有在该转座子中发现TpnA和TpnD的编码序列，也没有发现单独编码TpnA的序列。尽管在Tpn1转座子中没找到，但作者还是在基因组上分别找到了疑似TnpA和TnpD编码序列的同源拷贝，并且其转座酶的氨基酸序列十分保守，且与已知的En/Spm、金鱼草的Tam1转座酶拥有相同保守域。
5、基因预测和功能注释
一切按照精细图流程进行…于是重复序列预测之后，基因预测和功能注释登场...
利用来自叶子、花、胚胎、茎、根、种子，6个组织样的转录组数据进行转录组预测，同时以番茄为参考物种，利用Augustus进行从头预测。最终确定了42,783个基因，45,365个转录本。其中，44,916个转录本含有完整的ORF（开放阅读框，有确定的起始和终止密码子），并且95.54%的转录本存在于15条染色体中。17.52%的基因仅含单个外显子，2/3的转录本拥有的外显子数目小于等于5个。对预测基因进行功能注释，最终61.99%的基因注释到UniProt-Swiss-Prot数据库，未注释上的基因中又有16.93%注释到UniProt-Trembl数据库中。另外，61.92%的基因有Pfam结构域支持。最终共注释到79.12%的基因。
6、矮基因CONTRACTED的分析
常规分析之后，抓个重点细细研究~
隐性突变基因ct造成的矮杆植株拥有深绿色、厚且皱巴巴的叶片和子叶，并且花和种子也非常小。在遗传图谱上定位到ct基因在LG5位置，与编码花青素合成基因A3的位点仅相距1.2 cM。因为油菜素甾醇（BR）能促使ct和kobito（kbt）共同突变，并且kbt与star（s）是等位基因，而s与拟南芥中编码BR合成酶的DET2基因是同源基因，因此作者推测矮基因CONTRACTED（CT）是一个BR合成基因。众所周知，拟南芥中的ROT3基因编码P450蛋白，可以催化BR前体中的C-23羟基化。而作者就在牵牛花名为BDFN01000805的scaffold中，距离A3基因129 kb的位置，找到了ROT3的同源基因，将其命名为INIL05g09538，进一步证明了作者的推测。通过比较发现，在矮植株中杆19个ct突变基因的第一个外显子都有Tpn1转座子插入，但在正常植株中都没有该插入。另外，在牵牛花中检测到了3个ct等位基因：ct-1，ct-2，ct-w。其中ct-1和ct-w带有Tpn1家族转座子Tpn14和Tpn15，而ct-2由ct-1缺失5’端包括Tpn14部分序列在内的36 bp序列产生。最后利用RT-PCR比较BR合成酶中突变基因的转录本水平发现，转座子的插入会强烈抑制植物下胚轴基因中完整转录本的积累。也就是说，矮基因ct中Tpn1转座子的插入，可能就是形成矮杆植株的原因。
7、比较基因组学分析
最后，作者将牵牛花与近缘物种进行比较基因组学分析，来探索牵牛花特性。利用OrthoMCL对水稻、葡萄、猕猴桃（菊分支），以及茄科植物番茄、土豆、辣椒进行基因家族聚类，共找到1,353个单拷贝基因。4个茄目物种（牵牛花、番茄、土豆、辣椒），共同拥有10,549个基因家族。相对于茄科物种（番茄、土豆、辣椒），牵牛花拥有2,242个特有的基因家族。接着以水稻为单子叶植物外群，利用RaxML根据最大似然法构建进化树。然后利用BEAST评估得出牵牛花与其他茄目物种的分化时间在7,525万年前，这一结果在TTOL数据库中也得到了印证。
利用MCScanX分析共线性，在2,275个共线性区块中找到了17,376对同源基因存在与于染色体上。牵牛花与番茄仅拥有47.05%的同源基因，而两者与猕猴桃的共线性基因却极为相近：34.89%和36.01%。作者猜测造成这种现象，可能是由于猕猴桃发生了两次近期的全基因组复制事件（whole genome duplication，WGD）导致的。茄科物种在7,100±1940万年前共同发生了一次近期的WGD事件。根据Ks分布图发现，番茄自身的WGD事件发生在番茄和牵牛花的WGD事件之后，猜测这可能也是茄科特有的WGD事件，也就是说发生时间在7,525万年前。同时，在牵牛花中也发现了一个Ks峰，作者猜测这可能是旋花科特有的WGD事件。与茄目中直系同源基因相比，旁系同源基因存在基因家族的扩张现象。对牵牛花特有的基因家族进行注释，发现其中富含与授粉、生殖过程相关的基因。
结果讨论
实践证明，长片段数据在组装基因组中是非常有用的。已发表基因组的contigN50长度平均在50 kb左右，而本文组装得到的牵牛花基因组却高达1.87 Mb。本文中Tnp1转座子平均长度在7 kb左右，因此插入片段为7 kb的PacBio数据对于这一类重复序列区域的组装效果显著。
对基因组草图的研究能够极大帮助我们去了解牵牛花性状的遗传基础。作者在所有15条染色体中都检测到了Tpn1家族转座子的存在。TIRs和SRRs的存在意味着序列中可能会发生转座，而TpnA1，TpnA2，TpnA3和TpnA4可能是编码转座子的序列，因此作者推测这两种特征极有可能是突变个体中发生Tpn1转座子转座的诱导因素。进化分析表明，虽然牵牛花、番茄和猕猴桃都属于菊分支，但它们各有1次独立的WGD事件。
目前，牵牛花是旋花科中唯一组装到染色体水平的基因组，这不仅能促进未来对牵牛花及其相关物种的研究，还能为茄目的比较基因组学分析提供帮助。
参考文献
Hoshino A, Jayakumar V, Nitasaka E, et al. Genome sequence and analysis of the Japanese morning glory Ipomoea nil.[J]. Nature communications, 2016, 7.
