李文濱，吳 航，劉 冀，張翔超，郭志文，鄭立娜，趙 雪，韓英鵬，滕衛(wèi)麗

（東北農(nóng)業(yè)大學大豆研究所，大豆生物學教育部重點實驗室，農(nóng)業(yè)農(nóng)村部東北大豆生物學與遺傳育種重點實驗室，哈爾濱150030）

大豆是重要經(jīng)濟和糧食作物，其營養(yǎng)價值高，用途廣泛。倒伏是限制大豆高產(chǎn)、穩(wěn)產(chǎn)和優(yōu)質(zhì)重要因素之一，大豆品種抗倒伏性是獲得高產(chǎn)穩(wěn)產(chǎn)重要保證[1-2]。倒伏是受遺傳、栽培、氣候和生態(tài)等多種因素共同制約的復雜性狀。以往研究表明作物倒伏與植株地上部性狀（如株高、莖粗、莖稈強度）[3-4]和地下部性狀（如根重、根長、根體積等）有關(guān)[5-6]。在關(guān)于大豆倒伏及其相關(guān)性狀研究中，周蓉和鐘開珍等認為大豆倒伏性與莖稈性狀（株高、主莖節(jié)數(shù)、節(jié)間長度、分枝數(shù)等）和產(chǎn)量性狀（單株莢數(shù)、單株粒數(shù)、單株產(chǎn)量和百粒重）顯著相關(guān)[7-8]。由于大豆倒伏性與該性狀密切相關(guān)，且目前對于解析大豆抗倒性遺傳基礎(chǔ)研究報道相對較少，利用分子標記技術(shù)探索大豆倒伏性遺傳特性成為大豆科研工作者研究熱點。通過QTL定位方法已發(fā)掘出許多與大豆倒伏性相關(guān)位點，但利用QTL定位通常需構(gòu)建作圖群體，研究時間較長、定位精度低[9-10]，而基于自然群體的關(guān)聯(lián)分析具有研究時間短，檢測效率高、定位精度高等優(yōu)勢[11]。

諸多學者已對大豆主要農(nóng)藝性狀開展GWAS（Genome-wide association analysis）分析，例如，Jing等通過GWAS檢測到14個與株高相關(guān)SNP位點，預測6個可能與株高相關(guān)新基因[12]。Chang等對368份大豆種質(zhì)資源作GWAS分析，檢測到45個與株高相關(guān)位點和43個與主莖節(jié)數(shù)相關(guān)位點[13]。張軍等對190份大豆種質(zhì)多個農(nóng)藝性狀作GWAS分析，累計檢測到136個位點與目標性狀顯著關(guān)聯(lián)，其中14個位點與株高關(guān)聯(lián)，13個位點與倒伏性關(guān)聯(lián)[14]。張友誼利用224份大豆微核心種質(zhì)和1 749個SNP標記作GWAS分析，共檢測到36個SNPs位點，其中與株高關(guān)聯(lián)的13個位點，主莖節(jié)數(shù)4個，莖粗9個[15]。王自力對由573份種質(zhì)建立的江淮地區(qū)大豆育種群體作GWAS分析，兩年重復檢測到與倒伏性關(guān)聯(lián)SNP標記2個，株高130個，主莖節(jié)數(shù)92個[16]。關(guān)聯(lián)分析不僅可找到與大豆倒伏性狀有關(guān)遺傳位點，還可預測與目標性狀相關(guān)新基因。因此利用大豆自然群體作關(guān)聯(lián)分析研究抗倒伏性遺傳特性，對實現(xiàn)大豆高產(chǎn)穩(wěn)產(chǎn)優(yōu)質(zhì)等目標具有重要意義。

本研究利用150個大豆種質(zhì)資源組成的自然群體開展多年多點試驗，通過RTM-GWAS方法結(jié)合SNP位點對株高、主莖節(jié)數(shù)等性狀作全基因組關(guān)聯(lián)分析，挖掘候選基因，為分子標記輔助大豆抗倒伏育種提供依據(jù)。

1 材料與方法

1.1 試驗材料

供試材料選用150份不同大豆種質(zhì)資源構(gòu)成的自然群體。其中國外品種（系）3份、北京品種（系）8份、黑龍江品種（系）93份、吉林品種（系）20份、遼寧品種（系）3份、內(nèi)蒙古品種（系）15份、新疆品種（系）2份、河北品種（系）4份、山東品種（系）1份和甘肅品種（系）1份。

1.2 試驗設(shè)計

試驗于2018～2019年在向陽（A）、呼蘭（B）、阿城（C）3個試驗地點開展，采用隨機區(qū)組設(shè)計，3次重復，行長2 m、壟寬0.65 m，株距0.06 m。成熟期收獲并室內(nèi)考種，田間管理同大田。

1.3 方法

1.3.1 性狀調(diào)查

根據(jù)周蓉等方法于成熟期調(diào)查倒伏性[17]，田間倒伏分級標準：1級為小區(qū)植株全部直立；2級為輕度倒伏（植株傾斜＜15°）；3級為重倒伏（植株傾斜15°～44°）；4級為嚴重倒伏（植株傾斜＞45°）。從自然群體每個資源中選擇具有代表性連續(xù)5株測定株高、主莖節(jié)數(shù)、節(jié)間長度、莖粗、重心高度、莖稈強度、單株莢數(shù)、單株粒數(shù)、單株粒重、百粒重、莖稈抗倒指數(shù)等性狀，3次重復，以平均值為觀測值。其中重心高度和莖稈抗倒指數(shù)測定均參照陳喜鳳等方法[2]，其他室內(nèi)考種項目及標準均參照邱麗娟等方法[18]。

1.3.2 基因型檢測與分析

采用特異性位點擴增片段測序（SLAF-seq）技術(shù)，利用Illumina基因組分析儀Ⅱ獲得大豆染色體均勻分布高質(zhì)量單核苷酸多態(tài)性（SNPs）。提取150份大豆種質(zhì)資源基因型信息，使用Seqtk（https://github.com/lh3/seqtk）過濾，最終得到最小等位基因頻率（Minor allele frequency，MAF）＞0.05的23 309個標記用于進一步關(guān)聯(lián)分析。標記覆蓋所有染色體區(qū)段，標記之間平均距離為28 kb。

1.3.3 全基因組關(guān)聯(lián)分析

利用RTM-GWAS方法采用兩階段策略對大豆倒伏相關(guān)性狀表型值作關(guān)聯(lián)分析。

第一階段模型公式如下：

式中，yi表示個體i表型觀測值；μ表示總體平均數(shù)；J表示用于群體結(jié)構(gòu)矯正的特征向量個數(shù)，wij表示遺傳相似系數(shù)矩陣第j個特征向量在個體i上的系數(shù)，αj表示第j個特征向量效應；L表示測驗標記位點等位基因數(shù)目，xil表示測驗標記位點第l個等位基因?qū)τ趥€體i基因型指示變量，取值0或1；βl表示第l個等位基因效應；εi表示假定服從正態(tài)分布殘差效應。

第二階段模型公式如下：

式中，K表示總QTL數(shù)目；Lk表示第k個位點等位基因數(shù)目，xikl表示第k個位點第l個等位基因在個體i上基因型指示變量，取值0或1；βkl表示第k個位點第l個等位基因效應。其他符號含義與模型（1）相同。

1.4 數(shù)據(jù)處理與分析

采用Microsoft Office Excel 11.1.0.10314、SPSS 3.0軟件分析表型數(shù)據(jù)和作圖?；赑hytozome數(shù)據(jù)庫（https://phytozome.org/）和Soybase數(shù)據(jù)庫（https://www.Soybase.Org/）中基因功能注釋預測與大豆抗倒伏相關(guān)候選基因。在eFP數(shù)據(jù)庫（https://bar.Utoronto.ca）查詢候選基因表達量數(shù)據(jù)，利用R軟件CMplot程序包以及pheatmap程序包繪制曼哈頓圖和候選基因組織表達熱圖。

2 結(jié)果與分析

2.1 大豆倒伏相關(guān)性狀表型數(shù)據(jù)分析

分析多環(huán)境下大豆自然群體各農(nóng)藝性狀相關(guān)性（見表1）。結(jié)果表明，在不同環(huán)境下倒伏級別與各性狀相關(guān)性表現(xiàn)不同，除2019年呼蘭點外倒伏級別與株高均表現(xiàn)顯著相關(guān)，相關(guān)系數(shù)為0.243～0.527；兩年阿城點倒伏級別與主莖節(jié)數(shù)均表現(xiàn)顯著相關(guān)，相關(guān)系數(shù)為0.171～0.221；2018年呼蘭點倒伏級別與節(jié)間長度表現(xiàn)顯著相關(guān)，相關(guān)系數(shù)為0.275；2019年向陽點倒伏級別與莖粗均表現(xiàn)顯著相關(guān)，相關(guān)系數(shù)為0.290；除2018年阿城點外倒伏級別與重心高度均表現(xiàn)顯著相關(guān)，相關(guān)系數(shù)為-0.209～0.345；兩年呼蘭點倒伏級別與抗倒指數(shù)均表現(xiàn)顯著相關(guān)，相關(guān)系數(shù)為-0.254～0.178，表明倒伏受植株株高、主莖節(jié)數(shù)、節(jié)間長度、莖粗、重心高度以及抗倒指數(shù)等性狀影響顯著，該性狀均為倒伏顯著相關(guān)性狀。

分析多環(huán)境下大豆自然群體6個倒伏顯著相關(guān)性狀表型（見表2）。結(jié)果表明，不同環(huán)境下，大豆自然群體倒伏相關(guān)性狀表現(xiàn)均不同。各倒伏相關(guān)性狀在6個環(huán)境下偏度和峰度絕對值多數(shù)＜1，且具有廣泛變異，其中株高變異系數(shù)為0.13～0.29；主莖節(jié)數(shù)變異系數(shù)為0.11～0.23；節(jié)間長度變異系數(shù)為0.15～0.22；莖粗變異系數(shù)為0.14～0.23；重心高度變異系數(shù)為0.14～0.23；抗倒指數(shù)變異系數(shù)為0.18～0.44。

表1 大豆自然群體倒伏級別與不同農(nóng)藝性狀相關(guān)系數(shù)Table 1 Correlation coefficients between lodging grade and agronomic characters in soybean natural population

表2 大豆自然群體倒伏相關(guān)性狀表型分析Table 2 Phenotypic analysis of lodging traits in soybean natural population

大豆自然群體6個倒伏顯著相關(guān)性狀在6個環(huán)境下頻數(shù)分布情況見圖1。結(jié)果表明，該群體各倒伏相關(guān)性狀呈連續(xù)性變異，在群體中均呈正態(tài)分布，表明此性狀為多基因控制數(shù)量性狀，適合于關(guān)聯(lián)分析。

2.2 大豆倒伏相關(guān)性狀全基因組關(guān)聯(lián)分析

利用23 309個單核苷酸多態(tài)性位點基因分型信息，通過RTM-GWAS程序?qū)ψ匀蝗后w6個倒伏相關(guān)性狀作多環(huán)境關(guān)聯(lián)分析（見圖2）。

圖1 不同環(huán)境下大豆自然群體倒伏相關(guān)性狀分布Fig.1 Distribution of lodging traits of soybean natural population in different environments

結(jié)果表明，在P＜0.01顯著性水平上，共檢測到99個顯著關(guān)聯(lián)SNP位點，且這些顯著關(guān)聯(lián)SNP位點在大豆染色體上呈不均勻分布（見表3）。檢測到16個與株高顯著相關(guān)SNP位點，分布于1號、2號、3號、4號、6號、7號、8號、9號、13號、15號和18號染色體上，其中有9個效應顯著位點總表型變異解釋率為13.98%，14個與環(huán)境互作效應顯著位點總表型變異解釋率為19.47%，其中表型變異解釋率高于1%位點5個。檢測到19個與主莖節(jié)數(shù)顯著相關(guān)SNP位點，分布于3號、4號、7號、9號、12號、13號、15號、18號、19號和20號染色體上，其中有10個效應顯著位點總表型變異解釋率為14.54%，16個與環(huán)境互作效應顯著位點總表型變異解釋率為34.45%，表型變異解釋率高于1%位點7個。檢測到18個與節(jié)間長度顯著相關(guān)SNP位點，分布于2號、3號、4號、5號、10號、14號、15號、17號、18號、19號和20號染色體上，其中有7個效應顯著位點總表型變異解釋率為5.92%，18個與環(huán)境互作效應顯著位點總表型變異解釋率為37.81%，表型變異解釋率高于1%位點1個。檢測到20個與莖粗顯著相關(guān)SNP位點，分布于2號、3號、4號、6號、8號、10號、11號、15號、16號和19號染色體上，其中有10個效應顯著位點總表型變異解釋率為7.73%，18個與環(huán)境互作效應顯著位點總表型變異解釋率為19.45%，表型變異解釋率高于1%位點3個。檢測到17個與重心高度顯著相關(guān)SNP位點，分布于1號、4號、7號、8號、9號、10號、14號、15號、17號和20號染色體上，其中有8個效應顯著位點總表型變異解釋率為10.08%，14個與環(huán)境互作效應顯著位點總表型變異解釋率為24.00%，表型變異解釋率高于1%位點4個。檢測到9個與抗倒指數(shù)顯著相關(guān)SNP位 點，分 布 于2號、3號、4號、7號、9號、11號、12號和19號染色體上，其中有4個效應顯著位點總表型變異解釋率為2.67%，9個與環(huán)境互作效應顯著位點總變異解釋率為18.47%，表型變異解釋率高于1%位點1個。

2.3 候選基因預測

將以上檢測到的99個顯著SNP位點作候選基因挖掘分析，在Phytozome和Soybase數(shù)據(jù)庫中共檢索到22個有功能注釋基因（見表4）。根據(jù)注釋信息，該基因在大豆中主要與轉(zhuǎn)運蛋白、各種激酶、氧化還原酶等相關(guān)。

研究表明倒伏相關(guān)基因主要在作物莖、葉、根及頂芽等部位表達豐富[19-20]。本研究檢索Bar.utoronto.ca數(shù)據(jù)庫中已收錄大豆基因表達量數(shù)據(jù)，發(fā)現(xiàn)上述22個候選基因中20個基因在大豆10個組織中表達（見圖3），其中與株高關(guān)聯(lián)Glyma.15G124700基因在葉中高水平表達；與節(jié)間長度關(guān)聯(lián)Glyma.04G244100基因在根瘤中較高水平表達，Glyma.10G144600基因在莖尖、葉、根、根尖及根毛中高水平表達，Glyma.18G149700基因在莖尖和根尖中高水平表達，Glyma.18G299500基因在莖尖、根及根尖中高水平表達；與莖粗關(guān)聯(lián)Glyma.18G011400基因在根尖中表達水平相對較高；與抗倒指數(shù)關(guān)聯(lián)Glyma.09G199600基因在根瘤中高水平表達，Glyma.11G078600基因在大豆10個組織中均低水平表達。根據(jù)基因表達量和基因功能注釋推測，3個基因Glyma.04G244100、Glyma.18G149700、Glyma.18G011400為大豆倒伏性狀相關(guān)重要候選基因。

圖2 大豆倒伏相關(guān)性狀全基因組關(guān)聯(lián)分析曼哈頓圖Fig.2 Manhattan plots of genome wide association analysis of lodging-related traits in soybean

表3 與大豆倒伏性狀顯著相關(guān)位點Table 3 Loci significantly associated with soybean lodging traits(P＜0.01)

續(xù)表3

續(xù)表3

表4 大豆倒伏性狀相關(guān)聯(lián)候選基因信息Table 4 Information on candidate genes associated with lodging traits in soybean

圖3 候選基因組織表達熱圖Fig.3 Heatmap profiles of the candidate genes in tissues

3 討論

大豆具有極高營養(yǎng)和經(jīng)濟價值，從大豆種質(zhì)資源中挖掘與大豆倒伏相關(guān)功能標記，通過標記輔助選擇育種方法加快培育抗倒性強大豆品種是一種經(jīng)濟又有效的手段。近年來，隨大豆基因序列公布和基因組學發(fā)展[21-22]，促進大豆重要數(shù)量性狀遺傳解析進程[23]。RTM-GWAS是由He等提出的限制性兩階段多位點全基因組關(guān)聯(lián)分析方法[24]，與以往方法相比更適于自然群體（包括種質(zhì)資源群體），可降低假陽性，使檢測更為精準[25]。本研究利用RTM-GWAS對6個環(huán)境下大豆自然群體抗倒伏性狀作遺傳解析，以期發(fā)掘與大豆倒伏性相關(guān)重要標記和遺傳位點，預測候選基因，為通過分子育種技術(shù)改良大豆抗倒性狀奠定基礎(chǔ)。

本研究共檢測到99個與大豆倒伏相關(guān)性狀顯著關(guān)聯(lián)SNP位點，均是與莖稈性狀顯著關(guān)聯(lián)位點；共檢索出22個有功能注釋基因，其中2個基因與前人研究結(jié)果一致，如與重心高度相關(guān)聯(lián)Glyma.15G223700在Kim、張彥威等研究中被定位到[26-27]；與抗倒指數(shù)相關(guān)聯(lián)Glyma.11G078600在Espinosa等研究中被定位到[28]。其中3個基因在不同植物中出現(xiàn)相關(guān)報道，與節(jié)間長度相關(guān)聯(lián)UDP糖基轉(zhuǎn)移酶超家族蛋白Glyma.18G149700，糖基轉(zhuǎn)移酶基因家族主要功能是參與調(diào)節(jié)植物激素平衡、內(nèi)外源物質(zhì)解毒、防御反應和次生代謝產(chǎn)物修飾[29]，細胞壁合成主要由不同類型糖基轉(zhuǎn)移酶完成，糖基轉(zhuǎn)移酶影響細胞壁組成成分，進而改變細胞壁理化結(jié)構(gòu)與組成結(jié)構(gòu)[30]，以應對各種脅迫，根據(jù)Glyma.18G149700在莖尖和根尖中較高表達量，可將其作為與倒伏相關(guān)首要候選基因。如與節(jié)間長度相關(guān)細胞壁關(guān)聯(lián)蛋白激酶2Glyma.04G244100，細胞壁關(guān)聯(lián)蛋白激酶（WAK）具有信號轉(zhuǎn)導、細胞擴大與伸長、生殖生長、防御病原菌和機械傷害、金屬離子脅迫等生物學功能，在植物生長發(fā)育過程中發(fā)揮重要作用，Kanneganti和Wagner等研究推測WAK家族基因可能通過參與植物細胞壁元件相互作用調(diào)控細胞生長從而改變植株生長[31-32]，且Glyma.04G244100在根瘤中具有較高表達量，可作為與倒伏相關(guān)重要候選基因之一。與莖粗相關(guān)聯(lián)G蛋白偶聯(lián)受體Glyma.18G011400，G蛋白偶聯(lián)受體（GPCRs）主要負責生物細胞信號傳導[33]，參與植物激素和抗逆反應，在植物脅迫響應中發(fā)揮重要作用[34]。Dymock等發(fā)現(xiàn)擬南芥GCR1基因失活將導致植株根、莖對細胞分裂素不敏感，根、子葉和葉片伸展受到抑制[35]。Chakraborty等研究表明擬南芥GCR1突變體植株中脅迫相關(guān)基因表達量較野生型差異明顯，GCR1影響脅迫基因表達[36]。由此可知，G蛋白偶聯(lián)受體基因通過調(diào)控植物細胞對多種激素信號和外界脅迫響應過程抵抗逆境[34]；Glyma.18G011400在根尖中具有相對較高表達量，也可作為與倒伏相關(guān)重要候選基因之一。上述基因可在逆境中發(fā)揮作用，通過參與調(diào)控植物體細胞建成、信號傳導、生長發(fā)育等多種生理活動過程抵抗逆境脅迫，當植株發(fā)生倒伏時，基因具體調(diào)控機制尚不明確。

此外，本研究中還有4個分別在葉、根瘤、根及根毛、根尖中高表達基因（Glyma.15G124700、Glyma.09G199600、Glyma.10G144600、Glyma.18G299500）首次被定位到，但目前在這些基因相關(guān)研究中尚未發(fā)現(xiàn)其與倒伏之間聯(lián)系，還需進一步研究。

綜合基因表達量分析和基因功能注釋，可將Glyma.18G149700、Glyma.04G244100、Glyma.18G011400等3個基因依次作為研究與大豆倒伏性狀相關(guān)重要候選基因。

4 結(jié)論

a.在大豆自然群體中共檢測到99個顯著關(guān)聯(lián)SNP位點，在大豆染色體上呈不均勻分布，其中48個主效顯著位點，89個與環(huán)境互作顯著位點；其中21個位點表型變異解釋率高于1%。

b.檢測到與株高顯著關(guān)聯(lián)SNP位點16個，與主莖節(jié)數(shù)顯著關(guān)聯(lián)SNP位點19個，與節(jié)間長度顯著關(guān)聯(lián)SNP位點18個，與莖粗顯著關(guān)聯(lián)SNP位點20個，與重心高度顯著關(guān)聯(lián)SNP位點17個，與抗倒指數(shù)顯著關(guān)聯(lián)SNP位點9個。

c.共檢索出22個有功能注釋基因，通過20個候選基因在大豆10個組織中均表達，根據(jù)基因表達量分析和基因功能注釋，推測3個基因(Glyma.18G149700、Glyma.04G244100、Glyma.18G011400）可作為大豆倒伏性狀相關(guān)候選基因。