基于GBS-seq的青藏扁蓿豆SNP標記開發(fā)

周曉楠 徐金青 雷雨晴 王海慶

（1.中國科學院西北高原生物研究所 中國科學院高原生物適應與進化重點實驗室，西寧810001；2.中國科學院大學，北京100049；3.青海省作物分子育種重點實驗室，西寧810001）

青藏扁蓿豆（Medicago archiducis-nicolai）是豆科（Fabaceae）苜蓿屬（Medicago）多年生異花授粉二倍體植物［1］，分布于青藏高原及其毗鄰的高海拔高寒地區(qū)，對寒冷、干旱等極端環(huán)境具有極強的適應性［2-3］，是青藏高原天然高寒草地上具有馴化潛力的野生牧草種質(zhì)資源。遺傳多樣性是種內(nèi)個體間或一個群體內(nèi)不同個體遺傳變異的總和，研究物種的遺傳多樣性不僅可揭示該物種群體間的遺傳結(jié)構(gòu)、進化關(guān)系，以及與環(huán)境及地理分布之間的相關(guān)性［4］，也可以了解物種的遺傳背景及親緣關(guān)系，對種質(zhì)資源的評價利用具有重要意義。

目前，核基因組（ITS［3］）、葉綠體基因組（psbA-trnH［3］、trnL-trnF［5］）標記以及簡單重復序列（EST-SSR［6-7］）等分子標記技術(shù)已被應用于青藏扁蓿豆群體遺傳多樣性研究。利用ITS和psbA-trnH遺傳標記對不同海拔的青藏扁蓿豆野生群體進行遺傳分析，發(fā)現(xiàn)遺傳分化程度和采樣點海拔高度差之間存在顯著的相關(guān)性，表明海拔高度引起的遺傳隔離是產(chǎn)生群體間遺傳分化的主要原因?；贓ST-SSR標記進行的遺傳多樣性分析，發(fā)現(xiàn)青藏扁蓿豆群體間遺傳距離與地理距離和海拔差異之間均存在極顯著的相關(guān)性，同時也提示由海拔差異導致的局部氣候環(huán)境的異質(zhì)性可能對青藏扁蓿豆野生群體的遺傳多樣性產(chǎn)生影響，并使之獲得對局部極端生存環(huán)境的適應。但上述分子標記技術(shù)密度較低、多態(tài)性較低且操作復雜，具有一定的局限性。因此，為了進一步解析青藏扁蓿豆對極端環(huán)境的適應機制，評價青藏扁蓿豆種質(zhì)資源，需要對青藏扁蓿豆開發(fā)更高效的分子標記類型。

單核苷酸多態(tài)性標記（single nucleotide polymorphism，SNP）作為第三代分子標記，已在群體選擇分析［8］、重要性狀基因定位［9］、人類群體擴張研究［10］以及動植物遺傳連鎖圖譜構(gòu)建［11］等相關(guān)研究中廣泛應用，是研究物種遺傳變異的理想分子標記［12］。與以SSR為代表的二代分子標記相比，具有易于實現(xiàn)自動化分析、多態(tài)性高及密度高等優(yōu)點。但SNP標記由于前期測序階段的高成本，限制了其大規(guī)模的開發(fā)利用，GBS（genotyping-by-sequencing）是基于二代測序技術(shù)分型SNP的一種簡化基因組測序技術(shù)，是一種高效、簡單、低成本的基因分型方法，對了解種質(zhì)資源的遺傳背景和系統(tǒng)進化，研究群體遺傳多樣性及遺傳結(jié)構(gòu)具有重要意義［13］。

本研究以5個群體共80份野生青藏扁蓿豆為材料，采用GBS技術(shù)測序，進而開發(fā)SNP標記，基于開發(fā)的SNP標記，初步進行遺傳多樣性及遺傳結(jié)構(gòu)分析，為后續(xù)更加全面解析青藏扁蓿豆在青藏高原極端環(huán)境下的適應機制提供數(shù)據(jù)支持，同時也為青藏扁蓿豆種質(zhì)資源評價奠定基礎(chǔ)。

1 材料與方法

1.1 材料

80份供試材料于2017年采集自青海?。ū?），采樣群體海拔分布于2 462-3 311 m，每個群體采集16個個體，采集的單個個體放入信封袋用硅膠干燥保存?zhèn)溆谩?/p>

表1 青藏扁蓿豆材料來源Table 1 Sources of M.archiducis-nicolai materials

1.2 方法

1.2.1 基因組DNA的提取及純化 取約0.05 g干葉用組織研磨儀進行研磨，利用改良的CTAB法［14］提取總DNA，提取緩沖液中加入少量的聚乙烯吡咯烷酮（PVP）粉末和β-巰基乙醇（14.4 mol/L），水合時間延長至5 h，期間每隔10 min震蕩混勻一次，使干葉細胞充分裂解，異丙醇沉淀時間延長至90 min，提高DNA產(chǎn)率。DNA產(chǎn)物經(jīng)0.8%瓊脂糖凝膠電泳檢測，確認目的條帶清晰明亮且完整后，利用TIANGEN Purification Kit柱式膠回收試劑盒（增強型）進行純化，純化后的產(chǎn)物經(jīng)過濃度及純度檢測達到測序要求后，-20℃保存?zhèn)溆谩?/p>

1.2.2 GBS測序文庫的構(gòu)建 提取的DNA樣品送往北京奧維森基因科技有限公司，參照Elshire等［15］建庫技術(shù)（取樣-酶切-接頭連接-混樣-PCR擴增-純化）建庫，首先應用限制性內(nèi)切酶ApeKⅠ對基因組DNA進行酶切，加上帶有條形碼的接頭后，對每個樣品進行擴增，然后對樣品進行混合，電泳回收375-400 bp區(qū)間的DNA條帶，純化后產(chǎn)物用于測序，測序反應在Illumina HiSeqPE150測序平臺上進行雙末端150 bp的測序。

1.2.3 SNP開發(fā) 將測序得到的原始數(shù)據(jù)按照如下標準過濾：剔除帶接頭（adapter）的reads pair；單端read中含有N的含量超過該條read長度比例的10%時，去除此對paired reads；單端read中含有的低質(zhì)量（質(zhì)量值Q≤5）堿基數(shù)超過該條read長度比例的50%時，也需去除此對paired reads。最終得到Clean reads進行后續(xù)分析。因為青藏扁蓿豆尚未進行全基因組測序，所以首先選取堿基數(shù)量（176 814 300）最多的 HY_14基于 Stacks2.4［16］無參分析的流程構(gòu)建擬參考基因組，參數(shù)為-m=3，-M=2，-N=2；然后通過短序列比對軟件BWA（http://bio-bwa.sourceforge.net/）將每個樣本過濾后的高質(zhì)量序列與擬參考基因組比對，參數(shù)：MEM-t=4，-k=32，-M；最后根據(jù)比對結(jié)果，利用GATK3.8（https://software.broadinstitute.org/gatk/）開發(fā) SNP位點并使用Vcftools-v0.1.13（http://vcftools.sourceforge.net/）軟件對開發(fā)的SNP位點進行以下條件的過濾：對任意一個位點，群體里至少有80%的樣本有基因型、質(zhì)量值最低為20、次等位基因頻率最小為0.05并且次等位基因數(shù)量大于3，從而獲得高質(zhì)量的SNP位點。使用VCFTOOLS軟件對SNP分型結(jié)果進行轉(zhuǎn)換（transition）與顛換（transversion）統(tǒng)計并由EXCEL繪圖。

1.2.4 遺傳數(shù)據(jù)分析 利用ADMIXTURE（http://software.genetics.ucla.edu/admixture/download.html）軟件檢測群體的遺傳結(jié)構(gòu)，所得結(jié)果通過Excel進行可視化繪圖；采用GCTA（https://cnsgenomics.com/software/gcta/bin/gcta_1.93.2beta.zip）進行主成分分析（principle component analysis，PCA），并使用R 包 ggplot2［17］畫圖 ；利用 ARLEQUIN 3.5［18］軟件計算遺傳多樣性指數(shù)（Pi）、觀測和期望雜合度（HO和HE）以及遺傳分化指數(shù)（FST），進行遺傳多樣性和遺傳分化程度分析；利用GENALEX6.503［19］軟件中的Mantel test對成對的群體統(tǒng)計地理距離，并計算其與遺傳距離的相關(guān)性。

2 結(jié)果

2.1 測序質(zhì)量

對80份青藏扁蓿豆材料的測序數(shù)據(jù)進行統(tǒng)計，包括5個群體的平均測序片段數(shù)，平均堿基數(shù)，群體Q30與GC含量，以及酶捕獲率。結(jié)果（表2）顯示，測序共獲得60.79 Gb數(shù)據(jù)，平均每個樣本0.76 Gb，5個群體產(chǎn)生的測序片段數(shù)以及堿基數(shù)存在差異，其中，最多的為湟源（HY）群體（平均測序片段數(shù)為6 638 025條，平均堿基數(shù)為0.99 Gb），最少為西寧西山（XNXS）群體（平均測序片段數(shù)為4 360 916條，平均堿基數(shù)為0.65 Gb），平均Q30≥93.07%，平均GC含量≥42.29%，數(shù)據(jù)滿足分析要求。經(jīng)過低質(zhì)量數(shù)據(jù)過濾后得到的平均測序片段數(shù)為3 468 857（XNXS）-5 456 407（HY），群體測序平均酶捕獲率≥98.74%，總計保留下12 796個高質(zhì)量SNP位點用于后續(xù)分析。

表2 青藏扁蓿豆群體測序數(shù)據(jù)統(tǒng)計Table 2 Summary of sequencing data of M.archiducis-nicolai populations

2.2 SNP突變類型及數(shù)量

對SNP分型結(jié)果進行突變類型統(tǒng)計，結(jié)果（圖1）表明，C/T轉(zhuǎn)換類型最多（4 124個），占所有堿基突變類型的32.23%，A/G轉(zhuǎn)換類型（4 037個）占31.55%；G/C顛換類型最少（518個），占所有堿基突變類型的4.05%，A/C（1 164個）、A/T（1 777個）和G/T（1 176個）顛換類型分別占9.10%、13.89%和9.19%。轉(zhuǎn)換與顛換（Ts/Tv）之比為1.761。

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

基于SNP分型數(shù)據(jù)，對來自5個群體的80份青藏扁蓿豆材料進行群體遺傳結(jié)構(gòu)分析，提取K=2-5時的交叉驗證錯誤率（cross-validation error，CV error）（圖2-A），可見K從2-5交叉驗證錯誤率值逐漸增大，因此，可將K=2作為最佳K值。當K=2時，群體遺傳結(jié)構(gòu)圖將青海祁連（QLCD）群體與其他群體劃分開（圖2-C）。主成分分析以得到的SNP為基礎(chǔ)進行分析（圖2-B），結(jié)果顯示，個體的聚類情況能夠清楚地被反映出來，親緣關(guān)系的遠近也可以由群體之間的距離反映。第一主成分（PC1）將青海祁連（QLCD）群體與其他群體區(qū)分開，這與ADMIXTURE分析的結(jié)果相一致。

圖2 群體遺傳結(jié)構(gòu)分析Fig.2 Population genetic structure of populations

利用遺傳分化指數(shù)（FST）評估青藏扁蓿豆群體間的遺傳分化程度，發(fā)現(xiàn)群體間的遺傳分化指數(shù)為0.009 3-0.036 9，所有群體間均存在極顯著（P＜0.01）的遺傳分化（表3）?；贛antel test分析，結(jié)果表明，青藏扁蓿豆野生群體的地理距離與遺傳距離之間存在極顯著的正相關(guān)關(guān)系（R2=0.877 6，P=0.006，圖3）。

表3 群體間遺傳分化指數(shù)（FST）Table 3 Index of genetic differentiation（FST）between populations

圖3 群體遺傳距離與地理距離的相關(guān)性Fig.3 Correlation between population genetic distance and geographical distance

2.4 遺傳多樣性分析

基于開發(fā)的12 796個SNP位點進行遺傳多樣性分析，發(fā)現(xiàn)各群體的平均等位基因數(shù)（NA）均為2；觀測雜合度（HO）在0.187 68（RYS）-0.304 36（QLCD）；期望雜合度（HE）在0.201 97（RYS）-0.364 34（QLCD）；遺傳多樣性指數(shù)（Pi）為0.178 32（RYS）-0.241 34（QLCD）。 表 明 日 月 山（RYS）群體具有相對最低的遺傳多樣性水平，青海祁連（QLCD）群體具有相對最高的遺傳多樣性水平（表4）。

表4 青藏扁蓿豆野生群體遺傳參數(shù)統(tǒng)計Table 4 Statistics of genetic parameters in M.archiducis-nicolai wild populations

3 討論

簡化基因組測序技術(shù)適合于大樣本量的研究，可快速鑒定出高密度的變異，特別是SNPs變異，在物種進化、遺傳多樣性分析、全基因組關(guān)聯(lián)分析等領(lǐng)域中應用越來越廣泛。其中，GBS技術(shù)運用了甲基化敏感的限制性內(nèi)切酶，回避了基因組主要的重復區(qū)域，使得文庫構(gòu)建簡單，成本相對低廉，且有助于高通量SNPs分型技術(shù)的大規(guī)模應用。前人對煙草進行GBS測序，利用得到的SNP位點建立的聚類樹將供試材料進行了準確劃分［20］。也有研究表明，利用GBS技術(shù)對紫花苜蓿品系進行測序，并結(jié)合表型性狀進行全基因組關(guān)聯(lián)分析，可得到與黃萎病性狀顯著相關(guān)的SNP標記［21］。本研究利用GBS技術(shù)對80個青藏扁蓿豆樣本進行測序，得到了60.79 Gb的原始數(shù)據(jù)，經(jīng)過挖掘并過濾最終獲得了12 796個高質(zhì)量的SNP位點，用于后續(xù)遺傳多樣性和遺傳結(jié)構(gòu)分析。

堿基置換突變是形成生物個體多態(tài)性和推動種群進化的根本原因之一，大部分堿基置換突變只涉及單個堿基的替換，即轉(zhuǎn)換（transition）和顛換（transversion）［22］。轉(zhuǎn)換是嘧啶互換或嘌呤互換，如C/T和A/G型SNP；顛換則是嘧啶和嘌呤相互替換，如A/C、C/G、A/T和G/T型SNP。本研究中所獲得的SNP位點堿基轉(zhuǎn)換類型占63.78%，遠高于堿基顛換（36.22%）的占比，與嚴佳文等［23］研究的火龍果堿基變異類型（轉(zhuǎn)換占63.29%，顛換占36.71%）相同。由于轉(zhuǎn)換有4種可能性，顛換有8種可能性，因此，理論上發(fā)生轉(zhuǎn)換概率與發(fā)生顛換概率的比值（Ts/Tv）應該等于0.5。但實際上，Ts/Tv值往往大于0.5，這種差異性被稱為“轉(zhuǎn)換偏差”［24-26］。本研究挖掘的SNP突變類型存在明顯的轉(zhuǎn)換型偏差現(xiàn)象，Ts/Tv=1.761＞0.5，產(chǎn)生這種差異的原因可能與不同物種在進化中承受的選擇壓力有關(guān)［27］。

一個物種的進化潛力和抵御不良環(huán)境的能力既取決于種內(nèi)遺傳變異的高低［28］，也有賴于遺傳變異分布格局，即群體的遺傳結(jié)構(gòu)［29-30］。本研究對青藏扁蓿豆野生群體的群體遺傳結(jié)構(gòu)分析表明，5個青藏扁蓿豆群體劃分為2個大的類群，而之前研究并沒有將其分開［7］，因此，SNP標記與EST-SSR標記相比具有更高的分辨率。群體的遺傳分化指數(shù)（FST）是衡量群體間遺傳分化程度的重要參數(shù)，可以解釋影響種群發(fā)生遺傳分化的因素［31］。本研究5個群體中，群體分化指數(shù)為0.009 3-0.036 9，且所有群體間的群體分化指數(shù)均小于0.05，根據(jù)Wright［32］對遺傳分化指數(shù)的界定，F(xiàn)ST＜0.05表明群體遺傳分化較弱，0.05＜FST＜0.15表明群體中等程度遺傳分化，0.15＜FST＜0.25表明群體較大程度遺傳分化，F(xiàn)ST＞0.25表明群體極大程度遺傳分化。因此，本研究中群體間遺傳分化較弱。此前研究顯示，來自青海卓尼、青海樂都以及內(nèi)蒙古的扁蓿豆供試材料的群體分化指數(shù)為0.047 1-0.138 2，青藏扁蓿豆供試材料中的群體分化指數(shù)介于-0.055 2-0.062 4［5］，表明青藏扁蓿豆群體的傳分化程度較低。本研究中的5個青藏扁蓿豆群體間的分化指數(shù)也較低，表明各群體間遺傳分化程度均較弱、親緣關(guān)系相對均較近。對青藏扁蓿豆群體遺傳距離與地理距離的相關(guān)性分析表明，地理距離與群體的遺傳分化存在極顯著（P＜0.01）的相關(guān)性，與此前基于EST-SSR標記研究的結(jié)果相一致［7］，推測可能與冰期后居群由避難所向研究群體擴散過程中發(fā)生的奠基者效應有關(guān)。

高原抬升，亞洲季風、亞洲內(nèi)陸干旱化以及第四紀冰期等諸多因素共同作用造就了青藏高原復雜獨特的地貌環(huán)境和豐富的遺傳多樣性［33］。本研究發(fā)現(xiàn)，位于青藏高原北部的青海祁連（QLCD）群體具有相對較高的遺傳多樣性，而此前研究表明，來源于青藏高原東南邊緣的青藏扁蓿豆群體具有相對較高的遺傳多樣性［6-7］，推測高原東南邊緣可能存在青藏扁蓿豆冰期避難所。造成這種結(jié)果的差異可能是由于采樣群體數(shù)過少，且群體分布相對集中，依據(jù)并不充分。因此，下一步可擴大采樣群體范圍來進行后續(xù)研究分析。

4 結(jié)論

運用GBS測序技術(shù)在80份野生青藏扁蓿豆材料中開發(fā)了12 796個高質(zhì)量的SNPs。青海祁連（QLCD）群體的遺傳多樣性水平相對最高；5個群體被劃分為2個大的類群，此結(jié)果與主成分分析結(jié)果基本一致。另外，該青藏扁蓿豆群體間的遺傳分化程度較弱，地理距離與群體的遺傳分化程度呈現(xiàn)極顯著（P＜0.01）的相關(guān)性。