玉米种子活力相关性状的全基因组关联分析

刘文童，赵永锋，郭晋杰，黄亚群，陈景堂，祝丽英

（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 个，研究结果可为解析玉米种子活力遗传基础和分子辅助选择育种提供参考。