基于GBS技術的233份大蒜種質資源群體進化分析

管夢嬌，謝亞楠，顧啟玉，孫秀東，高榮廣

（山東農業大學園藝科學與工程學院，山東 泰安 271018）

大蒜（Allium sativum L.）為百合科（Liliacaea）蔥屬（Allium）兩年生草本植物，其幼苗（蒜苗）、花莖（蒜薹）和鱗莖（蒜頭）甚至青蒜和蒜黃都能作為蔬菜來食用，又因為大蒜擁有獨特的香氣，也是人們日常生活所必備的重要香辛料之一［1］。大蒜富含大蒜素、蛋白質、糖類、抗壞血酸和多種還原性硫化物等物質，具有非常高的烹飪價值、功能性價值和醫學價值［2，3］。大蒜原產自西亞地區的高原地帶，在中國已有2 000多年的種植歷史，我國已成為世界上大蒜種植面積和產量最多的國家。

大蒜是我國傳統的出口創匯蔬菜產品，屬于我國的名特蔬菜種類。由于我國大蒜種質資源保護工作起步晚、收集保存資源不足、遺傳多樣性及親緣關系缺乏系統研究，再加上環境和氣候變化使大蒜不斷發生變異及各地間頻繁引種，導致大蒜種質名稱混亂，種系不清，同名異物或異名同物的情況非常嚴重。因此，亟需對我國大蒜種質的遺傳進化特性進行鑒定，構建種質資源庫，為探索種質創新方法、研究配套關鍵技術以促進大蒜的豐產豐收奠定基礎。

簡化基因組測序（GBS）技術是指通過測序進行基因分型，是通過選取合適的限制性內切酶結合高通量群體測序來構建SNP分子標記的方法，可應用于超高密度遺傳圖譜構建［4，5］、群體GWAS分析［6，7］、群體遺傳分析［8］等領域。GBS避免了傳統方法的一些不足，雖不能達到覆蓋全基因組深度的程度，但由于其省時、價廉、有效等優點而得到大家的青睞［9］。GBS的發展有助于發現植物性狀及其相關變異來源，對今后品種選育、分子育種都有一定的指導意義。

群體進化分析是對同一物種的不同亞群或不同地理分布的品種進行全基因組測序或基于酶切的簡化基因組測序，通過與參考序列比對，獲得大量高準確性的SNP、InDel等變異信息，研究生物群體中的基因頻率和基因型頻率，探討群體的遺傳結構、遺傳平衡及影響群體遺傳平衡的因素，從分子層面揭示該物種的進化機制、環境適應性等系列問題［10］。現已有黃瓜基因組的馴化和多樣性研究［11］、黑脈金斑蝶遷徙和色變［12］、番茄的選育史［13］、大非洲稻的馴化［14］、豆的馴化［15］等模式生物的群體遺傳進化研究報道。

目前，已有多人使用不同的分子標記技術對大蒜進行了多樣性分析，如陳昕等利用RAPD和ISSR兩種分子標記技術對中國10個不同地區的大蒜進行種質資源多樣性研究［16］；韓曙等利用RAPD標記對18個大蒜品種進行遺傳多樣性和親緣關系分析［17］；陳書霞等利用SSR技術對40個大蒜品種進行聚類分析、主成分分析和遺傳多樣性評價［18］；王海平等用AFLP、SSR和InDel三種分子標記對212份大蒜資源進行分析［19］；孫亞麗等利用SSR分子標記技術對55份大蒜品種進行多樣性分析［20］。

本研究在前人研究的基礎上，應用GBS技術獲得了高質量的SNP位點信息，進而構建系統進化樹，并進行主成分分析，以期為大蒜種質資源的收集和保護提供依據，發掘優良種質、優良基因，為日后解析鱗莖等性狀的遺傳基礎、基因克隆奠定基礎，為大蒜性狀分子改良提供資源。

1 材料與方法

1.1 試驗材料

試驗于2018年9月至2019年6月在山東農業大學園藝作物生物學實驗室進行。所用233份大蒜種質材料及來源見表1，來自于不同時期的資源收集活動。其中，170份種質材料收集于國內，59份來源于國外，剩余4份來源不詳。

表1 試驗用大蒜種質材料及來源

表1（續）

表1（續）

1.2 群體進化分析

1.2.1 DNA提取和GBS文庫構建 DNA提取：2018年12月采集所有大蒜資源的嫩葉，利用植物基因組DNA提取試劑盒（Magen Hi Pure Plant DNA Mini Kit，Guangzhou，China）提取DNA。將提取的樣品DNA送諾禾致源進行DNA質檢、建庫和測序，使用Qubit?2.0熒光測定計（Invitrogen，Carlsbad，USA）檢測核酸濃度，使用Nanodrop檢測DNA的純度（OD260／280），使用Qubit對DNA濃度進行精確定量。選用高質量的DNA用于GBS文庫構建和測序。

GBS文庫構建：為了得到適合的marker密度，首先用限制性內切酶對0.1～1.0μg基因組DNA進行酶切；將P1和P2 Adapter接頭加在酶切得到的片段兩端后，再對樣品進行擴增，分別對含有P1和P2接頭的tag序列進行PCR兩端擴增，DNA片段pooling，然后混合樣品，電泳回收需要區間的DNA，選擇所需的片段進行GBS文庫構建。利用IlluminaHiSeq測序平臺，進行雙末端（Paired-End）150測序。

1.2.2 SNP檢測與注釋 用堿基識別（base calling）測序得到的原始圖像數據文件，再經過分析轉化為原始測序序列（sequenced reads）。分析前需將sequenced reads中可能包含的接頭信息、低質量堿基以及未測出的堿基（以N表示）等干擾信息去掉，最終得到的數據即為有效數據（clean data或clean reads）。原始數據過濾方法如下：首先過濾掉含有接頭序列的序列；然后檢查是否有單端測序序列中未測出堿基的含量超過該條序列長度10%的現象，如果有，也需要去除此paired reads；最后，若單端測序序列中含有的低質量（≤5）堿基數超過該條序列長度50%時，也需要去除此paired reads。經過嚴格過濾的測序數據就是高質量的有效數據。將高質量、有效的測序數據通過BWA軟件（參數：mem-t4-k32-M）比對到參考基因組（http：／／gigadb.org／dataset／view／id／100209／File_page／1）。

進行群體SNP檢測，使用的是SAMTOOLS等軟件。利用貝葉斯模型來檢測待測群體中的SNP位點。再利用ANNOVAR軟件對得到的高質量SNP進行群體SNP注釋，并對由多個基因組檢測出的基因變異進行功能注釋。

1.2.3 群體進化分析 SNP檢測之后，得到的個體SNP可以用于計算種群之間的距離（p-距離）。兩個體i和j之間的p-距離通過如下公式計算：

公式中dij表示兩個個體為不同基因型的賦值表，L為高質量SNPs區域長度。若在位置1的等位基因為A／C，那么：

運 用 TreeBest（http：／／treesoft.sourceforge.net／treebest.shtml）軟件計算距離矩陣，以此為基礎，通過鄰接法（neighbor-joining method）構建系統進化樹。引導值（bootstrap values）經過1 000次計算獲得。

1.2.4 主成分分析 主成分分析（PCA）是一種純數學的運算方法，可以將多個相關變量經過線性轉換選出較少個數的重要變量。PCA方法被廣泛運用于多學科，在遺傳學當中，主要用于聚類分析，它是基于個體基因組SNP差異程度，按照不同性狀特征將個體按主成分聚類成不同的亞群，同時用于和其它方法做相互驗證。本試驗通過GCTA（http：／／cnsgenomics.com／software／gcta／pca.html）軟件計算特征向量以及特征值，并利用R軟件繪制PCA分布圖。

2 結果與分析

2.1 測序質量

2.1.1 酶切效率 選用MseⅠ+Eco RⅠ限制性內切酶組合進行酶切，統計clean reads兩端為MseⅠ捕獲的reads數、酶捕獲率（捕獲的reads數占clean reads數的比率），用來評估酶切效率。酶切效率平均在99.0%左右，說明MseⅠ+Eco RⅠ限制性內切酶的組合酶切效率很高。

2.1.2 測序數據統計與質量評估 233個大蒜樣本總的測序數據量為1 486.6871 Gb，去除低質量序列后剩下的高質量有效數量為1 486.5896 Gb，平均每個樣本6.3802 Gb。測序質量高（Q20≥93.27%、Q30≥83.50%），GC分布正常，說明233個大蒜樣本都沒有被污染，此次建庫測序成功。

2.1.3 測序深度及覆蓋度 平均測序深度為7.68X～17.11X，1X覆蓋度（至少有一個堿基的覆蓋）在4.21%以上。各個樣本的比對結果顯示，它們與參考基因組的相似度達到重測序分析的要求，同時又有非常不錯的覆蓋深度和覆蓋度。

2.1.4 SNP的檢測 經SAMTOOLS軟件檢測，共獲得了67 817 466個SNPs位點。用個體深度DP＝4、缺失率Miss＝0.2、最小等位基因頻率MAF＝0.01等為條件進行過濾，最后獲得了高質量的SNPs位點2 036 116個（表2），可用于后續分析。

表2 SNPs統計信息

2.2 群體進化樹分析

系統進化樹是用來表示群體間進化關系、描述群體間進化順序的分支樹，可以根據群體物理或遺傳學特征等方面的共同點或差異推斷出其親緣關系遠近，即分析群體中個體間由于共同祖先而產生的相互關系。基于2 036 116個SNPs標記，構建系統進化樹（圖1），可以看出，233份大蒜自然群體被劃分為8個亞群，第一亞群為藍色部分，共包含10份材料，分別來源于中國的山東、四川、貴州、廣東以及埃及、美國、泰國和加拿大等；第二亞群為粉色部分，兩份分別來自于中國云南和中國臺灣，另一份來自歐洲；第三亞群為黃色部分，7份來源于加拿大，2份來源于烏克蘭；第四亞群為橙色部分，12份材料除1份來源于波斯外，其余全部來源于加拿大；第五亞群為棕色部分，18份材料中有6份來源于加拿大，其余12份則來源于中國的新疆、甘肅、青海、山西等地；第六亞群為紅色部分，56份材料中，有2份來源于日本，其余54份材料則全部來源于中國，主要來源于山東、河北、新疆、四川等地；第七亞群為綠色部分，12份材料中，3份來源于中國山東和陜西，其余9份則來源于美國、加拿大和澳洲；第八亞群為淡紫色部分，113份材料中，有4份來源不詳，其余則大部分來源于中國山東、新疆、河南、湖北等地。

2.3 主成分分析

本試驗利用經過濾得到的2 036 116個SNPs進行PCA分析，在忽略錯配數據和高于2個等位基因位點的條件下，選取貢獻率最大的前兩維數據繪制散點圖，見圖2，其中group1—group 8與進化樹依據顏色相對應。可知，PCA將233份大蒜群體分為4大類，將系統進化樹劃分的8個亞群用不同顏色呈現在PCA分析圖中，結果與進化樹結果大體一致。類群①對應群體進化樹中的藍色部分，大部分材料來自于中國；類群②對應群體進化樹中粉、黃、橙、棕、紅色部分，其中約有五分之三的材料來自于中國，另有約五分之一的材料來自于加拿大；類群③對應群體進化樹中的綠色部分，約有半數的材料主要來自于美國，其余材料則分別來自于美國周邊地區以及中國；類群④對應群體進化樹中的紫色部分，幾乎全部材料都來自于中國。

圖1 233份大蒜種質群體進化樹

圖2 233份大蒜種質的主成分分析結果

3 討論

本試驗得到高質量的有效數據量為1 486.5896 Gb，平均每個樣本6.3802 Gb，與參考基因組（大小為16 557 142 287 bp）相比，群體樣本平均比對率為82.94% ～99.23%，說明測序數據與參考基因組的相似程度較高。利用GBS簡化基因組測序技術，選用MseⅠ+Eco RⅠ酶切組合，測序深度為7.68X～17.11X，用個體深度DP＝4、缺失率Miss＝0.2、最小等位基因頻率MAF＝0.01等為條件進行過濾后，共獲得高質量SNPs位點2 036 116個。

233份大蒜資源的系統進化樹和PCA分析是基于SNP構建的，可以從中了解到大蒜不同種質之間的遺傳距離，進而了解他們之間的親緣關系遠近。在進化樹中，來源于同一地區的大蒜大都被聚集在一起，如第六亞群全部來自于國外，且幾乎都來自于加拿大，第七亞群大部分來自于北美洲，第五亞群主要來自于加拿大和中國西北部，第六亞群幾乎全部來自于中國，這與陳書霞［18］、王海平［19］和孫亞麗［20］等的試驗結果相同。不同地域來源的種質也可能被聚在一起，如第一亞群的9個種質來自三個大洲、五個國家、八個地區；第二亞群的3個種質分別來自于中國臺灣、中國云南和南歐，雖然三者地理位置相距較遠，但均擁有溫暖濕潤的環境條件，說明相似的地理環境可能使大蒜的遺傳多樣性趨于相同；第六亞群中僅有1／28的種質來自日本，其余均來自于中國，可能由于兩地在地理位置上距離較近、氣候相似；第三亞群大蒜種質來自加拿大和烏克蘭，這與孫亞麗等［20］的研究結果相同。

4 結論

本研究利用GBS技術對233份大蒜種質資源進行群體進化分析，共得到高質量的有效數據1 486.5896 Gb，平均每個樣本6.3802 Gb；共獲得高質量的SNPs位點2 036 116個，較均勻地分布在8條染色體上。基于獲得的SNPs標記進行PCA與系統進化樹分析，PCA將大蒜群體劃分為4大類，系統進化樹將大蒜群體劃分為8個亞群，將8個亞群用不同顏色在PCA圖中呈現，兩者結果基本一致，表明大蒜種質間的親緣關系主要受地理起源影響，另外還與各地大蒜種質資源交流頻繁、氣候環境等有一定關系。本研究結果可從分子水平上為大蒜親緣關系分析和系統分類等提供輔助作用，并為大蒜栽培和育種提供科學依據。