基于全基因組SNP分子標(biāo)記分析青海云杉遺傳結(jié)構(gòu)

張宏斌 呂 東 趙 明 趙 祜 趙興鵬 李 偉

（1. 甘肅省祁連山水源涵養(yǎng)林研究院，張掖 734000；2. 祁連山特有植物繁育及推廣國家地方聯(lián)合工程研究中心，張掖 734000；3. 北京林業(yè)大學(xué)生物科學(xué)與技術(shù)學(xué)院，北京 100083）

青海云杉（）是松科（Pinaceae）云杉屬（）常綠喬木，為我國特有樹種，具有涵養(yǎng)水源、保持水土的作用，是祁連山區(qū)主要的水源涵養(yǎng)林和造林更新樹種。具有極高的生態(tài)價值；木材質(zhì)量好，可供橋梁、家具及造紙等用材，是重要的用材樹種，具有重要的經(jīng)濟(jì)價值；四季常青、外觀優(yōu)美，被廣泛應(yīng)用于城市綠化、園林栽植中，具有觀賞價值。幾十年來，我國林業(yè)工作人員對青海云杉的進(jìn)化演變、空間分布、水分利用、氣候響應(yīng)等進(jìn)行了一系列研究；并開展了青海云杉母樹林與種子園的營建、表型變異、遺傳變異及子代測定等方面的遺傳改良工作，取得了一些成果。

青海云杉表型變異豐富，在分子水平上也具有較高的遺傳變異水平，在群體間和群體內(nèi)都存在著廣泛的變異，具有豐富的遺傳多樣性，因此，進(jìn)一步挖掘青海云杉群體變異信息，對青海云杉種質(zhì)資源保護(hù)和遺傳改良工作的科學(xué)開展具有重要意義。但青海云杉在改良過程中存在許多亟待解決的問題。第一，由于選育周期長，導(dǎo)致良種選育慢良種少，青海云杉遺傳改良進(jìn)程緩慢；已有許多針葉樹種相繼進(jìn)入第2 代種子園的建設(shè)當(dāng)中，但由于技術(shù)、科技投入不足，青海云杉的遺傳改良還停留在初級種子園階段。第二，由于全基因組較大，青海云杉遺傳改良基礎(chǔ)薄弱，在分子水平上的研究較少，缺少基因組信息，分子標(biāo)記多來自挪威云杉（），現(xiàn)有引物通用性差，數(shù)量少，導(dǎo)致標(biāo)記相關(guān)分析準(zhǔn)確性較低，缺少青海云杉特有的分子標(biāo)記。近年來，祁連山區(qū)生態(tài)環(huán)境逐步惡化，青海云杉作為主要造林綠化樹種，對其良種質(zhì)量和數(shù)量的要求逐漸升高，迫切需求青海云杉良種升級換代，亟待開展青海云杉的高世代親本選擇與遺傳評價等基礎(chǔ)研究，為青海云杉高遺傳改良層次提升奠定基礎(chǔ)。

傳統(tǒng)的選擇優(yōu)良單株和雜交的育種方式在青海云杉育種過程中所需周期太長，利用目標(biāo)性狀與分子標(biāo)記的關(guān)聯(lián)分析，篩選出與目標(biāo)性狀緊密關(guān)聯(lián)的標(biāo)記，為分子標(biāo)記輔助選擇提供依據(jù)，可以極大縮減青海云杉遺傳改良周期，提升青海云杉遺傳改良水平。因此，本研究采用SLAF-seq技術(shù)，開發(fā)全基因組水平上的青海云杉SNP 分子標(biāo)記，為青海云杉種質(zhì)資源的研究、遺傳連鎖圖譜構(gòu)建以及QTL定位等方面的研究奠定基礎(chǔ)。

1 材料與方法

1.1 試驗(yàn)材料

試驗(yàn)材料選取張掖龍渠青海云杉無性系種子園中的106 個青海云杉親本無性系。種子園建園材料來源于祁連山區(qū)的11 個林場（見表1），2019年采集每個無性系單株當(dāng)年生健康針葉于自封袋中，分別編號，及時存儲于液氮中，回到室內(nèi)，放冰箱里冷藏備用。

表1 試驗(yàn)材料信息Table 1 Experimental material information

1.2 基因組DNA提取與檢測

基因組DNA 的提取采用改良CTAB 法，提取的青海云杉基因組DNA 質(zhì)量濃度不低于20 ng·μL，總量不低于4 μg，保證質(zhì)量符合建庫要求。

1.3 SLAF文庫構(gòu)建和高通量測序

由于目前缺少青海云杉基因組序列的相關(guān)信息，選取挪威云杉基因組作為參考基因組（ftp：//plantgenie. org/Data/ConGenIE/Picea_abies/v1.0/），組裝出的基因組大小為12 G，GC 含量為37.88%。根據(jù)挪威云杉參考基因組，進(jìn)行電子酶切預(yù)測，確定最適酶切方案。

對質(zhì)量合格的各樣品基因組DNA 進(jìn)行酶切，得到的酶切片段（SLAF 標(biāo)簽）進(jìn)行3′端加A 處理、連接測序接頭后，進(jìn)行PCR 擴(kuò)增，通過純化、混樣、切膠等步驟選取目的片段，文庫質(zhì)檢合格后用Il?lumina HiSeq 進(jìn)行雙端測序。選用水稻日本晴（ssp.）作為試驗(yàn)建庫的對照，基因組數(shù)據(jù)參考來源于：http：//rapdb.dna.affrc.go.jp。由北京百邁克生物科技有限公司完成SLAF文庫的構(gòu)建以及測序。

1.4 SLAF測序數(shù)據(jù)分析和SNP位點(diǎn)篩選

對獲得的Raw Reads 進(jìn)行質(zhì)控，得到Clean Reads。將得到的reads利用BWA軟件比對到參考基因組上，統(tǒng)計(jì)SLAF 標(biāo)簽和多態(tài)性SLAF 標(biāo)簽的數(shù)量。SNP 分子標(biāo)記采用GATK 和samtools 兩種軟件共同開發(fā)，以次要基因型頻率（MAF）＞0.05、完整性＞0.8 為篩選標(biāo)準(zhǔn)，兩個軟件開發(fā)SNP標(biāo)記的交集作為最終可靠的SNP標(biāo)記。

1.5 遺傳結(jié)構(gòu)分析

根據(jù)開發(fā)的青海云杉高質(zhì)量SNP 標(biāo)記，采用MEGA X軟件，基于鄰接法（neighbor-joining），構(gòu)建各青海云杉樣品的系統(tǒng)發(fā)育樹。利用admixture軟件，研究青海云杉群體遺傳結(jié)構(gòu)，預(yù)先設(shè)定亞群數(shù)目（）為1～10，進(jìn)行聚類，并對聚類結(jié)果計(jì)算交叉驗(yàn)證錯誤率，依據(jù)其谷值確定最優(yōu)分群數(shù)。并且利用EIGENSOFT 軟件，進(jìn)行主成分分析，得到樣品的聚類情況。

2 結(jié)果與分析

2.1 試驗(yàn)建庫評估

根據(jù)參考基因組的電子酶切預(yù)測，最終確定限制性內(nèi)切酶為HinCII，長度在314～394 bp的酶切片段序列定義為SLAF標(biāo)簽，預(yù)測將獲得506 448個SLAF 標(biāo)簽。本研究中水稻日本晴數(shù)據(jù)的酶切效率高達(dá)98.12%，表明酶切反應(yīng)正常。通過SOAP軟件將測序的水稻日本晴reads 與其參考基因組進(jìn)行比對，顯示青海云杉基因組DNA構(gòu)建的SLAF文庫雙端比對效率為95.33%，說明青海云杉基因組DNA構(gòu)建的SLAF建庫正常。

2.2 測序數(shù)據(jù)評估

通過SLAF-seq 測序后，對照組水稻日本晴測序獲得1.48 Mb reads 的數(shù)據(jù)量，106 個青海云杉樣品共獲得1 375.57 Mb reads 數(shù)據(jù)。對照水稻日本晴reads 數(shù)為1 475 338，GC 含量為40.42%，Q30 為92.69%。各樣品的測序結(jié)果（見表2），所獲得的reads 數(shù)目范圍為4 411 192～30 427 534，其中，7 號所獲得的數(shù)據(jù)量最多，36 號獲得的數(shù)據(jù)量最少。GC 含量范圍在39.19%～43.62%，平均GC 含量為40.27%，測序質(zhì)量值Q30 范圍在90.97%～97.82%，平均Q30 為95.51%。本研究中所測序列的Q30 較高，表明堿基測序錯誤率低，說明測序質(zhì)量高，可以用作下一步的數(shù)據(jù)分析。

表2 青海云杉測序結(jié)果Table 2 Sequencing result of P.crassifolia

2.3 SLAF標(biāo)簽與SNP標(biāo)記開發(fā)

通過對測序數(shù)據(jù)的處理分析，本研究共開發(fā)4 058 883個SLAF標(biāo)簽，平均每個青海云杉樣本開發(fā)了389 459 個SLAF 標(biāo)簽，其中多態(tài)性SLAF 標(biāo)簽有1 573 899 個，占SLAF 標(biāo)簽總量的38.78%，標(biāo)簽的平均測序深度為21.21×（見表3）。共開發(fā)12 275 765個青海云杉SNP標(biāo)記，各青海云杉樣本的SNP數(shù)量為1 890 934～4 487 841，各樣品檢測到的SNP完整度為15.40%～36.56%，SNP的雜合率為5.41%～10.99%（見表4）。

表3 青海云杉SLAF 標(biāo)簽統(tǒng)計(jì)結(jié)果Table 3 SLAF tag statistics of P.crassifolia

表4 青海云杉SNP 信息統(tǒng)計(jì)結(jié)果Table 4 SNP information statistics of P.crassifolia

2.4 群體遺傳結(jié)構(gòu)

利用已開發(fā)的高質(zhì)量青海云杉SNP 標(biāo)記，構(gòu)建了106 個青海云杉的系統(tǒng)發(fā)育樹，結(jié)果如圖1所示，大致可分為6 小組，分別具有5、19、15、25、21、21個青海云杉無性系，來自多數(shù)種源的青海云杉在各組中分布比較均勻，不同種源的青海云杉多聚為一類，但每一類在不同種源的占比差別較大。

圖1 106個青海云杉材料的系統(tǒng)進(jìn)化樹Fig.1 The phylogenetic tree of 106 P.crassifolia materials

基于SNP 標(biāo)記，利用admixture 軟件分析了106 個青海云杉無性系群體結(jié)構(gòu)，結(jié)果表明當(dāng)=1時，交叉驗(yàn)證錯誤率最低（圖2：A～B），說明這些無性系來源于同一個祖先的可能性較大。

通過主成分分析，絕大多數(shù)來自不同種源的青海云杉形成一簇，只有小部分的青海云杉無性系分散分布，并未出現(xiàn)顯著的分化（圖2：C～E），表明無性系間親緣關(guān)系相近。

圖2 青海云杉遺傳結(jié)構(gòu)分析A.青海云杉樣品聚類結(jié)果；B.K值的交叉驗(yàn)證錯誤率；C～E.PCA聚類Fig.2 Analysis of genetic structure of P.crassifolia A.Clustering results of P.crassifolia；B.Cross validation error rates of K-values；C?E.Principal components analysis

3 討論

針葉樹種的基因組龐大復(fù)雜，目前僅有挪威云杉、白云杉、火炬松進(jìn)行了基因組測序，利用傳統(tǒng)的方法進(jìn)行分子標(biāo)記的開發(fā)，工作量大且復(fù)雜，花費(fèi)時間長，得到標(biāo)記的數(shù)量不能滿足試驗(yàn)的要求。袁行棟從已發(fā)表的挪威云杉60 對SSR 引物中篩選出了10 對具有多態(tài)性的SSR 引物，用作青海云杉遺傳多樣性分析，但分子標(biāo)記的數(shù)量少，引物之間的通用性較差，導(dǎo)致利用分子標(biāo)記的相關(guān)分析準(zhǔn)確性較低。

3.1 青海云杉全基因組SNP分子標(biāo)記開發(fā)

單核苷酸多態(tài)性（SNP）分子標(biāo)記是指DNA 序列發(fā)生多態(tài)性的變化，是由單個核苷酸發(fā)生缺失、插入、轉(zhuǎn)換等變異導(dǎo)致的，由于其在基因組中數(shù)量眾多，且存在范圍廣闊，在遺傳圖譜的構(gòu)建、數(shù)量性狀定位分析和全基因組關(guān)聯(lián)分析中被廣泛應(yīng)用。SLAF-seq 技術(shù)是一種快速開發(fā)大量SNP 標(biāo)記的技術(shù)，對于物種的要求較低，沒有其基因組信息也可以順利開發(fā)出大量SNP 標(biāo)記，對于未知基因組物種全基因組分子標(biāo)記的開發(fā)是非常有效的方法，在動植物中已經(jīng)得到了大量的使用。在針葉樹中，段紅靜通過SLAF-seq技術(shù)，首次在杉木群體中開發(fā)出166 646 個SNP 標(biāo)記。Bai 等采用SLAF-seq 技術(shù)對廣東馬尾松種質(zhì)資源進(jìn)行分析，在599 164個SLAF多態(tài)性標(biāo)簽中，共鑒定出471 660個SNP標(biāo)記。董明亮利用SLAF-seq技術(shù)在華北落葉松全基因組范圍內(nèi)挖掘SNP 標(biāo)記，共獲得6 323 943 個SLAF 位 點(diǎn)，檢 測 到324 352 個SNP 標(biāo)記，其中有122 785 個呈現(xiàn)多態(tài)性。王飛等利用SLAF-seq 技術(shù)開發(fā)華山松SNP 標(biāo)記，共獲得SLAF標(biāo)簽12 952 676 個，多態(tài)性SLAF 標(biāo)簽為1 456 486個，開發(fā)了3 469 074個群體SNP標(biāo)記。

在本研究中，首次將SLAF-seq 技術(shù)應(yīng)用于青海云杉群體，共獲得4 058 883 個SLAF 標(biāo)簽，開發(fā)了12 275 765個青海云杉SNP標(biāo)記，獲得了數(shù)量大且分布廣泛的分子標(biāo)記，并基于開發(fā)的分子標(biāo)記，分析了青海云杉無性系群體結(jié)構(gòu)，結(jié)果表明不同種源的青海云杉來源于同一個亞群的可能性較大。青海云杉SNP 標(biāo)記的開發(fā)為后續(xù)分子育種工作奠定了良好的基礎(chǔ)，可為遺傳評價、遺傳圖譜的構(gòu)建以及目標(biāo)性狀的關(guān)聯(lián)分析提供數(shù)據(jù)支撐。

3.2 青海云杉種子園遺傳結(jié)構(gòu)

種子園是按設(shè)計(jì)要求營建的，以優(yōu)樹的無性系或家系為材料，以生產(chǎn)優(yōu)良遺傳品質(zhì)的種子為目的營造的特種人工林。青海云杉種子園的研究包括林木生長，繁殖特性等。對于遺傳生物學(xué)形狀也有一些研究，呂東等以與本研究相同的材料進(jìn)行無性系結(jié)實(shí)性狀遺傳變異研究，表明青海云杉無性系間的結(jié)實(shí)能力差異明顯，其無性系重復(fù)力較高，球果主要集中在樹冠中上部的冠層分枝上。袁行棟2016年對其進(jìn)行了遺傳多樣性研究，表明青海云杉無性系較高的遺傳多樣性水平。本研究在前人的基礎(chǔ)上進(jìn)行，與之前遺傳多樣性水平較高的結(jié)論一致。

家系來源的11個林區(qū)基本涵蓋了青海云杉的自然分布區(qū)，結(jié)合地理分布分析。連城林區(qū)（LC）位于所有種源地的最東南方向，其15 個無性系分布在第Ⅱ、Ⅲ、Ⅴ、Ⅵ分組中，而僅Ⅴ組和Ⅵ組又多達(dá)11 個。大河口林區(qū)位于各林區(qū)的中間位置，其9 個無性系在5 個組分里面都有分布。隆暢河林區(qū)在最西北角上，10 個無性系中有7 個屬于第Ⅳ組分，其余Ⅲ、Ⅴ、Ⅵ組分各有1 個?？傮w來說，靠中間的種源地青海云杉在各組中分布相對比較均勻一般能聚為一類，靠邊緣的種源地雖然仍有多個組分的遺傳物質(zhì)來源，但都不能涵蓋所有組分且都有特定組分占大比例。青海云杉無性系群體結(jié)構(gòu)未出現(xiàn)分層情況，這些無性系來源于同一個祖先的可能性較大。

張掖龍渠青海云杉無性系種子園具有較高的遺傳多樣性水平，可為高世代種子園的建設(shè)提供豐富的原材料。該種子園建園初期包括了大量的無性系，擁有豐富的育種資源，但之后沒有連續(xù)補(bǔ)充新材料，良種基地發(fā)展的物質(zhì)基礎(chǔ)不夠充足。沒有充分重視遺傳測定工作，對收集的資源所做的作遺傳測定連續(xù)性不足，基地的更新?lián)Q代和發(fā)展的理論依據(jù)薄弱。建議今后應(yīng)不斷補(bǔ)充新的育種資源，對資源進(jìn)行有計(jì)劃的遺傳測定，同時可以采用多種管理技術(shù)措施；以充實(shí)基地發(fā)展的物質(zhì)基礎(chǔ)，為基地的更新?lián)Q代提供足夠的理論依據(jù)。

4 結(jié)論

通過SLAF-seq 技術(shù)，共開發(fā)了4 058 883 個SLAF標(biāo)簽，其中多態(tài)性SLAF標(biāo)簽有1 573 899個，占SLAF標(biāo)簽總量的38.78%，標(biāo)簽的平均測序深度為21.21×，共開發(fā)了12 275 765 個SNP 標(biāo)記，為今后遺傳多樣性的分析、遺傳圖譜的構(gòu)建等提供了基礎(chǔ)數(shù)據(jù)。

對青海云杉無性系遺傳結(jié)構(gòu)進(jìn)行分析，根據(jù)系統(tǒng)發(fā)育樹的結(jié)果，青海云杉可分為6 小組，來自不同種源的青海云杉在各組中分布整體比較均勻，但相距較遠(yuǎn)的種源之間存在一定差別。群體結(jié)構(gòu)分析中，當(dāng)=1 時，交叉驗(yàn)證錯誤率最低，說明這些無性系來源于同一個祖先的可能性較大，主成分分析顯示大部分來自不同種源的青海云杉無性系聚為一類，說明不同種源的青海云杉無性系間的親緣關(guān)系較近，青海云杉無性系群體結(jié)構(gòu)未出現(xiàn)分層情況，這些無性系來源于同一個祖先的可能性較大。