大豆種質耐堿性狀全基因組關聯分析

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

（東北農業大學，大豆生物學教育部重點實驗室，農業農村部東北大豆生物學與遺傳育種重點實驗室，哈爾濱 150030）

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

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

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

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

1 材料與方法

1.1 試驗材料

選用255份大豆種質構成的自然群體為試驗材料（見表1），試驗所用材料均由東北農業大學提供。

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

1.2 試驗設計

1.2.1 芽期試驗設計

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

1.2.2 芽期耐堿性狀鑒定

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

1.2.3 苗期試驗設計

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

1.2.4 苗期耐堿性狀鑒定

苗期耐堿性狀鑒定指標包括：莖粗（Stem diameter，SD）、葉綠素（Soil and plant analyzer develotrnent，SPAD）含量、主根長（Main root length，MSL）、地上部鮮重（Shoots fresh weight，SFW），相對堿害指標（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個標記用于進一步關聯分析。

1.2.6 全基因組關聯分析

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

1.2.7 數據處理與分析

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

2 結果與分析

2.1 大豆種質耐堿性狀表型和相關性分析

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

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

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

大豆種質耐堿性狀分布如圖2所示。相對吸脹率、相對葉綠素等芽期和苗期耐堿相關指標近似于正態分布，表明耐堿性狀由多基因控制，適合進行關聯分析。

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

2.2 大豆種質耐堿性狀全基因組關聯分析

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

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

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

表3 與大豆耐堿性狀相關聯位點Table 3 Loci associated with alkali tolerance traits in soybean

續表

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

2.3 大豆耐堿性候選基因預測

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

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

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

3 討 論

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

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

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

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

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

4 結 論

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

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

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