基于GBS測序的全基因組SNP揭示大涼山玉米地方品種資源的親緣關系與遺傳分化

冷益豐,羅 樊,陳從順,丁 鑫,蔡光澤

(1.西昌學院 農業(yè)科學學院,四川 西昌 615013; 2.攀西特色作物研究與利用四川省重點實驗室,四川 西昌 615013)

玉米(ZeamaysL.)是禾本科玉蜀黍屬一年生草本植物,又名苞谷、苞米棒子、玉蜀黍、珍珠米等,原產于中美洲和南美洲,是世界上重要的糧食作物,廣泛分布于美國、中國、巴西和其他國家。玉米是我國種植面積最大、產量最多的第一大作物,具有重要的糧用、飼用與工業(yè)用價值。西南山地玉米區(qū)作為我國五大玉米產區(qū)之一,玉米單產僅占全國平均水平的80%。2022年6月,農業(yè)農村部印發(fā)的《“十四五”全國種植業(yè)發(fā)展規(guī)劃》提出,加強糧食作物種質資源普查收集,開展種源關鍵核心技術攻關,加快突破性新品種培育推廣。近年來,隨著育種技術的發(fā)展,種質資源研究成為育種工作的重要任務。

基因分型測序(genotyping-by-sequencing, GBS)技術是在第二代測序技術基礎上發(fā)展起來的一種簡化基因組測序技術,該技術具有成本低廉、操作簡單、方便快捷和高通量等特點[1]。GBS技術通過酶切加標簽的方法,測序獲得酶切位點附近的基因序列信息,進而檢測大量高準確性的單核苷酸多態(tài)性(single nucleotide polymorphism, SNP)變異位點信息,是一種識別大規(guī)模單核苷酸多態(tài)性、降低基因組復雜性和基因分型最經濟有效的方法。SNP作為第三代分子標記技術,在植物群體進化、種群基因交流、構建遺傳圖譜和全基因組關聯分析(genome wide association study, GWAS)等方面有著廣泛的應用[2-6]。SuperGBS是一種改進型的GBS測序技術,主要利用甲基化敏感的限制性核酸內切酶將基因組DNA進行酶切,選擇性回收富含基因編碼區(qū)的酶切片段,然后對其進行高通量測序,通過分析獲得SNP信息[7],有或無參考基因組,都可以開發(fā)SNP標記[8]。

大涼山地處西南高原,光、熱資源豐富,其境內地形分布錯綜復雜,有河谷平壩、低山、中山、高山等。由于獨特的地理條件和生態(tài)環(huán)境,多數玉米雜交種難以適應高海拔玉米生態(tài)區(qū),導致當地玉米生產仍沿用白馬牙、白鶴、大黃等農家常規(guī)種,尤其是部分高山區(qū)域,老百姓的自留種成為玉米生產的主導品種。近年來,隨著雜交玉米的推廣,人們盲目追求高產,加上不斷變化的環(huán)境給玉米地方品種保存帶來了前所未有的威脅,導致玉米地方品種資源逐年減少,這嚴重影響優(yōu)異種質資源的挖掘與利用。到目前為止,尚未有人對大涼山玉米地方品種進行系統(tǒng)收集并對其群體結構、混雜程度、遺傳多樣性等開展相關研究。因此,本研究廣泛收集和整理了來自大涼山26個縣(市)的360份玉米地方品種,利用GBS簡化基因組測序技術分析品種間的親緣關系和群體結構,揭示不同品種間的進化關系,完善本區(qū)域的玉米分類系統(tǒng),以了解大涼山玉米地方品種資源的遺傳背景和系統(tǒng)演化,為下一步進行育種開發(fā)與利用提供科學理論依據。

1 材料與方法

1.1 試驗材料

2019—2020年廣泛收集并整理來自涼山州全境17個縣(市),甘孜州稻城、鄉(xiāng)城、得榮、巴塘,西藏察隅、芒康,云南巧家、永德,攀枝花米易的玉米地方品種360份(表1)。2021年種植于涼山州西昌市阿七鎮(zhèn)大田村,單行區(qū),行長5 m,株距0.25 m。每667 m2種植4 000株。

1.2 試驗方法

1.2.1 基因組DNA提取

苗期剪取各樣本植株幼嫩葉片200 g,放入事先標號且預冷的離心管,液氮速凍。待取完所有植株葉片后帶回實驗室,在-80 ℃超低溫冰箱存放。用天根生物科技有限公司的DNAsecure新型植物基因組DNA提取試劑盒(DNAsecure Plant Kit)提取DNA。通過1%瓊脂糖凝膠電泳和Nanodrop 2000分光光度計檢測提取樣品DNA的質量和濃度。

1.2.2 GBS文庫構建

質檢合格的樣品DNA由青島歐易生物科技有限公司利用SuperGBS技術[7]進行測序文庫構建。采用PstI-HF/MspI對DNA進行酶切,用T4連接酶對酶切后的片段兩端加接頭和樣品標簽(barcode),用磁珠回收系統(tǒng)通過調整磁珠溶液與連接產物的體積比回收300～700 bp的片段,使用高保真酶對回收片段進行PCR擴增,PCR產物濃度經Qubit測定(質量濃度需大于5 ng·μL-1),將混好的文庫上機(Illumina Nova,PE150)測序。

1.2.3 數據質控與比對

基于標簽(barcode)與酶切位點信息,使用Stacks軟件[9]對從測序儀下機的高通量測序數據進行拆分,得到每個樣本的raw reads。使用軟件fastp[10]對raw reads進行質控獲得clean reads。使用BWA軟件[11]將clean reads比對到B73參考基因組(B73_RefGen_v4)。

1.2.4 SNP檢測與注釋

基于樣品與參考基因組的比對結果,使用GATK軟件[12]的Haplotype Caller程序生成每個樣品中的g VCF文件,再通過GATK軟件的Genotype GVCFs程序進行群體SNP檢測,使用GATK軟件的Select Variants程序對得到的群體SNP預測結果進行篩選,得到初步SNP和InDel結果。使用VCFtools軟件[13]對獲得的SNP和InDel分型結果進行過濾,過濾條件為：(1)reads支持數不低于4;(2)剔除最小等位基因頻率(MAF)小于0.01的位點;(3)剔除SNP或者InDel分型缺失率高于20%的位點。使用SnpEff軟件[14]對得到的SNP進行注釋,以確定SNP在基因元件的位置、對氨基酸的變化影響等。

1.2.5 系統(tǒng)關系分析

采用鄰接法(neighbor-joining, NJ)[15]對樣本的SNP位點進行建樹,SNP位點中有缺失則用“-”代替。在treebest軟件[16]中用p-distance方法計算距離矩陣,進化樹的可靠性通過bootstrap法[17]進行重復1 000次檢驗。使用plink2[18]軟件對獲得的SNP標記進行主成分(principle component analysis, PCA)分析。利用ADMIXTURE軟件[19]按照K=2到K=10進行群體結構分析,選取10個不同的seed進行10次重復分析。根據交叉驗證誤差(cross validation error, CVE)確定最優(yōu)K值,對每個K值10次重復CVE值繪制盒形圖。

1.2.6 群體分化指數(fixation index,Fst)和核苷酸多樣性(π)的計算與選擇性清除分析

基于高質量的SNP,按照100 kb的窗口、10 kb的步長對染色體進行選擇性清除區(qū)域檢測。使用R語言的PopGenome軟件包,計算Fst、π和多樣性變化倍數(θπratio)等,并利用Fst值繪制曼哈頓圖。在篩選選擇性清除區(qū)域時,分別以Fst和π前5%分位數對應的數據作為閾值,取該區(qū)域的交集為選擇性清除區(qū)域。

1.2.7 候選基因的檢測與注釋功能富集分析

將鑒定的選擇性清除區(qū)域通過Interproscan軟件[20]進行GO注釋分析。在網站(http：//geneontology.org)上進行GO富集分析,并與Swiss-Prot[21]、GO[22]和KEGG[23]等數據庫進行對比得到注釋信息。參考NCBI數據庫(https：//www.ncbi.nlm.nih.gov/)和玉米數據庫(https：//www.maizegdb.org),對選擇區(qū)域的基因進行功能注釋,對候選基因進行GO功能富集分析。

2 結果與分析

2.1 基因組DNA提取效果

部分樣品的DNA電泳結果如圖1所示,條帶清晰明亮。同時,經Nanodrop 2000分光光度計檢測,樣品的DNA質量濃度均大于50 ng·μL-1,滿足GBS文庫構建要求。

M,DL2000 DNA分子量標準;1～24,不同樣本DNA。M, DL2000 DNA marker; 1-24, DNA samples from different samples.圖1 部分樣品的基因組DNA電泳條帶Fig.1 Genomic DNA electrophoresis of partial samples

2.2 GBS測序質量

360個樣本因低質量(low_quality)被過濾掉的reads占比0.945%～5.180%,因含N數量≥5(too_many_N)被過濾掉的reads占比0.003%～0.056%,因質控后R1長度短于142 bp或者R2長度短于150 bp(too_short_reads)被過濾掉的reads占比0.314%～3.139%,質控后得到的clean reads占比92%以上(92.90%～98.36%)。各樣本的單堿基平均錯誤率均在0.000 3左右,最高在0.000 5左右。所有樣本的A、T堿基比例均在30%左右,G、C堿基比例均大于20%,平均GC含量為48.20%(47.18%～49.57%)。

對360份玉米地方品種基因組DNA進行GBS測序,共產生250.99 GB有效數據,獲得1 659 033 712個clean reads,平均reads數量為4 608 427;平均每份樣品的raw base為0.70 GB,clean base為0.67 GB。各樣品的Q20平均值大于等于94.57%,Q30平均值大于等于87.14%。表明測序的質量高,可以用于后續(xù)的信息分析。

2.3 參考基因組比對

使用BWA軟件將clean reads比對到B73參考基因組上,通過比對率檢測樣品與參考基因組的相似程度,一般比對率越高,相似程度越高。根據比對結果統(tǒng)計樣品深度信息并計算測序數據對參考基因組的覆蓋度。全部樣本比對率為93.55%～98.74%,平均測序深度為14.79×,覆蓋度范圍為1.31%～3.02%。使用Qualimap 2[24]統(tǒng)計各樣本的插入片段長度分布,預估的插入片段長度符合正態(tài)分布(圖2),峰值位置接近建庫時DNA片段的平均長度。

圖2 丹紅_1的插入片段長度分布Fig.2 Distribution of insert size of Danhong_1

2.4 突變(SNP/InDel)檢測

經過突變位點檢測并過濾,共獲得124 342個SNP位點,32 063個InDel位點(其中插入15 738個位點、缺失16 325個位點),總共檢測到156 405個變異位點。發(fā)生堿基轉換位點12 955 634個、堿基顛換位點7 955 116個,堿基轉換/堿基顛換比例為1.628 591。在染色體上通過滑窗(窗口和步長均為500 kb)統(tǒng)計,得到各樣本SNP/InDel在染色體上的分布,如圖3所示。

2.5 SNP注釋

SNP注釋分析表明,124 342個SNP變異位點中對蛋白質具有高效(破壞性)影響,可能引起蛋白質截斷、功能喪失等的變異(high)有1 674個,占0.174 9%;低影響,可能不影響蛋白質變異的(low)有91 776個,占9.587 2%;中等影響,非破壞性變異,可能影響蛋白質功效的變異(moderate)有63 458個,占6.629 0%;修飾,通常為非編碼區(qū)變異,只影響非編碼基因(難以預測其影響程度)的變異(modifier)有800 364個,占83.608 8%(表2)。錯義突變有63 654個,占44.469 8%;無義突變有836個,占0.584 0%;同義突變,有78 650個,占54.946 2%(表2)。

表2 變異影響程度與功能級別影響統(tǒng)計

變異分布最多的是內含子區(qū)域的內含子變異,占32.793 6%,其次是基因下游區(qū)域的基因下游變異,占22.508 0%,第三是外顯子區(qū)域與基因上游變異,分別占14.878 9%和14.673 4%(表3、表4)。

表4 變異類型分布

2.6 系統(tǒng)關系分析

從系統(tǒng)進化樹(圖4)可以看出,整個群體被分成2大類群,類群之間分層明顯。各類群可進一步分離出若干亞群,且存在離群個體。

圖中的品種編號與表1一致。下同。The variety codes in the image are consistent with those in Table 1. The same as below.圖4 360份大涼山玉米地方品種的系統(tǒng)進化樹Fig.4 Phylogenetic tree of 360 maize landraces in Daliangshan Mountain area

2.7 主成分分析

對獲得的124 342個SNP標記進行PCA分析,得到影響最大的3個特征向量,3個主成分變異貢獻率分別為16.13%、2.57%和2.10%,對3個主成分作圖,結果見圖5。由圖5可知,360份大涼山玉米地方品種明顯分成兩部分,各形成一個緊密類群,說明2個類群遺傳背景相差較大;其中一個類群比較分散,而另一個類群相對集中,沒有明顯的遺傳差異,呈混雜在一起的趨勢。主成分分析結果與系統(tǒng)進化樹分析結果一致。

圖5 主成分分析Fig.5 PCA analysis of three-dimensional graph

2.8 群體遺傳結構分析

最佳K值分析結果表明(圖6),隨著K值的增大,CVE值順勢減小,當K=9時,CVE值最小,隨著K值再增大,CVE值也隨之增加,表明K=9時是最合適的分組,K=9時的分組情況見圖7。群體遺傳結構指遺傳變異在物種或群體中的一種非隨機分布。按照地理分布或親緣關系可將一個群體分為若干亞群,處于同一亞群內的不同個體親緣關系較近,而亞群與亞群之間則親緣關系稍遠。從圖7可以看出,360份大涼山玉米地方品種被分成9個亞群。

圖6 不同K值下群體結構的CVE值曲線Fig.6 The CVE curve of the population structure with different K values

圖7 360份大涼山玉米地方品種的群體遺傳結構Fig.7 Genetic structure of 360 maize landraces in Daliangshan Mountain area

2.9 群體遺傳分化分析

估算群體間的Fst可有效區(qū)分群體間相對遺傳變異大小,是解釋群體遺傳變異程度的主要來源。從圖8可見,類群A與類群B的Fst值大于0.462 2,說明類群間有很大的遺傳分化。類群A的核苷酸多樣性πA高于類群B的核苷酸多樣性πB?？赡苁且驗轭惾篈遺傳背景較為復雜,群體內玉米材料血緣豐富,導致類群A遺傳多樣性高于類群B,深層次原因需進一步研究。

2.10 差異化的候選區(qū)域與候選基因

以玉米B73基因組為參考,分析360份玉米差異基因組區(qū)域中的基因。從圖9可見,在A群體和B群體中共檢測到96個受選擇區(qū)域,其中418個基因具有強烈選擇信號。進一步對選擇性掃描的區(qū)域進行基因功能注釋,發(fā)現富集的基因組區(qū)域中包含寒冷響應基因Zm00001d027468(OJ2056_H01.33)、Zm00001d006508(YUP8H12.19)、Zm00001d020961(T24D18.5)、Zm00001d048168(F18A8.1),缺水響應基因Zm00001d033107(OsJ_16068)、Zm00001d013399(F2H15.17)、Zm00001d022002(OsJ_11395)、Zm00001d048165(F13H10.2),病菌響應基因Zm00001d028273(T9L24.32)、Zm00001d005727(T8O11.18)、Zm00001d006507(F8K4.6)、Zm00001d047690(MZN1.6)、Zm00001d025318(F1N19.11)等與逆境應答相關的基因。

lg(π ratio)和Fst分別對應上面的頻率分布圖和右側的頻率分布圖,中部的點圖則代表不同窗口內的相應的Fst和π比值。其中最上方綠色和紅色區(qū)域為π選擇出來的top 5%區(qū)域,右側橘色區(qū)域為Fst所選擇top 5%區(qū)域,中間紅色和紫色區(qū)域為Fst和π的交集,即為候選的位點。lg(π ratio) and Fst values, corresponding to the frequency distribution chart above and on the right, the point plot in the middle represents the corresponding ratio of Fst to π in different windows. The green and red areas at the top are the top 5% areas selected by π, and the orange area on the right is the top 5% area selected by Fst, and the red and purple areas in the middle are the intersection of Fst and π, which are the candidate loci.圖9 類群A和B的差異基因組區(qū)域Fig.9 Different genomic regions between group A and B

3 結論與討論

本研究采用GBS簡化基因組測序技術對收集自大涼山26個縣(市)的360份玉米地方品種進行了基因分型,獲得了大量的SNP分子標記,經過SNP calling并過濾,共獲得124 342個SNP位點、32 063個InDel位點(其中插入15 738個位點、缺失16 325個位點),總共檢測到156 405個變異位點?；讷@得的124 342個SNP位點,系統(tǒng)進化樹聚類分析和PCA主成分分析將群體分成A和B兩大類群,類群間遺傳分化明顯(Fst大于0.462 2),類群A的核苷酸多樣性高于類群B。玉米地方品種在經歷了長時間的人工選擇和對大涼山地區(qū)的氣候環(huán)境產生高度適應性后,類群間已經產生了明顯的遺傳差異,形成了獨立的遺傳背景,且在進化過程中受到的選擇強度較大,推測可能是由于長期的自然選擇和人工選育造成的。這些數據將為大涼山玉米地方品種的表型差異和復雜性狀的遺傳分析提供重要的遺傳基礎依據。研究結果將為本區(qū)域玉米重要目標性狀基因挖掘、農藝性狀QTL定位以及玉米新品系、新品種培育提供重要支持。

玉米種質類群的劃分問題一直存有爭議。我國玉米種質資源類群劃分研究始于20世紀50年代中后期,主要為瑞德群、蘭卡斯特群、塘四平頭群和旅大紅骨群4個類群。殷洪達等[25]根據1年2點產量和相關性狀的配合力與雜種優(yōu)勢表現,將40份中晚熟、耐密、抗病自交系分為Reid、Lancaster、塘四平頭和旅大紅骨4個雜種優(yōu)勢群。呂學高等[26]利用40個SSR分子標記將浙江省90份玉米地方種質資源劃分為塘四平頭、Reid、旅大紅骨和Lancaster等類群。吳承來等[27]利用112對SSR標記將我國97個玉米自交系劃分為PB、Reid、塘四平頭和旅大紅骨4個類群。類群多而復雜并不利于提高育種效率。張世煌等[28]認為雜種優(yōu)勢群的劃分應盡量向兩群靠攏,并最終歸納為兩個雜種優(yōu)勢群,故他提出了兩群劃分策略。本研究通過GBS測序技術進行系統(tǒng)發(fā)育與群體結構分析,把360份大涼山玉米地方品種分為兩大支系9個亞群,這與Mumm等[29]利用46個RFLP標記將148個美國常用自交系劃分為Lancaster Sure Crop(LSC)和Iowa Stiff Stalk Synthetic(BSSS)兩大類群11個亞群,孫友位等[30]利用70個SSR分子標記將85份早熟玉米自交系劃分A(Reid、熱帶種質、旅大紅骨)和B(Lancaster、塘四平頭、農家種)兩大類群6個亞群,張偉瑋等[31]利用43個SSR標記將56份玉米自交系劃分為SS(旅大紅骨、PA、Reid)和NSS(PB、四平頭、Lancaster)兩大種質5個亞群,王文斌等[32]利用2 846個高質量SNP標記將陜西31份玉米自交系劃分為A(丹340、掖478等)和B(Mo17、黃早四等)兩個類群,吳迅[33]利用56 110個SNP標記將367玉米自交系分為外引種質和本地種質兩大類群5個亞群(瑞德、蘭卡斯特、塘四平頭、溫熱I和P群)是基本一致的。本研究的系統(tǒng)進化樹結果支持玉米雜優(yōu)類群的兩群學說,但類群間是否具有育種上的雜交組配優(yōu)勢,還需要對具體材料的配合力進行分析。本研究的群體結構分析進一步將兩大類群分成9個亞群,與上述前人研究結果不同,這可能與研究材料本身的血緣有很大關系,導致分析所用的數據集不同而呈現差異。

大涼山地處西南高原,是玉米傳播到中國后較早的種植區(qū)之一。境內大山密布,交通極不方便。在海拔1 800～2 800 m的山地和山谷旱地之間,因雜交種難以適應其獨特的地方性氣候和生態(tài)條件,至今依然保留了較多的地方品種資源。這些品種由于長期所處氣候、生態(tài)條件的差異,形成了各自鮮明的特點,如抗旱、抗寒、耐貧瘠、抗病等,籽粒顏色豐富多彩(紅色、白色、紫色、黃色等),其中蘊含了復雜多樣的基因資源,使其成為我國玉米地方種質最為豐富的地區(qū)。此外,由于相對落后的文化、交通和經濟條件,使得該區(qū)域內種質資源交流較少,同時外來種質也難以適應局部小氣候,進一步導致種質交流和遺傳滲透遠遠落后于其他玉米生態(tài)區(qū),從而使得該區(qū)域玉米地方種質資源能夠得到很好的保存和利用,遺傳多樣性也遠遠高于其他玉米生態(tài)區(qū)。大涼山的玉米發(fā)展離不開地方特有玉米種質資源的挖掘,尤其是耐寒、耐旱和耐瘠群體、白粒適口玉米等自交系材料或群體的改良和擴增工作。本研究中,通過選擇信號分析在大涼山地方玉米品種群體中發(fā)現96個基因組區(qū)域受到選擇,檢測到受選擇區(qū)域的418個候選基因,部分基因已被證明與玉米的抗寒、抗旱和抗病等有關。楊宇昕等[34]利用群體分化指數和群體間擴展單倍型純合度(cross population extended haplotype homozygosity, XP-EHH)法分析溫帶和熱帶玉米群體間的選擇信號分布情況,選擇Fst和XP-EHH值的top 1%為閾值,分別鑒定到557和1 913個候選基因,多個候選基因與玉米的開花調控密切相關,包括ZmCCT9、COLl、GRMZM2G387528。周玲等[35]基于187份種質資源材料的全基因組重測序數據開發(fā)了120 583個高質量SNP標記,通過整合選擇信號檢測和GWAS分析結果,共識別出與1 153個候選位點緊密連鎖的324個候選基因,涉及氮化合物代謝、葉酸代謝、糖酵解、發(fā)育過程的負調控與轉錄調控等重要生物學途徑。

西南山地玉米區(qū)是我國三大玉米主產區(qū)之一,全區(qū)玉米種植面積500萬hm2左右,約占全國玉米生產面積的1/5[36]。該區(qū)地形地貌十分復雜,山地丘陵和高原土地總面積占90%以上,從海拔250 m河谷平壩直至2 500 m的高山均有玉米種植。生態(tài)條件區(qū)域差異很大,垂直分布明顯,農業(yè)立體性強。長期以來,西南地區(qū)玉米育種使用的自交系絕大多數為二環(huán)系,主要選育來自美國玉米帶雜交種和我國北方玉米優(yōu)良雜交種,種質基礎極為狹窄,導致所選育的雜交種產量沒有突破、穩(wěn)定性差,抗病性、適應性不強。這主要是由于急功近利的商業(yè)育種熱衷于抄近道,追求多(育種材料多)、短(組配時間短)、快(上市品種快),盲目跟風育種,卻未能靜下心來認真研究本土種質。西南地區(qū)類似大涼山這樣的山地玉米生態(tài)區(qū)還有很多,區(qū)內地形地貌復雜,生態(tài)條件多種多樣,種植玉米歷史悠久,加之長期以來部分區(qū)域交通不便造成相對閉塞,形成了許多適應性強、獨具特色的玉米地方種質資源。這部分地方玉米種質資源具有很大的開發(fā)利用潛力,目前存在的主要問題是鑒定、篩選和研究利用嚴重滯后。因此,對于西南玉米雜交育種,面對多種多樣的雜交種類型需求,必須把種質資源的拓展與創(chuàng)新作為今后重要的攻關方向;建議加強對本土地方玉米優(yōu)良種質的評價與遺傳多樣性研究,從特異抗逆基因的發(fā)掘與利用著手,簡化雜優(yōu)模式,分類合成與改良育種群體。