大豆抗菌核病的全基因組關(guān)聯(lián)研究

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

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

1 材料和方法

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

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

1.2 分析方法

1.2.1 表型值 用棉墊接種法接種核盤菌菌絲進(jìn)行表型鑒定[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é)果，不當(dāng)?shù)膮?shù)選擇會(huì)產(chǎn)生假陽性或假陰性結(jié)果。本研究位點(diǎn)間的關(guān)聯(lián)程度r2值0.8，LOD取值3，置信區(qū)間為95%，MAF值取0.01。由于大豆為自交作物，偏離哈迪-溫伯格平衡的SNPs很可能是與病害關(guān)聯(lián)的易感位點(diǎn)，因此，本研究不考慮哈迪-溫伯格平衡，得到30 125個(gè)SNPs(圖2)。126個(gè)材料共有3.7958E6個(gè)位點(diǎn)，其中雜合基因型位點(diǎn)有248 456個(gè)，占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)對(duì)126個(gè)供試材料進(jìn)行聚類(圖3)，聚為3大類，分別為67，45，14個(gè)材料。

圖3 基于NJ法的聚類圖

2.2 PC分析

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

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

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

圖4 Ward法聚類圖

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

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

2.4 連鎖不平衡分析

用PLINKv 1.07命令“--r2”運(yùn)行LD，過濾掉r2小于0.2的數(shù)據(jù)，共有82 911對(duì)，r2平均值為r0.73，r2大于等于0.8的有46 655對(duì)(占56.27%)，2等于1的有9 391對(duì)(占11.33%)。2個(gè)SNPs位點(diǎn)間相關(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。兩位點(diǎn)間距離≤500 kb(遺傳學(xué)上把位點(diǎn)間大于500 kb認(rèn)為沒有連鎖關(guān)系)的r2平均值為0.73。由此看出，LD隨著2個(gè)位點(diǎn)間的距離增大而衰減。從單條染色體和全基因組r2和2個(gè)SNPs位點(diǎn)間距離的關(guān)系圖也可以看出隨著物理距離的增加，2個(gè)位點(diǎn)間的r2越小，即連鎖強(qiáng)度越小，甚至沒有(圖5，6)。

圖5 一號(hào)染色體上的LD衰減圖

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

2.5 單體型分析

SNPs位點(diǎn)并不是獨(dú)立遺傳的，而是在染色體上傾向于以一個(gè)整體遺傳給后代，成組遺傳的SNPs位點(diǎn)在一代又一代的遺傳中絕少發(fā)生重組，于是，這樣的一組SNPs位點(diǎn)類型也就是單體型，也就是說相鄰位點(diǎn)是以單體型為單位傳遞的(前提是位點(diǎn)間的LD很強(qiáng))，當(dāng)然，單體型分析在統(tǒng)計(jì)學(xué)上的作用主要是減少多個(gè)位點(diǎn)組合分析時(shí)的自由度，從而增加研究的效能，例如，如果要分析3個(gè)二等位多態(tài)位點(diǎn)的組合和疾病/病害的關(guān)系，理論上有9種不同的組合，但是因?yàn)檫@3個(gè)位點(diǎn)間不是獨(dú)立的，在群體中可能只存在3種常見的單體型，如果3個(gè)位點(diǎn)間獨(dú)立，也就是連鎖平衡，那么群體中的單體型就會(huì)有9種，這時(shí)候做單體型分析是沒有意義的，因此，進(jìn)行單體型分析時(shí)最好以單體域Block為基礎(chǔ)。在PLINKv 1.07里有2種單體域推斷方法，Sliding windows和Haplotype list，Haplotype list要先計(jì)算Blocks，然后利用算出的Haplotype和Phenotype之間的關(guān)聯(lián)程度進(jìn)行分析，這種運(yùn)算比Sliding windows靈活，因?yàn)镾liding windows只能按SNP順序依次排列分析Haplotype，不能按實(shí)際跳躍分析。在比較了2種方法后，本研究采用“--blocks”命令推斷Blocks。20條染色體共有2 863個(gè)Haplotype，它們?cè)?0條染色體上的分布見圖7，其中18號(hào)染色體上的Haplotype最多(230個(gè))，11號(hào)染色體上的最少(57個(gè))。一個(gè)Haplotype最多的由36個(gè)SNPs組成，最少2個(gè)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號(hào)染色體上的4個(gè)SNPs組成的Haplotype解釋表型變異最大(17.56%)。其余依次為1，3，4，20號(hào)染色體上的Haplotype(表3)。

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

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

2.7 候選基因

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

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

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

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

表5 候選基因

3 討論

核盤菌侵染寄主的致病機(jī)理非常復(fù)雜，目前認(rèn)為核盤菌分泌的細(xì)胞壁降解酶(Cell wall degrading enzymes)和草酸在致病中起重要作用，由于草酸氧化酶(Oxalate oxidase，OXO)和草酸脫羧酶(Oxalate decarboxylase，OXDC)能把草酸降解成CO2和H2O2，從而減少草酸的危害，降低了病原菌的侵染，研究表明把從小麥和大麥中分離得到的OXO(一種類萌芽素蛋白 germin-like proteins)和OXDC轉(zhuǎn)化進(jìn)大豆中能提高大豆對(duì)菌核病的抗性[25-27]。研究發(fā)現(xiàn)，核盤菌致病僅有細(xì)胞壁降解酶或草酸還是不能完全致病，需要其他因子協(xié)同才能具有致病力。大豆抗菌核病的許多研究都是通過鑒定接種菌絲體后，不同時(shí)間段大豆植株抗逆酶譜/防御相關(guān)基因的表達(dá)譜的變化和細(xì)胞壁胼胝質(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號(hào)染色體上，在此區(qū)域內(nèi)的纖維素合酶基因(Glyma.17g072200)被鑒定出關(guān)聯(lián)。植物細(xì)胞壁主要由果膠質(zhì)、纖維素、木質(zhì)素等組成，是病原物侵入寄主植物的主要屏障。植物沒有強(qiáng)大的免疫系統(tǒng)，所以自身的防御機(jī)制就顯得尤為重要，細(xì)胞壁是植物抵御外來侵染的第一道屏障，在油菜中有研究報(bào)道，菌核病菌菌株的致病力與其纖維素酶活性呈正相關(guān)[29-30]。另外，Valera[31]從生理、解剖和分子技術(shù)對(duì)大豆菌核病的侵染過程進(jìn)行了研究，感染不同階段的大豆轉(zhuǎn)錄組測(cè)序研究證明感性品種的木質(zhì)素生物合成與JA/ET-related defense responses關(guān)閉有關(guān)，證明了細(xì)胞壁特性與木質(zhì)素酶在大豆抗菌核病中的重要作用。

基于單個(gè)SNP-Trait的最大關(guān)聯(lián)在3號(hào)染色體上，在此范圍內(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)。許多研究報(bào)道，在應(yīng)激引起的損害中，脯氨酸含量在植物體內(nèi)聚集[32-33]。植物中P450在寄主-病原物互作中的表達(dá)已經(jīng)有較多證據(jù)[34-35]。

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

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

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

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