洪森荣，曾清华，谭 鑫，陈永华，郑亚娇，徐迎昕，邱梦琴

（上饶师范学院 生命科学学院，江西 上饶334001）

上饶早梨Pyrus pyrifolia是江西省上饶县花厅、田墩、五府山等11个乡镇的国家地理标志农产品，被认为是药食兼优的夏令佳品［1-2］，栽培历史长，以 ‘黄皮消’ ‘Huangpixiao’和 ‘六月雪’ ‘Liuyuexue’2 个品种品质最为优良［3］。 对上饶早梨脱毒快繁［4-5］和品种鉴定［6-7］的研究表明： 来源上饶县花厅镇的上饶早梨道地性农家种 ‘六月雪’和 ‘黄皮消’经简单重复序列标记（SSR）和扩增片段长度多态性（AFLP）分子标记的聚类分析聚为1类，但生物学性状、农艺性状和果实品质均有一定差异［8］。因此，找出两者基因差异对地方品种的育种选种工作十分重要。自2000年模式植物拟南芥Arabidopsis thaliana的基因组被完全解析后，已有越来越多的种质被全基因组测序［9］；研究认为，野生种、农家种、栽培种的全基因组重测序是筛选重要性状关键基因的一个重要内容［10］。梨基因组测序的完成［11-13］和高通量测速技术的快速发展，使梨种质资源的全基因组变异分析成为可能，而全基因组单核苷酸多肽位点（SNP）和插入缺失位点（INDEL）相关标记的开发，对作物分子辅助育种和基因组学研究具有重要的作用。周贺等［14］对褐色砂梨 ‘黄花’Pyrus pyrifolia‘Huanghua’色泽形成期的果皮转录组数据进行SNP分子标记开发，共筛选到1 178个可能与果皮色泽形成相关的SNP标记。李节法［15］以6年生砂梨品种 ‘翠冠’P.pyrifolia‘Cuiguan’的果实膨大早期、中期和成熟期样品进行比较转录组学分析，鉴定到4 121个选择性剪切位点， 30 560 个 SNP 位点和 7 443 个 SSR 标记位点。 WU 等［16］、 MONTANARI等［17］和 TERAKAMI等［18］也通过SNP标记构建了梨遗传图谱。但梨的INDEL标记研究较少。本研究以上饶县花厅镇上饶早梨道地性农家种 ‘六月雪’和 ‘黄皮消’为材料，通过全基因组重测序，深度挖掘其基因组SNP，INDEL和结构变异位点（SV）等位点，为上饶早梨优质品种相关标记的开发、优异基因的挖掘提供参考。

1 材料和方法

1.1 材料

上饶县花厅镇上饶早梨道地性农家品种 ‘六月雪’和 ‘黄皮消’试管苗由上饶师范学院生命科学学院植物组织培养室提供。

1.2 方法

1.2.1 DNA提取 参照十六烷基三甲基溴化铵（CTAB）提取法提取样品DNA。

1.2.2 测序 以吸光度比值D（260）/D（280）和 Qubit 2.0（Life technologies,USA）， Bioanalyzer 2100（Agilent，Germany）软件分析完成总DNA样品的质量控制。称取1 μg基因组总DNA片段化处理至300～400 bp，进行末端修复、末端加 ‘A’处理、 接头（Adapters）连接反应，连接至Illumina公司测序平台的测序接头后进行聚合酶链式反应（PCR）扩增。连接产物用质量分数2%的琼脂糖凝胶电泳，选择400～500 bp的片段，随后用连接介导PCR（LM-PCR）进行扩增获得DNA文库。按照Illumina公司HiSeq 2500测序系统（Hiseq 2500）的操作说明对形成的DNA文库进行双端125 bp的高通量测序。

1.2.3 数据分析 鉴于Illlumina Hiseq 2500错误率对结果的影响，对原始数据进行质量控制（QC）。参考数据库为西洋梨基因组Pyrus communis Genome v1.0数据库。分别对每个样本使用bowtie 2软件进行测序短序列匹配（reads mapping）， 并用 UnifiedGenotyper软件进行 SNP和 INDEL的提取［19-21］。 采用 ANNOVAR软件对获得的SNP和INDEL进行功能注释［22］。SNP常见变异分析（common variation,CV）及差异INDEL分析首先获取2个样品相同位置的SNP/INDEL，再根据非同义SNP（nonsynonymous SNP，nsSNP）/INDEL获取相关基因［23］。将差异SNP和INDEL分别与转录组数据进行关联分析，分别考察差异SNP和INDEL相关的表达数据［24］。候选基因的富集分析递交至AgriGO软件用于富集基因本体术语（gene ontology terms）［25］。拷贝数变异位点（CNV）分析采用 VarScan 软件进行［26］。

2 结果与分析

2.1 测序数据预处理

以HiSeq 2500测序系统提供的起初测序数据为原始数据，即各样本测序得到的短序列（reads）数及碱基总数，共得到275 092 822个短序列和34 661 695 572个碱基，其中 ‘六月雪’中含短序列140 696 312个，碱基17 727 735 312个， ‘黄皮消’中短序列134 396 510个，碱基16 933 960 260个。为剔除Illlumina Hiseq 2500错误率对结果的影响，需对原始数据进行质量控制，包括去除低质量序列，去除接头，以进行后续工作。质量控制后得到232 434 654个短序列和29 286 766 152个碱基，总有效数据比例为84.5%；其中 ‘六月雪’含短序列119 056 332个，碱基150 010 978 322个，有效数据比例为84.6%； ‘黄皮消’含短序列113 378 320个，碱基14 285 668 320个，有效数据比例为84.4%。

2.2 短序列匹配（reads mapping）统计

经过测序将短序列匹配至参考基因组。 ‘六月雪’组中总匹配的短序列数为68 859 074个，占所有短序列数的57.8%，含唯一匹配的短序列数为31 889 494个，占总匹配数的26.8%；覆盖全基因的深度为24.35，覆盖全基因组的百分比为93.0%；当覆盖深度≥3时，覆盖全基因组的百分比为89.4%。‘黄皮消’组总匹配的短序列数为66 580 757个，占所有短序列数的58.7%，含唯一匹配的短序列数为32 165 247个，占总匹配数的28.4%；覆盖全基因的深度为23.89，覆盖全基因组的百分比为93.0%，覆盖深度≥3时，覆盖全基因组的百分比为89.5%。

2.3 SNP统计分析

‘六月雪’中共得到SNP 6 171 357个，在编码区的无义突变有335 659个，有义突变有281 871个；因SNP突变获得终止子6 001个，丢失终止子1 226个；在基因5′非翻译区（5′UTR内的SNP总数、在3′UTR内的SNP总数及位于5′UTR和3′UTR间的SNP总数均为0；在不同可变剪切的基因组区域内找到SNP 3 383个，内含子区域内找到1 298 966个，启动子区域内找到1 301 726个，基因间区域内找到2 942 525个。

‘黄皮消’中共得到SNP 6 140 603个，在编码区的无义突变有332 280个，有义突变有278 064个；因SNP突变获得终止子6 034个，丢失终止子1 210个；在基因5′非翻译区（5′UTR）内的SNP总数、在3′UTR内的SNP总数及位于5′UTR和3′UTR间的SNP总数均为0，在不同可变剪切的基因组区域内找到SNP 3 274个，内含子区域找到1 285 052个，启动子区域内找到1 291 598个，基因间区域内找到2 943 091个。

2.4 常见变异（common variant，CV）分析

对获得的335 659个（‘六月雪’）和332 280个（‘黄皮消’）nsSNP进行常见变异分析发现，2个品种共有2 282个nsSNP关联了2 067个基因，nsSNP对基因的平均关联频率超过了1。其中，烟草花叶病毒耐药蛋白N（PCP017781），假定的抗病 RPP13样蛋白1（PCP007357）， 可能的抗病 RPP8 样蛋白2（PCP030706）， 可能的LRR类受体丝氨酸/苏氨酸蛋白激酶 At3g47570（PCP021305），未注释蛋白1（PCP008176），烟草花叶病毒耐药蛋白N（PCP007457），烟草花叶病毒耐药蛋白N（PCP018815），未注释蛋白2（PCP021753），含重复锚蛋白的蛋白质 At5g0262 0（PCP022078），烟草花叶病毒耐药蛋白N（PCP030478），抗病蛋白RGA2（PCP014224），可能的LRR类受体丝氨酸/苏氨酸蛋白激酶At1g5342 0（PCP031574）， 假定的抗病蛋白RGA3（PCP023580）等蛋白编码的基因则关联了10个以上nsSNP。为进一步研究候选基因的选择压力，对得到的2 067个基因的nsSNP与同义SNP的比值（nsSNP/synonymous SNP，r）进行考察，发现r的对数呈现正态分布（图1），其值约为2，说明进化有一定的正向选择压力。

图1 候选基因的nsSNP/同义SNP比例分布曲线图Figure 1 Curve map of nsSNP/synonymous SNP ratio distribution of candidate genes

2.5 nsSNP关联的候选基因的基因本体论富集分析

由于生物学定义混乱的原因，不同的生物学数据库可能会使用不同的术语。为了解决这个问题，基因本体联合会（Gene Onotology Consortium）建立了 “基因本体论”（gene ontolog，GO）数据库，目的是通过利用统一化的、结构化的语言建立一个适用于不同物种的、对基因和蛋白质功能进行定义和描述，并且能够随着研究的不断深入而更新的语言词汇标准。GO数据库包含基因参与的生物过程、所处细胞位置及具有的分子功能3个方面信息，其注释信息可对基因功能进行预测。GO显著性富集分析以基因本体术语（GO term）为单位，确定差异表达基因行使的主要生物学功能。分析全局GO功能与候选基因所处的功能发现，刺激、结合反应等功能的基因相对于背景基因（1 306条）富集（图2）。通过差异基因富集得到的GO terms（图3及表1）可知，全部25 698条信息中，ADP结合、蛋白酪氨酸激酶活性、腺苷核糖核酸结合、嘌呤核苷结合、核苷结合、腺苷核苷酸结合、防御反应、蛋白丝氨酸/苏氨酸激酶活性、转移酶活性、转运含磷基团、嘌呤核糖裂解键、核糖核酸结合、嘌呤核苷酸结合、蛋白激酶活性、RNA导向的DNA聚合酶活性、磷酸转移酶活性、醇基作为受体、DNA聚合酶活性、RNA依赖性DNA复制、翻译后蛋白质修饰、脯氨酸氨基酸、磷酸化、蛋白质改性过程等GO terms具有显著性意义。

图2 nsSNP关联的候选基因总体GO分布图Figure 2 Overall GO distribution of candidate gene （nsSNP）

图3 nsSNP关联的候选基因GO富集柱状图Figure 3 GO enriched histogram of candidate gene （nsSNP）

2.6 INDEL统计分析

‘六月雪’样本共得到INDEL 800 388个，编码区内移码插入总数为6 092个，移码缺失总数8 884个；编码区内因INDEL突变获得终止子426个，丢失终止子107个；在基因5′UTR内的INDEL总数为31个，在基因3′UTR内的INDEL总数和位于5′UTR和3′UTR间的INDEL总数均为0，在不同可变剪切的基因组区域内找到INDEL 786个，内含子区域内找到201 635个，启动子区域内找到198 924个，基因间区域内找到383 503个。

‘黄皮消’样本共得到INDEL 799 603个，编码区内移码插入总数为6 021个，移码缺失总数8 708个；编码区内因INDEL突变获得终止子440个，丢失终止子105个；基因5′UTR内找到INDEL 26个，基因3′UTR内、不同基因的5′UTR和3′UTR间则未找到；在不同可变剪切的基因组区域内找到INDEL 758个，内含子区域内找到199 949个，启动子区域内找到198 089个，基因间区域内找到385 507个。

2.7 相关INDEL差异分析

2个样本共获得INDEL15 509和15 274个。对这些INDEL的差异分析发现，共有5 115个INDEL关联到了3 682个基因，其中24个终止子丢失（stop-loss）INDEL关联了24个基因，165个终止子获得（stop-gain）INDEL关联了160个基因，1 901个移码插入（frame-shift insertion）INDEL关联了1 629个基因，3 025个移码缺失（frame-shift deletion）INDEL关联了2 453个基因。分析发现1 276个基因内的INDEL数量超过1个；假定的抗病 RPP13样蛋白 1（PCP007357），烟草花叶病毒耐药蛋白 N（PPCP015254），未注释蛋白1（PCP015680），烟草花叶病毒耐药蛋白N（PCP030478），假定的酰胺酶C869.01（PCP023678），ATP依赖的 RNA解旋酶 DHX36（PCP003694），甘露糖寡糖 α-1,2-甘露糖苷酶MNS1（PCP005093）， 未注释蛋白 2（PCP017985）， 富有丝氨酸/精氨酸的分裂因子 SC35（PCP000011）， 来自转座子逆转录病毒相关的Pol聚蛋白TNT 1-94（PCP032808），未注释蛋白3（PCP031973）等蛋白编码的基因均具有7个以上INDEL。

2.8 INDEL关联的候选基因的GO富集分析

分析全局GO功能发现，INDEL关联的候选基因中，刺激响应、结合、催化活性等功能相对于背景基因（1 306条）来说更加富集；INDEL关联的差异基因的富集分析（图5及表2）结果与nsSNP关联的候选基因的GO富集分析一致，全部25 698条信息中。

表1 nsSNP关联的候选基因及其功能的GO富集分析Table 1 GO terms enriched by candidate gene （nsSNP）

2.9 拷贝数变异（copy number variations，CNV）分析

用VarScan软件对拷贝数变异（CNV）分析，共获得CNV 37 039个，其中缺失CNV（deletion CNV）20 629个， 中性 CNV（neutral CNV）7 577个， 扩增 CNV（amplification CNV）8 833个。

3 讨论

近年来，迅猛发展的基因组测序技术已被广泛应用于植物基因组重测序等研究［27］。SNP是基因组DNA序列上广泛存在的最基本的变异形式［28］，植物基因组上平均每数百bp就存在一个SNP［29］。将SNP和传统分子标记相结合用于分子辅助育种，通过全基因组重测序（WGR）得到测序数据，与参考基因组进行序列比对，可分析SNP遗传变异信息，开发出数量较为丰富的分子标记，实现遗传资源的高效利用［30］。

表2 INDEL关联的候选基因及其功能的GO富集分析Table 2 GO terms enriched by candidate gene （INDEL）

全基因组重测序为从基因组水平开发SNP标记提供了新的技术条件，将SNP识别、验证和基因型分析与传统分子标记相结合，能快速挖掘到候选基因和获得导致表型的SNP位点［9］。中国6个玉米Zea mays优良自交系的非重复区大约存在1 272 134个SNP和30 178个 INDEL，其中68 966个 SNP和571个INDEL位于功能基因内［31］。对梁Setaria italica‘SLX’品种的全基因组重测序发现， ‘SLX’基因组存在762 082个SNP，26 802个INDEL，10 109个SV［32］；对包括野生种、培育种及改良品种在内的360份番茄Lycopersicon esculentum进行全基因组重测序，发现番茄的驯化和改良是2个相对独立的过程，影响果实颜色的主效基因是SlMYB12［33］。对302株大豆Glycine max进行全基因组重测序，检测到了162个受选择的拷贝数变异（CNV），并发现植株进化、发育性状与受选择区域相关［34］。大豆品种 ‘齐黄34’‘Qihuang34’经全基因组重测序检测到1 519 494个SNP，357 549个INDEL，4 506个SV，17 748个基因变异，其中转录、复制、重组、修复、信号传导机制等 6个功能类序列存在较多的变异基因［35］。梨的全基因组重测序少见报道［11-13］，大多是通过转录组分析来发掘SNP位点并进行功能注释。编码区内的SNP一般可分为同义SNP（synonymous SNP）和非同义SNP（non-synonymous SNP），其中同义 SNP所致的编码序列的改变不会引起氨基酸序列变化，而非同义SNP会使氨基酸序列发生改变，最终影响蛋白质序列［36］，因此认为非同义SNP是导致生物性状改变的直接原因，开发并研究这类SNP标记往往具有更为重要的生物学意义。本研究中， ‘六月雪’和 ‘黄皮消’样本中总SNP数量分别为6 171 357和6 140 603个，编码区内无义突变内总数（nsSNP）分别为335 659和332 280个；对nsSNP的CV分析发现共有2 282个nsSNP关联了2 067个基因。2个样本的总INDEL数量分别为800 388和799 603个，位于编码区内的分别有15 509和15 274个；共5 115个差异INDEL关联到了3 682个基因，令烟草花叶病毒耐药蛋白N和抗病蛋白等关键基因发生变异。差异nsSNP基因和差异INDEL富集得到的GO terms一致。针对这些突变位点进行SNP和INDEL相关标记的开发、优异基因的挖掘将为分子标记辅助育种提供重要的标记资源，对上饶早梨 ‘六月雪’和 ‘黄皮消’育种研究具有重要的指导意义。

图4 INDEL关联的候选基因总体GO分布图Figure 4 Overall GO distribution of INDEL candidative gene

图5 INDEL关联的候选基因GO富集分布图Figure 5 GO enrichment distribution map of INDEL candidative gene