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

(1.陜西理工大學(xué) 生物科學(xué)與工程學(xué)院，陜西 漢中 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)完全免疫菌核病的大豆資源，前人研究報道，大豆種質(zhì)資源對菌核病的抗性存在著遺傳差異，有的大豆資源具有部分抗性，大豆菌核病是由多基因控制的數(shù)量性狀遺傳。數(shù)量性狀對環(huán)境條件敏感，在大田條件下目標性狀的選擇壓力難以保證一致，導(dǎo)致篩選抗菌核病種質(zhì)資源比較困難，因此，利用分子標記輔助選擇培育抗性品種顯得越來越重要。

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

1 材料和方法

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

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

1.2 分析方法

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

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

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

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

2 結(jié)果與分析

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份大豆品種進行主成分分析，結(jié)果顯示前2個主成分的方差貢獻率分別為6.24%和5.44%，Ward法聚類結(jié)果見圖4。126份供試材料聚為2類，分別為70，56個材料。

2.3 群體結(jié)構(gòu)分析

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

圖4 Ward法聚類圖

表1 群體結(jié)構(gòu)分析

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

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位點間相關(guān)(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(遺傳學(xué)上把位點間大于500 kb認為沒有連鎖關(guān)系)的r2平均值為0.73。由此看出，LD隨著2個位點間的距離增大而衰減。從單條染色體和全基因組r2和2個SNPs位點間距離的關(guān)系圖也可以看出隨著物理距離的增加，2個位點間的r2越小，即連鎖強度越小，甚至沒有(圖5，6)。

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

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

2.5 單體型分析

SNPs位點并不是獨立遺傳的，而是在染色體上傾向于以一個整體遺傳給后代，成組遺傳的SNPs位點在一代又一代的遺傳中絕少發(fā)生重組，于是，這樣的一組SNPs位點類型也就是單體型，也就是說相鄰位點是以單體型為單位傳遞的(前提是位點間的LD很強)，當然，單體型分析在統(tǒng)計學(xué)上的作用主要是減少多個位點組合分析時的自由度，從而增加研究的效能，例如，如果要分析3個二等位多態(tài)位點的組合和疾病/病害的關(guān)系，理論上有9種不同的組合，但是因為這3個位點間不是獨立的，在群體中可能只存在3種常見的單體型，如果3個位點間獨立，也就是連鎖平衡，那么群體中的單體型就會有9種，這時候做單體型分析是沒有意義的，因此，進行單體型分析時最好以單體域Block為基礎(chǔ)。在PLINKv 1.07里有2種單體域推斷方法，Sliding windows和Haplotype list，Haplotype list要先計算Blocks，然后利用算出的Haplotype和Phenotype之間的關(guān)聯(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范圍內(nèi)注釋)，最短跨越2 bp。

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

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

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

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

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

2.7 候選基因

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

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

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

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

表5 候選基因

3 討論

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

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

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

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

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

本研究引入Haplotype-phenotype的關(guān)聯(lián)可能提高真正與病害相關(guān)位點的檢出率，從Haplotype關(guān)聯(lián)的基因可以看出，與病害相關(guān)基因簇都包含在一個關(guān)聯(lián)的Haplotype中，這些基因簇在病害與寄主應(yīng)答中可能協(xié)同作用。通過對關(guān)聯(lián)出的基因簇的分析，可以加深對病害應(yīng)答通路的認識，Haplotype很可能作為一個功能單位在起作用。另外，在單個SNP和Haplotype與表型的關(guān)聯(lián)分析中，峰值SNP都包含在其相應(yīng)的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。目前，通過關(guān)聯(lián)分析在大豆的抗菌核病中也已經(jīng)定位了一些候選基因，如何在抗病育種中有效應(yīng)用這些QTLs，需要在分析的基礎(chǔ)上進一步結(jié)合蛋白組學(xué)進行功能驗證，直到成為穩(wěn)定的標記才能在抗病育種中應(yīng)用。