基于ddGBS的黑楊派楊樹SNP位點(diǎn)挖掘

劉粉粉，姜小龍，翁文源，龔細(xì)娟，徐剛標(biāo)

（1.中南林業(yè)科技大學(xué) 林學(xué)院，湖南 長(zhǎng)沙 410004；2.湖南省泰格林紙集團(tuán)有限公司，湖南 岳陽(yáng) 414000）

楊樹Populus是楊柳科Salicaceae 楊屬Populus的統(tǒng)稱，染色體組為2n=38，生長(zhǎng)迅速、適應(yīng)性強(qiáng)、無(wú)性繁殖和有性繁殖容易，是北半球中緯度地區(qū)重要的工業(yè)用材林、防護(hù)林及園林綠化樹種[1]。其中，黑楊派（SectionAigeiros）優(yōu)良無(wú)性系在全球楊樹人工林產(chǎn)業(yè)中占有重要地位[2]，世界上90%的楊樹優(yōu)良無(wú)性系是從黑楊派的美洲黑楊P.deltoides和歐洲黑楊P.nigra種內(nèi)或種間雜交子代中選育出來(lái)的。由于黑楊派無(wú)性系種質(zhì)遺傳基礎(chǔ)狹窄，親緣關(guān)系較近，很難鑒別形態(tài)特征差異[3]。

分子標(biāo)記揭示的DNA 多態(tài)性，是現(xiàn)代遺傳學(xué)研究的基礎(chǔ)，也是物種分類、種質(zhì)鑒定的強(qiáng)有力工具[4]。Rajora 等[5]采用SSR 和RAPD 標(biāo)記鑒定歐美楊無(wú)性系，Chauhan 等[6]利用AFLP 標(biāo)記鑒別緣毛楊P.ciliata與馬氏楊P.maximowiczii雜交子代，宋紅竹等[7]和李善文等[8]分別采用AFLP標(biāo)記分析楊樹無(wú)性系間、種間、派間及親本與子代間親緣關(guān)系，郭麗琴等[9]采用SRAP 標(biāo)記分別研究楊屬不同種親緣關(guān)系，賈會(huì)霞等[10]、鄭濤等[11]及劉超凡等[12]，先后采用SSR 標(biāo)記構(gòu)建楊樹無(wú)性系指紋圖譜。但是，這些標(biāo)記類型都有其局限性，不同學(xué)者使用的分子標(biāo)記類型及其引物不同，可能會(huì)導(dǎo)致研究結(jié)果不能共享。

單核苷酸多態(tài)性（single nucleotide polymorphism,SNP）是單堿基突變引起的DNA 序列多態(tài)性。雖然單個(gè)SNP 位點(diǎn)的多態(tài)信息含量較低，但數(shù)以百萬(wàn)計(jì)潛在的全基因組SNP 位點(diǎn)或小插入/缺失的信息，為高通量基因分型奠定了基礎(chǔ)[4,13]。SNP標(biāo)記摒棄了凝膠電泳檢測(cè)技術(shù)，位點(diǎn)多為雙等位、共顯性，基因分型結(jié)果便于記錄和分析解讀，適合自動(dòng)化分析，在不同的檢測(cè)技術(shù)間重現(xiàn)度高[13-14]。但是，挖掘海量SNP 位點(diǎn)的高通量基因測(cè)序分型，耗時(shí)且十分昂貴[4]。采用酶切技術(shù)針對(duì)物種基因組特定區(qū)域進(jìn)行簡(jiǎn)化基因組測(cè)序，降低了測(cè)序量和測(cè)序成本，已成為挖掘SNP 位點(diǎn)信息的主要方法[4,15]。

2011年，Elshire 等提出的基因分型測(cè)序（Genotyping by sequencing，GBS），相對(duì)于其它簡(jiǎn)化基因組測(cè)序技術(shù)，簡(jiǎn)化了文庫(kù)的構(gòu)建過(guò)程，不受參考基因組限制，無(wú)需進(jìn)行片段大小選擇，加快了DNA 序列數(shù)據(jù)獲取的進(jìn)度，尤其對(duì)多態(tài)性低、重復(fù)序列高的物種，具有獨(dú)特的優(yōu)勢(shì)[16]，已被廣泛應(yīng)用于SNP 標(biāo)記開發(fā)、物種分類、種質(zhì)鑒定、遺傳多樣性分析、基因組選擇、遺傳圖譜構(gòu)建及分子輔助選擇育種[4,14,18-20]。楊樹SNP 標(biāo)記研究報(bào)道較多[18-20]，測(cè)序方法多為從頭測(cè)序（De novo sequencing）[18]、酶切簡(jiǎn)化基因組測(cè)序RADseq (Restriction-site associated DNA sequencing)[19]和轉(zhuǎn)錄組測(cè)序[20-21]，基于GBS 技術(shù)挖掘楊樹基因組SNP 位點(diǎn)信息，鮮有報(bào)道[22]。本研究以46 份歐美楊無(wú)性系種質(zhì)為材料，利用ddGBS 技術(shù)挖掘高質(zhì)量SNP 位點(diǎn)信息，旨在為楊樹SNP 標(biāo)記開發(fā)、指紋圖譜構(gòu)建、遺傳多樣性分析以及重要經(jīng)濟(jì)性狀關(guān)聯(lián)分析等奠定前期研究基礎(chǔ)。

1 材料與方法

1.1 材 料

46 份黑楊派無(wú)性系樣本采自湖南省岳陽(yáng)市泰格林紙集團(tuán)有限責(zé)任公司國(guó)家楊樹種質(zhì)資源庫(kù)。采集的葉片樣品，用硅膠快速干燥后，保存在-4℃冰柜中。

1.2 方 法

1.2.1 DNA 提取

采用天根試劑盒提取基因組DNA，0.8%瓊脂糖凝膠電泳測(cè)定DNA 質(zhì)量，核酸蛋白儀測(cè)定DNA 的濃度與純度。

1.2.2 GBS 測(cè)序文庫(kù)構(gòu)建

采用MseI(T|TAA)和TaqaI(T|CGA)雙酶切基因組，T4DNA 連接酶連接其兩端的普通接頭和條形碼（barcode）接頭，全自動(dòng)磁珠純化系統(tǒng)回收純化連接后的片段。以回收片段為模板，進(jìn)行PCR擴(kuò)增，混合樣本，再次回收純化。利用 Illumina Hiseq2000 高通量測(cè)序平臺(tái)，進(jìn)行2×145 bp雙末端測(cè)序。GBS 文庫(kù)的構(gòu)建及測(cè)序，委托上海美吉生物醫(yī)藥科技有限公司完成。

1.2.3 測(cè)序數(shù)據(jù)質(zhì)控

采用Fastp v 0.19.6(https:// https://github.com/OpenGene/fastp)軟件，對(duì)原始數(shù)據(jù)進(jìn)行質(zhì)控。質(zhì)控標(biāo)準(zhǔn)為：去除雙端測(cè)序序列中接頭序列和含N比例大于10%的序列，5′端質(zhì)量評(píng)分低于Q20（堿基錯(cuò)誤識(shí)別概率為1%）、3′端質(zhì)量評(píng)分低于Q30（錯(cuò)誤識(shí)別概率為0.1%）的堿基，以及超過(guò)40%堿基質(zhì)量評(píng)分低于Q15（合格）的序列。

采用FastQC v 0.11.9(https://github.com/s-andrews/FastQC/archive/ v0.11.9.tar.gz)軟件，評(píng)估質(zhì)控后的數(shù)據(jù)。利 用Stacks v2.55(http://catchenlab.life.illinois.edu/stacks/)軟件中的process-radtags 程序，參數(shù)選擇默認(rèn)值，將質(zhì)控后的序列片段截成長(zhǎng)度為140 bp，用于后續(xù)SNP 變異位點(diǎn)開發(fā)。

1.2.4 SNP 位點(diǎn)開發(fā)

采用Bowtie2 v2.3.2(http://bowtie-bio.sourceforge.net/index.shtml)軟件，將剪切后的數(shù)據(jù)比對(duì)到毛果楊（P.trichocarpa）參考基因組（https://www.ncbi.nlm.nih.gov/assembly/GCF_0000027 75.4）上。采用Samtools v1.12(http://www.htslib.org/) 軟件，按照SNP 位點(diǎn)在基因組上的位置進(jìn)行排序。運(yùn)行Stacks v2.55 軟件中的Stacks-gstacks 程序，將測(cè)序片斷組裝為完整、無(wú)斷點(diǎn)的重疊克隆群，再比對(duì)到位點(diǎn)上，識(shí)別每個(gè)SNP 位點(diǎn)；運(yùn)行Stackspopulations 程序，對(duì)檢測(cè)到的 SNP 位點(diǎn)數(shù)據(jù)進(jìn)行過(guò)濾。SNP 位點(diǎn)過(guò)濾標(biāo)準(zhǔn)為：次要基因型頻率（Minor allele frequency,MAF）＞0.05，完整性＞0.8，最大雜合度為1，最小測(cè)序深度為 6，無(wú)缺失數(shù)據(jù)。

1.2.5 SNP 位點(diǎn)注釋與染色體定位

利用SnpEff v4.3(http://pcingola.github.io/SnpEff/)軟件，對(duì)SNP 位點(diǎn)進(jìn)行注釋。運(yùn)行R 語(yǔ)言中的“CMplot”包，繪制SNP 位點(diǎn)在染色體上的分布圖；運(yùn)行R 語(yǔ)言中的“ggplot2”包，統(tǒng)計(jì)SNP 位點(diǎn)數(shù)目，分析位點(diǎn)數(shù)目與染色體長(zhǎng)度的相關(guān)性并繪制相關(guān)性曲線圖。

2 結(jié)果與分析

2.1 測(cè)序質(zhì)量

46 份黑楊派無(wú)性系種質(zhì)經(jīng)ddGBS 測(cè)序后，得到的原始圖像數(shù)據(jù)經(jīng)堿基識(shí)別分析轉(zhuǎn)化的原始測(cè)序序列數(shù)據(jù)量為108 Gb，每個(gè)樣本平均數(shù)據(jù)量為2.35 Gb。去除低質(zhì)量的數(shù)據(jù)后，獲得的高質(zhì)量數(shù)據(jù)的基本信息見表1。共產(chǎn)生218 156 221 條序列片斷，總讀長(zhǎng)數(shù)據(jù)為62 832 509 890 bp；樣本序列片段數(shù)目為2 795 062(NL1383)～7 994 127(ZKY)，平均片斷數(shù)目為4 742 526條。各樣本讀長(zhǎng)數(shù)目在802 475 700(NL1383)～2 295 941 115 bp(BKY) 范圍內(nèi)，平均讀長(zhǎng)為 365 924 128 bp。Q30 測(cè)序質(zhì)量值為93.32%～94.44%，平均值為93.87%。所有樣本的Q30 值均在93%以上，說(shuō)明堿基測(cè)序錯(cuò)誤率較低，獲得的測(cè)序數(shù)據(jù)合格。GC 含量為41.83%(Y706)～44.96%(Z13)，平均值為42.65%。GC 含量普遍不高，分布正常，完全符合測(cè)序的要求。

表1 楊樹GBS 測(cè)序數(shù)據(jù)Table 1 GBS sequencing data of 46 Populus

續(xù)表1Continuation of table 1

FastQC 軟件對(duì)質(zhì)控后數(shù)據(jù)進(jìn)行評(píng)估的結(jié)果，見圖1。圖1中的前145 bp 和后145 bp 分別代表測(cè)序片段兩端的質(zhì)量值分布情況。由圖1A 可知，堿基G、C 含量波動(dòng)范圍為18%～25%；從圖1B可以看出，片斷堿基測(cè)序的平均錯(cuò)誤率低于0.04%。采用Bowtie2 軟件，將所有樣本序列比對(duì)到毛果楊基因組上，比對(duì)率為64.25%(NLC70_2)～73.53%(NL120)，平均為69.89%（表1），比對(duì)率較高，這表明參試樣本基因組與毛果楊參考基因組有較高的相似性。位點(diǎn)測(cè)序深度為12.4(Z7)～24.7(81S)，平均值為18.4。

圖1 雙端堿基含量及測(cè)序質(zhì)量值Fig.1 Distribution and sequence qualities of double-terminal bases

2.2 SNP 變異類型及其在基因組中的分布特征

46 個(gè)黑楊派無(wú)性系種質(zhì)材料中，共鑒定出386 558 個(gè)SNP 變異位點(diǎn)。運(yùn)行Stacks-populations程序，獲得131 194 個(gè)特異性SNP 位點(diǎn)。SnpEff軟件對(duì)SNP 位點(diǎn)注釋的結(jié)果，見圖2A。6 種可能的堿基突變中，C/T(38 056)和A/G(37 694)突變最多，分別占總突變數(shù)量的29.00%和28.73%。轉(zhuǎn)換與顛換的比例為1.37∶1，其中，同義突變19 543 個(gè)，錯(cuò)義突變20 368 個(gè)，無(wú)義突變225 個(gè)。

分析基因組中SNP位點(diǎn)分布的結(jié)果，見圖2B。SNP 變異位點(diǎn)主要位于基因下游、上游和轉(zhuǎn)錄區(qū)，分別占基因組總SNPs 數(shù)量的29.48%、25.51%和19.06%；分布于外顯子、內(nèi)含子和基因間隔區(qū)的SNPs 占比分別為8.09%、8.10%和5.95%；3′和5′非翻譯區(qū)(Untranslated region，UTR)分布的SNPs數(shù)目較少(3.38%)，位于剪切體上的SNP 位點(diǎn)更少(0.43%)。

圖2 SNPs 注釋結(jié)果圖Fig 2 The annotation of SNPs in the Populus genome

2.3 SNP 位點(diǎn)在染色體上分布特征

運(yùn)行R 語(yǔ)言中的“CMplot”包，繪制的染色體上SNP 位點(diǎn)分布圖（圖3A）表明，98.41%(129 104)的SNP 位點(diǎn)能成功定位到19 條染色體上，每條染色體分布的SNP 位點(diǎn)數(shù)目為3 402（染色體19）～15 566（染色體1），平均分布密度為1/3 188 bp。染色體長(zhǎng)度與SNP 數(shù)目相關(guān)性分析的結(jié)果顯示（圖3B），染色體長(zhǎng)度與SNP 數(shù)目呈極顯著線性正相關(guān)（r=0.84，p＝1.97e-08）。這表明，SNP 位點(diǎn)在染色體上分布均勻。

圖3 染色體上SNP 位點(diǎn)分布情況分析Fig.3 Distribution of SNP loci on the chromosomes

3 討 論

限制性內(nèi)切酶的酶切是各種簡(jiǎn)化基因組測(cè)序的共同點(diǎn)和起點(diǎn)[4,14-16]。通過(guò)對(duì)酶切片段測(cè)序，降低基因組的復(fù)雜度。因此，選擇合適的限制性內(nèi)切酶是構(gòu)建高質(zhì)量文庫(kù)的關(guān)鍵。不同物種，選擇的限制性內(nèi)切酶種類及其組合不同，其酶切效率存在差異[23]。五堿基酶Ape KⅠ(G|CWGC)是Ⅱ型限制性內(nèi)切酶，DNA 轉(zhuǎn)座子區(qū)幾乎沒有識(shí)別位點(diǎn)，且部分甲基化敏感，常用于簡(jiǎn)化基因組測(cè)序。四堿基酶因其酶切的片段最小，也是常用的內(nèi)切酶。采用ApeKI 單酶切美洲山楊（P.tremuloides）基因組[22]，四堿基限制性內(nèi)切酶Hae Ⅲ(GG|CC)單酶切金花茶（Camellia nitidssima）[24]和古茶樹（Camellia sinensis）[25]基因組，效果較為理想。與單酶切相比，雙酶切更大程度上降低了基因組復(fù)雜度，增加了構(gòu)建文庫(kù)的一致性，基因組上的位點(diǎn)分布更加均勻[14]。采用六堿基酶PstI(CTGCA|G)和四堿基酶MspI(C|CGG)組合，開發(fā)出較理想的桉樹（Eucalypts）基因組SNP 位點(diǎn)[26]。本研究采用四堿基酶MseI 和TaqaI 雙酶切楊樹基因組，構(gòu)建的文庫(kù)質(zhì)量較高。

本研究基于ddGBS 技術(shù)對(duì)46 份歐美楊無(wú)性系種質(zhì)進(jìn)行測(cè)序，數(shù)據(jù)過(guò)濾后，共獲得131 194個(gè)SNPs 變異位點(diǎn)，SNPs 分布密度為1/3 308 bp，接近于美洲山楊(160 183)[22]和擬南芥Arabidopsis thaliana(1 /3300 bp)[27]基因組中SNPs 分布密度，遠(yuǎn)低玉米Zea mays(1/120 bp～1/60 bp)[28]、水稻Oryza sativa(1 /268 bp)[29]、大豆Glycine max(1/185～1/260 bp)[30]和黃麻Corchorus capsularis(1/172 bp)[31]基因組中的SNPs 分布密度。這可能與物種基因組結(jié)構(gòu)與大小的差異有關(guān)。楊樹基因組相對(duì)較小(434 Mb，https://www.ncbi.nlm.nih.gov/genome/?term=Populus.+trichocarpa)，大部分基因組序列屬于編碼序列，存在于非編碼序列中的SNP 會(huì)提高SNP 出現(xiàn)的頻率[32]。也可能與檢測(cè)樣本大小及測(cè)序長(zhǎng)度有關(guān)[33]。本研究中，挖掘的潛在SNP 變異位點(diǎn)數(shù)目(386 558)僅為毛果楊轉(zhuǎn)錄組(561 302)[18]的3/5，這是由于ddGBS 技術(shù)獲得的測(cè)序片段數(shù)相對(duì)較少，對(duì)基因組覆蓋率低[14]。

不同物種間轉(zhuǎn)換與顛換的比值存在很大差異，大多數(shù)DNA 片段中，轉(zhuǎn)換發(fā)生頻率比顛換發(fā)生頻率要高[34]。甜瓜Cucumis melo基于2b-RAD 測(cè)序的結(jié)果，轉(zhuǎn)換與顛換比率為2.15[35]；GBS 技術(shù)檢測(cè)雪茄煙草Nicotiana tabacum基因組SNP 位點(diǎn)發(fā)現(xiàn)，轉(zhuǎn)換與顛換比率為1.24[36]；青楊Populus cathayana轉(zhuǎn)錄組測(cè)序的結(jié)果，轉(zhuǎn)換與顛換比率為1.57[21]。本研究中，轉(zhuǎn)換與顛換比率為1.37。與產(chǎn)生這種差異的原因，可能與不同物種在其長(zhǎng)期進(jìn)化過(guò)程中承受的進(jìn)化選擇壓力不同有關(guān)[37]。轉(zhuǎn)換類型中，C/T 突變最多（29.00%），這可能與CG序列上SNP 變異位點(diǎn)出現(xiàn)的最為頻繁，胞嘧啶（Cytosine,C）常以甲基化形式存在，脫氨后成為胸腺嘧啶（Thymine,T）有關(guān)[38]。

注釋后的楊樹基因組中的SNPs，主要位于基因上下游和轉(zhuǎn)錄區(qū)，部分位于基因間區(qū)和內(nèi)含子區(qū)域上，8.09%的SNP位點(diǎn)位于外顯子功能區(qū)域上，說(shuō)明部分變異位點(diǎn)對(duì)蛋白翻譯有影響，這可能是楊樹基因組相對(duì)較小、大部分基因組序列屬于編碼序列。本研究中，出現(xiàn)在基因間區(qū)和內(nèi)含子上的SNP 位點(diǎn)數(shù)目相接近，這與基于轉(zhuǎn)錄組測(cè)序的結(jié)果基本一致[18]，98.41%(129 104)的SNP 位點(diǎn)能成功定位到19 條染色體上，表明研究的材料與毛果楊親緣關(guān)系很近。SNPs 數(shù)目與染色體長(zhǎng)度呈極顯著的線性正相關(guān)，說(shuō)明這些SNP 位點(diǎn)在染色體上分布均勻。

GBS 技術(shù)最顯著的優(yōu)點(diǎn)是，DNA 用量少，建庫(kù)效率高效，成本低廉，發(fā)掘的SNP 位點(diǎn)在基因組中分布均勻，比其他各種基于高通量測(cè)序平臺(tái)的基因分型技術(shù)具有成本競(jìng)爭(zhēng)力，克服2b-RAD片段過(guò)短，易受重復(fù)干擾的不足，會(huì)使育種工作者在不必開發(fā)任何標(biāo)記的前提下進(jìn)行物種分類、種質(zhì)鑒定、基因組選擇或分子輔助選擇育種[14]。但是，ddGBS 獲得的測(cè)序片段數(shù)較少，對(duì)基因組覆蓋度相對(duì)降低。在今后的研究中，將楊樹轉(zhuǎn)錄組數(shù)據(jù)與ddGBS 測(cè)序數(shù)據(jù)相結(jié)合，將會(huì)提高SNP位點(diǎn)的分布密度和基因定位效率。

4 結(jié) 論

本研究表明，采用ddGBS 能有效地挖掘楊樹高質(zhì)量SNP 位點(diǎn)信息。ddGBS 具有高通量、低成本等特點(diǎn)，今后可用GBS 技術(shù)建立涵蓋楊樹國(guó)家種質(zhì)資源庫(kù)所有種質(zhì)材料的簡(jiǎn)化基因組基因分型數(shù)據(jù)，構(gòu)建較為完整的種質(zhì)無(wú)性系DNA 指紋圖譜，為楊樹無(wú)性系種質(zhì)鑒定和新品種知識(shí)產(chǎn)權(quán)保護(hù)提供理論依據(jù)，并為進(jìn)一步開展楊樹種質(zhì)遺傳多樣性分析及重要經(jīng)濟(jì)性狀關(guān)聯(lián)分析，乃至楊樹分子輔助選擇育種提供了重要的遺傳信息。