大豆抗旱遺傳位點及候選基因發掘

田蕊， 張華， 黃玫紅， 邵振啟， 李喜煥， 張彩英

（河北農業大學農學院，教育部華北作物種質資源研究與利用重點實驗室，河北 保定 071001）

大豆是重要的糧食和油料作物，是植物蛋白和食用油的主要來源[1]。干旱嚴重影響大豆生長發育，干旱條件下大豆植株生長緩慢，葉片萎蔫，生物量降低，以至減產甚至絕收。葉片是大豆植株重要光合器官，其葉綠素含量與光合速率及產量密切相關，易受干旱影響，是反映品種抗旱性的重要指標[2]。已有研究對作物干旱條件下的葉綠素含量進行定位分析，沈波等[3]利用水稻重組自交系（recombinant inbred lines， RIL）群體，在水分脅迫和正常條件下對葉綠素含量進行連鎖分析，定位到13 個控制葉綠素含量的數量性狀基因座（quantitative trait locus，QTL），其中水分脅迫條件下6 個，解釋表型變異4.74%～19.86%；胡頌平等[4]以水稻F9重組自交系（珍汕97B×IRAT109）為材料，在正常與水分脅迫下對葉綠素含量進行定位，檢測到13 個葉綠素含量QTLs，其中正常條件下檢測到7個，聯合貢獻率為56.19%；干旱處理下檢測到6 個，聯合貢獻率為47.39%。Kumar 等[5]利用C306×HUW206 構建的小麥RIL 群體在干旱條件下定位到1 個控制葉綠素含量QTL（QChl.ksu-3B），該QTL 解釋表型變異的14.2%。Peleg 等[6]以硬粒小麥和野生二粒小麥構建的RIL 群體為材料，在干旱及正常水分條件下，對葉綠素含量進行QTL定位，獲得8個葉綠素含量QTLs，解釋表型變異的2%～10%。

崔世友等[7]以‘波高’ב南農94-156’構建RIL群體為材料，在4 個不同生育期對大豆葉綠素含量進行連鎖分析，定位到與葉綠素含量相關的10 個位點，分別位于1、7、12、13、18 和19 號染色體，解釋表型變異的6.9%～23.4%。Fang 等[8]利用809 份大豆自然群體重測序，對葉綠素含量進行關聯分析，獲得9 個大豆葉綠素含量相關位點，分別位于11、13、18、20 號染色體。Kato 等[9]利用‘Minsoy’בT225H’構建的F2:3家系在14 及17 號染色體定位到2 個葉綠素相關位點。由此可見，相較其他農作物，目前大豆多就正常水分條件下葉綠素含量進行研究。因此，本研究以199 份大豆品種資源構建的自然群體為材料，在干旱和正常澆水條件下，測定其葉綠素含量，并計算抗旱指數，結合20×重測序基因型數據進行關聯分析，發掘基于葉綠素含量的抗旱遺傳位點和相關基因，為大豆品種抗旱分子遺傳改良提供依據。

1 材料與方法

1.1 供試材料

以本課題組前期構建的大豆自然群體為材料，包括152個育成品種和47個地方品種[10]。

1.2 大豆種植方法

將蛭石和土以體積比1∶1 進行混合，裝入花盆（直徑35 cm，高40 cm，花盆下鋪塑料布），每盆裝土20 kg，施復合肥5 g，澆水6 L；當基質含水量為20%時播種，每盆播種20 粒；待2 片真葉完全展開時定苗，每盆保留10株。試驗分別于2021年6月（E1）、7月（E2）和8月（E3）在河北農業大學作物育種中心試驗基地旱棚內進行，設置對照組和干旱脅迫處理組，2 次重復，其中干旱脅迫處理組于大豆V3 期開始進行干旱處理，對照組正常澆水，期間利用托普云農土壤水分含量測定儀對土壤水分含量進行測定，待干旱處理組的土壤含水量低于10%（植株受到干旱脅迫）時測定葉綠素含量[11-12]。

1.3 葉綠素含量測定

利用葉綠素儀SPAD-502Plus 分別測定干旱處理組和對照組的葉綠素含量。每盆隨機選擇3 個單株，分別測定其最上部完全展開三出復葉的葉綠素含量，其中每個葉片選擇非葉脈處測定3次。

1.4 葉綠素含量抗旱指數計算

利用干旱處理及正常澆水條件下的葉綠素含量計算抗旱指數（drought resistance index，DRI），公式如下。

式中，CCD為干旱條件下葉綠素含量；CCW為正常條件下葉綠素含量；CCDM為干旱條件下供試大豆自然群體葉綠素含量平均值。

1.5 葉綠素抗旱指數遺傳變異分析方法

利用SPSS 25.0 對基于葉綠素含量計算的抗旱指數進行描述統計及方差分析，利用公式計算其廣義遺傳力（broad-sense heritability，h2）[10]。

1.6 葉綠素抗旱指數關聯分析方法

利用GEMMA（genome-wide efficient mixed model association algorithm）軟件，采用混合線性模型（mixed linear model, MLM）對葉綠素抗旱指數進行全基因組關聯分析（genome-wide association study, GWAS），其中供試大豆自然群體的基因型數據包含566萬個單核苷酸多態性（single nucleotide polymorphism，SNP）標記[10]；此外，參照已有報道根據LD（linkage disequilibrium）衰減距離，同一LD區域顯著關聯SNP視為1個位點[10,13]。

1.7 候選基因的篩選

在顯著關聯SNP 上、下游100 kb 篩選候選基因，結合大豆公共數據庫Soybase（https://www.soybase.org/）注釋基因功能，初步篩選候選基因；然后利用已有大豆干旱脅迫條件下的轉錄組數據[14]，經標準化處理分析各基因表達量，根據基因差異表達結果確定候選基因。

2 結果與分析

2.1 大豆自然群體葉綠素含量抗旱指數遺傳變異分析

對199 份供試大豆品種資源進行葉綠素含量抗旱指數方差分析（表1）發現，每個環境條件下的供試品種資源抗旱指數均具有極顯著差異，說明干旱脅迫處理有效，品種資源抗旱能力存在較大差異。同時發現，不同環境下的葉綠素含量抗旱指數也存在差別，其平均值由大到小依次為E3、E2、E1，變化范圍為1.35～1.48。進一步分析不同品種資源葉綠素抗旱指數遺傳變異（表1）發現，變異系數為8.71%～11.64%，以E2 環境條件下的變異系數最大。分析葉綠素含量抗旱指數廣義遺傳力發現，其值為81.44%，表明主要由遺傳因素決定。另外，分析葉綠素含量抗旱指數的偏度和峰度系數發現，絕對值均趨近于0，說明葉綠素抗旱指數屬多基因控制數量性狀（圖1）。

表1 大豆自然群體抗旱指數遺傳變異分析Table 1 Genetic variation analysis of drought resistance index in soybean natural population

2.2 大豆葉綠素含量抗旱指數GWAS分析

對3 種環境條件下的葉綠素抗旱指數及其平均值、最佳線性無偏預測（best linear unbiased prediction，BLUP）值進行關聯分析（圖1）發現，關聯到237 個SNPs，分別位于1、2、3、7、10、16、17 及18 號染色體，其中，11 個SNPs 位于基因外顯子并引起非同義突變，28 個位于基因上游或下游，2 個位于基因3’UTR，這些SNPs位點可能影響基因功能或調控基因表達。此外，依據LD 衰減距離（100 kb）可將237 個SNPs 劃分為18 個位點，其中位點1、5、6、11、12、13、14、15、17、18 在E1、E2、E3、平均值及BLUP 值中均被關聯到，屬多環境關聯穩定位點（表2）。

進一步分析發現，在16號染色體關聯到181個SNPs（分布于Locus 10～Locus 14），占總SNPs 的76.37%，其中4 個位點在上述3 種環境以及平均值和BLUP 值同時關聯到，說明這4 個位點不僅是多環境穩定位點，而且是控制葉綠素抗旱指數熱點區域；并且，上述181 個SNPs 中，49 個位于基因內，其中9 個在內含子，17 個在外顯子（6 個屬于非同義突變，9 個屬于同義突變，2 個屬于提前終止），11 個在基因上游，10 個在下游，2 個在3’UTR。此外，在10 號染色體關聯到27 個SNPs（位于Locus 5～Locus 9），其中Locus 5 和Locus 6屬多環境穩定關聯位點，且在上述27 個SNPs 中，2 個位于基因上游或下游，2 個在內含子區域。

2.3 葉綠素抗旱指數顯著關聯SNP驗證

為驗證顯著關聯SNPs及其位點準確性，選取18 個位點的最顯著P值關聯SNPs，對參試品種含有的葉綠素抗旱指數優異等位SNPs 數目及其抗旱指數進行回歸分析，結果（圖2）表明，隨優異等位變異數目增加，葉綠素抗旱指數逐漸增加（R2=0.76）。依據參試品種葉綠素抗旱指數，分別選取高、低極端品種各5個，對每個品種含有的18個位點中的11 個純合SNPs 數目進行分析，結果發現（表3），5個抗旱品種中，3個具有所有的11個優異等位變異，另2 個分別含有9 和10 個優異等位變異，而5個敏感品種則不含任何優異等位變異。

圖2 大豆葉綠素抗旱指數與優異等位變異數量回歸分析Fig. 2 Regression analysis of drought resistance index and elite allele numbers of soybean accessions

表3 大豆抗旱特異種質顯著關聯SNP分析Table 3 Analysis of associated SNPs in different soybean accessions with different drought resistance index

2.4 大豆抗旱相關基因發掘

在顯著關聯SNPs上、下游100 kb范圍內尋找到57 個候選基因（表4），其中10 個基因外顯子存在等位變異，10 個基因內含子存在等位變異，15 個基因上游或下游存在等位變異，2 個基因3’UTR存在等位變異。基于此，根據Soybase網站（https://www.soybase.org）Wm82 基因組注釋信息，從57 個候選基因中初步篩選出16 個相關基因，分別參與植物激素信號轉導、蛋白磷酸化、抗旱反應等；利用已有大豆轉錄組數據[14]分析候選基因表達量（圖3），最終篩選出3 個抗旱相關基因，即Glyma.16G063600、Glyma.10G007000和Glyma.17G143900。

圖3 大豆抗旱候選基因表達分析Fig. 3 Expressions of candidate genes for soybean drought tolerance

表4 大豆抗旱指數候選基因Table 4 Candidate genes of drought resistance index in soybean

Glyma.16G063600存在1 個非同義突變（C/A，天冬氨酸/谷氨酸）、3 個同義突變、2 個下游SNPs、17 個基因間SNPs（圖4）。該基因編碼F-box 蛋白，可與小熱激蛋白結合，進而調控植物耐旱性[15]。依據已發表轉錄組數據[14]發現，該基因受干旱脅迫誘導顯著上調表達。利用該基因的非同義突變分析供試群體發現，含有優異等位變異A與含有非優異等位變異C的品種間的葉綠素抗旱指數差異極顯著。

圖4 候選基因Glyma.16G063600分析Fig. 4 Identification of candidate gene Glyma.16G063600

Glyma.10G007000編碼AP2-EREBP 轉錄因子。AP2/EREBP 轉錄因子含有1 個或多個AP2/ERF 結構域，該類轉錄因子可響應植物鹽、旱、冷等多種非生物脅迫[16]。轉錄組數據分析發現，干旱脅迫條件下，Glyma.10G007000表達顯著上調；利用其上游SNP 分析群體發現，含有優異等位變異的品種其葉綠素抗旱指數顯著高于含有非優異等位變異的品種（圖5）。

圖5 候選基因Glyma.10G007000分析Fig. 5 Identification of candidate gene Glyma.10G007000

Glyma.17G143900編碼AP2/ERF 結構域蛋白。利用基因間最顯著SNP 分析群體發現，含有優異等位變異T的品種其葉綠素抗旱指數顯著高于含有非優異等位變異C的品種（圖6）。

圖6 候選基因Glyma.17G143900分析Fig. 6 Identification of candidate gene Glyma.17G143900

2.5 大豆抗旱優異基因型篩選

利用Glyma.10G007000、Glyma.16G063600及Glyma.17G143900等位變異，對供試品種資源進行分析，結果（表5）發現，可將其分為8種類型，其中基因型Ⅶ（AAT）含3 個優異等位變異，其葉綠素抗旱指數最高；基因型Ⅳ（CAT）、基因型Ⅵ（ACT）及基因型Ⅷ（AAC）含有2 個優異等位變異，葉綠素抗旱指數次之；而不含任何優異等位變異的基因型Ⅰ（CCC）葉綠素抗旱指數最低，說明3 個候選基因與干旱脅迫下的抗旱指數相關，為大豆抗旱育種提供了基礎材料。

表5 供試大豆自然群體基于入選抗旱基因等位變異分析Table 5 Analysis of soybean natural population based on the allele variations of three selected causal genes

3 討論

葉綠素是光能轉化為化學能的關鍵色素[17]。對干旱條件下的葉綠素含量進行研究發現，隨干旱脅迫增加，葉綠素含量呈現先增加后降低的趨勢，說明葉綠素含量受干旱脅迫的影響[18-19]。抗旱指數是指性狀抗旱系數與干旱條件下該性狀指數的乘積[20]。劉桂茹等[21]對小麥抗旱系數、干旱敏感系數、單株產量降低指數分析發現，抗旱指數能反應品種的抗旱性，可作為小麥抗旱篩選的重要指標。楊玉敏等[22]以100 份小麥資源為材料，對抗旱系數、抗旱指數、耐旱指數進行分析發現，抗旱指數綜合考慮了各品種在干旱條件下的敏感度和干旱脅迫下目標性狀在整個群體中的地位，比較適合作為作物抗旱性評價的標準。基于此，本研究利用干旱處理和正常澆水處理下的葉綠素含量計算抗旱指數，用于大豆品種資源抗旱種質篩選以及遺傳位點與候選基因挖掘；并依據群體供試品種資源的葉綠素抗旱指數，篩選出優異種質5 份，分別是‘冀10B5’‘滄0627’‘邯15-685’‘冀11B9’和‘HN1030’，其抗旱指數為1.70～1.83；進一步分析發現，其含有9～11 個優異等位變異，且均含有3 個入選相關基因的優異等位變異，為大豆抗旱育種奠定了重要物質基礎。

目前，有關水稻、小麥等作物的干旱脅迫葉綠素含量遺傳位點發掘研究已有報道，但大豆中開展相關研究較少[23-26]。本研究利用SPAD-502Plus測定干旱及正常條件下供試大豆自然群體葉綠素含量，計算抗旱指數，并依據本團隊前期重測序數據對群體的葉綠素含量抗旱指數進行關聯分析，結果獲得18 個位點，通過與已有研究進行比較發現，Locus3 與李廣軍等[27]定位到的葉綠素含量QTLLeaflet chlorophyll 1-4重合，其加性效應值為0.97，解釋表型變異5%。同時發現，本研究定位到的Locus15與Du等[28]定位到的干旱敏感指數相關QTLDrought index 1-9距離700 kb，該QTL 加性效應值為0.128 1，可解釋表型變異6.42%；且QTLDrought index 1-5與本研究中的Locus10 相距4.2 Mb，該QTL 解釋表型變異9.36%。另外，本研究獲得的Locus3 與Carpentieri 等[29]定位到的耐旱相關QTLDrought tolerance 6-1及Shi 等[30]定位到的葉綠素相關QTLLeaflet chlorophyll 3-1相距2.2 Mb。上述結果證實了本研究結果的可靠性，為大豆抗旱遺傳改良提供了選擇標記。

為篩選大豆抗旱相關基因，依據關聯區域內的基因SNP突變、基因注釋及其表達量篩選出3個抗旱相關基因，其中候選基因Glyma.16G063600編碼F-box 蛋白，該類蛋白既可調控植物生長發育，也可對植物多種生物及非生物脅迫作出響應[31-32]。Qu 等[33]發現，與擬南芥F-box 基因突變體fof2相比，過表達FOF2株系表現為ABA 含量增加及耐旱性提高；Lim 等[34]發現，辣椒F-box 類基因CaDIF1沉默植株表現出對干旱脅迫的敏感性，而過表達CaDIF1植株則表現為耐旱；Xu 等[15]發現，在干旱條件下，過表達大豆F-box 類基因GmFBL144可提高擬南芥過氧化氫及丙二醛含量。本研究篩選的F-box蛋白基因Glyma.16G063600與上述已報道的GmFBL144屬同一亞族，說明其可能在大豆抵御干旱脅迫逆境中具有重要作用。

本研究還篩選出AP2/EREBP（APETALA2/ethylene-responsive element binding proteins）轉錄因子家族基因Glyma.10G007000，該類轉錄因子通常包含1 或2 個AP2/ERF 結構域，依據結構域數目，可將其分為2 個亞族，即EREBP 亞族（包含1 個AP2/ERF 結構域）和AP2 亞族（包含2 個AP2/ERF 結構域）；AP2 亞族可調控花、胚珠及種子發育過程，EREBP 亞族（包括DREB 和ERF 類）參與植物激素調節，并響應非生物脅迫[16]。研究發現，過表達蘋果AP2/EREBP 基因MdSHINE2可增加擬南芥植株莖和葉表面蠟質，并提高其抗旱性[35]。大豆DREB 同源基因GmDREB2在干旱、高鹽及低溫等非生物脅迫條件下被誘導表達，且過表達該基因會激活下游轉錄因子進而提高抗旱性[36]。因此，推斷AP2/EREBP 轉錄因子基因Glyma.10G007000可能參與大豆抗旱反應。

另外，本研究還獲得了諸如MYB 轉錄因子、Dof （DNA binding with one finger）結構域蛋白、CCR-4（carbon catabolite repressor 4）、LRR（Leucinerich repeat）候選基因，已有研究發現這些蛋白均可響應植物干旱脅迫，Cominelli 等[37]發現，AtMYB60通過調控葉片保衛細胞的開關影響植物抗旱；Ma 等[38]在干旱脅迫下對白菜9 個Dof 亞族基因表達模式進行分析發現，不同亞族Dof基因均可響應干旱脅迫；Li 等[39]發現，大豆LRR-RLK基因GmSARK可通過調控葉綠素降解基因Gmlls1及葉綠素合成基因Gmgtr1調控葉片衰老，說明該類基因位于Gmlls1和Gmgtr1上游，可間接調控葉綠素合成或降解。因此，本研究篩選出的這些候選基因可作為今后開展大豆抗旱分子育種研究的重要目標基因。