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

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

倒伏是農(nóng)作物的莖稈從相對的直立形態(tài)發(fā)生不可逆錯(cuò)位的現(xiàn)象[1]，作物的倒伏將會(huì)對群體的結(jié)構(gòu)和莖稈的輸導(dǎo)能力造成一定程度的破壞，并且能夠顯著地降低光合效率和光合作用產(chǎn)物向穗部運(yùn)輸?shù)哪芰?，最終造成作物大量減產(chǎn)[2]。谷子作為一種重要的雜糧作物，具有光合效率高、耐旱耐瘠、營養(yǎng)均衡等特點(diǎn)，近些年廣受人們的喜愛，然而，由于莖稈偏軟、種植密度較大，在遭遇風(fēng)雨氣候等因素時(shí)，谷子容易發(fā)生倒伏，對產(chǎn)量和商品價(jià)值造成極大的影響[3]。造成谷子倒伏的因素有外因與內(nèi)因兩個(gè)方面：內(nèi)因是指品種的區(qū)別，不同谷子品種間抗倒伏性存在巨大的差異；外因則主要包含具有災(zāi)害性的氣候(例如暴雨、冰雹、大風(fēng)等)、施肥情況、病蟲害、栽培密度[4]。施肥和栽培密度配合可以一定程度上改善水稻莖稈理化特性，有效提高抗倒伏性[5]。通過分析施肥對谷子抗倒伏性和產(chǎn)量的影響，發(fā)現(xiàn)不同品種對施肥處理的反應(yīng)不同，有些品種抗倒伏能力提高而有些卻降低，這主要取決于品種對肥料的需求與否[6]。也有研究發(fā)現(xiàn)矮化處理后能有效防止谷子倒伏[7]。

目前有關(guān)谷子抗倒伏性評價(jià)的研究較多，如袁立新[2]研究發(fā)現(xiàn)，谷子的株型與莖稈的倒伏性狀關(guān)系十分密切，莖部性狀中的株高、節(jié)間長度和穗部性狀中的穗長、穗柄長、穗彎曲度與谷子的抗倒伏性呈現(xiàn)出顯著負(fù)相關(guān)，而莖部性狀莖粗和谷子根量則與谷子的抗倒伏性呈現(xiàn)出正相關(guān)。研究結(jié)果表明，株型相關(guān)性狀的選擇對于選育出抗倒伏性強(qiáng)的谷子品種具有重要意義；袁志華等[8]建立的力學(xué)模型可用來評估谷子的抗倒伏能力；劉艷麗等[9]應(yīng)用倒伏系數(shù)研究表明，倒伏性狀對于谷子的產(chǎn)量、品質(zhì)和機(jī)械化收獲等有著明顯的影響，同時(shí)確定了對谷子倒伏性有著重要影響的一些關(guān)鍵性狀指標(biāo)；賈小平等[10]進(jìn)行了不同谷子品種(系)的抗倒性分析，發(fā)現(xiàn)谷子的倒數(shù)第2莖節(jié)倒伏指數(shù)對谷子植株的倒伏性起到關(guān)鍵作用。

當(dāng)前對谷子倒伏性的研究主要集中在形態(tài)學(xué)特征、理化特性方面，有關(guān)抗倒伏性的QTL定位等分子水平的研究尚未見報(bào)道，在模式作物水稻中相關(guān)研究報(bào)道較多[11-13]。近些年谷子基因組學(xué)研究取得了重大進(jìn)展，全基因組測序工作已經(jīng)分別由中國和美國完成[14-15]，而基于重測序的全基因組關(guān)聯(lián)分析方法已經(jīng)廣泛應(yīng)用于植物QTL定位研究[16-17]。Jia等[18]首次對916份谷子材料進(jìn)行了重測序，開發(fā)了250多萬個(gè)SNP標(biāo)記，并用其中80萬個(gè)SNP標(biāo)記構(gòu)建了谷子基因組單倍型圖譜，對47個(gè)性狀進(jìn)行了全基因組關(guān)聯(lián)分析，為開展谷子全基因組關(guān)聯(lián)分析研究提供了借鑒。本研究選擇98份來自國內(nèi)外不同地區(qū)的谷子材料，在洛陽、吉林兩地調(diào)查了這98份谷子材料的倒伏率，并對這些谷子材料進(jìn)行基因組重測序，開發(fā)大量的SNP標(biāo)記，進(jìn)行突變位點(diǎn)與抗倒伏性的全基因組關(guān)聯(lián)分析，發(fā)掘與抗倒伏性顯著相關(guān)的SNP位點(diǎn)，篩選抗倒伏性候選基因，為進(jìn)一步定位、克隆控制谷子抗倒伏性的基因奠定基礎(chǔ)。

1 材料與方法

1.1 試驗(yàn)材料及種植方法

本研究所用的98份谷子材料由中國農(nóng)業(yè)科學(xué)院作物科學(xué)研究所(國家作物遺傳育種改良中心)提供(表1)，98份谷子材料于2016年5月中旬分別播種于河南科技大學(xué)試驗(yàn)田(洛陽)、吉林市農(nóng)業(yè)科學(xué)院試驗(yàn)田(吉林)，株距5 cm，行距40 cm，每品種種植兩行，行長2 m。

1.2 谷子抗倒伏性的測定

在谷子成熟期調(diào)查倒伏率，每個(gè)品種除去行兩端各2株，統(tǒng)計(jì)每一品種的總株數(shù)，并仔細(xì)調(diào)查出每個(gè)品種的倒伏株數(shù)，倒伏的標(biāo)準(zhǔn)為莖稈與地面的傾角小于45°，計(jì)算倒伏率。

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

在洛陽試驗(yàn)田當(dāng)98份谷子材料長到3葉期時(shí)，采集幼嫩葉片保存于液氮中，隨后用2×CTAB法提取基因組DNA。質(zhì)量合格的DNA送上海美吉生物工程有限公司進(jìn)行基因組文庫構(gòu)建，用Illumina HiSeqTM進(jìn)行測序。

1.4 突變位點(diǎn)檢測和注釋

運(yùn)用BWA軟件將質(zhì)控后的測序片段與谷子參考基因組序列進(jìn)行比對，利用GATK進(jìn)行SNP檢測。對于檢測到的SNP位點(diǎn)，用Samtools提供的vcfutils工具對其過濾。測序深度大于等于3，即支持該變異的測序片段數(shù)目要大于等于3，質(zhì)量值Qual>50.0。然后利用ANNOVAR軟件對由多個(gè)基因組檢測出的基因變異進(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(cross validation error)最低點(diǎn)對應(yīng)的K值來確定最佳分群數(shù)。用VCFtools軟件計(jì)算標(biāo)記間連鎖不平衡系數(shù)(r2)，進(jìn)行連鎖不平衡分析，r2越接近1，代表連鎖關(guān)系越強(qiáng)。

1.6 抗倒伏性的關(guān)聯(lián)分析及候選基因預(yù)測

利用Tassel5.0軟件進(jìn)行關(guān)聯(lián)分析，采用混合線性模型(mixed linear model，MLM)，同時(shí)對群體結(jié)構(gòu)和個(gè)體親緣關(guān)系進(jìn)行校正，以減少計(jì)算誤差。以群體結(jié)構(gòu)得到的Q值和TASSEL軟件計(jì)算各材料之間的K系數(shù)矩陣作為協(xié)變量，以P Marker值來衡量關(guān)聯(lián)標(biāo)記的顯著性。用Manhattan圖來顯示關(guān)聯(lián)分析結(jié)果，在檢測到與抗倒伏性顯著相關(guān)的標(biāo)記位點(diǎn)兩側(cè)200 kb范圍內(nèi)尋找候選基因，并在GO數(shù)據(jù)庫對候選基因進(jìn)行功能注釋。

2 結(jié)果與分析

2.1 洛陽、吉林兩地谷子倒伏率調(diào)查

在洛陽地區(qū)，谷子材料倒伏率變異范圍在0～42.5%，倒伏率為0的品種最多，占有效調(diào)查品種數(shù)目的50%以上，其次為倒伏率在10%以內(nèi)的品種，約占有效調(diào)查品種數(shù)的39%，倒伏率在15%以上的品種最少，約占有效調(diào)查品種數(shù)的5%。在吉林地區(qū)，谷子材料倒伏率變異范圍在0～30%，與洛陽調(diào)查結(jié)果相似，倒伏率為0的品種最多，也達(dá)到有效調(diào)查品種數(shù)的50%以上，其次是倒伏率在10%～15%的品種，占有效調(diào)查品種數(shù)的29%，倒伏率在15%以上的品種最少，約占有效調(diào)查品種數(shù)的14%(圖1)。總體來看，洛陽、吉林兩地谷子倒伏率分布規(guī)律表現(xiàn)一定的一致性。

2.2 基于重測序谷子SNP標(biāo)記的開發(fā)

通過對98份谷子材料進(jìn)行基因組重測序檢測到9 740 706個(gè)SNP，經(jīng)過濾獲得4 482 208個(gè)高質(zhì)量的SNP標(biāo)記，這些SNP標(biāo)記中有4 037 579個(gè)位于基因間隔區(qū)，約占SNP總數(shù)的90%；其次為分布在內(nèi)含子區(qū)域、轉(zhuǎn)錄起始位點(diǎn)的上游區(qū)域以及基因下游1 kb區(qū)域的SNP位點(diǎn)，其相應(yīng)數(shù)目為122 698、102 355和81 070；剩下的SNP主要分布于外顯子區(qū)、5′非翻譯區(qū)、3′非翻譯區(qū)以及非編碼RNA的外顯子、內(nèi)含子區(qū)域(表1)。

圖1 洛陽(A)、吉林(B)兩地倒伏率分布情況Fig.1 Distribution of lodging percentage in Luoyang(A) and Jilin (B)

表1 重測序開發(fā)的SNP位點(diǎn)統(tǒng)計(jì)

2.3 98份谷子品種群體結(jié)構(gòu)分析

在SNP標(biāo)記開發(fā)的基礎(chǔ)上，通過STRUCTURE軟件分析98個(gè)樣品的群體結(jié)構(gòu)，分別假設(shè)98個(gè)樣品的分群數(shù)，即K值為1～7，進(jìn)行聚類分析。根據(jù)CV error最低點(diǎn)對應(yīng)的K值來確定最佳分群數(shù)為3(圖2-A)。以各品種的遺傳成分值Q>0.5作為標(biāo)準(zhǔn)進(jìn)行品種歸類：亞群1包括了44個(gè)品種，占總品種的44.9%，主要源自河南、河北、吉林、甘肅，還包括部分黑龍江、遼寧、陜西、山東、內(nèi)蒙古、四川、山西的品種，2個(gè)北京品系和2個(gè)日本引進(jìn)品系、1個(gè)美國引進(jìn)品系；亞群2包括了10個(gè)品種，占總品種的10.2%，主要源自國內(nèi)，包括陜西、新疆、寧夏、湖南、吉林、甘肅、河南各地區(qū)的品種；亞群3包括了38個(gè)品種，占總品種的38.8%，主要源自河北、河南、山西，還包括了黑龍江、遼寧、吉林、山東、新疆、寧夏、甘肅、內(nèi)蒙古的部分品種，以及朝鮮、美國、日本、南非各一個(gè)引進(jìn)品系?；旌闲蛠喨喊?個(gè)品種，占總品種的6.1%，主要源自黑龍江，還包括了河南、吉林、河北，以及一個(gè)法國引進(jìn)品系。98份材料的群體結(jié)構(gòu)較簡單，有效減弱了群體結(jié)構(gòu)對關(guān)聯(lián)分析所造成的干擾(圖2-B)。

2.4 連鎖不平衡分析

通過計(jì)算覆蓋谷子全基因組的4 482 208個(gè)SNP標(biāo)記間的連鎖不平衡系數(shù)(r2)來分析谷子群體全基因組的LD水平，而常用r2衰退到一半時(shí)所對應(yīng)的距離作為連鎖不平衡的衰減的數(shù)值。由圖3可知，伴隨著SNP間距的不斷增加，谷子基因組的r2在逐漸下降，r2衰退到一半時(shí)LD-decay長度為47.5 kb，說明谷子基因組平均LD衰減距離是47.5 kb，因此可以在47.5 kb范圍對應(yīng)的LD區(qū)域內(nèi)掃描抗倒伏性相關(guān)候選基因。

2.5 谷子抗倒伏性的全基因組關(guān)聯(lián)分析

在洛陽生態(tài)區(qū)，檢測到與倒伏率存在顯著關(guān)聯(lián)的SNP位點(diǎn)有764個(gè)(P<10-4)，這些位點(diǎn)分布在谷子所有9個(gè)染色體上，其中8號染色體分布的SNP位點(diǎn)最多，為221個(gè)；其次為9號染色體，達(dá)到141個(gè)；3號、5號染色體分布的SNP位點(diǎn)數(shù)也達(dá)到104、102個(gè)；6號染色體分布的SNP數(shù)目最少，為15個(gè)。在吉林生態(tài)區(qū)，檢測到24個(gè)SNP位點(diǎn)與倒伏率顯著關(guān)聯(lián)(P<10-4)，這些位點(diǎn)分布在谷子1～8號染色體上，其中5號、8號染色體分布的SNP位點(diǎn)最多，均為6個(gè)；其次為2號染色體和1號染色體，分布的SNP位點(diǎn)數(shù)分別為4個(gè)和3個(gè)；其余染色體分布的SNP數(shù)目為1～2個(gè)(圖4)。

A，K值分布曲線；B，98份谷子材料群體結(jié)構(gòu)。A，K value distribution curve; B，Population structure of 98 foxtail millet materials.圖2 九十八份谷子材料群體結(jié)構(gòu)分析Fig.2 Population structure analysis of 98 foxtail millet materials

橫坐標(biāo)代表SNP 的間距，縱坐標(biāo)代表r2值。r2衰退到一半時(shí)對應(yīng)的SNP 間距即為LD-decay。The abscissa represented the spacing of SNP; the ordinate represented the r2 value.The corresponding SNP spacing was LD-decay when r2 fell to half.圖3 LD衰減分布圖Fig.3 LD decay distribution diagram

在關(guān)聯(lián)SNP位點(diǎn)兩側(cè)區(qū)域200 kb的范圍搜索候選基因，洛陽生態(tài)區(qū)在谷子1～9號染色體上均檢測到候選基因，其中激酶類基因分布在谷子的1～8號染色體上，且數(shù)量最多，包括LRR受體類絲氨酸/蘇氨酸蛋白激酶、L-型受體激酶、亮氨酸重復(fù)受體類絲氨酸/蘇氨酸蛋白激酶、蛋白激酶G11A、細(xì)胞壁關(guān)聯(lián)受體激酶等，抗病蛋白數(shù)量也較多，分布在除1號、4號、5號、6號染色體之外的其他5條染色體上。此外，三角狀五肽蛋白、糖類分解、轉(zhuǎn)移、運(yùn)輸相關(guān)的蛋白，如糖轉(zhuǎn)運(yùn)蛋白、蔗糖運(yùn)輸?shù)鞍住?，4-α-葡聚糖分支酶、EGF結(jié)構(gòu)域特異性O(shè)-連接N-乙酰葡糖胺轉(zhuǎn)移酶、葡萄糖醛酸轉(zhuǎn)移酶、葡聚糖內(nèi)切-1，3-β-葡糖苷酶、β-1，4-甘露糖基糖蛋白4-β-N-乙酰葡糖胺基轉(zhuǎn)移酶、多聚半乳糖醛酸酶、UDP糖基轉(zhuǎn)移酶也在各染色體廣泛分布，關(guān)聯(lián)到數(shù)量較少的基因包括GA合成途徑相關(guān)酶類、E3泛素-蛋白連接酶、液泡蛋白分選相關(guān)蛋白、含有NAC結(jié)構(gòu)域的蛋白質(zhì)及細(xì)胞色素P450相關(guān)蛋白，此外還關(guān)聯(lián)到一些微管蛋白、肌動(dòng)蛋白、甲基轉(zhuǎn)移酶及生長素合成途徑的一些酶相關(guān)基因。吉林生態(tài)區(qū)在谷子1～8號染色體上均關(guān)聯(lián)到抗倒伏性相關(guān)候選基因，但是數(shù)量遠(yuǎn)少于洛陽生態(tài)區(qū)，主要集中在5號、6號、8號染色體上。吉林關(guān)聯(lián)到的候選基因同樣是激酶類最多，包括絲氨酸/蘇氨酸蛋白激酶、受體類蛋白激酶、LRR受體類絲氨酸/蘇氨酸蛋白激酶和細(xì)胞壁關(guān)聯(lián)受體。此外關(guān)聯(lián)到的候選基因還有UDP糖基轉(zhuǎn)移酶、細(xì)胞色素P450相關(guān)、甲基轉(zhuǎn)移酶、驅(qū)動(dòng)蛋白、幾丁質(zhì)酶類蛋白質(zhì)等少數(shù)基因。

洛陽、吉林兩地都關(guān)聯(lián)到較多的激酶類基因，特別是6號染色體上兩地都關(guān)聯(lián)到同一個(gè)候選基因——LRR受體類絲氨酸/蘇氨酸蛋白激酶基因(LOC101776261)。此外，在5號、6號、8號染色體上還關(guān)聯(lián)到一些位點(diǎn)不同但功能接近的蛋白激酶基因，如5號染色體上洛陽關(guān)聯(lián)到基因(LOC101754749)和吉林關(guān)聯(lián)到的基因(LOC101756923)均為受體類蛋白激酶，6號、8號染色體上兩地均關(guān)聯(lián)到一些細(xì)胞壁相關(guān)的受體激酶基因(LOC105914550、LOC101786429、LOC105913474、LOC101766879)，還有一個(gè)8號染色體上的激酶基因(LOC101760794)，在洛陽地區(qū)檢測到的也是LRR受體類絲氨酸/蘇氨酸蛋白激酶基因，與6號染色體洛陽、吉林兩地共同檢測到的激酶基因類似。還有一些其他類型的基因同樣雖然在兩地位于不同的染色體上，但是功能很接近，如甲基轉(zhuǎn)移酶LOC101766453和LOC101771920、UDP糖基轉(zhuǎn)移酶LOC101764556和LOC101753540等(表2)。

圖4 洛陽(A)、吉林(B)關(guān)聯(lián)分析結(jié)果的曼哈頓圖Fig.4 Manhattan map of correlation analysis results in Luoyang (A) and Jilin (B)

表2 洛陽、吉林關(guān)聯(lián)到的抗倒伏性候選基因

續(xù)表2

續(xù)表2

續(xù)表2

3 討論

倒伏是影響作物產(chǎn)量和品質(zhì)的重要因素，定位控制倒伏性的基因位點(diǎn)可以揭示抗倒伏性的遺傳機(jī)制，為利用分子育種技術(shù)培育抗倒伏作物新品種奠定基礎(chǔ)。目前有關(guān)抗倒伏性的QTL定位、基因克隆方面的研究在水稻中報(bào)道較多。研究表明，控制水稻株型、莖稈性狀、細(xì)胞壁性狀的相關(guān)基因和抗倒伏性關(guān)系密切，如和株高相關(guān)的赤霉素合成途徑基因SD1突變可以導(dǎo)致植株矮化，提高抗倒伏能力[19]；DEP1為控制水稻直立穗型的主要基因，具有直立穗型的水稻品種抗倒伏性明顯增加[20]。與水稻莖稈物理特性相關(guān)的基因SMOS1編碼一個(gè)轉(zhuǎn)錄因子，通過生長素依賴形式調(diào)節(jié)細(xì)胞伸長，該基因突變能增加稈壁厚度，提高抗倒伏性[21]；控制水稻莖稈硅含量的硅轉(zhuǎn)運(yùn)蛋白基因突變會(huì)導(dǎo)致抗倒伏性、抗病性降低[22]。與水稻次生細(xì)胞壁相關(guān)的基因與抗倒伏性密切相關(guān)，如對香豆酸輔酶A連接酶是木質(zhì)素合成的關(guān)鍵酶，編碼該酶的基因表達(dá)量下調(diào)會(huì)導(dǎo)致植株變矮，莖稈強(qiáng)度變低[23]。

本研究通過全基因組關(guān)聯(lián)分析，在洛陽、吉林兩地均檢測到較多的蛋白激酶基因，特別是在谷子6號染色體，洛陽、吉林都檢測到同一個(gè)LRR受體類絲氨酸/蘇氨酸蛋白激酶基因(LOC101776261)，在6號、8號染色體還檢測到一些細(xì)胞壁關(guān)聯(lián)的類受體激酶基因(LOC105914550、LOC101786429、LOC105913474、LOC101766879)。目前報(bào)道的類受體蛋白激酶多數(shù)集中在植物抗病防御反應(yīng)，有關(guān)在抗倒伏性方面的報(bào)道極少[24]。近期，華中農(nóng)業(yè)大學(xué)王石平教授領(lǐng)導(dǎo)的團(tuán)隊(duì)發(fā)現(xiàn)水稻抗白葉枯病主效基因Xa4編碼一個(gè)細(xì)胞壁相關(guān)的激酶，該激酶通過促進(jìn)纖維素的合成增強(qiáng)了細(xì)胞壁的強(qiáng)度，不僅防御了白葉枯病菌的侵染，而且增加了莖稈機(jī)械強(qiáng)度，提高了水稻抗倒伏性[25]。據(jù)此推測，本研究關(guān)聯(lián)到的LRR受體類絲氨酸/蘇氨酸蛋白激酶基因、細(xì)胞壁相關(guān)的受體激酶基因可能與谷子抗倒伏性有密切關(guān)系。此外，一些雖然在洛陽、吉林兩地檢測到的位點(diǎn)不同但功能接近的基因，如細(xì)胞色素P450相關(guān)蛋白、UDP糖基轉(zhuǎn)移酶基因、驅(qū)動(dòng)蛋白基因、甲基轉(zhuǎn)移酶基因也可能與抗倒伏性有一定關(guān)系，這些基因仍需要進(jìn)一步進(jìn)行功能驗(yàn)證。