大豆種質(zhì)耐堿性狀全基因組關(guān)聯(lián)分析

滕衛(wèi)麗，史飛飛，林 峰，劉晨煦，高 鵬，趙慧艷，趙 雪，韓英鵬，李文濱

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

大豆作為我國五大作物之一，是重要的高蛋白糧飼兼用及油料作物。當(dāng)前大豆產(chǎn)量不能滿足人們?nèi)找嬖鲩L的生活需求，土壤鹽堿化是影響大豆產(chǎn)量重要因素之一。土壤中存在的鹽類影響植物正常生長發(fā)育，降低產(chǎn)量[1]，研究認(rèn)為由Na2CO3、NaHCO3等堿性鹽比NaCl、Na2SO4等中性鹽對作物造成的脅迫危害更嚴(yán)重[2]。我國是受土壤鹽漬化危害影響比較嚴(yán)重的國家，有逐年擴大趨勢[3-4]。篩選耐堿大豆種質(zhì)，找到與耐堿有關(guān)且穩(wěn)定表達的位點，對培育耐堿大豆新品種、有效利用鹽堿化土壤具有重要意義。

種子芽期和苗期是作物在鹽漬環(huán)境下生長的關(guān)鍵及敏感階段，避開敏感期，作物在各發(fā)育時期均可較好生長發(fā)育。大豆耐堿相關(guān)性狀是復(fù)雜的數(shù)量性狀，受多基因控制[5]。隨著現(xiàn)代分子生物技術(shù)發(fā)展，耐堿基因發(fā)掘和使用成為近年植物抗逆研究熱點。目前，許多學(xué)者已檢測到多個芽期和苗期耐堿相關(guān)位點和QTL。Zhang 等對257 份大豆品種在芽期作Na2CO3溶液處理，通過測定主根長度、幼苗生物量等芽期耐堿指標(biāo)評估品種耐受性，結(jié)合135個SSR標(biāo)記進行互作關(guān)聯(lián)分析，檢測到86個耐堿性相關(guān)QTL[6]。張新草利用混鹽堿溶液模擬鹽堿脅迫，對281份大豆品種進行芽期脅迫鑒定，調(diào)查發(fā)芽率、根長、幼苗長等8個芽期耐鹽堿指標(biāo)，經(jīng)全基因組關(guān)聯(lián)分析，關(guān)聯(lián)到268個SNP標(biāo)記，其中優(yōu)異增效等位變異有44 個，占SNP 標(biāo)記總數(shù)16.42%[7]。Tuyen等利用F6RIL群體和F2群體進行QTL分析，在17號染色體上檢測到耐堿性顯著的QTL，占F6和F2群體STR（耐鹽堿等級）總變異的50.20%和13.00%[8]。Tuyen 等進一步以剩余雜合系（RHL46）為材料，利用表型評估和SSR 標(biāo)記分析，在標(biāo)記Satt447和GM17-11.6之間3.33 cM處定位到耐堿QTL[9]。張金龍以NJRIZM6 和NJRILM6 兩個RIL群體為試驗材料，利用NaHCO3溶液進行堿脅迫處理，對STR、Na+濃度、K+濃度等7個苗期耐堿指標(biāo)進行QTL定位，在兩個群體中分別定位到22 個、19個與耐堿性狀相關(guān)的QTL位點[5]。

目前對玉米、水稻和大豆等作物耐鹽堿鑒定時期雖集中在芽期和苗期，但主要采用單一時期和單一成分的鹽或堿開展試驗。黑龍江省作為大豆主產(chǎn)地之一，鹽堿土類型主要為蘇打鹽堿土，且鹽堿土面積有逐年增加趨勢[3]，因此利用混合成分模擬實際堿脅迫環(huán)境，開展多時期大豆種質(zhì)堿脅迫鑒定研究，對實際生產(chǎn)意義重大。

本研究針對黑龍江省大慶市蘇打鹽堿地，以混合堿模擬其堿脅迫環(huán)境，評價255份大豆種質(zhì)芽期和苗期耐堿相關(guān)性狀，結(jié)合全基因組關(guān)聯(lián)分析，獲得關(guān)聯(lián)位點，挖掘候選基因，以期為分子標(biāo)記輔助選育大豆耐堿品種提供依據(jù)。

1 材料與方法

1.1 試驗材料

選用255份大豆種質(zhì)構(gòu)成的自然群體為試驗材料（見表1），試驗所用材料均由東北農(nóng)業(yè)大學(xué)提供。

表1 255份大豆供試材料Table 1 Total 255 soybean test materials

1.2 試驗設(shè)計

1.2.1 芽期試驗設(shè)計

針對黑龍江省蘇打鹽堿地，以混合堿模擬李憲偉研究中采用的大慶地區(qū)中度蘇打鹽堿土壤基本情況[10]，按其土壤中主要成分HCO3-與CO32-質(zhì)量比為8.37∶1，計算HCO3-與CO32-物質(zhì)的量比為10.56∶1，配置堿溶液，pH 為8.9～9.1；借鑒Tuyen等試驗方法[8]，采用濃度梯度方式在人工氣候箱25 ℃、相對濕度60%條件下開展試驗，隔天更換相應(yīng)25 mL溶液，試驗中每天調(diào)查發(fā)芽種子數(shù)，用于計算活力指數(shù)，種子發(fā)芽標(biāo)準(zhǔn)參考徐玲秀[11]。

1.2.2 芽期耐堿性狀鑒定

芽期耐堿性狀鑒定指標(biāo)包括：吸脹率（Imbibition rate-48h，IR-48 h）、胚根干重（Radicle dry weight，RDW）、子葉下部干重（Dry weight below cotyledon，DWBC）、活力指數(shù)（Vigour index，VI）[12]，相對堿害指標(biāo)（Alkali tolerance，AT）包括：相對吸脹率（Alkali tolerance-imbibition rate，AT-IR）、相對胚根干重（Alkali tolerance-radicle dry weight，ATRDW）、相對子葉下部干重（Alkali tolerance-dry weight below cotyledon，AT-DWBC）、相對活力指數(shù)（Alkali tolerance-vigour index，AT-VI）[13]。

1.2.3 苗期試驗設(shè)計

苗期同樣以混合堿模擬堿脅迫環(huán)境，借鑒Tuyen 等試驗方法[8]，將大小均勻的大豆種子播種于裝有1.15 kg土的盆中，待所有品種達到V2期進行定苗，對長勢一致的大豆植株進行堿脅迫處理。根據(jù)李憲偉描述的中度蘇打鹽堿土壤情況[10]，計算每次每盆需加入堿溶液或自來水318.5 mL，待最大濃度180 mmol·L-1處理2 d后進行表型測定。

1.2.4 苗期耐堿性狀鑒定

苗期耐堿性狀鑒定指標(biāo)包括：莖粗（Stem diameter，SD）、葉綠素（Soil and plant analyzer develotrnent，SPAD）含量、主根長（Main root length，MSL）、地上部鮮重（Shoots fresh weight，SFW），相對堿害指標(biāo)（Alkali tolerance，AT）包括：相對莖粗（Alkali tolerance-stem diameter，AT-SD）、相對葉綠素（Alkali tolerance-soil and plant analyzer develotrnent，AT-SPAD）、相對主根長（Alkali tolerancemain root length，AT-MSL）、相對地上部鮮重（Alkali tolerance-shoots fresh weight，AT-SFW）[5]。

1.2.5 基因型鑒定與分析

利用Illumina 基因組分析儀Ⅱ進行測序，完成基因型鑒定，將最小等位基因頻率MAF＜0.05 的SNP 通過BLINK 和FarmCPU 軟件過濾淘汰，最終得到1 355 928個標(biāo)記用于進一步關(guān)聯(lián)分析。

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

利用BLINK 和FarmCPU 兩種方法對大豆種質(zhì)芽期和苗期堿脅迫下相關(guān)性狀表型值進行關(guān)聯(lián)分析，為更好降低假陽性率，識別真正的標(biāo)記關(guān)聯(lián)性狀，利用Bonferroni 校正方法，將顯著P值設(shè)置為P＜3.6875e-8（P=0.05/n，n為總標(biāo)記數(shù)），繪制曼哈頓圖顯示SNP 與表型性狀之間相關(guān)性。

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

采用Microsoft Excel 2021 處理表型數(shù)據(jù)和作圖，利用Soybase 和Phytozome 數(shù)據(jù)庫獲得基因功能注釋和基因在組織中表達數(shù)據(jù)，利用CMplot R包和Prism 8 分別繪制曼哈頓圖和基因組織表達熱圖。

2 結(jié)果與分析

2.1 大豆種質(zhì)耐堿性狀表型和相關(guān)性分析

大豆自然群體耐堿相關(guān)性狀均值、峰度、偏度等見表2。與堿脅迫相比，無堿條件下芽期和苗期相關(guān)性狀表型值更高。芽期和苗期多數(shù)耐堿性狀偏度和峰度絕對值小于1，且在8 個性狀中，群體變異系數(shù)也不同，其中，葉綠素含量對照組群體變異系數(shù)最小，為7.06%，活力指數(shù)處理組變異系數(shù)最大，為123.70%。

表2 大豆種質(zhì)耐堿性狀表型分析Table 2 Phenotypic analysis of alkali tolerance traits in soybean germplasm

為消除不同大豆種質(zhì)間基因型背景差異，以芽期和苗期性狀相對值為相關(guān)耐堿指標(biāo)分析相關(guān)性（見圖1）。在芽期階段，相對胚根干重與相對下部干重、相對活力指數(shù)呈極顯著正相關(guān)，其中相對胚根干重與相對下部干重相關(guān)程度最高，相關(guān)系數(shù)為0.90；相對活力指數(shù)與相對下部干重呈極顯著正相關(guān)，相關(guān)系數(shù)為0.22。在苗期階段，相對葉綠素與相對主根長呈顯著負(fù)相關(guān)；相對莖粗與相對地上部鮮重呈顯著正相關(guān)，相關(guān)系數(shù)為0.16。在芽期測得的耐堿性狀與苗期測得的耐堿性狀之間未觀察到顯著相關(guān)性，表明大豆兩個發(fā)育階段耐堿相關(guān)性狀在很大程度上受獨立遺傳控制。

大豆種質(zhì)耐堿性狀分布如圖2所示。相對吸脹率、相對葉綠素等芽期和苗期耐堿相關(guān)指標(biāo)近似于正態(tài)分布，表明耐堿性狀由多基因控制，適合進行關(guān)聯(lián)分析。

圖2 大豆種質(zhì)耐堿性狀分布Fig.2 Distribution of alkali tolerance traits of soybean germplasm

2.2 大豆種質(zhì)耐堿性狀全基因組關(guān)聯(lián)分析

本研究利用Bonferroni 校正方法識別真正的標(biāo)記關(guān)聯(lián)性狀，因Bonferroni 校正方法較嚴(yán)格，借鑒Li等GWAS分析果實大小性狀時適當(dāng)調(diào)整閾值的方法，將閾值由3e-8降至1e-5，以消除Bonferroni方法校正過于嚴(yán)格導(dǎo)致的假陰性[14]。利用BLINK 和FarmCPU 兩種方法對大豆種質(zhì)芽期和苗期耐堿性狀進行全基因組關(guān)聯(lián)分析（見圖3 和4），兩種方法共檢測到64 個與大豆耐堿性狀關(guān)聯(lián)的SNP 位點（見表3）。

圖3 大豆芽期耐堿相關(guān)性狀全基因組關(guān)聯(lián)分析曼哈頓圖Fig.3 Manhattan plots of genome wide association analysis of alkali tolerance related traits in soybean at germination stage

圖4 大豆苗期耐堿相關(guān)性狀全基因組關(guān)聯(lián)分析曼哈頓圖Fig.4 Manhattan plots of genome wide association analysis of alkali tolerance related traits in soybean at seedling stage

表3 與大豆耐堿性狀相關(guān)聯(lián)位點Table 3 Loci associated with alkali tolerance traits in soybean

續(xù)表

這些位點主要分布在除12、13、16、20 號染色體外其他16 條染色體上，呈不均勻分布。其中，有2個SNP位點與芽期相對吸脹率相關(guān)聯(lián)，分布在2、14 號染色體上，表型貢獻率為0.36%～0.44%；9 個SNP 位點與相對胚根干重相關(guān)聯(lián)，分布在3 號和7 號染色體上，表型貢獻率為0.09%～9.62%，其中6 個SNP 位點貢獻率高于1%，位點7:8915107 貢獻率最高。BLINK 和FarmCPU 兩種方法在芽期下部干重共檢測出8個關(guān)聯(lián)SNP位點，其中，BLINK方法單獨檢測出1個顯著關(guān)聯(lián)位點，為9:37033917，貢獻率為2.41%，還有7 個關(guān)聯(lián)SNP位點由兩種方法同時檢測出，分布在3、4、8、9號染色體上，表型貢獻率為0.11%～8.16%，其中3:37525453 貢獻率最高，該位點與相對下部干重顯著關(guān)聯(lián)。兩種方法共檢測出9個與相對活力指數(shù)相關(guān)聯(lián)SNP位點，其中FarmCPU方法單獨檢測出8個SNP 位點，分布在1、6、7、8、17、18、19 號染色體上，還有1 個相同位點由BLINK 和FarmCPU方法同時檢測出，分布在6號染色體上。兩種方法共檢測出11 個與苗期相對葉綠素含量相關(guān)聯(lián)的SNP位點，其中有2個SNP位點由FarmCPU方法單獨檢測出，均位于5號染色體上，還有9個關(guān)聯(lián)位點由BLINK 和FarmCPU 方法同時檢測出，分布在2、5、15、18 號染色體上，貢獻率為0.10%～1.70%。兩種方法共檢測出16個與苗期相對莖粗相關(guān)聯(lián)SNP位點，其中有8個位點由FarmCPU方法單獨檢測出，分布在2、3、5、8、9、10、18號染色體上，還有8個位點由兩種方法同時檢測出，分布在4 號和8 號染色體上，貢獻率為0.15%～1.35%。兩種方法同時檢測出3 個與相對主根長相關(guān)聯(lián)的SNP 位點，分布在6、9、14 號染色體上，貢獻率為0.24%～1.65%。兩種方法共檢測出6個與苗期相對地上鮮重相關(guān)聯(lián)SNP 位點，其中，F(xiàn)armCPU 方法單獨檢測出1個與相對地上鮮重關(guān)聯(lián)位點，位于11號染色體上，還有5個關(guān)聯(lián)SNP位點由兩種方法同時檢測出，位于11號和17號染色體上，貢獻率為0.47%～1.91%。

2.3 大豆耐堿性候選基因預(yù)測

將檢測到的64 個SNP 位點在上下游100 kb 區(qū)間內(nèi)挖掘基因，共挖掘到26 個有功能注釋基因。根據(jù)Soybase和Phytozome數(shù)據(jù)庫中基因功能注釋數(shù)據(jù)（見表4），發(fā)現(xiàn)這些基因主要在編輯蛋白激酶超家族蛋白、可逆糖基化多肽方面發(fā)揮作用。結(jié)合Phytozome 數(shù)據(jù)庫中已收錄基因表達量數(shù)據(jù)，發(fā)現(xiàn)上述基因主要在花、葉、根等9 個組織中表達（見圖5）。

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

分析上述基因功能，發(fā)現(xiàn)有7個基因與耐逆性相關(guān)。與芽期相對下部干重相關(guān)聯(lián)基因Glyma.03G160100、Glyma.03G160200 均具有編碼脂肪酸ω-羥化酶的功能，脂肪酸ω-羥化酶參與角質(zhì)合成，角質(zhì)與植物抵御外界脅迫有關(guān)[15]，其中Glyma.03G160200 在根瘤、根、根毛中表達量較高。與相對下部干重相關(guān)聯(lián)基因Glyma.09G149900在根中表達水平最高，該基因與bHLH 有關(guān)，bHLH 轉(zhuǎn)錄因子在抗非生物脅迫方面發(fā)揮重要作用[16-17]。與相對活力指數(shù)相關(guān)聯(lián)基因Glyma.19G067400 在花、葉等9個組織中均有較高水平表達，具有編碼Arf-GTPase 激活蛋白的功能，ArfGTPase 激活蛋白與植物應(yīng)答非生物脅迫有關(guān)[18]。與芽期相對活力指數(shù)關(guān)聯(lián)基因Glyma.18G272200、Glyma.01G057100 分別在（根瘤、根、根毛）和（豆莢、根）中表達量較高，兩個基因分別具有編碼可逆糖基化多肽2 和UDP-葡糖6-脫氫酶家族蛋白的功能，可逆糖基化多肽2 可能直接參與植物防御反應(yīng)[19]，UDP-葡糖6-脫氫酶家族蛋白與細胞壁合成有關(guān)[20]，細胞壁能提高植物耐鹽性[21]。與苗期地上部鮮重相關(guān)聯(lián)基因Glyma.11G238400 編碼硫酸鹽轉(zhuǎn)運蛋白，在葉和莖尖分生組織中表達量較高，硫酸鹽轉(zhuǎn)運蛋白可能參與相應(yīng)非生物逆境[22]。

3 討 論

芽期和苗期是大豆在鹽漬環(huán)境下生長關(guān)鍵和敏感階段，高濃度堿延遲種子萌發(fā)和出苗，導(dǎo)致植株莖稈變細，干重降低，老葉逐漸枯黃發(fā)干甚至脫落[11]。本研究在大豆芽期和苗期進行混堿鑒定，對耐堿相關(guān)性狀進行GWAS 分析，以期發(fā)掘重要標(biāo)記和遺傳位點，挖掘候選基因，為分子標(biāo)記輔助選育大豆耐堿品種提供依據(jù)。

目前，BLINK 和FarmCPU 等方法在挖掘全基因組遺傳變異方面發(fā)揮重要作用。本研究利用這兩種方法檢測到64 個與大豆耐堿性狀相關(guān)聯(lián)SNP 位點，其中7 個位點與已知QTL 位點物理距離小于5 Mb[23]。與相對胚根值和相對下部干重值關(guān)聯(lián)位點3:37525453，與相對胚根值關(guān)聯(lián)位點3:37513903、3:37490361，這3 個關(guān)聯(lián)位點與耐鹽相關(guān)的ss715585865 分別相距377.48、389.03、412.57 kb[24]，本研究發(fā)現(xiàn)位點3:37525453與相對下部干重顯著關(guān)聯(lián)，貢獻率較高，為8.16%，因此該位點可靠性較高，可用于進一步基因挖掘。與相對下部干重顯著關(guān)聯(lián)位點9:37033917 與苗期耐堿性狀相關(guān)的qCa09.1 相距240.70 kb[5]，相對活力指數(shù)關(guān)聯(lián)位點17:1589196 與鹽脅迫相關(guān)BARC-031827-07220 相距196.70 kb[25]，相對主根長關(guān)聯(lián)位點14:9834716 在堿脅迫性狀相關(guān)的9 806 588-9 879 926 區(qū)間內(nèi)[5]，相對莖粗關(guān)聯(lián)位點9:43068852與堿脅迫下苗期性狀相關(guān)qCa09.3相距887.82 kb[5]。由此可知，3:37513903、3:37490361、9:37033917等7個位點相對穩(wěn)定，可靠性較高，具有重要參考價值。

在GWAS 分析過程中發(fā)現(xiàn)，與FarmCPU 方法相比，BLINK方法的運算速度更快，在檢測功效方面，F(xiàn)armCPU方法更具明顯優(yōu)勢，可挖掘更多與耐堿相關(guān)的位點，如19:19200072、18:55518213 等，與Liu 等對相關(guān)性狀進行GWAS 分析時發(fā)現(xiàn)Farm-CPU方法的檢測功效較高的觀點相似[26]。

本研究通過功能分析，發(fā)現(xiàn)7個基因與植物耐逆性相關(guān)。其中與相對下部干重相關(guān)聯(lián)的Glyma.03G160100、Glyma.03G160200 具有相同的擬南芥同源基因，功能注釋為脂肪酸ω-羥化酶。脂肪酸ω-羥化酶，屬于細胞色素P450超家族，與角質(zhì)合成有關(guān)，角質(zhì)在保護植物免受脅迫具有重要作用[15]。Mao 等研究發(fā)現(xiàn)細胞色素P450 單加氧酶基因突變體過表達可增強轉(zhuǎn)基因擬南芥鹽脅迫耐受性[27]。Xia 等以蒺藜苜蓿為材料再次證實MtP450 樞紐基因提高植物鹽脅迫下耐受性[28]。Glyma.03G160200在根瘤、根、根毛中表達量較高，是由重要位點3:37525453 挖掘出的基因，可將其作為與耐堿相關(guān)的候選基因。與相對下部干重顯著關(guān)聯(lián)的基因Glyma.09G149900，功能注釋為堿性螺旋-環(huán)-螺旋（bHLH）DNA 結(jié)合超家族蛋白。bHLH 轉(zhuǎn)錄因子在調(diào)控植物的耐鹽性方面，分為正負(fù)調(diào)控兩種，如Zhou 等研究發(fā)現(xiàn)OrbHLH2 基因超表達會提高植株在高鹽脅迫下的種子萌發(fā)率[16]，Verma 等研究發(fā)現(xiàn)，高鹽脅迫下bHLH 轉(zhuǎn)錄因子編碼基因AtMYC2負(fù)調(diào)控脯氨酸的合成，進而降低轉(zhuǎn)基因植株耐鹽性[17]。基因Glyma.09G149900 在根中表達水平很高，是由重要位點9:37033917 挖掘出的基因，可將其作為與耐堿相關(guān)的候選基因之一。與相對活力指數(shù)相關(guān)聯(lián)的基因Glyma.19G067400，功能注釋為編碼Arf GTPase 激活蛋白，劉芳等以美洲商陸為材料，發(fā)現(xiàn)ArfGTPase 激活蛋白基因PaAGAP 參與植物逆境反應(yīng)[18]。根據(jù)Glyma.19G067400 在花、葉、根瘤、豆莢等9個組織中均有較高水平表達，且是由顯著關(guān)聯(lián)位點19:19200072 挖掘出的基因，可將此基因也作為耐堿候選基因。與相對活力指數(shù)相關(guān)聯(lián)基因Glyma.01G057100 和Glyma.18G272200、與苗期地上部鮮重相關(guān)聯(lián)基因Glyma.11G238400的功能注釋參與植物耐逆性[19-22]，但尚未發(fā)現(xiàn)挖掘這3個基因的位點與其他耐鹽堿位點重疊或相近。在功能分析方面，上述7 個基因與擬南芥、蒺藜苜蓿、美洲商陸等作物的非生物脅迫相關(guān)，但目前尚未在大豆堿脅迫方面發(fā)現(xiàn)這些基因的相關(guān)報道，可能是由于所用材料類型不同，遺傳背景存在差異導(dǎo)致。為進一步探究這些基因在大豆耐堿方面的作用，后期將驗證候選基因功能，為大豆在耐堿方面分子育種提供有效依據(jù)。

綜上，根據(jù)挖掘基因位點情況，查詢基因功能注釋及分析基因表達量等，可將Glyma.03G160200、Glyma.09G149900 和Glyma.19G067400 等3 個 基 因作為大豆堿脅迫性狀相關(guān)重要候選基因。

4 結(jié) 論

a.BLINK和FarmCPU兩種方法共同檢測出64個與大豆芽期和苗期耐堿性狀關(guān)聯(lián)SNP 位點，其中3:37513903、3:37490361、9:37033917 等7 個 位 點相對穩(wěn)定，可靠性較高。

b. 在上下游區(qū)間100 kb 進行基因挖掘，根據(jù)Soybase和Phytozome篩選到26個有功能注釋基因。

c. 根據(jù)挖掘基因位點情況、功能注釋和基因表達量，推測3 個基因與耐堿性相關(guān)，分別為Glyma.03G160200、Glyma.09G149900 和Glyma.19G067400。