川芎基因组survey 测序及其特征分析

毛常清，沙秀芬, ，黄 静，陶 珊，彭 芳，李 群，张 超，袁 灿*

1.四川省农业科学院 经济作物育种栽培研究所，四川 成都 610300

2.四川师范大学 生命科学学院，四川 成都 610101

川芎Ligusticum chuanxiongHort.是伞形科藁本属草本植物，以干燥的根茎入药，始载于《神农本草经》，已有1500 多年的种植历史，是我国传统的大宗中药材。川芎广泛分布于我国四川彭州、什邡、眉山等地，在云南、贵州、广西、湖北、江西等省也少量引进种植。川芎性温，味辛、微苦，具有活血行气、祛风止痛的功效，其主要化学成分包含挥发油、生物碱、有机酸和多糖等，在临床上广泛地用于治疗心脏病、脑梗死及尿路结石等疾病[1-3]。

基因组序列是研究一个物种遗传背景的基础，随着高通量测技术的逐渐成熟，许多植物基因组序列相继被发表，包括大量药用植物基因组。如通过高通量测序技术成功测算出灵芝、丹参、人参、三七、天麻、穿心莲、黄花蒿、广藿香、铁皮石斛等几十种重要药用植物的基因组大小和特征[4-5]。但我国药用植物种类丰富，约占中药材资源总数的87%[6]，同大多数药用植物一样，现报道川芎的研究主要集中在化学成分[2]和药理药效机制上[7]，分子生物学方面仅开展了川芎转录组分析和利用通用引物分析其遗传多样性[8-10]，在分子遗传学系统研究上存在较大空白。虽然药用植物基因测序技术的应用为川芎全基因组测序提供了技术基础，但由于川芎基因组结构庞大，遗传背景复杂，直接进行全基因组测序存在一定困难，因此在进行全基因组测序之前，有必要对川芎基因组大小进行调研。

本研究采用流式细胞术和Illumina Hiseq 2500高通量测序技术相结合的方式对川芎基因组大小进行估算，对所得的基因组数据进行K-mer 分析、基因组预测及注释，关注阿魏酸等成分合成途径的基因注释，为进一步全基因组精细测序和药效成分合成分子机制研究提供参考依据和基因资源。

1 材料

分别采集2 个月左右的绿豆Vigna radiata(Linn.) Wilczek.、陆地棉Gossypium hirsutumL.幼嫩叶片和1 个月左右川芎幼嫩叶片（种植于四川省农科院经济作物育种栽培研究所基地）用于流式细胞分析，并采集川芎叶片用于基因组测序。绿豆由四川省农业科学院经济作物育种栽培研究所叶鹏盛研究员鉴定为豆科豇豆属植物绿豆，陆地棉由中国农业科学院棉花研究所杜雄明研究员鉴定为锦葵科棉属植物陆地棉，川芎由四川农业大学陈兴福教授鉴定为伞形科藁本属植物川芎。

2 方法

2.1 样品的制备

分别选取绿豆、棉花、川芎幼嫩叶片各2 份，每份120 mg，洗净置于预冷的培养皿中，向培养皿中加入1 mL 预冷的OttoI 细胞裂解液，快速切碎叶片后，用移液枪上下吹打混匀（避免气泡），所得的提取液用42 μm 尼龙膜滤过到离心管中，低速冷冻离心后，弃上清液，向沉淀中加入1 mL冰浴的OttoII缓冲液重悬细胞，放置4 ℃备用。

采用改良CTAB 法提取川芎基因组DNA，使用NanoDrop 2000C 超微量分光光度计和1%琼脂糖凝胶电泳检测DNA 浓度及完整性。

2.2 流式细胞术测定川芎基因组大小

将上述制备好的细胞悬浮液样品中加50 μL 1 mg/mL RNase，50 μL、50 μg/mL PI（DNA 荧光染料碘化丙啶，预先经0.22 μm 微孔滤膜滤过，-20 ℃保存），混匀，4 ℃避光染色 10 min。随后用FACSCalibur 流式细胞仪检测PI 在488 nm 激发光下发出的荧光，CellQuest 软件捕捉荧光信号数据，ModFit 软件分析结果。测定基因组大小。

待测样品DNA 量=对照DNA 量×待测样品的荧光强度/对照品的荧光强度

2.3 基因组测序和质量评估

构建270 bp 和500 bp 的小片段文库，利用Hiseq2500 测序技术对文库进行双端测序。从文库中随机取10 000 条单端read 与NCBI 数据库中的核苷酸数据库（NT）进行BLAST[11]比对，判断样本是否被外源物种污染。数据测序由北京百迈客生物科技有限公司完成。

2.4 基因组大小、重复序列和杂合率预测

选取K值为21 对基因组大小、重复序列、杂合率进行预测。用程序Jellyfish 计算K-mer 分布，并通过K-mer 分布曲线初步评估基因组重复序列含量和杂合度。

基因组大小=K-mer 总数/K-mer 期望深度值

2.5 基因组初步组装和GC-Depth 分布分析

利用SOAPdenovo 进行组装得到contig，利用双末端信息进行gap 填充，将无overlap 关系的contig 拼接组装成scaffold，获得含有N（重复序列）的初级基因组序列。过滤后的read 比对到已组装好的基因序列上，获得碱基深度，以10 kb 为窗口，在序列上无重复前进，计算每个窗口的平均深度与GC 含量，做出GC_depth 图[12]。

2.6 重复序列分析

对测序所得的数据进行重复序列分析。使用4个互补程序 LTR_FINDER[13]、MITE-Hunter[14]，RepeatScout[15]和PILER-DF[16]构建川芎重复序列文库，随后由PASTEClassifier[17]分类，并与Repbase[18]转座因子库结合起来作为最终的库，然后运行软件Repeat-Masker[19]在最终文库中找到同源重复序列。

2.7 基因预测和注释

基因从头预测，在滤过掉小于1000 bp 大小的scaffold 后，用Genscan 和Augustus 软件通过拟南芥的训练集预测川芎基因。将预测到的基因比对到非冗余蛋白序列（non-redundant，Nr）、真核生物蛋白直系同源簇（clusters of euKaryotic orthologous groups，KOG）、基因本体（gene ontology，GO）、swiss-prot 和京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes，KEGG）等数据库进行BLAST 分析来对预测基因进行注释。然后，通过KEGG 在线网站（http://www.genome.jp/kegg/）检索KEGG 途径，使用Blast2GO 软件处理得到的Nr 注释结果进行 GO 分类，使用在线网站（http://www.ncbi.nlm.nih.gov/COG/）处理KOG 注释结果进行KOG 分类。

2.8 基因家族鉴定及系统发育树构建

使用ORTHOMCL v2.0.9 将预测到的川芎蛋白序列与丹参Salvia miltiorrhizaBunge、胡萝卜Daucus carotavar.sativaHoffm.、葡萄Vitis viniferaL.和拟南芥Arabidopsis thalianaL.等4 种植物中氨基酸数目大于50 的序列汇集到一个蛋白数据库中，通过blastp 比对获得所有物种蛋白序列之间的相似性关系，E值为1×10-5，去掉序列一致度小于30%，覆盖率小于30%的序列；并对比对结果进行聚类，默认膨胀系数为1.5。在ORTHOMCL 结果中检索单拷贝基因家族，利用单拷贝同源基因基于最大似然法（maximum likelihood，ML）进行进化树构建。

3 结果与分析

3.1 流式细胞术检测基因组大小

采用绿豆和陆地棉作为内标植物，用流式细胞仪检测其混合样品的荧光强度。已知绿豆基因组大小为579 Mb[20]，陆地棉基因组为2173 Mb[21]，通过2 个内标植物计算出川芎的基因组大小分别为525/78×579=3 897.12 Mb，变异系数（coefficient of variation，CV）为3.40%（图1）；613/439×2173=3 034.28 Mb，CV 为2.01%（图2）。

图1 绿豆与川芎混合样品流式细胞术测定Fig.1 FCM determination for mixed samples of V.radiata and L.chuanxiong

图2 陆地棉花与川芎混合样品流式细胞术测定Fig.2 FCM determination for mixed samples of G.hirsutum and L.chuanxiong

3.2 川芎基因组测序数据统计及质量评估

使用川芎基因组DNA 构建270 bp 的文库，通过Illumina Hiseq2500 测序平台测序并过滤得到222.04 Gb 高质量的数据，测序深度为72.59 X，测序数据Q20 比例均在97.59%以上，Q30 比例均在94.74%以上。随机筛选的1000 条单端read 能够比对上NT 核酸数据库的read 占总read 的7.62%，其中比对上野胡萝Daucus carotaL.、细叶藁本Ligusticum tenuissimum(Nakai) Kitagawa 的read 数分别占比对上NT 库reads数的61.81%、3.01%，且未发现动物、微生物等异常比对，说明样本不存在污染。

3.3 川芎基因组特征

通过270 bp 文库数据构建K=21 的K-mer 分布图（图3），K-mer 深度62 X 为主峰（由于杂合度较高，本研究对主峰判断参考流式细胞结果），测序得到K-mer 的总数为189 618 848 178，估算出川芎基因组大小为3 058.37 Mb，与流式细胞实验中以陆地棉为内标植物测定的结果相近。在主峰对应深度的1/2 处出现明显的杂合峰，深度为31 X，说明川芎基因组具有较高的杂合度。进一步根据变异数目占基因组大小的比例即杂合度估算该基因组杂合率，在所有有效数据中检测到每139 个碱基对中就有3 个SNP，估算出川芎基因组具有较高杂合率，约为2.16%。K-mer 分布存在较长拖尾，暗示川芎基因组存在较高的重复率，估算重复序列的基因组大小估计为 2 440.13 Mb，约为川芎基因组的79.79%。说明川芎属于高重复、高杂合、大基因组等基因组特征的复杂物种。

图3 K-mer 分布估算基因组大小Fig.3 K-mer distribution to estimate genome size

3.4 川芎基因组组装和GC 含量分析

使用222.04 Gb 高质量数据，基于K=41 组装产生12 973 787 个contig 和8 119 089 个scaffold。contig N50 为286 bp，N90 为131 bp，最大长度达到26 842 bp，总长度为3 198 598 874 bp；scaffold N50 为493 bp，N90 为191 bp，最大长度为29 779 bp，总长度为3 284 536 081 bp。其中，contig 和scaffold 的N50 值相对较短，可能是由于川芎高杂合率引起的。通过对组装的contig 进行GC 含量的统计，结果显示川芎基因组的GC 含量约为36.32%（图4），说明测序不具有明显的GC 偏向性，不影响测序分析的准确性。

图4 川芎基因组GC 含量和平均测序深度Fig.4 GC content and average sequencing depth of L.chuanxiong genome

3.5 川芎基因组重复序列分析

重复序列检测显示其总长度为2 250 901 770 bp，约为基因组大小的73.60%，低于K-mer 分析估算的重复序列含量，原因可能是组装效果的限制，导致组装过程中重复序列损失6.19%。注释上，能够找到明确重复序列元件的总长度约为2 026.25 Mb。其中，长末端重复序列（long terminal repeated，LTR）是最丰富的重复元件，占基因组的14.14%，其次是长散在重复元件（long interspersed nuclear elements，LINE），占基因组的0.49%。SSR 重复序列总长度约49.41 Mb，占基因组的1.62%，占重复序列的2.44%。

3.6 基因预测和注释

总共预测到34 864 个基因，总长度为27 532 625 bp。在预测的基因中，查找到79 130 个外显子，总长度为18 557 406 bp；79 129 个内含子，总长度为8 975 219 bp。有30 737 个基因在功能数据库中能比对到注释信息，Nr 具有最高的注释率（30 584，87.72%）KEGG 具有最低的注释率（9623，27.60%）。在预测的基因中，15 184 个基因被分类为GO 功能类别；15 598 个基因被分类为KOG 功能类别；9623个基因注释到125 个KEGG 代谢途径。在GO 功能类别中，包含分子功能、细胞组成和生物过程（图5）；在KOG 功能类别中，通用功能预测的基因最多，其次是翻译后修饰、蛋白质周转、分子伴侣和信号转导机制（图6）；在参与基因数目最多的前10条KEGG 代谢途径中，注释基因分别参与核糖体代谢、植物激素信号转导、内质网蛋白质、剪接体、碳代谢、氨基酸生物合成、RNA 转运、氧化磷酸化、植物-病原菌互作、淀粉和蔗糖代谢等途径。

图5 川芎基因的GO 注释Fig.5 GO annotations of L.chuanxiong genes

图6 川芎基因功能注释KOG 功能分类Fig.6 KOG functional classification of L.chuanxiong genes

3.7 基因家族鉴定及系统发育分析

通过与胡萝卜、丹参、葡萄和拟南芥等物种蛋白序列比对（图7-A、B），川芎中常见的基因家族中的基因数量小于同科的胡萝卜，每个基因家族的平均基因数量与其他物种相当，但川芎独特基因家族的数量比其他物种的独特基因家族中的基因数量要大得多，共计2058 个基因家族。所有5 个物种共有的基因家族为7112 个，其中2001 个基因是单拷贝基因，即每个基因家族中只存在一个直系同源基因，可用于系统发育推断和发散时间估计。对2001 个单拷贝基因利用ML 构建系统发育树（图7-C）。从川芎系统发育树分析，川芎、胡萝卜、丹参、葡萄和拟南芥来自于共同祖先。其中，川芎与胡萝卜最晚与其他物种发生分歧，二者进化分支长度差异最小，分歧时间更短，亲缘关系更近。从遗传变异度上来看，胡萝卜所在的分支最短，遗传变异度最小，进化距离最近，川芎的遗传变异度和进化距离仅次于胡萝卜，丹参的遗传变异度最高，进化距离最远。

图7 川芎基因家族鉴定及系统发育树Fig.7 Gene family identification and phylogenetic tree of L.chuanxiong

3.8 参与阿魏酸生物合成的基因

阿魏酸是川芎的主要药效成分，属于苯丙烷生物合成途径，是木质素合成的中间体。近年来研究者证实了在川芎、当归等伞形科物种中，阿魏酸的合成途径主要包括COMT 途径和CCoAMT 途径[9,22-24]，主要参与的酶有苯丙氨酸解氨酶、肉桂酸-4-羟化酶、香豆酸-3-羟化酶、咖啡酸-O-甲基转移酶、4-香豆酸辅酶A 连接酶、莽草酸O-羟基肉桂酰转移酶、奎宁O-羟基肉桂酰转移酶、咖啡酰辅酶A-O-甲基转移酶、肉桂酰辅酶A 还原酶、醛脱氢酶家族成员C4。为挖掘川芎中参与阿魏酸生物合成的主要基因，本研究从注释的KEGG 代谢通路中提取苯丙烷生物合成参考途径（map00940）相关数据，对阿魏酸合成过程的相关基因进行了检索，共有53 个功能基因覆盖到阿魏酸合成途径中（表1）。其中，编码咖啡酸-O-甲基转移酶、4-香豆酸辅酶A 连接酶、莽草酸O-羟基肉桂酰转移酶、咖啡酰辅酶A-O-甲基转移酶等酶的基因均有多个同源拷贝存在，预测川芎基因组在进化的过程中某一时间点发生过基因扩张。

4 讨论

基因组包含了一个生物体所有基因的总和，了解生物体基因组信息有助于深入了解生物体遗传、进化、生物合成、次生代谢等全部过程。通过基因组测序技术可以对特定物种基因组进行测序，利用生物信息学方法对测序序列进行拼接和组装，最终获得该物种基因组序列，进而了解其基因组信息[25]。2009 年，陈士林团队[6]首次提出本草基因组计划，此后，越来越多的学者开始在药用植物基因组学研究上投入大量精力。药用植物全基因组水平研究，有助于阐明药用植物活性成分生物合成和代谢调控途径之间的关系，为具有相似有效成分和药理活性的近缘物种间的系统发育关系研究奠定了基础，也为药用植物的遗传育种和基因资源保护提供了重要依据[26]。

本研究利用流式细胞术和高通量测序相结合的方式对川芎基因组进行测算。通过流式细胞术估算出川芎基因组大小分别为3 897.12 Mb和3 034.28 Mb，通过高通量测序的K-mer 分析，综合2 个分析结果估算川芎基因组大小为3 058.37 Mb，属于基因组较大的物种。本研究测得川芎基因组GC 含量为36.32%，处于植物基因组GC 含量应介于25%～65%的合理范围[27]，说明川芎基因组测序的结果和组装是正确可靠的。通过K-mer 分析，估算出川芎重复序列含量为79.79%，杂合率约为2.16%，与地黄[28]、黄芪[29]等药用植物类似，呈现高重复、高杂合的基因组特征，进一步说明川芎基因属于高重复、高杂合、基因组庞大的物种。

在基因预测、注释和基因家族的鉴定中，本研究共预测到川芎编码蛋白基因34 864 个，远高于其伞形科亲缘关系较近的胡萝卜（32 113 个）[30]，可能是因为组装的都是短片段测序文库，川芎中的基因数量可能被高估了。此外，本研究中测序完成后从头组装中产生的contig N50 为286 bp，scaffold N50 为493 bp，明显较预期短，这与广藿香[4]、罗汉果[12]等药用植物的全基因组调研结果一致。提示对于川芎这种具有复杂基因组特征的物种来说，利用二代高通量测序对其全基因组进行精确测序仍然存在技术难度。因此，提高川芎基因组的测序深度和组装质量，建议后续的研究可采用二代和三代测序技术相结合，并利用全基因组染色体构象捕获技术（high-through chromosome conformation capture，Hi-C），解析全基因组范围内整个染色质DNA 在空间位置上的关系，获得完整准确的全基因组图谱[31]，得到高质量的川芎基因组序列。

本研究首次利用流式细胞技术和基因组survey分析，初步获得川芎基因组大小和结构特征，即基因组庞大、序列重复率高、序列杂合度高，为下一步进行全基因组精细测序奠定基础。本研究中组装产生的大量川芎基因组序列和注释基因为后续分子标记的开发和基因功能研究提供了大量的资源。同时，本研究挖掘了阿魏酸合成途径中的参与基因，对川芎阿魏酸生物合成途径潜在分子机制的初步研究，为进一步研究川芎生物学和选育具有优良药用性状的品种奠定了基础。

利益冲突所有作者均声明不存在利益冲突