单叶省藤和黄藤基因组组装研究

王思宁 赵韩生 高志民

(国际竹藤中心 北京 100102 )

棕榈藤是世界上最重要的非木材森林植物资源之一，属于棕榈科(Arecaceae)，是自然进化中出现的攀爬植物的代表。近年来的研究表明，在棕榈科省藤族省藤亚科中棕榈藤植物可分为11属，共有631种，其中省藤属(NCBITaxonID:4711)和黄藤属(NCBITaxonID:93268)种类最多，分别占藤种的约65%和20%。这2个属也是最重要的棕榈藤藤材来源，为工业生产提供了95%以上的藤条。全球有超过500万人在经济上依赖棕榈藤，藤工业产值每年约70亿美元，包括国内工业生产、国际贸易，以及藤编织的篮子、座椅和家具等。同时，对棕榈藤遗传育种的关注在持续增长，并且预计将来在几年内棕榈藤的种植面积将逐渐超过天然棕榈藤面积。

单叶省藤(Calamussimplicifolius)(NCBITaxonID:746888)是中国原产、种植面积最广的藤种之一，茎长约50m，直径约15mm。单叶省藤是海南岛特有的棕榈藤，可生产高品质、中等直径的藤条，用于捆扎和编织。黄藤(Daemonoropsjenkinsiana)(NCBITaxonID:1510057)，是常绿攀援棕榈藤的代表，自然生长在海拔1000m以下的低地雨林，分布于孟加拉国、不丹、柬埔寨、印度、老挝、缅甸、尼泊尔、泰国、越南和中国的东南部。黄藤丛生生长，茎可达50m, 节间长40cm，直径30mm。这2种产值最高的藤种现种植面积超过1000hm2，栽培在中国小于北纬23°30′的地区，例如海南岛、广东、广西、云南、福建等华南地区。

单叶省藤和黄藤由于其藤条韧性强和耐久度好而具有广泛的用途和巨大的开发潜力。需要应用分子育种技术开展棕榈藤育种，来满足对棕榈藤质量和数量不断增长的需求。然而，棕榈藤遗传背景不清，严重阻碍了对其分子生物学的研究和实际生产的应用，更影响了相关物种之间比较基因组分析的深入开展。因此，我们使用最新的高通量测序技术，完成了单叶省藤和黄藤的测序和染色体水平基因组组装。随着这2个棕榈藤基因组的完成，使许多比较基因组分析和其他下游分析成为可能，包括分子标记的开发、功能基因的鉴定和分子辅助育种等。此外，高质量的基因组将有助于对棕榈藤生物性状进行基因组学、转录组学和代谢组学分析。本文选取并鉴定了棕榈藤中改良藤条性能的一些重要候选基因——木质素生物合成基因家族，为今后利用木质素生物合成基因来改良藤条的品质以及研究物种多样性奠定了基础。

1 材料方法与数据处理

1.1 DNA分离、文库构建和高通量测序

单叶省藤和黄藤材料取自中国林业科学研究院热带林业研究所(北纬23°11′29″，东经113°22′40″，海拔87m)，于2015年春季采集处于营养生长阶段的幼嫩叶片组织用于总DNA提取。使用DNeasyPlantMiniKits(Qiagen)分离和提取总DNA。根据高分子量核DNA分离法纯化基因组DNA获得多个DNA文库，并在IlluminaHiSeq4000和PacBioSequel平台上进行测序(原文表1)。即：在基于Illumina的方法上构建了3个不同插入片段大小的文库(270bp、500bp和800bp)，用于PE测序，以及4个不同插入片段大小的文库(2kb、5kb、10kb和20kb)用于MP测序。另外，按照PacBio方法，构建了5个插入片段大小为20kb的PacBioSequel文库，进行测序。通过对测序数据清理和预处理，最终得到约494.08Gb的单叶省藤高质量数据(322.3GbPEreads,93.4GbMPreads,78.38GbPacBiodata)，覆盖率为单叶省藤基因组的252倍， 约426.17Gb的黄藤高质量数据(244.58GbPEreads,103.21GbMPreads,and78.38GbPacBiodata)，覆盖率为黄藤基因组的266倍。

作为与Illumina和PacBio文库构建测序，进行了另一平行分析，即使用相同组织材料为单叶省藤和黄藤建立了2个高通量测序文库。使用MboI限制性内切酶在甲醛固定后消化DNA，然后用生物素标记5个残基。在原位连接平末端片段后，对分离出的DNA进行反向交联、纯化和过滤，得到包含生物素标记的片段。随后，按此顺序进行DNA片段端修复、适配体连接和聚合酶链式反应。然后，根据BGISEQ-500的使用说明，对其进行100PE的测序。在分别获得了大约148Gb和154Gb的原始数据后，经过低质量数据的过滤和其他预处理，最终获得约6.7Gb和13.1Gb的高质量数据，并使用HiC-Pro(version2.8.0devel)软件分别对单叶省藤和黄藤进行了评价和分析(原文表1)。

1.2 基因组调查

对给定物种基因组大小、杂合性等基因组特征的了解解，将有助于开发一套适合的测序方案和组装策略。因此，本研究采用4种独立的方法来估算基因组大小，即KmerSpectrumPlot.plinALLPATHS-LG(版本r52488)、GCE(GenomeCharacteristicsEstimation,released7Jan.2015)、JELLYFISH(版本2.0)和流式细胞术(原文表S1、S2，图S1、S2)。在调研中，我们测序获得单叶省藤约为98Gb，黄藤约为60Gb的二代原始测序数据作为调研数据。经过数据预处理和组装，结果显示(原文表S1)，最终预测单叶省藤和黄藤基因组大小分别约为1.98Gb和1.61Gb，相关杂合度分别为1.32%～1.52%和1.19%～1.31%。因此，基因组分析结果表明，这2种棕榈藤的基因组比较适用于Illumina和PacBio测序数据的联合组装。

1.3 利用Illumina和PacBio高通量测序数据联合组装基因组

在Illumina数据的预处理过程中，过滤掉了低质量的读长和接头序列，获得单叶省藤和黄藤的有效数据分别为416Gb和348Gb。对于PacBio数据，使用了MECAT(发布于2017年6月27日)软件，得到单叶省藤和黄藤分别为52Gb和32Gb的高质量数据。随后，使用FALCON(版本0.3)软件进行第1次contig组装。如原文表S3所示，使用不同参数为单叶省藤基因组生成了2个组装结果，一个ContigN50为67.2kb，基因组约为1.59Gb(约为预测基因组大小的80%)，另一个ContigN50为66.7kb，基因组约为1.53Gb(约为预测基因组大小的77%)。此外，获得一个黄藤组装结果ContigN50为81.5kb，基因组约为1.27Gb(约为预测基因组大小的79%)。通过以上结果显示，仅用三代数据和MECAT软件不能获得对2种高质量棕榈藤基因组，可能是由于2种棕榈藤的高杂度(单叶省藤1.32%～1.52%，黄藤1.19%～1.31%)、高比例重复序列(单叶省藤为54.15%，黄藤为70%；有关详细信息请参阅原文中的后续分析)和使用了测序深度不足的三代数据(PacBio数据经错误校正后分别为26×和20×)造成的。因此结合以上结果，我们尝试利用Illumina和PacBio测序数据，对单叶省藤和黄藤进行了联合基因组组装。首先，利用Platanus(版本1.2.4，是一种用于高杂合数据的全新基因组汇编器，可以通过构建具有自动优化的k-mer大小的DeBruijn图)将PE片段整合到contig序列中。其次，用DBG2OLC(2015年7月11日发布，参照以下参数：DBG2OLCContigscontig.faLD0K17KmerCovTh4MinOverlap25AdaptiveTh0.007RemoveChimera1fscaffold.fa.) 矫正PacBio测序数据并组装scaffolds，以分别获得单叶省藤和黄藤约1.92Gb和1.56Gb的初始组装结果。再次，在进行SSPACE分析前，参照DBG2OLC(原文表S4)进行了平行分析，这一步骤有助于提高基因组组装的质量，减少SSPACE过程中的错误。然后，使用SSPACE(版本3.0)和MP大片度测序数据延长scaffolds，并通过使用GapCloser(版本1.12)和PBJelly(2015年8月24日发布) 进行补洞。最终，获得了1.96Gb的单叶省藤组装数据，包含5116个Contigs，ContigN50长度为107kb，ScaffoldN50长度为803kb，同时获得了一个1.60Gb的黄藤组装数据，包含ContigN50和ScaffoldN50长度分别为108kb和784kb(原文表2)。

随后，通过3D-DNA分析流程(版本170123)对有效的高通量测序数据与上述装配数据一起进行处理，得到了一个合理的连接模式，准确的染色体组装结果。如原文图2所示，使用Juicebox(版本1.5.2)显示连接模式图谱。原文表S5展示了单叶省藤12个最长的染色体水平Scaffold和黄藤13个最长的染色体水平Scaffold。预测染色体总长度分别占单叶省藤和黄藤基因组总长度的92.08%和92.01%，ScaffoldN50分别为169Mb和119Mb。

1.4 基因组评估

采用3种独立的方法评估了单叶省藤和黄藤基因组组装的准确性和完整性。首先，对基因组中不确定碱基的百分比(Ns)和鸟嘌呤、胞嘧啶含量(GC)进行了分析。结果显示在每个基因组中低丰度的Ns(单叶省藤约为0.6%和黄藤约为0.7%)和总GC含量(单叶省藤为41.78%和黄藤为41.07%)与相关转录组数据相似(单叶省藤为41.68%和黄藤为41.89%)。此外，利用BasicLocalAlignmentSearchTool(BLAST)-like比对工具(版本1.0)(默认参数)进行RNA测序(RNAseq)数据所组装出来的功能基因序列与组装得到的基因组序列比对，结果显示，在一个大片段上90%以上的序列可以与基因组的组装序列匹配(单叶省藤为92.89%和黄藤为81.81%)(原文表S6)。最后，利用BUSCO(版本3.0)对2种棕榈藤基因组的完整性进行了评估，这是一种从进化的角度定量对单拷贝直系同源基因含量进行计算，以评估基因组完整性的方法。BUSCO结果显示，在单叶省藤基因组组装结果中检测到保守性为96.4%的BUSCO蛋白(基于embryophytaodb9库文件)(包括3.8%的BUSCO蛋白片段)。同时，在黄藤基因组组装结果中检测到保守性为87.3%的BUSCO蛋白(包括4.0%的BUSCO蛋白片段)(原文表S7)。

1.5 重复注释

本研究在蛋白质编码基因模型预测之前，完成了单叶省藤和黄藤基因组中转座因子(TEs)和串联重复序列的鉴定。采用2种独立的方法来预测序列重复：同源注释和从头预测(denovo)法。同源注释结果显示，利用RepeatMasker(版本4.0.5)和RepeatProteinMasker(版本4.0.5)对Repbase文库(发布于2017年12月01日)进行搜索，查找转座因子序列。从头预测结果显示，利用Repeat-Modeler(版本1.0.8)和LTRFINDER方法在去除污染和多拷贝基因后，可以构建一个重复序列文库。然后，使用RepeatMasker对基因组序列进行了分类，并与重复序列文库进行了对比。此外，串联重复序列由TandemRepeatFinder(版本4.09)识别，其参数为：Match=2,Mismatch=7,Delta=7,PM=80,PI=10,Minscore=50andMaxPeriod=2000。最终结果表明，长末端重复序列(longterminalrepeat,LTR)是最丰富的重复序列类型，2种非LTR逆转座子(短穿插核元件和长穿插核元件)在2种棕榈藤基因组中占比例最少(原文表S8)。TEs在单叶省藤和黄藤基因组中分别占54.15%和70%。TEs序列分歧度结果表明，从头预测的重复序列比Repbase预测的重复序列更多(原文图3)。

1.6 RNA样本收集、文库构建和转录组组装

分别选取单叶省藤和黄藤3个不同发育阶段的纤鞭样本，每个样本3个生物学重复(原文表S9)。为保证与2种棕榈藤基因组研究的一致性，RNA取样的位置和DNA取样的位置一致。利用RNA提取试剂盒(TaKaRa, 日本)提取各样本的总RNA，并使用NanoDrop2000分光光度计测定RNA的质量和浓度。使用RecombinantDNaseI于37℃孵育30min去除残余DNA，并使用反转录试剂盒(Promega, 美国)合成cDNA第一条链，操作参照相关试剂盒说明书进行。然后，使用BGISEQ-500平台对文库池中的短100PE片段测序。在对转录组数据进行预处理时，使用SOAPnuke(版本1.5.6)对匹配序列和低质量的读长序列进行过滤，过滤参数为：-n0.001-l20-q0.4-q2。所有样本均使用Trinity(版本2.0.6)按以下参数组装：(1)grouppairsdistance500, (2)mincontiglength200, (3)minkmercov2, (4)minglue2, (5)bflyopts-V5, (6)edge-thr=0.1, (7)stderr,and(8)SSlibtypeRF。然后，对Trinity的输出结果进行聚类分析，使用TGI聚类工具(版本v2.0.6)生成一组非冗余集合，参数如下：(1)aminimumof95%identitybetweenthecontigs, (2)aminimumof35overlappingbases, (3)aminimumscoreof35,and(4)amaximumof20unmatchedoverhangingbasesatthesequenceends。最后，根据序列相似度将整合的转录本分为2类：成簇转录本和单一转录本。在成簇转录本中，转录本之间的序列相似区都超过了70%，由这些转录本被剪接成同源基因或旁系同源基因。此外，所有鉴定的基因将被用于后续的功能分析。

1.7 基因建模与预测

使用3种独立的方法对完整的蛋白质编码基因模型进行了综合预测，包括从头预测、同源蛋白预测和转录组数据预测。使用AUGUSTUS(版本3.3)的默认参数对重复区段进行第一次注释，该程序是一个基于自我调试模型的从头预测软件。通过进行多次调试优化数据，在单叶省藤和黄藤中分别预测获得了85246和87613个基因模型。在基于同源性的预测中，以油棕(Elaeisguineensis)、长叶刺葵(Phoenixdactylifera)、二穗短柄草(Brachypodiumdistachyon)、水稻(Oryzasativa)、粟(Setariaitalica)、高粱(Sorghumbicolor)和玉米(Zeamays)等7个物种作为参考基因组。从ENSEMBL数据库下载它们的蛋白序列并应用TBLASTN(版本2.2.26)对单叶省藤和黄藤中E-value不小于1e-5的组装序列进行比对。然后，由GeneWise(版本2.0)生成拼接结果。转录组数据预测结果显示，使用HISAT2(版本2.0.2)来识别外显子-内含子剪接位置，并进一步进行转录组序列与基因组序列的比对。使用Cufflinks(版本2.2.1)分别在单叶省藤和黄藤中分别鉴定了56024和58134个蛋白质编码基因模型(原文表S10)。最后，使用MAKER(版本2)整合了上述3个独立预测的数据结果，随着相同基因序列在单叶省藤和黄藤基因组中的出现，最终共分别鉴定出51235和53342个完整的蛋白编码基因模型。

1.8 基因注释和基因功能预测

使用2种独立的方法对预测的基因集合进行评估：基因功能评估和BUSCO完整性评估。在基因功能评估中，比照基因在近缘物种中同源蛋白序列和预测的蛋白集合，评估了预测的注释和蛋白质序列的一致性。参照NCBI非冗余蛋白质数据库(2018年3月13日发布)、SWISS-PROT(2018年1月1日发布)、GeneOntology(GO)(2013年10月30日发布)、KyotoEncyclopediaofGenesandGenomes(KEGG)(数据集v81)和InterPro(数据集v.53)等5个权威蛋白数据库进行注释比对，结果表明(原文表S11)在单叶省藤和黄藤预测的基因模型中分别有5.34%和2.89%的基因为未注释的基因。此外，BUSCO的评估显示，在单叶省藤和黄藤中分别有88.7%和91.3%的保守BUSCO蛋白(embryophytaodb9)存在。在这些保守的BUSCO蛋白中，分别有76.2%和81.2%是完整的。此外，tRNA，rRNA，miRNA和snRNA4种非编码RNA也被预测出来(原文表S12)。

1.9 基因家族结构和棕榈藤特异性基因家族分析

本研究中，通过逐一序列比对的方法来预测基因组水平的直系同源基因。这种方法因速度快而常用于处理大数据量分析。采用BLAST-based的OrthoMCL(版本2.0.9)来识别单叶省藤和黄藤中E-value≤1e-5且匹配系数≥80的同源基因。Markovchain聚类(默认膨胀参数)也共同用于与其他8个物种的全面BLASTP分析，包括无油樟(Amborellatrichopoda)、油棕、拟南芥、二穗短柄草、水稻、紫萍(Spirodelapolyrhiza)、长叶刺葵和高粱。在所有10个物种中鉴定得到30936个基因家族，其中单叶省藤和黄藤基因组中分别检测到44700和44537个同源基因。10个物种中共有的基因家族近6132个(占19.8%)，在单叶省藤和黄藤基因组中别检测到的2366和2707个特有的基因家族(原文图4b)。此外，家族分析结果显示其中637个基因家族是棕榈藤所特有的。这些特异的基因家族大多富集在膜组成成分、转录因子活性相关的基因簇上(原文表S13)，且聚集在植物与病原体相互作用、植物激素信号转导相关的KEGG通路中(原文表S14)。

1.10 系统发育进化和分化时间分析

利用获得的962个种间高度保守的单拷贝直系同源基因，进行棕榈藤和其他物种的进化关系分析。首先，通过MUSCLE(版本3.8.31)进行蛋白质序列的多序列比对；在此基础上，基于蛋白质序列进行了DNA编码序列(CDS)比对。随后，将所有比对的CDS序列连接起来，使用Perl脚本为每个物种生成一个超级基因。然后，提取了每个密码子第2个位置上(相位1)的核苷酸使用RAxML(版本8.2.3)构建系统发育树，模型为“GTRGAMMA”。结果表明，4种棕榈科植物分布在一支上，由2个独立的子分支组成，其中一个分支上包含单叶省藤和黄藤，另一个分支上包含油棕和长叶刺葵(原文图4a)。

此外，使用了PAML(版本4.5)的MCMCTree程序并按照参数“-nsample200000-burnin40000”估计单叶省藤和黄藤与其他8个物种的分化时间。时间标准是根据参考物种发生分化的时间推算得出的。结果表明，单叶省藤和黄藤之间的分化时间约在1930万年前，而棕榈科的油棕和长叶刺葵的分化时间约在4080万年前(原文图4c)。

1.11 木质素生物合成通路基因家族全基因组鉴定分析

木质素是一种由木质素单体组成的复杂芳香族聚合物，与纤维素和半纤维素基相互作用共同构成次生细胞壁。木质素通常由3种木质素单体对-羟基苯基(p-hydroxyphenyl，H)、香草醛(vanillin，G)和丁香醛(syringaldehyde，S)组成。本研究利用拟南芥、二穗短柄草、水稻、高梁、毛竹(Phyllostachysedulis)、毛果杨(Populustrichocarpa)、单叶省藤和黄藤等8个物种的基因组进行木质素生物合成通路13个基因家族的全基因鉴定。多数基因组序列(拟南芥、二穗短柄草、水稻、高梁和毛果杨)从ENSEMBL数据库下载获得，从毛竹基因组数据库(BambooGenomeDatabase)中下载了毛竹的基因组序列。在前人研究的基础上，收集得到了经过实验验证的140个木质素生物合成通路基因(原文表S15)，这些基因作为查询序列进行进一步的基因功能鉴定。如上所述，在全基因组基因鉴定过程中使用了BLAST比对和保守结构域分析方法。简单地说，使用已知基因的编码序列，对所有基因组序列(包括两种棕榈藤的基因组序列)进行了目标蛋白的BLAST比对(版本2.2.26)，选取E-value<1e-10，相似度>40%且序列的覆盖率>95%的序列。随后，利用Pfam-A.hmm数据库(2017年3月31日发布)对已过滤的序列进行hmmsearch(版本3.1b2)分析，通过人工校正的方法去除结果中结构域不完整的序列。结果表明，单叶省藤和黄藤中木质素相关基因数量均有明显增多(原文表3)。单叶省藤和黄藤中每个木质素生物合成通路基因家族平均成员数量分别约为15和13个，总数分别为193和172个。在2种棕榈藤中过氧化物酶基因家族是最常见的基因，单叶省藤的苯丙氨酸解氨酶基因和黄藤中的香豆酸3-羟化酶基因与肉桂酸4-羟化酶基因是最不常见的。棕榈藤中木质素生物合成基因的数目扩增可能是由于全基因组复制(WGD)事件造成的，因为全基因组复制可以提供更多的基因拷贝，从而促进具有新功能基因的进化。

2 结论

本研究采用多种类型的测序数据，首次完成了单叶省藤和黄藤染色体水平基因组的组装，2种棕榈藤的基因组作为重要的参考基因组，将有助于促进其他藤种的从头基因组测序组装和重测序研究，同时通过与不同物种进行比较研究，为物种进化研究提供证据。借助两种棕榈藤高质量基因组数据，使木质素生物合成通路关键基因的鉴定更加便捷，这些候选基因对棕榈藤的生长发育非常重要。本研究为进一步对棕榈藤及相关物种基因组的研究奠定了基础。