單建偉，索海翠，王 麗，安 康，劉計濤，李成晨，白建明，李小波

（1.廣東省農(nóng)業(yè)科學(xué)院作物研究所/廣東省農(nóng)作物遺傳改良重點實驗室，廣東 廣州 510640；2.云南省農(nóng)業(yè)科學(xué)院經(jīng)濟(jì)作物研究所，云南 昆明 650000）

【研究意義】馬鈴薯為茄科（Solanaceae）茄屬（Solanum）馬鈴薯組（Petota）草本植物，原產(chǎn)于南美洲安第斯高原［1］。馬鈴薯的植物分類學(xué)是一個不斷更新和完善的過程，Spooner 等綜合前人有關(guān)馬鈴薯表型、分子標(biāo)記、系統(tǒng)發(fā)生、生殖生物學(xué)、倍性、核型等方面的研究將馬鈴薯組分類為107 個野生種和4 個栽培種［2］。對于栽培作物而言，通常在其起源地具有最豐富的遺傳資源和遺傳多樣性，而在其傳播過程中遺傳多樣性逐漸降低［3］。20 世紀(jì)20 年代，前蘇聯(lián)、英國、德國和美國學(xué)者先后赴南美洲考察和收集馬鈴薯野生種及栽培種種質(zhì)資源，為馬鈴薯的育種和遺傳學(xué)研究開創(chuàng)了新紀(jì)元［4］。當(dāng)今世界范圍內(nèi)廣泛種植的現(xiàn)代馬鈴薯品種（modern cultivar）是以馬鈴薯普通栽培種為基礎(chǔ)，通過與其他地方栽培種或野生種雜交利用了有利基因的產(chǎn)物，使其能夠適應(yīng)更寬泛的生態(tài)條件和地理區(qū)域，但目前廣泛種植的馬鈴薯品種遺傳多樣性依然很低［5］。廣泛收集國內(nèi)外馬鈴薯品種（系）資源，對其群體結(jié)構(gòu)和遺傳多樣性進(jìn)行分析評價，將為我國馬鈴薯遺傳育種提供理論指導(dǎo)。

【前人研究進(jìn)展】我國馬鈴薯的栽培歷史并不長［6］，對其育種起步也較晚，始于20 世紀(jì)30年代中末期，早期的馬鈴薯育種工作以引種為主，先后從英國、美國和前蘇聯(lián)等引進(jìn)一批育成品種和雜交組合。由于馬鈴薯適應(yīng)性廣、產(chǎn)量高，新中國成立以后為了盡快恢復(fù)農(nóng)業(yè)生產(chǎn)，國家有計劃地組織調(diào)運(yùn)種薯，擴(kuò)大馬鈴薯種植面積，全國開展了馬鈴薯育種協(xié)作，真正開始馬鈴薯育種工作，促進(jìn)了馬鈴薯在我國的迅猛發(fā)展，并相繼育成了一批有代表性的馬鈴薯品種，如東農(nóng)303、克新系列、隴薯系列等［7-8］。當(dāng)前，我國馬鈴薯種植面積和總產(chǎn)量均超過全球的1/4，居世界首位。由于我國馬鈴薯育種工作起步較晚，所用親本多引自歐洲和北美等國，對種質(zhì)資源特別是野生種和地方栽培種的收集、研究、創(chuàng)新不夠，因此我國栽培的馬鈴薯品種遺傳多樣性更低，許多品種之間同質(zhì)性高，部分品種存在同種異名現(xiàn)象，僅靠表型差異難以分辨。

【本研究切入點】為了指導(dǎo)馬鈴薯遺傳育種，有必要收集國內(nèi)外馬鈴薯品種（系）資源，并對其遺傳多樣性進(jìn)行評價?！緮M解決的關(guān)鍵問題】本研究對收集的國內(nèi)外馬鈴薯品種（系）資源，利用限制性酶切位點關(guān)聯(lián)DNA 測序（Restrictionsite Associated DNA sequencing，RADseq）中 的ddRADseq 技術(shù)對其進(jìn)行簡化基因組測序，結(jié)合群體結(jié)構(gòu)分析和遺傳多樣性分析對其進(jìn)行群體結(jié)構(gòu)和遺傳多樣性評價，開發(fā)馬鈴薯SNP-Panel，為馬鈴薯育種、品種鑒定提供理論依據(jù)。

1 材料與方法

1.1 試驗材料

供試材料為廣東省農(nóng)業(yè)科學(xué)院作物研究所馬鈴薯研究室收集的185 份馬鈴薯品種（系），包含育成品種以及一些高代品系（表1）。這些材料分別來自我國的華南、華北、西南、華中地區(qū)以及國外的秘魯和意大利。

表1 馬鈴薯樣品來源Table1 Sources of potato samples

1.2 試驗方法

1.2.1 基因組DNA 制備 取每個馬鈴薯品種（系）的嫩葉，利用液氮速凍后儲存于-80℃超低溫冰箱，分離樣品基因組DNA 后，利用瓊脂糖凝膠電泳檢測DNA 樣品完整性及有無RNA 污染，并利用分光光度計對DNA 樣品的質(zhì)量和濃度進(jìn)行定量檢測。

1.2.2 文庫構(gòu)建、質(zhì)控及測序 取約200 ng 基因組DNA 用限制性內(nèi)切酶MseⅠ和SacⅠ（紐英倫生物技術(shù)有限公司）按照產(chǎn)品說明書進(jìn)行完全酶切；酶切產(chǎn)物連接特異性接頭，然后對連接產(chǎn)物進(jìn)行純化；用高保真聚合酶KOD-Plus-Neo（東洋紡生物科技有限公司）按照產(chǎn)品說明書進(jìn)行PCR 擴(kuò)增富集；將所有擴(kuò)增產(chǎn)物進(jìn)行低壓電泳，切取300～400 bp 大小的電泳片段，用瓊脂糖凝膠回收試劑盒（QIAGEN，德國）進(jìn)行純化；用分光光度計對純化產(chǎn)物進(jìn)行初步定量，用安捷倫2100 生物分析儀對文庫插入片段大小進(jìn)行檢測；采用實時熒光定量PCR 手段對文庫的有效濃度進(jìn)行準(zhǔn)確定量，有效濃度＞2 nmol/L 為合格；最后用Illumina HiSeq 測序平臺進(jìn)行雙末端測序。

1.2.3 測序數(shù)據(jù)過濾 測序得到的原始數(shù)據(jù)被稱為raw read，含有測序接頭序列或低質(zhì)量的堿基。使用Cutadapt（Version 1.13）軟件去除raw read中的接頭序列，用Trimmomatic（Version 0.36）去除測序數(shù)據(jù)中的低質(zhì)量堿基，得到clean read 用于后續(xù)分析；測序read 的長度必須大于50 bp。

1.2.4 參考基因組比對 使用BWA（Version 0.7.15-r1140）軟件中的MEM 算法將測序數(shù)據(jù)與參考基因組（http://solanaceae.plantbiology.msu.edu/pgsc_download.shtml）進(jìn)行比對，得到SAM 格式的比對結(jié)果，然后使用 Samtools（Version 1.3.1）將其轉(zhuǎn)換為BAM 格式，使用Picard（Version 1.91）軟件中的SortSam 對BAM 文件中的序列進(jìn)行排序，將處理后的BAM 文件用于覆蓋度和覆蓋深度統(tǒng)計以及變異位點識別。

1.2.5 變異位點檢測及注釋 利用GATK（Version 3.7）軟件包中的HaplotypeCaller 模塊對所有樣品生成gvcf 格式文件，然后用GenotypeGVCFs 模塊對所有樣品進(jìn)行SNP 和 InDel 變異檢測。并按照測序深度≥5、所有樣品基因型缺失率≤50%、稀有等位基因頻率≥10%、雜合率≤50%、相對雜合率≤67%等標(biāo)準(zhǔn)對變異位點進(jìn)行篩選，得到的變異位點為有效變異位點。用ANNOVAR（Version 2016Feb1）對變異位點進(jìn)行注釋。依據(jù)變異位點在參考基因組上的位置信息以及參考基因組中的基因位置信息，可得到變異位點在基因組中的區(qū)域（基因間區(qū)、基因區(qū)或CDS 區(qū)等）信息。

1.2.6 馬鈴薯群體結(jié)構(gòu)分析 通過Structure 軟件對群體結(jié)構(gòu)進(jìn)行分析，結(jié)果以ΔK值變化圖和祖先堆疊圖展示，可直觀反映個體間的分類關(guān)系以及每個樣品的“混血”程度。采用PLINK（Version v1.90p）軟件對群體結(jié)構(gòu)進(jìn)行主成分分析（Principal Component Analysis，PCA）。PCA 散點圖中，兩個樣品距離越遠(yuǎn)，表明遺傳背景差異越大，而遺傳背景相似的個體在PCA 散點圖中聚類在一起。采用MEGA7（Version 7.0）軟件中的鄰接法（Neighbor-joining methods）構(gòu)建系統(tǒng)發(fā)生樹，并通過ggtree（Version 1.7.10）進(jìn)行可視化處理。基于分子標(biāo)記信息，利用Cervus 3.0 軟件計算多態(tài)性信息含量（Polymorphism Information Content，PIC）、觀測雜合度（Observed Heterozygosit，Ho）、期望雜合度（Expected Heterozygosity，He）等。

1.2.7 SNP 分子標(biāo)記開發(fā) 在綜合區(qū)分度和檢測成本的基礎(chǔ)上，結(jié)合SNP 標(biāo)記在染色體上的位置信息以及多態(tài)性，選擇120 個在染色體上分布相對均勻的SNP 標(biāo)記，分別在SNP 位點上下游開發(fā)PCR 引物，PCR 產(chǎn)物長度介于150～350 bp 之間。

2 結(jié)果與分析

2.1 變異位點檢測及注釋

通過測序最終獲得clean read 7.50×108條、2.04×1011個堿基，平均每個樣品4.05×106條clean read，1.10×109個堿基。為了檢測變異位點，將測序得到的clean read 與馬鈴薯基因組（http://solanaceae.plantbiology.msu.edu/pgsc_download.shtml）進(jìn)行比對，并進(jìn)行群體結(jié)構(gòu)和遺傳多樣性進(jìn)行分析。本研究在綜合考慮可行性和經(jīng)濟(jì)性的基礎(chǔ)上對測序覆蓋度進(jìn)行了控制，樣品平均測序覆蓋度為2.62%。利用GATK 軟件包共檢測到39 038 個有效變異位點，其中SNP 位點36 267個，InDel 位點2 771 個。SNP 位點在馬鈴薯12條染色體上的平均分布密度為50.02 個/Mb SNP位點，但在不同染色體上的的分布密度不均，其中10 號染色體上的SNP 位點密度最低（26.22 個/Mb），5 號染色體上SNP 位點密度最高（83.50個/Mb）；InDel 變異位點在12 條染色體上的分布密度介于1.97～6.05 個/Mb 之間，同樣10 號染色體分布密度最低，5 號染色體分布密度最高（表2）。同一條染色體不同區(qū)域的變異位點密度也存在較大差異，如1 號染色體10～20 Mb 范圍內(nèi)的變異位點密度最高，每100 kb 存在約90 個變異位點，而在70～90 Mb 范圍內(nèi)變異位點密度較低（圖1）。

表2 變異位點檢測統(tǒng)計結(jié)果Table 2 Detection and statistic results of variation sites

圖1 變異位點在染色體上的分布Fig.1 Distribution of variation sites on chromosomes

利用ANNOVAR 軟件結(jié)合變異位點在參考基因組和自身基因組上的位置信息對變異位點進(jìn)行注釋，結(jié)果（表3）顯示，共有26 697 個SNP 位點位于基因間區(qū)，占所有SNP 位點的73.65%；1 392 個SNP 位點位于基因轉(zhuǎn)錄起始位點上游1 kb 范圍內(nèi)，1 314 個SNP 位點位于轉(zhuǎn)錄終止位點下游1 kb 范圍內(nèi)，116 個SNP 位點位于一個基因轉(zhuǎn)錄起始位點上游1 kb 范圍內(nèi)并同時在另一個基因的轉(zhuǎn)錄終止位點下游1 kb 范圍內(nèi)；6 729個SNP 位點位于基因內(nèi)部，占所有SNP 位點的18.56%，其中3 172 個SNP 位點位于內(nèi)含子區(qū)域，2 605 個SNP 位點位于外顯子區(qū)；352 個SNP 位點位 于5′ UTR 區(qū)，597 個位 點位于3′ UTR 區(qū)；3 個SNP 位點位于可變轉(zhuǎn)錄本上。此外，共有1 761 個（63.57%）InDel 位點位于基因間區(qū)；308個InDel 位點位于基因編碼區(qū)上下游1 kb 范圍內(nèi)；25.31%的InDel 位點位于基因內(nèi)部，其中406 個位于內(nèi)含子內(nèi)，134 個位于外顯子內(nèi)；3′ UTR 和5′UTR 區(qū)域內(nèi)的InDel 位點數(shù)分別為101、56 個；4個InDel 位點位于可變轉(zhuǎn)錄本內(nèi)。

表3 注釋的變異位點分布Table 3 Distribution of annotated variation sites

2.3 馬鈴薯群體結(jié)構(gòu)分析

為了明確本研究所用馬鈴薯的群體結(jié)構(gòu)，通過PCA 法對篩選得到的有效SNP 標(biāo)記進(jìn)行降維處理，從中提取關(guān)鍵信息將分子標(biāo)記差異信息反映在二維坐標(biāo)圖上，坐標(biāo)軸為對樣品差異性解釋度最高的特征向量，兩個樣品遺傳背景差異越小，則樣品分子標(biāo)記位點信息越相似，反映在PCA 圖中的距離越近。PCA 散點圖顯示，第一特征向量（PC1）和第二特征向量（PC2）可以將本研究所用群體區(qū)分為G1 和G2 兩個亞群（圖2A）。利用Structure 軟件對群體結(jié)構(gòu)進(jìn)行分析，群體結(jié)構(gòu)分析的ΔK值變化趨勢圖顯示，當(dāng)K=2 時對應(yīng)的ΔK值最大，因而本研究群體的最佳亞群數(shù)為2，這一結(jié)果與PCA 分析結(jié)果一致，即將本研究群體分為兩個亞群較為合適（圖2B）。G2 亞群包含29 個品種（系）全部來自云南省農(nóng)業(yè)科學(xué)院，分別為S17-192-12、S17-114-4、S17-188-13、S17-119-2、S17-188-4、S17-114-2、S17-166-11、S17-147-12、S17-114-10、S17-1-1、S17-176-1、S17-173-97、S17-4-13、S17-201-1、S17-173-4、S17-147-23、S17-187-4、S17-1-9、S17-4-21、S17-173-30、S17-1-11、S17-129-3、S17-27-1、S17-8-1、S17-1-15、S17-181-12、S17-172-2、S17-181-1、紫云一號-1；剩余的156 個品種（系）構(gòu)成G1 亞群，包括來自我國華南、華北、華中、西南、西北等地的品種（系），以及引自秘魯、意大利等國的品種（系）。祖先堆疊圖（圖2C）顯示，兩個亞群間存在明顯基因交流，部分樣本存在混血。群體系統(tǒng)進(jìn)化分析結(jié)果（圖2D）表明，所有G2 亞群中的29 個品種（系）聚類在一起且其分枝離圓心的距離最長，表明G2 亞群中的個體親緣關(guān)系較近，同時G2 亞群中的個體在進(jìn)化上積累了更多變異。

圖2 馬鈴薯群體結(jié)構(gòu)分析Fig.2 Analysis of potato population structure

PIC＞0.5 為高度多態(tài)性［9］，而本研究群體的平均PIC=0.3107，表明群體遺傳多樣性中等偏低。同時，群體的期望雜合度He=0.3932，觀測雜合度Ho=0.1852，群體自交系數(shù)Fis=0.553，也表明群體遺傳多樣性偏低。

2.4 SNP 分子標(biāo)記開發(fā)

靶向測序（target-seq）技術(shù)通過設(shè)計多重PCR 引物對目標(biāo)區(qū)域DNA 片段進(jìn)行特異性擴(kuò)增，再通過雙端標(biāo)簽區(qū)分不同樣本，混樣建庫，進(jìn)行多樣本多位點高通量測序分析，具有信息量大、針對性強(qiáng)、低成本、效率極高、精準(zhǔn)度高等特點，被廣泛應(yīng)用于分子標(biāo)記輔助選擇育種（Marker-Assisted Selection，MAS）、數(shù)量性狀位點（Quantitative Trait Locus，QTL）精細(xì)定位、品種指紋分析、品種鑒定、基因編輯位點檢測等領(lǐng)域。為了便于后續(xù)開展馬鈴薯親緣關(guān)系鑒定、品種鑒別、群體遺傳分析以及分子標(biāo)記輔助選擇育種，我們開發(fā)了包含120 個馬鈴薯SNP 標(biāo)記的SNP-Panel，這些SNP 標(biāo)記均勻分布在馬鈴薯的12 條染色體上（圖3），并針對這120 個SNP 標(biāo)記設(shè)計了分離目標(biāo)DNA 區(qū)域的特異引物用于后續(xù)多樣本高通量SNP標(biāo)記檢測。為了驗證SNP 標(biāo)記的有效性，選取了60 個SNP 標(biāo)記，利用與所選SNP 標(biāo)記對應(yīng)的引物對46 份馬鈴薯品種（系）進(jìn)行了靶向高通量測序，結(jié)果表明，60 個SNP 標(biāo)記能夠準(zhǔn)確將46 份馬鈴薯品種（系）區(qū)分開。理論上，所選標(biāo)記越多，越能區(qū)分更多的品種（系），因此需要在區(qū)分度和檢測成本之間尋求平衡。因此在綜合區(qū)分度和檢測成本的基礎(chǔ)上，本研究確定了SNP 標(biāo)記的個數(shù)為120。

圖3 SNP-Panel 在馬鈴薯染色體上的分布Fig.3 Distribution of SNP-Panel on potato chromosomes

3 討論

RADseq 是在第二代測序技術(shù)基礎(chǔ)上發(fā)展而來的一項基于全基因組酶切位點的簡化基因組測序技術(shù)，該技術(shù)具有檢測速度快、高通量、成本低等優(yōu)點，可以檢測出上千萬的SNP 位點，能夠滿足群體結(jié)構(gòu)分析、遺傳圖譜構(gòu)建和QTL 定位等研究需求［10-11］。依據(jù)所使用的限制性內(nèi)切酶的數(shù)量和種類，RADseq 技術(shù)可細(xì)分為Original RADseq、GBS、2b-RADseq、ezRAD、ddRADseq等技術(shù)方法。ddRADseq 技術(shù)由于檢測準(zhǔn)確、成本低、SNP 位點豐富等優(yōu)點得到廣泛應(yīng)用［12-13］。本研究采用ddRADseq 技術(shù)對185 個馬鈴薯品種（系）進(jìn)行了簡化基因組測序分析。基于變異位點在各樣品中的測序深度均不小于5、在所有樣品中的缺失率不超過50%、低頻等位基因頻率不低于10%、雜合率不超過50%、相對雜合率不超過67%等標(biāo)準(zhǔn)共獲得39 038 個高質(zhì)量變異位點，其中SNP 位點36 267 個，InDel 位點2 771 個。高質(zhì)量變異位點為后續(xù)馬鈴薯群體結(jié)構(gòu)分析、SNP 標(biāo)記Panel 開發(fā)、分子標(biāo)記輔助選擇育種、QTL 位點定位奠定了基礎(chǔ)。

物種的遺傳多樣性越高對生存環(huán)境的適應(yīng)性越強(qiáng)，抵御生物脅迫、非生物脅迫的能力越強(qiáng)，易于拓展其生存分布范圍。本研究通過對來自不同地方的185 份馬鈴薯品種（系）的遺傳多樣性分析發(fā)現(xiàn)，多態(tài)性信息含量PIC=0.3107，處于中等偏低水平；群體期望雜合度He=0.3932，而觀測雜合度Ho=0.1852，即Ho＜He。這一結(jié)果表明這些馬鈴薯品種（系）的遺傳多樣性較低，即不同品種（系）間的同質(zhì)性較高，這與前人研究結(jié)果較為一致［14-18］。為了提升馬鈴薯品種的遺傳多樣性，在未來的馬鈴薯育種工作中應(yīng)加強(qiáng)對馬鈴薯種質(zhì)資源的收集和創(chuàng)新，特別是對馬鈴薯野生種的利用。馬鈴薯種質(zhì)資源豐富，包含107個野生種和4 個栽培種。原始栽培種和野生種含有眾多調(diào)控馬鈴薯晚疫病、瘡痂病、青枯病、病毒病、黑脛病抗性的等位基因，也含有能夠調(diào)控馬鈴薯高淀粉和蛋白質(zhì)含量、低還原糖含量、耐霜凍、休眠期、早熟等品質(zhì)和農(nóng)藝性狀的優(yōu)良等位基因［8］。加速對馬鈴薯種質(zhì)資源的收集、評價、研究、創(chuàng)新和利用，有利于突破當(dāng)前馬鈴薯育種的瓶頸［19-22］。

本研究結(jié)果表明，供試馬鈴薯種（系）間的遺傳多樣性較低，工作中也發(fā)現(xiàn)一些品種間表型差異很小，如果不借助分子標(biāo)記技術(shù)，僅從表型差異很難區(qū)分。鑒于此，本研究還開發(fā)了包含120 個馬鈴薯SNP-Panel，這些SNP 標(biāo)記均勻覆蓋馬鈴薯12 條染色體，并針對這些SNP 標(biāo)記設(shè)計了分離目標(biāo)DNA 區(qū)域的特異引物，以用于后續(xù)進(jìn)行多樣本高通量靶向SNP 標(biāo)記檢測，為后續(xù)開展馬鈴薯親緣關(guān)系鑒定、品種鑒別、群體遺傳分析以及分子標(biāo)記輔助選擇育種奠定基礎(chǔ)。

4 結(jié)論

本研究對收集的185 份馬鈴薯品種（系）種質(zhì)資源進(jìn)行簡化基因組測序，鑒定了39 038個變異位點，其中SNP 位點36 267 個，InDel位點2 771 個；PCA 分析和Structure 分析均將供試的185 份馬鈴薯品種（系）劃分為兩個亞類；遺傳多樣性分析表明，群體多態(tài)性信息含量PIC=0.3107，期望雜合度He=0.3932，觀測雜合度Ho=0.1852，自交系數(shù)Fis=0.553，群體遺傳多樣性較低；開發(fā)了覆蓋馬鈴薯12 條染色體且包含120 個SNP 標(biāo)記的SNP-Panel。本研究為后續(xù)進(jìn)行馬鈴薯遺傳多樣性研究、分子標(biāo)記輔助選擇育種、分子身份證開發(fā)、品種鑒定、遺傳圖譜構(gòu)建奠定了一定基礎(chǔ)。