《Plants》：Multi-Season Multi-Model GWAS Prioritizes Stable Genomic Loci and Candidate Genes for Six Agronomic Traits in Soybean Under Conditions in Southeastern Kazakhstan
编辑推荐：
：在可变生长季节中可靠地定位控制复杂农艺性状的基因组区域，仍是大豆遗传育种面临的主要挑战。本研究利用一个由252份大豆种质组成的多样化群体，在2018–2023年连续六个生长季节中对开花时间、成熟期、株高、单株粒数、单株产量和千粒重进行了表型评价。通过全基因组
：在可变生长季节中可靠地定位控制复杂农艺性状的基因组区域，仍是大豆遗传育种面临的主要挑战。本研究利用一个由252份大豆种质组成的多样化群体，在2018–2023年连续六个生长季节中对开花时间、成熟期、株高、单株粒数、单株产量和千粒重进行了表型评价。通过全基因组重测序和变异过滤，获得2,019,772个高质量单核苷酸多态性（SNP），并采用整合式包容性输入多位点随机SNP效应混合线性模型（IIIVmrMLM）、贝叶斯信息与连锁不平衡迭代嵌套关键通路法（BLINK）以及多位点混合模型（MLMM）评估关联信号，同时结合基于连锁不平衡（LD）的位点整合。跨模型优先化保留了21个由三种全基因组关联分析（GWAS）模型共同支持的高置信度位点，分布于10条染色体上。在鉴定到的位点中，18个与先前报道的SoyBase基因和数量性状位点（QTL）重叠或共定位，而3个位点（q.VER2.13-1、q.YP.01-1和q.TSW.15-1）与已知基因和QTL无位置重叠，因此被认为可能为新位点。这3个位点分别与开花时间、单株产量和千粒重相关，分别解释1.60%、2.13%和5.53%的表型变异。10个位点与包含已确立的大豆调控因子的基因组区域共定位，包括E2、E3、GmDt2和POWR1，支持了关联结果的生物学合理性。整合基因组位置、功能注释和组织表达证据，在18个位点中优先筛选出112个候选基因。这些发现为在大陆性多变生长条件下大豆适应性和产量形成的遗传基础提供了聚焦的基因组位点和候选基因集，以供独立验证和进一步研究。
论文解读
一、研究背景与意义
大豆是全球最重要的豆科作物之一，为食品、饲料和工业应用提供植物蛋白和植物油的重要来源。在气候条件日益多变的情况下，维持大豆生产力已成为育种领域的重大挑战，因为温度和降水量的年际波动会强烈影响植株发育、物候和籽粒产量。因此，鉴定能够在不同生长季节中赋予稳定农艺表现的遗传因子，是现代大豆改良的关键目标。
在育种所关注的大豆农艺性状中，开花时间、生理成熟期、株高、粒数、粒重和籽粒产量尤为重要，因为它们共同决定了环境适应性和生产力。开花和成熟调控营养生长与生殖生长的持续时间，进而影响生物量积累、冠层结构、同化物分配和籽粒充实。株高影响抗倒伏性和光截获能力，而单株粒数和千粒重则是决定最终籽粒产量的主要组分。这些性状通过复杂的生理和遗传网络在发育上相互关联，其表达受环境条件的强烈调控。
复杂性状的遗传结构以众多效应大小各异的位点、广泛的互作效应以及显著的基因型与环境互作（G×E）为特征。尽管包括E、J、Tof和FT基因家族成员在内的多个控制开花和成熟的主效基因已被功能鉴定，但它们仅解释了所观察表型变异的一部分。大多数农艺性状仍呈高度多基因遗传，众多位点贡献小至中等效应的等位变异，且常受环境条件影响。因此，解析这些性状的稳定遗传结构仍是大豆遗传育种的主要目标。
高通量基因分型技术的快速发展和大豆基因组资源的持续完善，使全基因组关联分析（GWAS）成为解析复杂数量性状最有力的方法之一。通过利用多样化种质资源中的历史重组，GWAS比传统双亲连锁作图提供更高的定位分辨率，同时捕获更广泛的自发等位变异。然而，GWAS发现向实际育种的转化仍受到许多关联信号跨环境可重复性有限的制约。大量在单一年份或单一环境实验中鉴定的位点在气候条件变化时无法维持显著效应，因为表型表达受环境变异强烈影响。这一局限性在产量组分性状中尤为突出，G×E互作常占总表型方差的相当比例。因此，仅鉴定统计学显著的位点是不够的；育种应用需要那些在年际环境波动中仍表现稳定遗传效应的基因组区域。多环境GWAS因此成为联合建模环境异质性、提高统计功效、改进效应估计精度以及促进鉴定具有更高育种价值的环境稳定位点的有效策略。
方法学进展进一步提高了GWAS解析复杂遗传结构的能力。传统的单位点混合线性模型能有效减少群体结构和亲缘关系导致的假阳性关联，但常因严格的多重检验校正而统计功效有限。为克服这些局限，多种多位点方法被提出，包括MLMM、BLINK以及新近发展的IIIVmrMLM框架。通过同时建模多个数量性状核苷酸，这些方法显著提高了小效应和中效应位点的检测能力，同时保持对假发现的适当控制。特别是IIIVmrMLM通过明确纳入多环境和数量性状核苷酸与环境互作，扩展了传统的多位点关联分析，从而能够更全面地刻画环境响应性和组成性表达的位点。
尽管过去十年报道了大量大豆GWAS研究，但将长期表型评价、多环境多位点关联分析和鉴定位点的系统生物学解释相结合的研究相对较少。大多数先前研究侧重于大群体或相对短期试验。相比之下，在单一试验点对遗传多样化种质进行连续多个生长季节的重复评价，为评估遗传效应在年际环境变异下的一致性提供了互补框架。这一设计最小化了地点间地理差异带来的变异，同时捕获了气候条件的年际变化。然而，从这一设计推断的稳定性仅针对试验地点的生长季节，并不能证明跨地理上不同环境的广泛适应性。
这一策略对大陆性气候条件下的大豆生产尤为重要。在哈萨克斯坦，大豆生产集中于东南部地区，该地区温度、降水和季节长度的显著年际波动强烈影响作物表现。在该环境中，开花和成熟决定生长与有利环境条件的同步性，而株高和产量组分受季节变化强烈塑造。因此，利用多年GWAS补充现有遗传作图数据，对于发现能够推动气候适应型大豆品种培育的韧性基因组区域至关重要。
对统计显著关联信号的生物学解释同样重要。大豆参考基因组、功能注释、转录组学以及SoyBase等精选基因组资源的发展，显著增强了GWAS位点下生物学合理候选基因的鉴定能力。将关联作图与连锁不平衡信息、先前报道的数量性状位点、比较功能注释和基因本体分析相结合，为候选基因优先化提供了更强证据，并促进了从统计关联向生物学机制和分子育种应用的转化。
本研究的目的是：(i) 表征在连续六个生长季节中评价的全球多样化大豆群体的六个主要农艺性状的表型变异、方差组分和广义遗传力；(ii) 通过整合长期表型评价、最佳线性无偏估计（BLUE）表型调整、全基因组重测序和互补GWAS方法，鉴定物候、形态和产量相关性状的稳定基因组位点；(iii) 通过基于连锁不平衡的位点整合和独立跨模型验证优先化稳健基因组区域；(iv) 利用基因组注释、基因本体分类和组织特异性表达分析鉴定并功能表征优先化稳定位点下的候选基因；(v) 评价所鉴定稳定基因组位点和候选基因作为大陆性环境条件下分子育种和大豆改良基因组资源的潜在应用。
二、主要技术方法
研究人员构建了一个包含252份大豆种质的关联群体，这些种质来自欧洲、亚洲和北美洲20个国家，代表广泛的遗传多样性。田间试验于2018至2023年连续六个生长季节在哈萨克斯坦东南部的哈萨克农业与植物栽培研究所（KRIAPG）试验田进行，采用随机完全区组设计（RCBD），设两次田间重复。对六个性状进行表型评价：开花时间（VER2）、成熟期（VER8）、株高（PH）、单株粒数（NSP）、单株产量（YP）和千粒重（TSW）。利用混合线性模型计算BLUE值作为GWAS表型数据。全基因组重测序数据比对至中黄13（ZH13 v2.0）参考基因组，经严格质控后获得2,019,772个高质量SNP。群体结构通过主成分分析（PCA）和基因组亲缘关系矩阵评估。采用三种互补的多位点GWAS方法——IIIVmrMLM、BLINK和MLMM进行分析，初始筛选阈值为p < 1.0 × 10?4。基于全基因组LD衰减距离（2.196 Mb）进行位点整合，并通过跨模型一致性验证确定高置信度位点。候选基因鉴定基于ZH13 v2.0基因组注释，功能注释和表达分析利用SoyMD数据库。
三、研究结果
3.1 六个生长季节的表型变异
六个农艺性状均表现出显著的表型变异，各性状的年际变异幅度存在差异。所有性状均呈现数量遗传特征的连续表型分布。物候性状（VER2和VER8）的变异系数低于产量相关性状。VER2的平均值从2018年的43.6天逐渐降至2023年的36.0天，VER8的平均值介于112至118天之间。PH、NSP和YP表现出较宽的表型范围和较高的变异系数，其中YP的年际变异最大，而TSW的季节波动相对较小。基于六个生长季节计算的BLUE值相对于单年观测呈中间分布。基于BLUE值的皮尔逊相关分析显示农艺性状间存在强正相关，VER2与VER8相关性最强（r = 0.87），VER8与PH（r = 0.87）、NSP与YP（r = 0.84）次之，PH与NSP（r = 0.79）和YP（r = 0.71）正相关，TSW与其余性状仅呈弱至中等相关。
3.2 方差分析与遗传力
方差分析（ANOVA）显示基因型、生长季节和基因型×生长季节互作对全部六个性状均有显著效应。广义遗传力（H2）范围为0.235至0.913，VER8最高（0.913），其次为PH（0.758）和VER2（0.743），NSP（0.348）、YP（0.313）和TSW（0.235）较低。基因型效应在所有性状中均显著（p < 2 × 10?16），F值范围为NSP的14.64至VER8的427.98。基因型×生长季节互作效应也均显著，F值范围为PH的2.60至TSW的55.83。总体而言，H2按VER8 > PH > VER2 > NSP > YP > TSW的顺序递减。
3.3 基因分型数据与连锁不平衡
高质量SNP标记密集分布于全部20条大豆染色体，提供了近乎全基因组覆盖的标记。全基因组LD分析显示，随着SNP间物理距离增加，成对LD逐渐下降。基于平滑的全基因组平均LD曲线，LD衰减阈值（r2 = 0.20）在2.196 Mb处达到。各染色体间LD衰减模式相似，但LD程度存在染色体间差异。
3.4 群体结构
PCA和基因组关系矩阵显示该大豆关联群体存在中等程度的群体结构。前两个主成分解释了11.44%的分子变异，表明尽管种质地理来源广泛，群体分层相对较弱。来自相似地理区域的种质倾向于形成部分重叠的聚类。基因组亲缘矩阵显示252份种质间成对遗传关系普遍较弱，仅有限数量的基因型密切相关。
3.5 全基因组关联分析与稳健关联位点优先化
在六个农艺性状中，初始多位点GWAS筛选在p < 1.0 × 10?4下保留了8,537个标记-性状关联（MTA），其中IIIVmrMLM鉴定124个，BLINK鉴定6,679个，MLMM鉴定1,734个。基于全基因组LD结构，将LD定义区间（2.196 Mb）内的SNP整合为单一QTL，由统计支持最强的先导SNP代表。该程序将关联信号从8,537个SNP减少至198个非冗余基因组位点。其中118个仅由单一GWAS框架检测，80个由两个或以上模型检测。经跨模型验证，80个位点至少由两个GWAS模型支持而被保留。应用更严格的三模型共同支持标准后，最终保留21个高置信度稳定位点。
3.6 高置信度稳定基因组位点及其与已知大豆基因和QTL的关系
21个高置信度稳定位点分布于10条染色体上，VER2（5个）和PH（5个）最多，其次为NSP（4个）、VER8（3个）、TSW（2个）和YP（2个）。基于2,019,772个SNP的Bonferroni校正全基因组显著阈值为p < 2.48 × 10?8，18个先导关联超过该阈值，其余3个基于跨模型一致性保留。与已知大豆基因的位置比较发现10个位点与主要基因邻近：10号染色体上的q.VER2.10-1、q.NSP.10-1和q.PH.10-1与E2（GmGIa）区域共定位；19号染色体上的多个位点与E3（GmPhyA3）共定位；q.PH.18-1与GmDt2邻近；q.NSP.20-1和q.YP.20-1位于POWR1/GmCCT基因组区域内。与SoyBase数据库比较显示，21个位点中17个与至少一个先前报道的QTL存在位置重叠，而q.VER2.13-1、q.NSP.19-1、q.YP.01-1和q.TSW.15-1未检测到重叠QTL。其中三个位点（q.VER2.13-1、q.YP.01-1和q.TSW.15-1）被认为可能为新位点，分别解释1.60%、2.13%和5.53%的表型变异。
3.7 候选基因鉴定、功能注释与表达谱分析
利用中黄13（ZH13 v2.0）参考基因组注释在21个高置信度稳定位点的LD区间内鉴定候选基因。基因本体（GO）分析显示候选基因在生物学过程类别中广泛分布，物候性状相关基因在蛋白质生物合成（33个基因）、胁迫与防御（22个）和初级代谢（21个）类别中占主导。分子功能分析中，结合（33个基因）和结构分子活性（31个）为最大功能组。细胞组分分析显示细胞质包含最多候选蛋白。组织特异性表达分析显示物候相关基因在所有15个组织中最为丰富，产量相关基因在豆荚和茎组织中表达较高。在21个稳定QTL中，18个包含转录本丰度超过200 TPM的高表达候选基因，最终优先化出112个候选基因。这些基因被分为12类主要生物学机制，其中激素调控（18个基因）、细胞壁重塑（15个）、初级代谢（12个）和胁迫与物候（11个）为最大功能类别。
四、讨论与结论
讨论部分指出，六年表型评价和BLUE估计为GWAS提供了稳健的表型基础。VER2、VER8和PH具有较高遗传力和较大基因型贡献，而NSP、YP和TSW受生长季节和互作效应影响更大。TSW的遗传力最低（0.235）且基因型×生长季节互作最强（F = 55.83），这可能解释了为何仅有有限数量的稳定TSW位点被保留。研究强调该设计为多季节单地点实验，稳定性推断限于该地点的年际环境变异，跨地理环境的验证仍属必要。
LD整合和跨模型验证显著减少了初始关联空间，8,537个MTA经整合后减少至198个位点，跨模型比较保留80个，最终21个高置信度位点仅占初始MTA的0.25%。三种GWAS方法提供互补统计框架，一致性检测作为稳健性标准。
已知发育基因的恢复为关联作图框架提供了生物学支持。E2和E3区域与多个物候和形态位点共定位，GmDt2和POWR1/GmCCT区域与株高和产量相关位点共定位，支持了关联结果的生物学合理性。
三个可能新位点（q.VER2.13-1、q.YP.01-1和q.TSW.15-1）与开花时间、单株产量和千粒重相关，其中q.TSW.15-1因效应较大而成为独立验证的有前景靶标。这些区域应被视为候选基因组资源而非可立即部署的育种位点。
112个优先化候选基因代表多样生物学过程，组织表达谱为优先化提供了补充信息。研究强调这些证据是优先化标准而非功能验证，精细定位、单倍型分析和功能实验仍是区分因果基因与连锁候选所必需的。
结论：连续六年田间评价结合全基因组重测序和互补GWAS模型，鉴定了21个与大豆物候、株型结构和产量组分相关的高置信度基因组位点。10个位点与包含E2、E3、GmDt2和POWR1等已确立调控因子的区域共定位，支持关联结果的生物学合理性。3个位点（q.VER2.13-1、q.YP.01-1和q.TSW.15-1）与SoyBase比较中纳入的基因和QTL无位置重叠，被分类为可能新位点。整合基因组位置、功能注释和组织表达证据进一步优先化了18个位点中的112个候选基因。三个可能新位点及其相关候选基因代表独立验证、精细定位、单倍型分析和功能表征的优先靶标。由于本研究在单一地点进行，需要在额外大豆群体和地理上不同环境中进行验证，以确定这些位点更广泛的环境稳定性和育种价值。
