張 羽，F(xiàn)ran?ois Belzile

(1.陜西理工大學 生物科學與工程學院，陜西 漢中 723000；2.Département de Phytologie，Université Laval，Québec，QC，Canada G1V 0A6)

大豆菌核病是20世紀80年代中期開始出現(xiàn)的世界范圍大豆病害，主要分布在美國、加拿大、中國、阿根廷、匈牙利、日本和印度等國，特別是在北美和加拿大的寒涼大豆種植區(qū)已經(jīng)成為影響大豆產(chǎn)量的主要病蟲害之一[1-2]。它主要由核盤菌(Sclerotiniasclerotiorum(Lib.)de Bary)侵染寄主植株所引起。1986年，加拿大的豆類作物菌核病發(fā)病率達25%[3]。近年來，大豆菌核病有繼續(xù)擴大蔓延、加重危害的趨勢，嚴重威脅大豆的生產(chǎn)[4]。生產(chǎn)實踐證明，培育抗(耐)病大豆品種是防治大豆菌核病最經(jīng)濟有效的途徑。目前為止，還沒有發(fā)現(xiàn)完全免疫菌核病的大豆資源，前人研究報道，大豆種質資源對菌核病的抗性存在著遺傳差異，有的大豆資源具有部分抗性，大豆菌核病是由多基因控制的數(shù)量性狀遺傳。數(shù)量性狀對環(huán)境條件敏感，在大田條件下目標性狀的選擇壓力難以保證一致，導致篩選抗菌核病種質資源比較困難，因此，利用分子標記輔助選擇培育抗性品種顯得越來越重要。

大豆抗菌核病的QTLs定位自2000年Kim等[5]利用分子標記RFLP、SSR和RAPD進行研究以來，很多研究都是基于SSR分子標記對大豆抗菌核病進行QTLs定位，但由于一代和二代分子標記在基因組中的多態(tài)性有限，使其在QTLs定位中的利用不太理想，隨著二代測序技術的發(fā)展，基于SNPs的全基因組關聯(lián)研究在生物領域成為熱點，自從全基因組關聯(lián)研究在人類疾病研究中取得進展后，植物上通過全基因組關聯(lián)研究進行病害候選基因定位也越來越受到重視。由于全基因組測序面臨費用高和較多重復序列的緣故，為了提高測序性價比，基于限制性核酸內切酶識別位點相關DNA序列(RAD tags)的簡化基因組測序技術得到了普遍應用，其中GBS(Genotyping by sequencing)技術由于能大幅降低基因組的復雜度，尤其適合于大樣本量的研究，可快速鑒定出高密度的SNPs位點，其結果在QTLs定位中應用越來越廣[6-12]。目前為止，在大豆上多選擇Apek Ⅰ、MspⅠ和PstⅠ酶進行簡化基因組后測序，基于SNPs的全基因組關聯(lián)研究進行大豆菌核病QTLs定位的研究報道較少，由于定位樣本量大小、表型鑒定方法、分析算法及標記多態(tài)性位點數(shù)目等因素，影響了QTLs定位的一致性。且由于定位材料及標記數(shù)目不同，研究結果也不一致[13-16]。分子標記輔助選擇培育大豆抗菌核病育種還未見報道，抗性種質資源的選擇很多還是用草酸浸泡葉、根、莖后，根據(jù)植株萎蔫程度和莖中可溶性色素水平高低評價抗性強弱[17]，因此，大豆菌核病QTLs的定位對大豆抗性育種及抗性機理的研究進展相對緩慢。本研究通過對126個加拿大大豆品種的30 125個SNPs，利用PLINKv 1.07分析軟件，進行SNP-phenotype和Haplotype-phenotype的關聯(lián)分析，找到大豆與抗菌核病強關聯(lián)的位點/候選基因，為大豆菌核病QTLs分析及大豆抗菌核病育種提供參考。

1 材料和方法

1.1 DNA提取、文庫準備、測序、數(shù)據(jù)質量控制

126個大豆品種的DNA提取、文庫準備、測序、質量控制等步驟參考文獻[13]。

1.2 分析方法

1.2.1 表型值 用棉墊接種法接種核盤菌菌絲進行表型鑒定[18]，莖稈損傷長度介于28.6～192.4 mm，供試材料的表型值分布見圖1。

圖1 供試材料主莖損傷長度表型值分布

1.2.2 數(shù)據(jù)質量控制與分析 數(shù)據(jù)質量控制是全基因組關聯(lián)分析的關鍵之一，不同參數(shù)的選擇代表了不同的遺傳學意義，產(chǎn)生不同的分析結果，不當?shù)膮?shù)選擇會產(chǎn)生假陽性或假陰性結果。本研究位點間的關聯(lián)程度r2值0.8，LOD取值3，置信區(qū)間為95%，MAF值取0.01。由于大豆為自交作物，偏離哈迪-溫伯格平衡的SNPs很可能是與病害關聯(lián)的易感位點，因此，本研究不考慮哈迪-溫伯格平衡，得到30 125個SNPs(圖2)。126個材料共有3.7958E6個位點，其中雜合基因型位點有248 456個，占6.546%，雜合體的質量控制在10%以下，說明了數(shù)據(jù)可信度高[19-21]。

圖2 每條染色體上的SNP數(shù)

2 結果與分析

2.1 126份大豆材料的NJ法聚類分析

用TASSEL 5.0的鄰接法(The neighbor-joining algorithm，NJ)對126個供試材料進行聚類(圖3)，聚為3大類，分別為67，45，14個材料。

圖3 基于NJ法的聚類圖

2.2 PC分析

用SPSS對126份大豆品種進行主成分分析，結果顯示前2個主成分的方差貢獻率分別為6.24%和5.44%，Ward法聚類結果見圖4。126份供試材料聚為2類，分別為70，56個材料。

2.3 群體結構分析

很多研究表明，群體結構分析是對目標性狀進行關聯(lián)分析的前提，因為復雜的群體結構會導致基因型與表型之間的假陽性。研究發(fā)現(xiàn)，考慮群體結構和不考慮群體結構時檢測到的與表型變異相關的多態(tài)性位點數(shù)不同[22]。本研究采用Structure 2.3.4軟件模擬群體遺傳結構[23]，其基本原理是：假設樣本存在K個(K值未知)等位變異頻率特征類型數(shù)，將供試材料歸到K個群體，使得該群體位點頻率都遵循同一個Hardy-Weinberg 平衡。設定群體數(shù)K的取值為1～10，將MCMC(Markov chain Monte Carlo)開始時的不作數(shù)迭代(Length of burn-in period)設置為10 000 次，再將不作數(shù)迭代后的MCMC設置為10 000 次，每個K運行20次，然后依據(jù)ΔK值最大的原則選取一個合適的K 值。本研究的ΔK值最大時的K值為2，見表1。將126個供試材料劃分為2個亞群，分別為68，58個材料。

圖4 Ward法聚類圖

表1 群體結構分析

表2 用Structure、PCA 和UPGME方法的聚類結果比較

2.4 連鎖不平衡分析

用PLINKv 1.07命令“--r2”運行LD，過濾掉r2小于0.2的數(shù)據(jù)，共有82 911對，r2平均值為r0.73，r2大于等于0.8的有46 655對(占56.27%)，2等于1的有9 391對(占11.33%)。2個SNPs位點間相關(r2大于0.2)的平均距離為114 979.2 bp，r2小于0.8的平均距離為124 596.9 bp，r2大于等于0.8的平均距離為107 505.3 bp，r2等于1的平均距離為80 397.08 bp。兩位點間距離≤500 kb(遺傳學上把位點間大于500 kb認為沒有連鎖關系)的r2平均值為0.73。由此看出，LD隨著2個位點間的距離增大而衰減。從單條染色體和全基因組r2和2個SNPs位點間距離的關系圖也可以看出隨著物理距離的增加，2個位點間的r2越小，即連鎖強度越小，甚至沒有(圖5，6)。

圖5 一號染色體上的LD衰減圖

圖6 全基因組上的LD衰減圖

2.5 單體型分析

SNPs位點并不是獨立遺傳的，而是在染色體上傾向于以一個整體遺傳給后代，成組遺傳的SNPs位點在一代又一代的遺傳中絕少發(fā)生重組，于是，這樣的一組SNPs位點類型也就是單體型，也就是說相鄰位點是以單體型為單位傳遞的(前提是位點間的LD很強)，當然，單體型分析在統(tǒng)計學上的作用主要是減少多個位點組合分析時的自由度，從而增加研究的效能，例如，如果要分析3個二等位多態(tài)位點的組合和疾病/病害的關系，理論上有9種不同的組合，但是因為這3個位點間不是獨立的，在群體中可能只存在3種常見的單體型，如果3個位點間獨立，也就是連鎖平衡，那么群體中的單體型就會有9種，這時候做單體型分析是沒有意義的，因此，進行單體型分析時最好以單體域Block為基礎。在PLINKv 1.07里有2種單體域推斷方法，Sliding windows和Haplotype list，Haplotype list要先計算Blocks，然后利用算出的Haplotype和Phenotype之間的關聯(lián)程度進行分析，這種運算比Sliding windows靈活，因為Sliding windows只能按SNP順序依次排列分析Haplotype，不能按實際跳躍分析。在比較了2種方法后，本研究采用“--blocks”命令推斷Blocks。20條染色體共有2 863個Haplotype，它們在20條染色體上的分布見圖7，其中18號染色體上的Haplotype最多(230個)，11號染色體上的最少(57個)。一個Haplotype最多的由36個SNPs組成，最少2個SNPs組成。最長Haplotype跨越200 kb(因此，后續(xù)候選基因在200 kb范圍內注釋)，最短跨越2 bp。

2.6 基因型與表型關聯(lián)分析

在基因型-表型的關聯(lián)中，為了發(fā)現(xiàn)較多的關聯(lián)信息，把P小于1E-5的分析結果都列表顯示。

2.6.1 基于Haplotype-phenotype關聯(lián)分析 基于Haplotype-phenotype關聯(lián)分析中，17號染色體上的4個SNPs組成的Haplotype解釋表型變異最大(17.56%)。其余依次為1，3，4，20號染色體上的Haplotype(表3)。

圖7 每條染色體上的Haplotype數(shù)

2.6.2 基于單個SNP-phenotype關聯(lián)分析 在單個SNP的基因型-表型關聯(lián)分析中，α≤0.05(P≤0.05/30125=1.660E-6)時，最強關聯(lián)在3號染色體的34387780位置和其臨近的幾個SNP位點(P值為8.669E-7)，可解釋表型變異的17.80%。在單個SNP-phenotype的關聯(lián)分析中發(fā)現(xiàn)，往往多個SNPs組成了一個Haplotype，如在3號染色體上，從SNP3818到SNP3823一共6個SNP組成了一個Haplotype，這種情況同時也在其他染色體上看到。但由于126個材料中某些材料在一個Haplotype的幾個SNPs中的個別SNP位點為雜合子，所以導致單個SNP解釋的表型變異不一致(表4)，但峰值SNP都在相應的Haplotype中。其次關聯(lián)在20，1，4，17號染色體上(圖8)。

2.7 候選基因

由于峰值SNP都出現(xiàn)在相應的Haplotype中，因此結合單個SNP和Haplotype的基因型-表型關聯(lián)分析位點，本研究將較強關聯(lián)的位點在Soybeanbase數(shù)據(jù)庫(http：//soybase.org/tools.php)中比對分析，可能參與大豆菌核病防御機制的蛋白質見表5，大致可分為參與生物過程、細胞組成、分子功能及未知功能蛋白。

圖8 基于單個SNP的基因型-表型關聯(lián)

表3 用PLINYKv 1.07進行的基于Haplotype-表型關聯(lián)分析

表4 用PLINYKv 1.07進行的基于單個SNP-表型關聯(lián)分析

表5 候選基因

3 討論

核盤菌侵染寄主的致病機理非常復雜，目前認為核盤菌分泌的細胞壁降解酶(Cell wall degrading enzymes)和草酸在致病中起重要作用，由于草酸氧化酶(Oxalate oxidase，OXO)和草酸脫羧酶(Oxalate decarboxylase，OXDC)能把草酸降解成CO2和H2O2，從而減少草酸的危害，降低了病原菌的侵染，研究表明把從小麥和大麥中分離得到的OXO(一種類萌芽素蛋白 germin-like proteins)和OXDC轉化進大豆中能提高大豆對菌核病的抗性[25-27]。研究發(fā)現(xiàn)，核盤菌致病僅有細胞壁降解酶或草酸還是不能完全致病，需要其他因子協(xié)同才能具有致病力。大豆抗菌核病的許多研究都是通過鑒定接種菌絲體后，不同時間段大豆植株抗逆酶譜/防御相關基因的表達譜的變化和細胞壁胼胝質沉積來鑒定有哪些抗逆酶參與了此過程，前人用蛋白組學技術鑒定出超氧化物歧化酶(SOD)、多酚氧化酶(PPO)、過氧化物酶(POD)、苯丙氨酸解氨酶(PAL)等在菌核病侵染過程中抗感材料的活性有所變化[28]。本研究中沒有關聯(lián)出OXO、OXDC、SOD、PPO、POD、PAL等基因。

本研究基于Haplotype-Traits的最大關聯(lián)在17號染色體上，在此區(qū)域內的纖維素合酶基因(Glyma.17g072200)被鑒定出關聯(lián)。植物細胞壁主要由果膠質、纖維素、木質素等組成，是病原物侵入寄主植物的主要屏障。植物沒有強大的免疫系統(tǒng)，所以自身的防御機制就顯得尤為重要，細胞壁是植物抵御外來侵染的第一道屏障，在油菜中有研究報道，菌核病菌菌株的致病力與其纖維素酶活性呈正相關[29-30]。另外，Valera[31]從生理、解剖和分子技術對大豆菌核病的侵染過程進行了研究，感染不同階段的大豆轉錄組測序研究證明感性品種的木質素生物合成與JA/ET-related defense responses關閉有關，證明了細胞壁特性與木質素酶在大豆抗菌核病中的重要作用。

基于單個SNP-Trait的最大關聯(lián)在3號染色體上，在此范圍內的候選基因有pyrroline-5-carboxylate(P5C)reductase(Glyma.03g129100)和cytochrome P450、family 86、subfamily A、polypeptide 1(Glyma.03g129200)。P5CR是廣泛存在于原核和真核生物中的一種重要的管家蛋白，其主要功能是催化脯氨酸生物合成的最后一步反應，即在NAD(P)H的作用下，將吡咯啉-5-羧酸(P5C)轉化為脯氨酸的反應。許多研究報道，在應激引起的損害中，脯氨酸含量在植物體內聚集[32-33]。植物中P450在寄主-病原物互作中的表達已經(jīng)有較多證據(jù)[34-35]。

用全基因組SNP關聯(lián)研究進行大豆菌核病QTLs定位顯示出諸多優(yōu)點，但QTLs定位工具的選擇是QTLs定位研究中需要考慮的主要問題之一，在本研究中利用幾種模型分析后發(fā)現(xiàn)，并不能斷言某個軟件或某種方法最適合QTLs定位，因此，在關聯(lián)分析中盡可能多用幾種軟件或多種方法分析比較，然后對共同定位的位點進行進一步研究驗證。

在與復雜性狀的關聯(lián)分析中，樣本數(shù)和標記數(shù)多少影響定位效果，從理論上講樣本數(shù)目越多且標記覆蓋基因組范圍越大，其分析結果越客觀。

本研究引入Haplotype-phenotype的關聯(lián)可能提高真正與病害相關位點的檢出率，從Haplotype關聯(lián)的基因可以看出，與病害相關基因簇都包含在一個關聯(lián)的Haplotype中，這些基因簇在病害與寄主應答中可能協(xié)同作用。通過對關聯(lián)出的基因簇的分析，可以加深對病害應答通路的認識，Haplotype很可能作為一個功能單位在起作用。另外，在單個SNP和Haplotype與表型的關聯(lián)分析中，峰值SNP都包含在其相應的Haplotype中，但由于Haplotype包含許多SNPs的遺傳信息，采用單體型比單個SNP具有更好的統(tǒng)計分析效果，能降低自由度。

以α≤0.05為閾值，3號染色體上的候選基因為Glma.03g129100、Glma.03g129200、Glma.03g129300、Glma.03g129500、Glma.03g129800、Glma.03g129900。17號染色體上為Glma.17g071300、Glma.17g072200、Glma.17g073300。目前，通過關聯(lián)分析在大豆的抗菌核病中也已經(jīng)定位了一些候選基因，如何在抗病育種中有效應用這些QTLs，需要在分析的基礎上進一步結合蛋白組學進行功能驗證，直到成為穩(wěn)定的標記才能在抗病育種中應用。