牛場辣椒的全基因組SNP標(biāo)記分析

黃冬福， 何建文， 江葉莎， 付文婷， 范高領(lǐng)， 吳 迪， 詹永發(fā), 石燕金， 王楠藝

(1.貴州省農(nóng)業(yè)科學(xué)院辣椒研究所， 貴陽 550009； 2.遵義市農(nóng)業(yè)農(nóng)村局， 貴州 遵義 563000)

辣椒為茄科辣椒屬植物，有極大的利用價值，營養(yǎng)價值高可鮮食，是一種重要的調(diào)味品，富含的辣椒堿具有一定的藥用價值，辣椒紅素可用于食品及化妝品的著色。貴州六枝特區(qū)的牛場辣椒于2014年被認(rèn)定為國家地理標(biāo)志產(chǎn)品，具有果色深紅，肉厚，辣味適中，香氣濃郁等優(yōu)異品質(zhì)[1]。辣椒基因組上含有多種分子標(biāo)記。在眾多的分子標(biāo)記中，SNP(單核苷酸多態(tài)性)最重要且最具吸引力，其高水平的多態(tài)性、共顯性、高通量、豐富的信息量[2]廣泛用于作物育種中的遺傳多樣性分析、基因組關(guān)聯(lián)分析及比較基因組學(xué)[3-6]。

國內(nèi)外研究者獲得了辣椒的大量SNP標(biāo)記。Siddique等[7]對188個辣椒重組自交系個體和352個辣椒材料進(jìn)行重測序，開發(fā)了666 405個SNP標(biāo)記，結(jié)合經(jīng)典的QTL定位和全基因組關(guān)聯(lián)分析，獲得了3個賦予疫病廣譜抗性的主效QTL。Wu等[8]對287個辣椒材料進(jìn)行重測序，獲得了9 557 790個SNP，通過全基因組關(guān)聯(lián)分析，發(fā)現(xiàn)調(diào)控26個辣椒農(nóng)藝性狀的2 126個候選基因。Wang等[9]對辣椒不育系和保持系的線粒體基因組測序，獲得了兩者間的112個SNP，結(jié)合已知的CMS(細(xì)胞質(zhì)雄性不育)基因特征，篩選出2個最有可能決定CMS的ORF。Han等[10]對208份辣椒材料進(jìn)行重測序，開發(fā)了109 610個SNP標(biāo)記，發(fā)現(xiàn)99個SNP與辣椒素顯著關(guān)聯(lián)。孫茜[11]對辣椒抗感黃瓜花葉病毒(CMV)的基因池進(jìn)行重測序，獲得了51 969 152個SNP標(biāo)記，結(jié)合關(guān)聯(lián)分析和經(jīng)典的QTL定位，發(fā)現(xiàn)了抗CMV的1個主效QTL和2個微效QTL。趙紅[12]對349份國內(nèi)辣椒核心種質(zhì)進(jìn)行重測序，平均每份種質(zhì)獲得了7 425 498個SNP，通過全基因組關(guān)聯(lián)分析，發(fā)現(xiàn)94個SNP與果實(shí)辣味等20個性狀顯著關(guān)聯(lián)。

就目前的研究來看，各研究者開發(fā)了大量辣椒SNP標(biāo)記，但是所用品種沒有涉及牛場辣椒，而利用已有的SNP標(biāo)記重新篩選牛場辣椒特異的SNP費(fèi)時費(fèi)力且無法保證數(shù)量與質(zhì)量。另外，各研究者開發(fā)SNP采用的是簡化基因組測序法，基于此法開發(fā)的SNP無法覆蓋全基因組，SNP的數(shù)量及密度遠(yuǎn)不如全基因組測序法。因此，本研究利用全基因組重測序分析牛場辣椒的SNP標(biāo)記，為牛場辣椒遺傳圖譜構(gòu)建、重要農(nóng)藝性狀基因挖掘、遺傳改良、品種鑒定與保護(hù)奠定基礎(chǔ)。

1 材料與方法

1.1 材 料

牛場辣椒(CapsicumannuumL.)，2014年被農(nóng)業(yè)部認(rèn)定為地理標(biāo)志農(nóng)產(chǎn)品，由貴州省辣椒研究所保存。取苗期的葉片用于全基因組重測序。

1.2 基因組DNA提取與DNA質(zhì)量檢測

采用DNA secure Plant Kit(TIANGEN)試劑盒提取基因組DNA；1%瓊脂糖凝膠電泳檢測DNA是否降解和污染；NanoPhotometer?spectrophotometer (IMPLEN, CA, USA)檢查 DNA 純度；使用2.0 Flurometer (Life Technologies, CA, USA)檢測DNA 濃度。

1.3 測序文庫構(gòu)建

分別取檢測合格的DNA 樣品700 ng，通過Covaris破碎機(jī)打斷成長度為350 bp的片段，使用NEB Next?Ultra DNA Library Prep Kit(NEB, USA)構(gòu)建文庫，并將index codes添加到每個測序樣本中。使用AMPure XP系統(tǒng)(Beckman Coulter，Beverly，USA)純化DNA，DNA片段的3′末端腺苷酸化后，連接具有發(fā)夾環(huán)結(jié)構(gòu)的NEB下游銜接子以準(zhǔn)備雜交；PCR反應(yīng)前使用電泳來選擇指定長度的DNA片段，在USER酶(NEB，USA)的作用下，37 ℃ 15 min、95 ℃ 5 min；然后用Phusion High-Fidelity DNA聚合酶、Universal PCR引物和Index(X)引物進(jìn)行PCR， PCR產(chǎn)物用AMPure XP系統(tǒng)進(jìn)行純化。文庫構(gòu)建完成后，先使用Qubit 2.0軟件進(jìn)行初步定量，稀釋文庫至1 ng/μL，隨后使用Agilent Bioanalyzer 2100軟件對文庫的insert size進(jìn)行檢測，insert size符合預(yù)期后，使用Q-PCR方法對文庫的有效濃度進(jìn)行準(zhǔn)確定量(文庫有效濃度>2 nmol/L)，保證文庫質(zhì)量。

1.4 全基因組測序和質(zhì)量控制

利用Illumina HiSeq 2000平臺進(jìn)行全基因組測序，測序生成的原始圖像數(shù)據(jù)文件經(jīng)堿基識別轉(zhuǎn)化為原始測序序列(Illumina pipeline CASAVA v 1.8.2)，然后通過質(zhì)量控制去除以下不能用的reads，且兩端reads均去除：

1) 帶文庫構(gòu)建接頭的reads；

2) 未知堿基超過10%的reads；

3) 低質(zhì)量堿基(測序質(zhì)量值≤5)超過50%的reads。

1.5 序列比對

質(zhì)控后的有效測序數(shù)據(jù)經(jīng)BWA軟件比對到參考基因組(https://www.ncbi.nlm.nih.gov/genome/10896，v1.0)，比對結(jié)果經(jīng)SAMTOOLS和PICARD軟件(http://picard.sourceforge.net)去除重復(fù)。

1.6 變異檢測及注釋

設(shè)定參數(shù)(-q 1-C 50-m 2-F 0.002-d 1000)，用samtools軟件檢測原始的SNP，然后用如下標(biāo)準(zhǔn)進(jìn)行過濾：

1) 變異位置的測序深度>4；

2) 質(zhì)量值>20。

用ANNOVAR軟件對SNP進(jìn)行變異注釋。

2 結(jié)果與分析

2.1 測序數(shù)據(jù)產(chǎn)出及與參考基因組比對分析

由表1可知，經(jīng)過全基因組重測序，共獲得785 408 260條Raw reads，平均讀長150 bp，堿基總長為117.8 Gb；通過質(zhì)量控制去除帶接頭的reads、未知堿基超過10%的reads以及低質(zhì)量(Q值≤5)堿基數(shù)超過50%的reads，得到高質(zhì)量的clean reads；clean reads共有783 349 390條，覆蓋98.23%的基因組；去除clean reads中的非特異reads，獲得771 948 477條有效reads，其中758 748 158條reads能錨定到“遵辣1號”參考基因組上，配對率為98.29%，平均測序深度為36.35×。

表1 牛場辣椒測序數(shù)據(jù)產(chǎn)出及與參考基因組比對情況Table 1 The sequencing data of Niuchang pepper and the result compared with the reference genome

2.2 染色體上SNP的分布規(guī)律和特征

根據(jù)與參考基因組的比較，SNP分為純合和雜合類型，分別占SNP總數(shù)的59.12%、40.88%。每條染色體上的SNP總數(shù)不同，10號染色體上SNP總數(shù)最多(1 371 387個)，4號染色體上SNP總數(shù)最少(349 894個)。每條染色體上純合和雜合SNP的數(shù)量也不同， 10號染色體上的純合SNP最多(889 193個)，9號染色體上的純合SNP最少(169 010個)，9號染色體上的雜合SNP最多(779 326個)，4號染色體上的雜合SNP最少(93 614個)，具體見表2。

每條染色體上密度最高區(qū)域SNP出現(xiàn)頻率不同，10號染色體上密度最高區(qū)域SNP出現(xiàn)頻率最高(1 904個/100 kb)，4號染色體上密度最高區(qū)域SNP出現(xiàn)頻率最低(1 162個/100 kb)，具體見表3和圖1。

表2 每條染色體上SNP的數(shù)量Table 2 The number of SNP on each chromosome

表3 每條染色體上SNP密度最高區(qū)域及其出現(xiàn)頻率Table 3 The frequency of SNP in highest SNP densityarea of chromosomes

2.3 基因組不同位置SNP分布特征

牛場辣椒基因組中的SNP分布在5個不同位置：基因上游、基因內(nèi)、基因下游、基因上游/下游、基因間。(基因上游是指基因上游1 kb區(qū)域；基因內(nèi)指基因內(nèi)部；基因下游指基因下游1 kb區(qū)域；基因上游/下游指基因上游1 kb區(qū)域，同時也在另一基因的下游1 kb區(qū)域；基因間指兩個基因間區(qū))。5個不同位置的SNP數(shù)量不同且差異顯著，SNP數(shù)量從多到少依次為基因間>基因內(nèi)>基因上游>基因下游>基因上游/下游(表4)，基因間、基因內(nèi)、基因上游、基因下游的SNP占比依次為94.68%、3.64%、0.9%、0.74%(圖2)。

基因內(nèi)不同位置所含SNP數(shù)量也不同，基因內(nèi)包含外顯子、內(nèi)含子、剪接位點(diǎn)3個位置，所含SNP數(shù)量分別為51 242、281 002、288個，SNP數(shù)量從多到少依次為內(nèi)含子、外顯子、剪接位點(diǎn)。針對外顯子區(qū)域，根據(jù)SNP變異引起的密碼子變化， 可將SNP分為4種類型：終止子獲得、終止子缺失、同義突變、非同義突變，數(shù)量分別為710、188、19 079和31 265，SNP數(shù)量依次為非同義突變>同義突變>終止子獲得>終止子缺失。

圖1 每條染色體上的SNP密度熱圖 Fig.1 The density heat map of SNP on each chromosome

圖2 牛場辣椒基因組中不同位置的SNP數(shù)量差異Fig.2 The number difference of SNP on the different position in the genome of Niuchang pepper

表4 牛場辣椒基因組中SNP的位置及相應(yīng)數(shù)量Table 4 The position and corresponding number of SNP in the genome of Niuchang pepper

2.4 SNP突變頻譜

全基因組SNP突變可分為6類：T∶A>G∶C，T∶A>C∶G，T∶A>A∶T，C∶G>T∶A，C∶G>G∶C和C∶G>A∶T。以T∶A>C∶G為例，此種類型SNP突變包括T>C和A>G。由于測序數(shù)據(jù)既可比對到參考基因組的正鏈，也可比對到參考基因組的負(fù)鏈，當(dāng)T>C類型突變出現(xiàn)在參考基因組正鏈上，A>G類型突變即在參考基因組負(fù)鏈的相同位置，所以將T>C和A>G劃分成一類。C∶G>T∶A的數(shù)量最多(3 109 688個)，C∶G>G∶C的數(shù)量最少(474 542個)。T∶A>G∶C、T∶A>A∶T、C∶G>G∶C、C∶G>A∶T為顛換，總數(shù)為3 057 408個，T∶A>C∶G及C∶G>T∶A為轉(zhuǎn)換，總數(shù)為6 094 214個，發(fā)生轉(zhuǎn)換的數(shù)量是顛換的1.99倍(圖3)。

圖3 SNP突變頻譜Fig.3 The mutation frequency and type of SNP

3 討論與結(jié)論

3.1 基于WGRS技術(shù)開發(fā)辣椒SNP標(biāo)記的高效性

辣椒全基因組序列的公布為其分子育種帶來前所未有的機(jī)遇。想要開展辣椒分子育種，就必須對群體中所有個體進(jìn)行基因分型。利用傳統(tǒng)方法對辣椒進(jìn)行基因分型費(fèi)用高、耗時耗力，低水平的分子標(biāo)記也是基因分型的重要挑戰(zhàn)?；谛乱淮鷾y序技術(shù)(NGS)的基因分型通量高、成本低、分子標(biāo)記密度高。

全基因組重測序(WGRS)是新一代測序技術(shù)(NGS)的一種。利用WGRS技術(shù)開發(fā)分子標(biāo)記具有標(biāo)記密度大、有效標(biāo)記多、準(zhǔn)確率和特異性高、穩(wěn)定性好的優(yōu)勢。本研究利用WGRS技術(shù)共鑒定出9 141 358個SNP，SNP的出現(xiàn)頻率為1個/366 bp，其中51 242個SNP位于外顯子。Ahn利用WGRS技術(shù)獲得了6 840 889個辣椒SNP，其中39 955個SNP位于外顯子[13]。

相比而言，利用簡化基因組測序開發(fā)SNP標(biāo)記的效率較低。Nimmakayala等[14]利用簡化基因組測序僅獲得66 960個辣椒SNP，SNP的出現(xiàn)頻率為1個/40.7 kb，其中僅有2 521個SNP位于外顯子。Nimmakayala等[15]采用簡化基因組測序獲得77 407個辣椒SNP，SNP的出現(xiàn)頻率為1個/35.6 kb，其中26 697個SNP位于外顯子。由此可見，利用WGRS獲得的SNP標(biāo)記數(shù)量、出現(xiàn)頻率、外顯子上的SNP數(shù)量遠(yuǎn)高于簡化基因組重測序[14-15]。

功能標(biāo)記是分子標(biāo)記的一種，基于功能基因內(nèi)的多態(tài)序列開發(fā)，與常規(guī)的分子標(biāo)記相比，與表型的連鎖程度更緊密[16]，外顯子上的SNP最可能作為功能標(biāo)記。因此，通過檢測功能標(biāo)記能更準(zhǔn)確地預(yù)測表型，在加速育種進(jìn)程的同時極大地提高了標(biāo)記輔助選擇的準(zhǔn)確性。本研究獲得的51 242個外顯子SNP作為潛在的功能標(biāo)記，將成為辣椒分子育種的有力武器。

3.2 牛場辣椒SNP標(biāo)記分布特征

SNP在牛場辣椒基因組上的分布表明，基因間的SNP數(shù)量比基因內(nèi)多，是基因內(nèi)的26倍，內(nèi)含子區(qū)域的SNP數(shù)量比外顯子多，是外顯子區(qū)域的5.5倍，這種分布規(guī)律與Kim等[17]的研究結(jié)果相似。外顯子上，SNP引起的終止子獲得有710處，終止子獲得是指堿基突變導(dǎo)致終止密碼子提前出現(xiàn)。其產(chǎn)生截短的蛋白質(zhì)，從而使基因散失原來的功能，并進(jìn)一步引發(fā)作物表型變異。因此，終止子獲得對基因功能研究具有重要意義。AFS1基因發(fā)生5 bp缺失及G/A轉(zhuǎn)換，翻譯提前終止，導(dǎo)致水稻小穗異常，小穗上多出一個外稃狀器官，內(nèi)稃發(fā)生不同程度的退化[18]。GmSG基因發(fā)生A/G轉(zhuǎn)換，翻譯提前終止，導(dǎo)致大豆種皮顏色由黃色轉(zhuǎn)變?yōu)辄S/綠色[19]。SNFL1基因單內(nèi)含子上最后一個堿基發(fā)生單核苷酸突變，翻譯提前終止，導(dǎo)致水稻旗葉變短變窄[20]。OsCUL 3 a蛋白翻譯提前終止，導(dǎo)致水稻flg 22、幾丁質(zhì)誘導(dǎo)的活性氧、發(fā)病相關(guān)基因的表達(dá)量明顯增加，進(jìn)而產(chǎn)生類病斑[21]。

本研究發(fā)現(xiàn),牛場辣椒10號染色體上的SNP數(shù)量最多，且18.17%的純合SNP位于該染色體上，4號染色體上的SNP數(shù)量最少，與Nimmakayala[15]的研究結(jié)果不同，Nimmakayala的研究顯示辣椒3號染色體上SNP數(shù)量最多，9號染色體上SNP數(shù)量最少，與Ahn[13]的研究結(jié)果也不同，Ahn的研究顯示辣椒品種PRH 11號染色體上SNP數(shù)量最多，8號染色體上SNP數(shù)量最少。造成不同染色體上SNP分布頻率差異的原因可能是不同研究者使用的辣椒品種不同。

SNP引起的點(diǎn)突變可分為轉(zhuǎn)換和顛換?，F(xiàn)有研究表明，人類基因組CpG中的胞嘧啶C能自發(fā)脫氨基變成胸腺嘧啶T，導(dǎo)致C∶G>T∶A的轉(zhuǎn)換比T∶A>C∶G的轉(zhuǎn)換多，同時，轉(zhuǎn)換比顛換更容易發(fā)生，且轉(zhuǎn)換的SNP是顛換的2倍[22]。牛場辣椒C∶G>T∶A的轉(zhuǎn)換比T∶A>C∶G的轉(zhuǎn)換多，且轉(zhuǎn)換的SNP數(shù)量明顯高于顛換，是顛換的1.99倍，與前人的研究結(jié)果基本一致。

本研究對牛場辣椒進(jìn)行全基因組重測序，獲得了9 141 358個SNP，SNP的出現(xiàn)頻率為1個/366 bp，展示出WGRS技術(shù)開發(fā)辣椒SNP標(biāo)記的高效性；其中51 242個SNP位于外顯子上，外顯子上的SNP具有開發(fā)成功能標(biāo)記的巨大潛力，功能標(biāo)記與表型的連鎖程度更緊密，能極大地提高標(biāo)記輔助選擇的準(zhǔn)確性，是辣椒分子育種的有力武器；外顯子上的SNP產(chǎn)生了710處終止子，終止子獲得會產(chǎn)生截短的蛋白質(zhì)，使基因散失原來的功能，并進(jìn)一步引發(fā)作物表型變異，對基因功能研究具有重要意義。