海南短日照條件下谷子穗部性狀的全基因組關(guān)聯(lián)分析

賈小平,張 博，董志平，全建章，王永芳，張小梅，袁璽壘，李劍峰，戴凌峰

(1.河南科技大學(xué) 農(nóng)學(xué)院，河南 洛陽 471203； 2.河北農(nóng)林科學(xué)院 谷子研究所/國家谷子改良中心，河北 石家莊 050031)

谷子(Setariaitalica)是我國的傳統(tǒng)作物，因營養(yǎng)組成均衡而受到人們的喜愛。同時(shí)，其較小的基因組(約510 Mb)、嚴(yán)格自花授粉、抗旱耐瘠、高光合的特點(diǎn)符合理想的旱作農(nóng)業(yè)C4作物模型的要求，目前被作為能源作物的模型受到國內(nèi)外重視[1-4]。谷子產(chǎn)量相對較低一直是影響其在我國實(shí)現(xiàn)產(chǎn)業(yè)化發(fā)展的重要因素。改良谷子的穗部性狀一直是谷子高產(chǎn)育種的主要目標(biāo)之一。對谷子穗部性狀的大量田間調(diào)查鑒定研究表明，谷子的莖稈性狀、穗部性狀與產(chǎn)量有著較為密切的關(guān)系[5]；穗部性狀存在豐富的變異，性狀間存在相互制約依從的關(guān)系，通過良種選育和改進(jìn)栽培措施等方法改善這些性狀，使其達(dá)到一定的平衡，才能培育出高產(chǎn)品種[6]；谷子穗部性狀與抗倒伏能力存在一定的關(guān)系，特別是穗長與谷子倒數(shù)第2節(jié)倒伏指數(shù)存在顯著相關(guān)[7]；通過穗體積和穗側(cè)面積可以精準(zhǔn)地預(yù)測出谷子的穗粒質(zhì)量[8]。產(chǎn)量與穗粗、穗長、單穗質(zhì)量、千粒質(zhì)量間存在顯著正相關(guān)[9]，與谷子穗粒質(zhì)量的關(guān)聯(lián)度由高到低依次為穗質(zhì)量、碼粒數(shù)、小區(qū)產(chǎn)量、穗粗、穗長、千粒質(zhì)量、穗碼數(shù)和株高[10]。

穗部性狀是構(gòu)成作物產(chǎn)量的重要因子，了解穗部性狀的遺傳基礎(chǔ)有利于通過分子設(shè)計(jì)育種的方法將各性狀有利等位基因聚合在一起，發(fā)揮作物的產(chǎn)量潛能。近些年，基于經(jīng)典的雙親作圖群體的谷子穗部性狀QTL定位研究已有報(bào)道[11-14]，但這種方法產(chǎn)生的變異局限于2個(gè)親本間，信息量有限，而且由于重組事件比較少，定位的區(qū)間一般都比較大。全基因組關(guān)聯(lián)分析采用高密度的DNA標(biāo)記對自然群體產(chǎn)生的性狀變異進(jìn)行全基因組掃描，發(fā)掘控制性狀的遺傳位點(diǎn)，克服了經(jīng)典遺傳作圖的不足，是挖掘有利等位基因的重要手段，在作物遺傳育種中發(fā)揮著越來越重要的作用[15]。目前，全基因組關(guān)聯(lián)分析已經(jīng)廣泛應(yīng)用于大田作物、蔬菜等重要性狀的QTL定位[16-20]。谷子主要農(nóng)藝性狀關(guān)聯(lián)分析方面的報(bào)道較少，Gupta等[21]調(diào)查了184份谷子材料的產(chǎn)量相關(guān)性狀，用50對SSR標(biāo)記進(jìn)行了性狀與標(biāo)記的關(guān)聯(lián)分析，結(jié)果發(fā)現(xiàn)，有2個(gè)標(biāo)記與旗葉寬和籽粒產(chǎn)量有較強(qiáng)關(guān)聯(lián)，且位于5號染色體。谷子全基因組測序的完成為開展基于重測序的全基因組關(guān)聯(lián)分析提供了方便[22-23]，Jia等[24]對不同的谷子資源進(jìn)行了重測序，鑒定了258萬個(gè)SNP，并使用80萬個(gè)高質(zhì)量的SNP標(biāo)記構(gòu)建谷子基因組的單倍型圖譜，獲得了多個(gè)地理環(huán)境下穩(wěn)定表達(dá)的QTL位點(diǎn)。谷子屬于短日照作物，對光周期比較敏感，隨著日照的延長，抽穗延遲，營養(yǎng)生長增強(qiáng)，株高及穗部性狀發(fā)生較大變化，而在短日照條件下可以排除光周期的影響，使一些性狀得以穩(wěn)定表達(dá)，進(jìn)而充分展示品種間的遺傳差異。鑒于此，在海南短日照條件下調(diào)查了98份不同生態(tài)區(qū)谷子品種的穗部性狀，并對這些谷子進(jìn)行基因組重測序，利用重測序開發(fā)的SNP標(biāo)記對谷子5個(gè)穗部性狀進(jìn)行全基因組關(guān)聯(lián)分析，發(fā)掘短日照條件下與穗部性狀相關(guān)的SNP突變位點(diǎn)，為揭示谷子穗部性狀遺傳機(jī)制、開展谷子高產(chǎn)分子育種奠定基礎(chǔ)。

1 材料和方法

1.1 材料

供試98份谷子材料中，89份來自河南、河北、山東、山西、陜西、甘肅、寧夏、內(nèi)蒙古、吉林、遼寧、黑龍江、新疆、湖南、四川等地區(qū)，9份來自美國、法國、日本、朝鮮、南非及國際半干旱研究所(表1)。98份谷子材料種子均由中國農(nóng)業(yè)科學(xué)院作物科學(xué)研究所/國家作物遺傳改良中心、河北農(nóng)林科學(xué)院谷子研究所提供。

表1 98份谷子材料來源信息

1.2 谷子材料的種植及穗部性狀的調(diào)查

將98份谷子材料于2015年10月中旬播種于海南省樂東縣九所鎮(zhèn)試驗(yàn)田，株距3～5 cm，行距50 cm，每品種2行，行長2 m。2016年2月谷子成熟后將谷穗自然風(fēng)干，每個(gè)谷子材料取中部10株調(diào)查穗長、穗粗、穗碼數(shù)、千粒質(zhì)量、穗粒質(zhì)量5個(gè)性狀，取平均值作為最終測定值。

1.3 谷子基因組DNA的提取及重測序

谷子材料長到3葉期時(shí)，采集幼嫩葉片利用2×CTAB法提取基因組DNA，提取的DNA滿足OD260/OD280=1.8～2.0，且沒有明顯降解，送上海美吉生物技術(shù)公司進(jìn)行重測序。

1.4 SNP標(biāo)記的開發(fā)

運(yùn)用BWA軟件將測序片段與參考基因組序列進(jìn)行比對，利用GATK軟件進(jìn)行SNP檢測。對檢測出的SNP位點(diǎn)用Samtools提供的vcfutils工具進(jìn)行過濾，過濾的條件為：測序深度≥3，質(zhì)量值Qual>50.0。過濾后的SNP用ANNOVAR軟件進(jìn)行注釋。

1.5 群體結(jié)構(gòu)分析及連鎖不平衡分析

選擇均勻分布于谷子9條染色體上的27 000個(gè)SNP標(biāo)記，用STRUCTURE Version 1.8軟件對98份谷子材料進(jìn)行群體結(jié)構(gòu)分析，首先假設(shè)98份谷子材料的分群數(shù)(k值)為7，進(jìn)行聚類，根據(jù)CV error最低點(diǎn)對應(yīng)的k值來確定最佳分群數(shù)。SNP間的連鎖不平衡分析用VCFtools軟件來計(jì)算連鎖不平衡系數(shù)(r2)，r2越接近1代表連鎖關(guān)系越強(qiáng)。

1.6 5個(gè)穗部性狀與SNP標(biāo)記間全基因組關(guān)聯(lián)分析

采用GAPIT軟件中的混合線性模型(Mixed linear model,MLM)，以群體結(jié)構(gòu)為固定效應(yīng)，個(gè)體親緣關(guān)系為隨機(jī)效應(yīng)，校正群體結(jié)構(gòu)和個(gè)體間親緣關(guān)系的影響，進(jìn)行性狀與標(biāo)記間的關(guān)聯(lián)分析。對篩選出的關(guān)聯(lián)SNP標(biāo)記所在基因與GO、NR、SwissProt、COG、KEGG等功能數(shù)據(jù)庫比對，尋找性狀候選區(qū)域內(nèi)的基因。

2 結(jié)果與分析

2.1 群體穗部性狀表現(xiàn)及性狀間的相關(guān)性分析

在海南短日照條件下，所調(diào)查的5個(gè)穗部性狀中有3個(gè)在98份谷子品種中存在較豐富的變異，分別是穗長、穗碼數(shù)、穗粒質(zhì)量，變幅分別為3.8～19.8 cm、28.9～97.1個(gè)、1.47～13.45 g；穗粗和千粒質(zhì)量的變幅較窄，分別在3.3～9.0 cm、1.96～4.56 g。對5個(gè)穗部性狀進(jìn)行相關(guān)性分析(表2)可見，穗粒質(zhì)量與穗長、穗粗、穗碼數(shù)均存在極顯著正相關(guān)(P<0.01)；穗長與穗碼數(shù)也存在極顯著正相關(guān)(P<0.01)。

表2 谷子穗部性狀的相關(guān)性分析

注：**表示在0.01水平上極顯著相關(guān)。

2.2 利用重測序開發(fā)SNP標(biāo)記

通過對98份谷子材料基因組重測序，開發(fā)了4 482 208個(gè)有效的SNP標(biāo)記，這些標(biāo)記主要分布在基因間隔區(qū)(4 037 579個(gè))、內(nèi)含子區(qū)(122 698個(gè))、轉(zhuǎn)錄起始位點(diǎn)上游區(qū)域(102 355個(gè))，少量位于剪切位點(diǎn)區(qū)域(222個(gè))。

2.3 群體結(jié)構(gòu)分析及SNP標(biāo)記間的連鎖不平衡分析

從圖1可以看出，98份谷子材料被分為3個(gè)群，第1群主要來自河南、山東、河北3個(gè)地區(qū)，還有少數(shù)來自遼寧以及國外，共31個(gè)；第2群主要來自東北、內(nèi)蒙古、甘肅等春谷地區(qū)，還有一些來自河北，共31個(gè)；第3群主要來自國外，還有一些來自東北、湖南、四川，共12個(gè)；剩余的24個(gè)材料來源較復(fù)雜。將基因組上SNP間距與r2進(jìn)行擬合，一般用r2衰退到1/2時(shí)對應(yīng)的位點(diǎn)間的距離作為連鎖不平衡(LD)的衰退距離。當(dāng)r2衰退到1/2時(shí)LD衰退距離為47.5 kb，說明谷子的基因組LD衰退距離比較短，SNP間發(fā)生連鎖的概率小，適合進(jìn)行關(guān)聯(lián)分析研究。

谷子材料編號

2.4 穗部性狀與SNP標(biāo)記關(guān)聯(lián)分析

全基因組關(guān)聯(lián)分析過程中，個(gè)體親緣關(guān)系和群體分層是造成假關(guān)聯(lián)的主要因素，因此，采用混合線性模型進(jìn)行性狀關(guān)聯(lián)分析，以群體遺傳結(jié)構(gòu)作為固定效應(yīng)，個(gè)體親緣關(guān)系作為隨機(jī)效應(yīng)，矯正群體結(jié)構(gòu)和個(gè)體親緣關(guān)系的影響。關(guān)聯(lián)分析顯著性檢驗(yàn)值P<10-4,即-lgP>4時(shí)認(rèn)為性狀與標(biāo)記位點(diǎn)存在顯著關(guān)聯(lián)，圖2中橫線以上的點(diǎn)代表與對應(yīng)性狀存在顯著關(guān)聯(lián)的SNP位點(diǎn)。結(jié)果獲得與穗長關(guān)聯(lián)的SNP位點(diǎn)27個(gè)，這些位點(diǎn)分布于谷子所有9條染色體上；與穗粗關(guān)聯(lián)的位點(diǎn)有95個(gè)，這些位點(diǎn)分布于谷子的Chr1、Chr2、Chr3、Chr5、Chr6、Chr7、Chr8等7條染色體上；與穗碼數(shù)關(guān)聯(lián)的位點(diǎn)有18個(gè)，這些位點(diǎn)分布于Chr2、Chr3、Chr4、Chr6、Chr8、Chr9等6條染色體上；與穗粒質(zhì)量關(guān)聯(lián)的位點(diǎn)有22個(gè)，這些位點(diǎn)分布于Chr1、Chr2、Chr3、Chr5、Chr6、Chr7、Chr8等7條染色體上；與千粒質(zhì)量關(guān)聯(lián)的位點(diǎn)有28個(gè)，這些位點(diǎn)分布于Chr1、Chr2、Chr3、Chr4、Chr5、Chr7、Chr8等7條染色體上。

A.穗長； B.穗粗； C.穗粒質(zhì)量； D.千粒質(zhì)量； E.穗碼數(shù)

2.5 候選基因的預(yù)測

從與性狀顯著關(guān)聯(lián)SNP位點(diǎn)所在的連鎖不平衡區(qū)域檢索到與穗長相關(guān)聯(lián)的基因有24個(gè)，與穗粗相關(guān)聯(lián)的基因有33個(gè)，與穗碼數(shù)相關(guān)聯(lián)的基因有18個(gè)，與穗粒質(zhì)量相關(guān)聯(lián)的基因有18個(gè)，與千粒質(zhì)量相關(guān)聯(lián)的基因有22個(gè)。表3為這些基因在GO、NR、SwissProt、COG、KEGG等功能數(shù)據(jù)庫的注釋情況。根據(jù)注釋情況，從候選基因中進(jìn)一步篩選出一些關(guān)鍵的目標(biāo)基因，其中選出與穗長關(guān)聯(lián)的目標(biāo)基因4個(gè)，分布在Chr3、Chr5、Chr7、Chr8等4條染色體上，這些基因主要功能包括信號傳遞、泛素-蛋白酶系統(tǒng)介導(dǎo)的蛋白質(zhì)降解、蛋白質(zhì)運(yùn)輸及防御反應(yīng)；與穗粗關(guān)聯(lián)的基因有9個(gè)，分布在Chr1、Chr2、Chr3、Chr5、Chr6、Chr8等6條染色體上，這些基因主要功能包括信號轉(zhuǎn)導(dǎo)、木質(zhì)素合成、RNA加工、赤霉素合成、賴氨酸組氨酸轉(zhuǎn)運(yùn)等；與穗碼數(shù)關(guān)聯(lián)的基因有5個(gè)，分布在Chr2、Chr3、Chr4、Chr8等4條染色體，主要功能包括泛素-蛋白酶系統(tǒng)介導(dǎo)的蛋白質(zhì)降解、光系統(tǒng)Ⅱ合成、信號傳遞、防御反應(yīng)、脂肪酸合成；與穗粒質(zhì)量關(guān)聯(lián)的基因有8個(gè)，分布在Chr1、Chr2、Chr3、Chr5、Chr6、Chr7等6條染色體上，主要功能包括木質(zhì)素合成、信號轉(zhuǎn)導(dǎo)、防御反應(yīng)、細(xì)胞運(yùn)動(dòng)、泛素-蛋白酶系統(tǒng)介導(dǎo)的蛋白質(zhì)降解；與千粒質(zhì)量關(guān)聯(lián)的基因有7個(gè)，分布在Chr1、Chr2、Chr3、Chr4等4條染色體，主要功能有信號轉(zhuǎn)導(dǎo)、防御反應(yīng)、脂肪酸合成、鋅離子運(yùn)輸、泛素-蛋白酶系統(tǒng)介導(dǎo)的蛋白質(zhì)降解等(表4)。

表3 關(guān)聯(lián)區(qū)域檢測候選基因數(shù)目及其在各數(shù)據(jù)庫的注釋情況 個(gè)

從表4可以看出，在所篩選出的候選基因中，位于6號染色體上的與穗粒質(zhì)量關(guān)聯(lián)的基因XB3是一個(gè)重要的候選基因，該基因與水稻中報(bào)道的環(huán)型蛋白同源，具有E3泛素連接酶活性，參與泛素-蛋白酶系統(tǒng)介導(dǎo)的蛋白質(zhì)降解過程。在水稻研究中，已證明該基因可以調(diào)控籽粒質(zhì)量和寬度，缺失該基因?qū)е伦蚜＜?xì)胞數(shù)目增加，加速籽粒灌漿[25]。本研究中，與XB3關(guān)聯(lián)的性狀是穗粒質(zhì)量，與單個(gè)籽粒的平均質(zhì)量密切相關(guān)，因此推測，XB3基因與水稻中報(bào)道的環(huán)狀E3泛素連接酶具有相同的功能，可以調(diào)控谷子的籽粒質(zhì)量。而另一個(gè)位于3號染色體的PP2C基因與谷子千粒質(zhì)量關(guān)聯(lián)，在大豆的研究中發(fā)現(xiàn)，PP2C是與百粒質(zhì)量密切相關(guān)的基因，該基因等位變異PP2C-1能與油菜素內(nèi)酯信號通路的轉(zhuǎn)錄因子相互作用，通過去磷酸化激活這些轉(zhuǎn)錄因子，促進(jìn)下游控制種子大小的基因表達(dá)進(jìn)而提高籽粒質(zhì)量[26]。推測PP2C基因在谷子籽粒質(zhì)量調(diào)控中也起到重要作用。

表4 篩選出的與5個(gè)穗部性狀關(guān)聯(lián)的候選基因

3 結(jié)論與討論

隨著測序技術(shù)成本的下降，利用重測序開發(fā)SNP標(biāo)記進(jìn)行重要性狀的全基因組關(guān)聯(lián)分析已經(jīng)廣泛應(yīng)用于農(nóng)作物QTL定位研究。谷子全基因組測序已經(jīng)完成，為開展基于重測序的全基因組關(guān)聯(lián)分析提供了條件。Jia等[24]通過對960份谷子材料進(jìn)行重測序，開發(fā)了250多萬個(gè)SNP標(biāo)記，并用其中80萬個(gè)SNP標(biāo)記構(gòu)建了谷子基因組單倍型圖譜。本研究開發(fā)的標(biāo)記數(shù)量高于Jia等[24]報(bào)道的數(shù)量，表明本研究所篩選的谷子材料雖然數(shù)量有限，但代表了較高的遺傳多樣性水平，而且本研究測序深度平均覆蓋谷子5倍基因組，高于Jia等[24]報(bào)道的3倍測序深度，這也是獲得較高SNP檢出率的一個(gè)原因。高密度的SNP標(biāo)記使谷子基因組LD衰退距離變小，本研究中谷子基因組的LD衰退距離為47.5 kb，低于Jia等[24]報(bào)道的結(jié)果(約100 kb)，因此性狀定位區(qū)間較窄，精度相對提高。

目前，有關(guān)谷子穗部性狀QTL定位的研究主要是利用雜交分離群體進(jìn)行的。利用SSR標(biāo)記對F2群體進(jìn)行QTL定位，將控制穗長的QTL定位于谷子1、2、7號染色體上，將控制主穗粗的QTL定位在1、4、5號染色體上，將控制單株穗粒質(zhì)量的QTL定位到谷子1、4號染色體上，將控制千粒質(zhì)量的QTL定位到4、5號染色體上[11-12]。近期，有學(xué)者利用測序產(chǎn)生的SNP標(biāo)記對雜交分離群體進(jìn)行農(nóng)藝性狀的QTL定位分析，將控制穗長的QTL定位到1、2、4、5、8、9號染色體，將控制穗粗的QTL定位到1、2、4、5、7、9號染色體，將控制穗碼數(shù)的QTL定位到2、4、9號染色體，將控制千粒質(zhì)量的QTL定位到3、4、5、7號染色體[13-14]。本研究利用全基因組關(guān)聯(lián)分析定位到與穗長、穗粗、穗碼數(shù)、穗粒質(zhì)量、千粒質(zhì)量關(guān)聯(lián)的位點(diǎn)數(shù)分別為27、95、18、22、28個(gè)；控制穗長的QTL除在前人報(bào)道的所有染色體都檢測到，還在3、6號染色體檢測到關(guān)聯(lián)位點(diǎn)；控制穗粗的QTL位點(diǎn)沒有在前人報(bào)道的4、9號染色體檢測到，但是在前人沒有報(bào)道的3、6、8號染色體上檢測到新的QTL位點(diǎn)；控制穗碼數(shù)的QTL除了在前人報(bào)道的染色體都檢測到，還在3、6、8號染色體上檢測到新的QTL；控制穗粒質(zhì)量的QTL在前人報(bào)道的4號染色體上沒有檢測到，但是在前人沒有報(bào)道的2、3、5、6、7、8號染色體上檢測到新的QTL位點(diǎn)；控制千粒質(zhì)量的QTL沒有在前人報(bào)道的7號染色體檢測出，但是在前人沒有報(bào)道的1、2、8號染色體上檢測到新的QTL位點(diǎn)。總體上，本研究檢測到的QTL位點(diǎn)多于前人基于雜交分離群體定位到的QTL位點(diǎn)，基于自然群體各性狀的變異遠(yuǎn)多于來自2個(gè)親本的變異。而對于本研究沒有檢測出的QTL，有可能與生長環(huán)境有關(guān)，本研究是在海南短日照條件下定位各穗部性狀，而前人研究有些是在中、長日照條件下進(jìn)行的，因此，一些對光周期敏感的位點(diǎn)可能在本研究中未檢測到。

通過對關(guān)聯(lián)位點(diǎn)所在區(qū)域進(jìn)行候選基因查找、注釋，每個(gè)穗部性狀篩選出4～9個(gè)候選基因，這些候選基因中參與泛素-蛋白酶系統(tǒng)介導(dǎo)的蛋白質(zhì)降解過程、防御反應(yīng)、木質(zhì)素合成的比較多，還有一些參與脂肪酸、赤霉素合成及信號傳遞，少數(shù)參與氨基酸、鋅離子的運(yùn)輸。此外，本研究在谷子3號染色體發(fā)現(xiàn)一個(gè)與千粒質(zhì)量關(guān)聯(lián)的候選基因PP2C，該基因編碼蛋白磷酸酶，最近的研究表明，PP2C基因是控制大豆百粒質(zhì)量的主效基因[26]，因此，推測該基因也可能對谷子千粒質(zhì)量有重要作用，但仍需要進(jìn)一步進(jìn)行基因的表達(dá)分析和轉(zhuǎn)化驗(yàn)證。在水稻研究中也發(fā)現(xiàn)一些控制粒質(zhì)量的基因參與泛素介導(dǎo)的蛋白質(zhì)降解過程[25，27]，因此，推測泛素-蛋白酶系統(tǒng)介導(dǎo)的蛋白質(zhì)降解過程可能廣泛參與作物穗部性狀的調(diào)控。

綜上，對谷子材料重測序開發(fā)了4 482 208 個(gè)高質(zhì)量的SNP標(biāo)記，LD分析表明谷子LD衰退距離較短(47.5 kb)。于海南短日照條件下調(diào)查了谷子的穗長、穗粗、穗碼數(shù)、穗粒質(zhì)量、千粒質(zhì)量5個(gè)性狀，對5個(gè)穗部性狀進(jìn)行全基因組關(guān)聯(lián)分析，每個(gè)性狀篩選出4～9個(gè)候選基因，這些基因主要參與泛素-蛋白酶系統(tǒng)介導(dǎo)的蛋白質(zhì)降解過程、防御反應(yīng)、信號傳導(dǎo)、木質(zhì)素合成，還有一些參與脂肪酸合成及氨基酸、鋅離子運(yùn)輸。本研究為進(jìn)一步克隆谷子穗部性狀相關(guān)基因奠定了基礎(chǔ)。