玉米種子活力相關性狀的全基因組關聯分析

劉文童，趙永鋒，郭晉杰，黃亞群，陳景堂，祝麗英

（1. 河北農業大學 農學院 / 國家玉米改良中心河北分中心 / 河北省作物種質資源實驗室， 河北 保定 071001；2. 青島農業大學 農學院，山東 青島 266109）

玉米是我國播種面積最大的作物，在農業生產中有舉足輕重的地位。種子活力是衡量種子質量的重要指標，對田間發芽、幼苗建成和品種增產性能的發揮有重要的影響。高活力種子發芽早，出苗迅速且整齊一致，對環境的適應性和抵抗力強，有明顯的生長優勢［1-2］。選育和應用高種子活力玉米品種，是我國玉米生產全程機械化發展的必然趨勢。開展種子活力的遺傳研究，對選育高種子活力玉米品種具有重要意義。

種子活力是一個抽象且復雜的綜合概念，1977年國際種子檢驗協會將種子活力定義為：決定種子或種子批在發芽和出苗期間的活性水平和行為的那些種子特性的綜合表現［3］。它與種子發芽、幼苗生長速度和均勻性、廣泛環境條件下種子出苗能力等多方面的表現有關，并非一個單一可測量的性狀，涉及形態學指標、生理生化指標、耐逆性指標，這些相關性狀多為多基因調控的數量性狀，會受到環境和遺傳因素的影響［1-2］。隨著分子標記、測序技術等相關研究理論和技術的發展，連鎖分 析［4］和全基因組關聯分析（GWAS，Genome-wide association studies）［5］被廣泛應用，極大地促進了植物領域性狀的遺傳解析研究。對種子活力相關性狀的遺傳解析，已經在水稻［6］、小麥［7］、大 豆［8］、油菜［9］等多種作物中開展。Hu 等［10］檢測到15 個QTL 與玉米低溫發芽性狀有關，單個QTL貢獻率在3.39%～11.29%之間。Wang 等［11］對重組自交系（RIL，Recombinant inbred line）群體和永久F2群體種子進行不同天數老化處理，基于田間試驗對種子活力進行研究，共檢測到49 個QTL。Ku 等［12］使用2 組相關玉米RIL 群體和加速老化處理方法對種子相關性狀進行研究，連鎖分析共檢測到74 個QTL，Meta-QTL 分析將其整合至20 個MQTL 之中，其中存在32 個候選基因。Han 等［13］采用類似方法研究得到65 個種子活力相關QTL，Meta-QTL 分析將其中貢獻率大于10%的初始QTL整合進18 個MQTL 內，得到23 個候選基因。Shi等［14］使用2 組相關玉米RIL 群體和低溫發芽試驗方法對種子活力開展研究，共檢測到26 個QTL，Meta-QTL 分析將其整合進5 個MQTL 內。田潤苗等［15］對玉米種子萌發相關性狀進行全基因組關聯分析，共得到15 個顯著關聯的SNP，推測可能的候選基因GRMZM2G148411，該基因可能參與種子休眠和萌發的信號調控。Huang 等［16］對玉米發芽期和苗期耐冷性進行全基因組關聯分析，檢測到43 個顯著關聯位點，得到40 個候選基因，并未發現同時與發芽期和苗期耐冷性顯著關聯的位點。Hu 等［17］對玉米發芽期耐冷性進行全基因組關聯分析，得到17 個顯著關聯位點，檢索到18 個候選基因，其中10 個候選基因得到前人QTL 定位研究支持。目前，玉米種子活力的連鎖分析相對較多，全基因組關聯分析相對較少，并且各研究得到的QTL 或顯著關聯位點在數量、位置及效應上存在差異，候選基因也不盡相同，因此需要更多的研究進一步解析玉米種子活力的遺傳基礎。

本研究使用201 份自交系組成的自然群體，通過標準發芽試驗和加速老化試驗測定種子活力相關性狀，采用全基因組關聯分析檢測顯著關聯位點，并搜索相關候選基因，為解析玉米種子活力的遺傳基礎和玉米種子活力分子輔助育種提供理論支持。

1 材料與方法

1.1 試驗材料

201 份玉米自交系組成的自然群體，其中包括我國玉米育種中的骨干自交系及其衍生系、部分國外玉米自交系和部分自選系。2018 年12 月在海南三亞種植自然群體，開花前套袋隔離，人工自交授粉。成熟后統一收獲晾曬，按系混合脫粒，挑選整齊一致的種子保存備用。

1.2 發芽試驗方法

將玉米自交系種子放入托普LH-150 型種子老化箱中進行加速老化處理，設置溫度為45 ℃，相對濕度為90%，處理時間為72 h。

發芽試驗于人工氣候室中進行，采用完全隨機區組試驗設計，3 次重復。標準發芽試驗（SGT，Standard germination test）用正常種子，加速老化試驗（AAT，Accelerated aging test）用加速老化處理的種子。使用種子發芽盒，BP 紙間置床方法進行發芽試驗，每個處理50 粒種子。人工氣候室溫度設置為26±1 ℃，相對濕度為65%±5%，光周期為 14/10 h（晝長/夜長），每天按時加水、通氣。

1.3 相關性狀的測定

參考張紅生和胡晉［1］的方法，發芽培養第4 天開始逐天記錄發芽數至第7 天。調查結束后，每系選5 株長勢均勻的幼苗，測量苗長和根長，稱量苗鮮重和根鮮重。然后放入烘箱中，65 ℃烘干48 h，稱量苗干重和根干重。以5 株幼苗的平均值作 為表型值。

根據以下公式計算發芽勢、發芽率、發芽指數、活力指數：

發芽勢 =（第4 天正常發芽數）/（供試種子總數）× 100%

發芽率 =（第7 天正常發芽數）/（供試種子總數）× 100%

發芽指數 = ∑（Gt/Dt）

活力指數 = ∑（Gt/Dt）×S

以上公式中，Dt為發芽天數，Gt代表對應當天Dt的正常發芽數，S為一定時期的生長量，本研究選用幼苗干重。

使用Microsoft Excel 2010 對數據進行整理，對發芽勢和發芽率進行反正弦轉換用于方差分析，使用IBM SPSS 23 進行描述統計、方差分析。

參考孔繁玲［18］的方法，利用方差分析結果計算廣義遺傳率，計算公式如下：

1.4 全基因組關聯分析

采用CTAB 法［19］提取自交系基因組DNA，GBS（Genotype by sequencing）分型法［20］構建文庫，采用雙末端測序技術、使用第2 代測序儀進行基因測序。測序數據與B73 AGPv2 參考基因組比對，后續數據經質控后，得到83 057 個標記用于全基因組關聯分析。SNP 標記按“染色體編號_位置”規則命名。

將各性狀3 次重復的平均值作為表型值用于全基因組關聯分析。全基因組關聯分析使用R 平臺GAPIT 軟件［21］，分析模型選擇SUPER［22］，參數設置為“CV=CV, sangwich.top=”MLM”, sangwich. bottom=”SUPER”, LD=0.1”，其中CV變量為PCA 矩陣，由TASSEL 軟件［23］計算得到。顯著關聯位點的閾值使用GEC 軟件［24］建議的P值“2.46E-5”，對應LOD值為4.61。

1.5 候選基因搜索

用PLINK 軟件［25］計算基因組上200 kb 內兩兩SNP 間的r2，參考Hill 和Weir 的方法［26］取r2等于 0.1 時的衰減距離120 kb 為該群體的LD 衰減距離。在B73 AGPv2 參考基因組上搜索顯著位點上下游 120 kb 范圍內的基因，基因功能等信息參考UniProt等相關數據庫網站。

2 結果與分析

2.1 種子活力相關性狀的描述統計分析

分別對自然群體在標準發芽試驗和加速老化試驗得到的10 個種子活力相關性狀進行描述統計分析，結果見表1。發芽勢、發芽率、發芽指數、苗長、根長、苗鮮重、苗干重、活力指數的平均值在加速老化試驗條件下略低于標準發芽試驗。標準發芽試驗條件下，各性狀的變異系數范圍為11.20% ～35.01%；加速老化試驗條件下，變異系數為14.69% ～32.66%，表明不同自交系種子活力相關性狀存在較大差異，該群體的種子活力相關性狀有廣泛的遺傳變異，適合于全基因組關聯分析研究。苗長、根長、苗鮮重、苗干重、根干重、活力指數在兩種試驗條件下峰度和偏度的絕對值小于1，表明這些性狀數據呈正態分布，符合數量性狀的特征。

2.2 種子活力相關性狀的方差分析

分別對標準發芽試驗和加速老化試驗測得的種子活力相關性狀進行方差分析，結果見表2。各性狀的區組和基因型效應在2 種試驗條件下均達到極顯著水平，表明兩種試驗條件下各自交系種子活力相關性狀的表現均同時受到環境和基因型效應的影響。各相關性狀的廣義遺傳率在2 種試驗條件下均較高，標準發芽試驗下為83.33%～93.10%，加速老化試驗下為77.88%～91.67%，說明遺傳因素是影響種子活力相關性狀的主要因素。

表1 種子活力相關性狀的描述統計Table 1 Descriptive statistics of seed vigor related traits

表2 種子活力相關性狀的方差分析Table 2 Analysis of variance of seed vigor related traits

2.3 種子活力相關性狀的全基因組關聯分析

分別對自然群體在標準發芽試驗和加速老化試驗條件下的種子活力相關性狀進行全基因組關聯分析，結合測序質控后的83 057 個SNP 分子標記，使用SUPER 模型結合PCA 矩陣方法，在LOD 值大于或等于4.61 判別條件下，10 個性狀共檢測到32個性狀SNP 關聯（見表3），其中標準發芽試驗條件下16 個，加速老化試驗條件下16 個。這些性狀SNP 關聯涉及26 個位點，分布在第1、2、3、4、5、6、7、9 染色體的22 個Bin 區段上，對應染色體的性狀SNP 關聯存在的個數分別為7、10、5、2、3、1、1、3；單個表型貢獻率在0.01%～12.64%之間。

2 個位點與發芽勢關聯，僅在標準發芽試驗中檢測到，分布在第2、4 染色體上，可分別解釋4.80%和1.66%的表型變異。5 個位點與發芽率關聯，分布在第1、2、3、9 染色體上，可解釋0.11%～8.08%的表型變異。7 個位點與發芽指數關聯，分布在第1、2、6、9 染色體上，可解釋0.14%～7.89%的表型變異。2 個位點與苗長關聯，僅在標準發芽試驗中檢測到，分布在第3、4 染色體上，可分別解釋4.34%和12.64%的表型變異。5 個位點與根長關聯，分布在第1、3、5、7 染色體上，可解釋0.30%～3.59%的表型變異。3 個位點與苗鮮重關聯，分布在第1、2 染色體上，可解釋1.09%～8.39%的表型變異，其中位點2_180176428 在2 種試驗條件中同時檢測到。僅1 個位點與根鮮重關聯，在加速老化試驗中檢測到，位于第3 染色體上，可解釋0.35%的表型變異。5 個位點與苗干重關聯，僅在加速老化試驗中檢測到，分布于第1、2、5、9 染色體上，可解釋0.01%～10.85%的表型變異。僅1 個位點與活力指數關聯，在標準發芽試驗中檢測到，位于第5 染色體上，可解釋0.75%的表型變異。

有5 個位點在不同性狀中同時檢測到，其中位點1_127505950、1_288806019、2_3108015、 和9_97280218 同時與發芽率和發芽指數顯著關聯；位點2_180176428 同時與苗鮮重和苗干重顯著關聯，可能存在一因多效的現象。

表3 種子活力相關性狀的全基因組關聯分析結果Table 3 Results of seed vigor related traits in genome-wide association studies

續表：

2.4 種子活力相關性狀的候選基因

基 于B73 AGPv2 參 考 基 因 組， 在 顯 著 關 聯SNP 位點上下游120 kb 的范圍內搜索候選基 因，綜合文獻和數據庫信息推測可能的候選基因（見表4），種子特異蛋白GRMZM2G163912、組蛋白GRMZM2G164020、WRKY轉錄因子GRMZM2G164082、 果糖呋喃糖苷酶GRMZM2G018716、MYB 基因家族GRMZM2G035370、Ⅱ型內含子成熟酶家族蛋白GRMZM2G375999 與發芽率、發芽指數有關，酸性內切幾丁質酶GRMZM2G023650、LG2 轉錄因子GRMZM2G060216、蔗糖合酶GRMZM2G045171與苗長有關，胞質醛脫氫酶GRMZM2G071021、醛 脫 氫 酶GRMZM2G097706 與 根 長 有 關，TIFY 蛋 白GRMZM2G064775、ABA 誘 導 的 蛋 白GRMZM2G063693 與苗干重有關。

表4 種子活力相關性狀可能的候選基因Table 4 Possible candidate genes of seed vigor related traits

3 討論

3.1 種子活力相關性狀的定位結果比較

前人對于種子活力相關性狀基因定位的研究多基于雙親群體進行QTL 定位，劉海英［27］在第1、2、3、4、5、8、9、10 染色體和韓贊平［2］在10 條染色體上均檢測到玉米種子活力相關性狀的QTL，本研究使用自然群體201 份自交系材料對玉米種子活力相關性狀進行全基因組關聯分析，共獲得26 個顯著關聯位點，分布在8 條染色體上22 個Bin 區段中，并未在第8 和10 染色體上檢測到顯著關聯位點。

由于玉米種子活力相關性狀定位研究較少以及標記類型和參考基因組間存在差異，本研究定位結果位點與其他相關研究定位結果很難進行清晰比較。本研究僅發現檢測到的部分關聯位點與前人研究位于同一Bin 區段，如劉海英［27］研究中Bin1.02（qGEF1，umc1568-bnlg1007），Bin1.05（qNGE1，umc1323），Bin1.08（qNGI1，phi039-umc2029），Bin3.05（qSDW3，bnlg1957），Bin9.03（qLGE9，umc1691）；Wang 等［11］研究中Bin1.03（umc1044- phi339017），Bin1.08（phi038-dupssr12），Bin6.04（umc2006-umc1979）。另外，本研究還有部分顯著位點位于相關研究中MQTL 區域內。如與Ku 等［12］研究相比較，位點2_109576143 位于MQTL2-1（Chr2:33590242-144073650） 內， 位點3_176783445 位 于MQTL3-3（Chr3:175554472-184720933） 內， 位 點5_36667115 和5_44370621位 于MQTL5-1（Chr5:17380699-59293558） 內；與Han 等［13］研究相比較，位點3_176783445 位于MQTL3-4（Chr3:175554472-184720933）內，位點6_106468212 位 于MQTL6-1（Chr6:86257528-123864232）內。目前，限于玉米種子活力相關性狀基因定位研究數量較少，暫未發現玉米種子活力相關性狀定位的熱點染色體區域。

3.2 種子活力相關性狀的候選基因分析

本試驗搜索到的GRMZM2G163912、GRMZM2G164020、GRMZM2G164082、GRMZM2G018716、GRMZM2G035370、GRMZM2G375999與發芽率、發芽指數有關。GRMZM2G163912編碼種子特異蛋白，Lopez-Malvar 等［28］研究認為該基因模型是對香豆酸的候選基因，對香豆酸與細胞壁加強和強化有關、與植物防御害蟲和病原體有關；?ukalovi? 等［29］認為在植物細胞發育過程中對香豆酸多轉化為木質素，可能參與防御機制，對香豆酸和阿魏酸產生化感作用對種子萌發產生影響。GRMZM2G164020編碼組蛋白，Teoh 等［30］發現其與DNA 合成有關并在胚成熟期間特異大量表達。GRMZM2G164082屬于WRKY 基因家族，WRKY 轉錄因子作用于植物生長發育，并對生物和非生物脅迫做出反應，WRKY 蛋白被認為與發育進程調控有關，參與種子發育、胚胎發生、葉片衰老、休眠、植物生長和代謝途徑等發育過程［31］。GRMZM2G018716編碼果糖呋喃糖苷酶，Zhou 等［32］認為其與棉子糖代謝有關，而植物細胞棉子糖的積累與冷、熱、旱等環境脅迫因素響應有關。GRMZM2G035370屬于MYB 基因家族，MYB基因家族參與次生代謝調控、細胞形態發生調控、分生組織形成調控、花器官和種子發育調控、細胞周期調控等多種生理生化過程，有些還參與各種防御和應激反應以及光和激素信號通路［33］。GRMZM2G375999與Ⅱ型內含子成熟酶家族蛋白有關，可能參與種子萌發等相關生理過程。

GRMZM2G023650、GRMZM2G060216、GRMZM2G045171與苗長有關。GRMZM2G023650編碼酸性內切幾丁質酶，Hawkins 等［34］認為該基因模型屬于糖苷水解酶家族，其他相關研究結果表明植物誘導表達的植物幾丁質酶蛋白的主要功能是單獨或結合β-1,3-葡聚糖酶聯合作用防御真菌病原體，組成型表達的幾丁質酶在胚胎發生、成花、衰老、種子萌發等一系列生理和形態學過程中發揮作用。GRMZM2G060216與LG2 轉錄因子有關，可能參與葉片發育正調控，Zhang 等［35］認為其可能為bZIP轉錄因子，bZIP 轉錄因子在調控種子成熟和發芽方面發揮重要的作用。GRMZM2G045171編碼蔗糖合酶，Schlüter 等［36］認為其與纖維素合成轉錄因子存在相關性，可能參與細胞壁代謝；Trucillo Silva 等［37］認為其與氮代謝有關，參與蔗糖和二磷酸尿苷葡糖UDPG 間的可逆反應，動員蔗糖參與如細胞壁建成、淀粉形成等多種利用活性糖的途徑。

GRMZM2G071021、GRMZM2G097706與根 長有關。GRMZM2G071021編碼胞質醛脫氫酶，GRMZM2G097706編碼醛脫氫酶，可以代謝各種內源性和外源性醛［38］。

GRMZM2G064775、GRMZM2G063693與 苗 干重有關。GRMZM2G064775編碼TIFY 蛋白，Zhou等［39］認為其是耐旱性和生長調控的重要基因，在JA、ABA 和GA 信號途徑處于交叉樞紐位置。GRMZM2G063693編碼ABA 誘導的蛋白，Hu 等［40］認為其可能屬于油體蛋白基因，可能與一系列包括應激反應、植物氧脂生物合成、細胞質信號轉導、脂質運輸、膜擴張和油體生物發生等細胞和生理過程有關，并發現該基因在種子和雄蕊中特異表達。本研究發現的多數候選基因與抗逆性研究有關，這些候選基因對玉米種子活力相關性狀研究提供了參考，對種子活力的遺傳基礎解析仍需大量相關后續研究來驗證。

4 結論

種子活力相關性狀的表現受到環境和基因型效應的共同影響，但基因型效應的影響較大。標準發芽試驗和加速老化試驗條件下，全基因組關聯分析共檢測到32 個性狀SNP 關聯，涉及26 個SNP 位點，推測可能的候選基因13 個，研究結果可為解析玉米種子活力遺傳基礎和分子輔助選擇育種提供參考。