全基因组重测序分析4个毛竹变型的秆形和秆色变异

牟少华，李雪平，李 娟，高 健

(国际竹藤中心，国家林业和草原局/北京市共建竹藤科学与技术重点实验室，北京 100102)

全基因组重测序是对已完成基因组序列的同一物种或者近缘种的其他个体进行全基因组测序，并在此基础上对个体或群体进行差异性分析。新个体基因组序列通过映射参考基因组获得同源数据集，然后检测基因组范围内的单核苷酸多态性(Single nucleotide polymorphism，SNP)、插入缺失(Insertion/deletion，InDel)及结构变异(Structural variation，SV)等多态位点[1]。运用高通量的重测序技术进行遗传变异分析已经在谷子[2]、大豆[3]、水稻[4]、花生[5]、鼠尾草[6]、葡萄[7]、木兰[8]、烟草[9]、竹黄菌[10]和蘑菇[11]等植物上得到了广泛应用。

毛竹(Phyllostachysedulis)是我国特有的、利用历史悠久的重要经济竹种，也是栽培面积最广的竹种，其生长速度快，材性优良，繁殖能力强，一次造林可以永续利用。毛竹笋是一种全天然营养食品，含有丰富的蛋白质、脂肪、糖类、胡萝卜素、维生素、钙、磷、铁以及18种氨基酸等成分。因此，毛竹是优良的笋材两用竹种，有着巨大的开发和应用潜力[12]。毛竹广泛分布于我国南方16个省区，跨越北、中、南3种亚热带气候，在长期适应不同的生长环境和栽培经营措施下，加上自身发生遗传漂移等因素，发生了一系列的形态变异。目前，在天然毛竹林中已经发现21个毛竹变型[13]，这些变型在形态上表现出丰富的多态性，尤其是秆形、秆色等形态性状差异显著。为研究这些形态差异究竟是因为基因变异而导致的可遗传变异还是由于生境不同而产生的表型可塑性，非常有必要在基因水平上揭示毛竹的遗传多样性和变异程度。

近年来，随着分子生物学技术的突飞猛进以及组学技术在竹类植物上的应用，特别是毛竹全基因组测序工作的完成[14]，竹类植物基因组研究进入后基因组时代。目前，毛竹已得到一个完整的可变剪接图谱[15]，对P450s[16]、AP2/ERF[17]、NAC[18]、SAUR[19]、AQP[20]、IQD[21]、SBP-like[22]、WRKY[23]、HD-Zip[24]和TCP[25]、Hsp[26]、CO-Like[27]等基因家族已采用生信方法进行了鉴定和功能验证，但尚未发现有关毛竹变型基因组序列变异的研究。本研究采用成熟且经济的二代测序技术，以毛竹全基因组为参考，对4个有代表性的毛竹变型进行重测序，对其SNP、Indel、结构变异SV等进行深度挖掘，从而揭示其全基因组的突变类型，分析其广适性状的变异基因，为从全基因组尺度分析毛竹的“种性”和变型种质“品性”获得海量基因组序列数据。

1 材料和方法

1.1 试验材料

试验材料为毛竹的4个变型(表1)：圣音竹、黄槽毛竹、黄皮花毛竹和厚壁毛竹，编号分别为R01、R02、R03和R04，均采自国际竹藤中心安徽太平试验中心。选取刚长出且尚未展开的嫩叶，液氮速冻后超低温保存。

表1 4个毛竹变型形态特征Tab.1 Morphological characteristics of 4 variations of Moso bamboo

1.2 试验方法

1.2.1 DNA提取与基因组测序 基因组DNA采用改良CTAB法提取[28]。将检测合格的样品基因组DNA，用超声波方法将DNA进行片段化，然后进行DNA片段的纯化、末端修复、在3′端加A、连测序接头，进行琼脂糖凝胶电泳完成片段大小选择，通过PCR扩增建成测序文库，进行文库的质检，合格文库采用Illunima Hiseq 2500平台进行测序，得到原始测序序列(Raw reads)，对其进行数据过滤，去掉带接头的reads，过滤N含量高(超过10%)和质量值低(低于10的碱基数超过50%)的reads，得到Clean reads，用作后续的信息分析。

1.2.2 与参考基因组进行比对统计 将重测序得到的测序reads重新定位到毛竹参考基因组上，进行后续的变异分析[14]。运用Bwa软件将测序获得的短序列与毛竹的参考基因组(毛竹基因组数据库 BambooGDB下载网址：http://www.bamboogdb.org/page/download.jsp)进行比对，通过Samtools flagstat命令对Clean reads比对定位到参考基因组上的位置，然后统计4个样品的测序深度和基因组覆盖度等方面的信息，进行序列变异的检测。

1.2.3 检测基因组变异 使用GATK软件工具包对样品的SNP和InDel进行检测[29]。根据参考基因组上Clean reads的定位结果，使用Picard软件去重复、GATK软件进行局部重比对和碱基质量值的校正等预处理，然后进行SNP和InDel的变异检测。对SNP进行严格的过滤，包括：SNP cluster过滤(5 bp内若有2个SNP过滤掉)，InDel附近5 bp内的SNP过滤掉；相邻2个InDel距离小于10 bp过滤掉，从而获得SNP位点集。使用BreakDancer软件检测SV，主要是基于序列的比对结果，获得测序文库的插入片段大小、方差，再查找序列和参考基因组间比对的异常结果，寻找可能的结构变异。

1.2.4 注释SNP、InDel和SV 通过SnpEff软件进行SNP、InDel和SV的注释[30]。根据检测获得的变异位点比对到参考基因组的位置基因、CDS位置等，获得变异位点在基因组的发生区域(基因区、基因间区或CDS区等)，以及产生的同义非同义突变等。

1.2.5 挖掘变异基因和功能注释 通过寻找样品间与参考基因组发生非同义突变SNP基因、CDS区发生的InDel基因与SV基因，挖掘可能存在的功能变异基因。通过Blast将变异基因与NR[31]、GO[32]、COG[33]、KEGG[34]等功能数据库进行比对，获得基因的注释，以便分析基因的功能。

2 结果与分析

2.1 与参考基因组比对统计

4个竹种通过高通量测序(Illunima Hiseq 2500等测序平台)得到原始图像数据文件，经碱基识别分析转化为原始测序序列(Raw reads)，然后去除带接头的reads，过滤得到Clean reads，分别为180 506 436，184 142 480，213 567 146，188 998 942 bp，定位到毛竹参考基因组的占所有Clean reads数的百分比分别为96.13%，96.32%，96.33%和97.44%，双端均定位到毛竹参考基因组上并且距离符合测序片段的长度分布占所有Clean reads数的百分比分别为88.18%，88.42%，88.57%和89.65%。样品平均覆盖深度在10×以上，各深度对应的基因组覆盖比例如表2所示。

表2 4个毛竹变型的测序数据统计Tab.2 Summary of sequencing data of 4 variants of Moso bamboo

2.2 SNP的检测与注释

2.2.1 SNP检测 对4个毛竹变型的基因组DNA进行SNP检测，获得SNP位点集(表3)，4个毛竹变型样品分别检测到的SNP数量为1 479 567～1 616 020，SNP总数为2 035 983，Ti/Tv值约为2.90。样品之间的SNP数统计结果见表4。

表3 4个毛竹变型SNP数据统计Tab.3 Summary of SNPs detected in 4 variants of Moso bamboo

表4 毛竹变型间的SNP统计Tab.4 Summary of SNPs detected between variants of Moso bamboo

2.2.2 SNP注释 运用SnpEff软件对4个毛竹变型全基因组SNP进行注释，得到其变异位点在基因组发生的区域或类型(表5)。圣音竹发生在CDS区域内的SNP共计24 426个，其中，同义突变占比为47.28%，非同义突变占比为51.56%。黄槽毛竹发生在CDS区域内的SNP共计25 319个，其中，同义突变占47.26%，非同义突变占51.50%。黄皮花毛竹发生在CDS区域内的SNP共计25 628个，其中，同义突变占47.26%，非同义突变占51.52%。厚壁毛竹发生在CDS区域内的SNP共计24 738个，其中，同义突变占47.49%，非同义突变占51.31%。

表5 4个毛竹变型SNP注释结果Tab.5 SNP annotations in 4 variants of Moso bamboo

2.3 InDel检测与注释

2.3.1 InDel检测 从4个样品InDel检测结果(表6)可以看出，圣音竹全基因组范围与编码区分别检测出174 657，3 173个InDel，CDS区有插入突变2 296个和缺失突变877个，其中纯合突变为578个。黄槽毛竹全基因组范围与CDS区分别检测出185 096，3 262个InDel，CDS区有插入突变2 323个和缺失突变939个，其中629个为纯合突变。黄皮花毛竹全基因组范围与CDS区分别检测出195 518，3 296个InDel，CDS区有插入突变2 371个和缺失突变925个，其中纯合突变为639个。厚壁毛竹全基因组范围与CDS区分别检测出176 393，3 121个InDel，CDS区有插入突变2 246个和缺失突变875个，其中纯合突变为611个。对样品在全基因组范围和CDS区的InDel长度进行统计(图1)，基因组范围存在较多的+1、-1、+2、-2类型突变，在CDS区域存在较多的+1、-1、+3、-3类型突变。

表6 4个毛竹变型全基因组和编码区InDel统计Tab.6 Summary of InDels detected in 4 variants of Moso bamboo

横坐标小于0为缺失，大于0为插入。

2.3.2 InDel注释 根据4个样品与参考基因组的InDel检测结果，进行两两比较，获得样品间的InDel变异位点(表7)。黄皮花毛竹与厚壁毛竹间的InDel数最少，为45 967；圣音竹与黄槽毛竹间InDel数最多，为48 551。

表7 毛竹变型间的InDel统计Tab.7 Summary of InDels detected between variants of Moso bambooo

2.4 SV检测与注释

2.4.1 SV检测 使用BreakDancer进行SV的检测(表8)，圣音竹共检测到148 682个SV，其中插入3 701个，占总变异2.49%；缺失31 439个，占总变异21.15%；染色体间易位107 945个，占总变异72.60%；染色体内部易位4 991个，反转321个，其他285个。黄槽毛竹共检测到176 960个SV，其中插入17 377个，占总变异9.82%；缺失36 496个，占总变异20.62%；染色体间易位116 941个，占总变异66.08%。黄皮花毛竹共检测到181 687个SV，其中插入9 498个，占总变异5.23%；缺失37 880个，占总变异20.85%；染色体间易位127 646个，占总变异70.26%。厚壁毛竹共检测到148 968个SV，其中插入2 481个，占总变异1.67%；缺失30 971个，占总变异20.79%；染色体间易位110 054个，占总变异73.88%。

表8 4个毛竹变型结构变异数量统计Tab.8 Summary of SVs detected in 4 variants of Moso bamboo

2.4.2 SV注释 对结构变异深入分析表明(表9)，圣音竹外显子区变异总数为1 894个，内含子区变异总数为789个，基因间区变异总数为32 778。黄槽毛竹外显子区变异总数为2 732个，内含子区变异总数为1 209个，基因间区变异总数为50 241个。黄皮花毛竹外显子区存在2 498个变异，内含子区存在1 064个变异，基因间区存在44 150个变异。厚壁毛竹外显子区分别存在1 731，内含子区存在773个，基因间区存在31 237个。

表9 4个毛竹变型结构变异注释结果Tab.9 SV annotations in 4 variants of Moso bamboo

2.5 DNA水平的变异基因挖掘与功能注释

2.5.1 变异基因挖掘 分别提取毛竹4个变型与参考基因组间发生非同义突变的SNP以及CDS区发生InDel与SV的基因，各样品的变异如表10所示。秆形变异的2个竹种圣音竹和厚壁毛竹基因组分别存在9 307，9 148个基因变异，其中SNP突变基因分别为4 973，4 981个，InDel基因分别为2 735，2 699个，SV突变的基因分别为1 599,1 468个，多个基因同时存在多种类型突变。秆色变异的2个竹种黄槽毛竹和黄皮花毛竹基因组分别存在10 194,9 977个基因变异，其中，SNP突变基因分别为5 094,5 095个，InDel基因分别为2 799,2 816个，SV突变基因2 301,2 066个，且多个基因同时存在多种类型突变。

表10 4个毛竹变型的变异基因统计Tab.10 Summary of gene variations in 4 variants

2.5.2 变异基因的功能注释 运用Blast软件，将挖掘的变异基因分别与GO、KEGG和COG等功能数据库进行比对，获得这些变异基因的注释，进而分析基因功能。4个毛竹变型注释到数据库的变异基因数分别为8 351，9 195，8 966，8 197个。

通过GO数据库，获得变异基因GO分类统计结果(图2)。从图中可以看出，细胞组件、分子功能和生物过程3个基因功能分类体系中，GO各分类内容对应的基因数以及基因数目所占百分比。对GO数据进一步分析发现，在细胞组分方面，叶绿素合成相关基因有1 848个，细胞壁相关基因750个，细胞骨架相关基因194个；在生物过程方面，参与类胡萝卜素合成过程的基因有60个，参与花青素合成过程中的调控以及紫外光下组织中花青素积累的相关基因有55个；在功能基因方面，参与信号转导的基因283个，转录因子475个，激素相关基因75个。

图2 变异基因的GO注释分类图Fig.2 Classification of gene variations compared with GO database by Blast

COG数据库是对基因产物进行直系同源分类的数据库。从变异基因COG分类统计结果(图3)中可以看出，在不同的功能分类中，基因数目相差很大，其中参与功能预测、转录、复制重组修复和信号转导机制的基因数目较多。基因所占比例多少，反映对应时期和环境下代谢或者生理偏向等内容。COG注释获得的参与功能预测的基因1 131个，参与糖代谢途径的基因355个，细胞壁、细胞膜合成相关的基因有276个，信号转导机制的基因数为228个，细胞周期调控、细胞分裂和染色的相关基因145个。

通过KEGG数据库，变异基因注释到114个KEGG代谢通路中，分析共获得了1 310个转录因子及813个信号转导相关单基因；在功能基因方面，共获得618个与细胞壁构建相关的代谢基因。KEGG数据库能够系统地分析基因产物在细胞中的代谢途径和功能，有助于把基因及表达信息作为一个整体的网络进行研究。从变异基因的糖代谢通路(图4)可以看出，整个通路有46种不同的酶通过复杂的生化反应形成，框内的数字分别代表酶的号码。其中的4种酶：1.2.4.1(丙酮酸脱氢酶)、1.8.1.4(二氢硫辛酸脱氢酶)、5.3.1.1(磷酸丙糖异构酶)和5.4.2.2(葡萄糖磷酸变位酶)对应的框代表变异基因中与此通路相关。

图3 变异基因的COG注释分类图Fig.3 Classification of gene variations compared with COG database

将变异基因与GO、COG、KEGG等3个功能数据库进行比对、注释，共获得1 739个参与细胞壁和细胞膜合成的基因、501个细胞骨架构建相关基因、1 785个转录因子、1 324个信号转导、104个赤霉素及76个激素代谢相关基因、1 938个叶绿素合成相关的基因、73个类胡萝卜素合成代谢相关基因和68个花青素合成调控基因。

3 结论与讨论

3.1 4个毛竹变型的秆形、秆色变异表征

经过长期的自然演化和人工栽培后，毛竹种内产生丰富的遗传变异，从而形成不同的结构形态。圣音竹、黄皮花毛竹、黄槽毛竹和厚竹是毛竹种的4个变型，具有明显的秆形或秆色特征。圣音竹竹秆从上面向基部逐渐膨大呈喇叭状，同时节间缩短，且上粗下细或下部紧靠节处有不规则凹陷，节部呈波浪状歪斜。厚竹秆略呈四方形或椭圆形，壁厚，上部近实心，分枝以下节隆起较高。黄皮花毛竹秆黄色或黄绿色基本各半，有宽窄不等的绿色纵条纹。黄槽毛竹秆主要为绿色，但节间分枝一侧纵沟槽为黄色或淡黄色。这些在竹秆形态和颜色等方面的不同变异使其具有很高的观赏价值。

3.2 4个毛竹变型的基因组变异分析

通过对4个毛竹变型进行重测序，初步统计分析了这些形态变异发生的分子数据。样品检测到的SNP总数为2 035 983，Ti/Tv值约为2.90，说明同种类型的碱基之间突变比不同类型间更容易发生。分别提取毛竹4个变型与参考基因组间发生非同义突变的SNP、CDS区发生InDel与SV的基因，可以发现其基因组基因的变异数目和变异类型呈现一定的规律，其中，2个秆形变异竹种圣音竹和厚壁毛竹基因组基因的变异数目和变异类型相近，而2个秆色变异竹种黄槽毛竹和黄皮花毛竹基因组的变异数目和变异类型相近，但秆形和秆色变异的竹种两两之间差异较大。

3.3 变异基因的功能注释与变异性状的关联分析

DNA水平的变异基因挖掘和功能注释，可以分析基因产物在细胞中的代谢途径以及这些基因产物的功能。将变异基因通过多个数据库注释，共获得1 739个基因参与细胞壁、细胞膜合成，501个细胞骨架构建相关基因，以及1 785个转录因子、1 324个信号转导、104个赤霉素和76个激素代谢相关基因。参考矢竹的研究结果，其地下茎生长发育受到从各类激素、转录因子到其下游功能基因等组成的复杂分子网络的协同调控，细胞壁、细胞骨架下游功能相关基因在矢竹地下茎节间快速生长阶段表达最高，而赤霉素合成及信号传导相关基因则在节间伸长起始阶段最高[35]。因此，进一步对毛竹变型细胞组分、生物过程和分子功能相关变异基因的研究，有利于推测毛竹秆形变异的分子调控机制。

竹类植物中的色素主要有叶绿素、类胡萝卜素和花青素3大类，色素的变化会影响植株的呈色变异。通过数据库注释，发现的叶绿素合成相关基因、类胡萝卜素合成代谢相关基因和花青素合成调控基因，深入研究这些差异基因的调控途径利于从DNA水平上解释秆色的变异。研究表明[36]，毛竹不同变异类型在净光合速率、叶绿素含量、β胡萝卜素含量、硝酸还原酶活性、游离氨基酸含量这5个生理指标中存在着显著性差异，其中，黄皮花毛竹最高，其次为厚壁毛竹，然后是圣音毛竹、黄槽毛竹。对毛竹的栽培变种进行ISSR和AFLP分子标记分析，分子遗传的分析结果与阮晓赛[37]进行的形态学研究结果相似，变型之间的遗传变异程度较小。进一步分析这4个毛竹变型基因产物的功能及其代谢途径，也有利于从基因水平上对其生物学和生理特性加以解释和分析。数据的挖掘与深入分析工作正在进行，后续分析将解析毛竹众多种质的分子遗传基础，阐析不同变异类型的分子品性、基因家族、功能基因等遗传根基。