上饒早梨 ‘六月雪’和 ‘黃皮消’全基因組重測序分析

洪森榮，曾清華，譚 鑫，陳永華，鄭亞嬌，徐迎昕，邱夢琴

（上饒師范學院 生命科學學院，江西 上饒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