上位性在复杂性状遗传剖析、物种进化机制和杂种优势遗传基础解析等方面具有重要作用, 因此有必要研究上位性检测方法^[1-3]。基于双亲分离群体数量和多歧性状上位性检测方法, 已有惩罚最大似然^[4]、经验贝叶斯^[5]、LASSO^[6]、分层广义线性模型^[7-8]和EBLASSO^[9]等, 它们同时考虑主效QTL、QTL间互作和QTL与环境互作效应, 能更全面解析复杂性状的遗传基础。目前, 在人类和植物数量性状关联分析中上位性检测方法有长足的发展, 包括参数^[10-13]和非参数方法^[14-15]两大类。

作物品种是多次杂交重组的结果, 由它们组成的品种资源群体遗传变异大, 蕴藏可供作物遗传育种利用的优异等位基因, 为育种设计和优异亲本组合预测提供素材。在探索数量性状上位性方面, Wang等^[16]的自适应混合LASSO和Lü等^[17]的经验Bayes可处理模型变量个数少于样本容量10倍的数据集, Wen等^[13]的EBLASSO可处理更大的数据集。若要提高运算速度, 可采用Zhang等^[12]的基于GPU计算和混合模型的方法。抗性性状是主要的育种目标, 在统计遗传中称为多歧性状。为了挖掘这类性状的优异等位基因, 有必要探索品种资源群体多歧性状上位性检测方法。然而, 多歧性状的非孟德尔遗传和信息量少等特点, 使得其关联分析研究的工作更加困难。所以, 品种群体多歧性状QTL上位性检测工作相对较少。目前, 混合线性模型和Logistic回归模型、贝叶斯和惩罚似然等参数估计方法以及函数逼近和迭代等技术的应用, 为多歧性状上位性关联分析方法的研究打下坚实基础^{[9, 18-21]}。

Feng等^[22]基于伪数据正态逼近算法和经验Bayes算法^[23]提出的分层广义线性混合模型方法可检测多歧性状的主效QTL, 本研究在此基础上拓展了多歧性状QTL与环境互作和QTL间互作检测的新方法, 并通过Monte Carlo模拟试验和大豆幼苗主根长耐盐碱性指数分级数据分析, 证实该方法的有效性。

1 材料与方法 1.1 试验材料

257份大豆品种的135个SSR分子标记数据及其2009与2010两年耐盐碱性指数分级表型数据^[24]用于新方法的试验验证。其中将大豆品种的耐盐碱性指数^[24]按照0~20%、20%~40%、40%~60%、60%~80%和80%~100%将其耐盐碱性级别划分为高耐、耐、中耐、感和高感5种, 依次记为1到5, 作为本研究的表型数据。ATI代表耐碱性性状, STI代表耐盐性性状。

用STRUCTURE软件^[25]分析135个SSR标记基因型数据。随着亚群数增加, 对数似然函数值增加。因此, 采用ΔK方法和育种经验, 确定亚群数K=4为适宜亚群数^[24], 并获得各品种属于各亚群后验概率的Q矩阵, 以控制群体遗传结构。

1.2 广义混合线性模型

设多歧性状表型观测向量为Y={y_j}_j=1ⁿ, 其中y_j是第j品种表型观测值, 分成C类, 即y_j的取值属于1~C类。为建立数学模型, 引入一个连续的潜在变量l_j与多歧性状表型y_j相对应, 即:

式中:t_c为C+1个假定阈值, 且t₀=-∞ < t₁ < … < t_C=+∞。由此, l_j的混合线性模型为:

(1)

式中:β₀为群体平均值; β₀₀为群体结构效应向量; (β₁, …, β_q-1)^T为环境效应, q为环境个数; 向量γ_k(k=1, …, p) 是第k个位点的g_k-1维QTL效应, g_k是位点k的等位基因个数; 向量ζ_mk是第m个位点和第k个位点的 (g_m-1)(g_k-1) 维上位性效应; 向量ξ_mk是 (q-1)(g_k-1) 维的QTL与环境互作效应, 与γ_k、ζ_mk均为随机效应; ε_j是随机误差, 服从标准正态分布。(z₁^T z₂^T … z_n^T)^T是群体结构的Q矩阵; x_jm和u_jk分别是第j品种的参数β_m和γ_k的哑变量; 向量w_mk是参数ζ_mk的哑变量, 其值由u_jm和u_jk确定。

将群体平均值、环境效应和群体结构等固定效应合并记为向量β, 其设计阵记为X_j; 将QTL效应、上位性效应和QTL与环境互作效应等随机效应合并为随机向量γ, 其设计阵记为Z_j, 则将模型 (1) 简化为混合线性模型:

(2)

当参数β、γ和t=(t₀, t₁, …, t_C) 已知时, 表型y_j条件独立, 其对数似然函数为:

(3)

式中:Φ(·) 为标准正态分布的累积分布函数; θ={t, β, γ}为参数向量; I_yj=c为示性函数。

1.3 参数估计 1.3.1 先验分布

在模型 (2) 中, 参数γ为随机效应。采用分层广义线性模型^[7], 以解决随机效应个数远远大于样本容量的超饱和线性模型问题。假定每个随机QTL效应γ_k(k=1, …, vq) 的方差σ_k²不同, 并通过假定如下的先验分布建立了分层模型:

式中:τ和a是常数。通过模拟研究发现, (τ, a)=(1, 0.5) 的效果较好。

1.3.2 参数估计

根据正态逼近似然方法^[26-27], 可将式 (3) 的似然函数Pr (y_j|X_j, Z_j, θ) 转化为正态似然函数N(w_j|μ_0j, σ_0j², t)。其中, , μ_0j=X_jβ+Z_jγ, σ_0j²=; L_j=log{Φ(u_c)-Φ(u_c-1)}为Pr (y_j|X_j, Z_j, θ) 的对数似然函数, L′_j和L″_j分别为L_j的一阶和二阶微分, u_c=t_c-X_jβ-Z_jγ。

由牛顿迭代、后验分布和经验Bayes可得到β、γ_k、σ_k²、ω_k、b_k和阈值t_c的估计量^[22]。

1.3.3 迭代过程

1) 初始值设置γ_k⁽⁰⁾、t_c⁽⁰⁾、ω_k⁽⁰⁾、σ_k²⁽⁰⁾和b_k⁽⁰⁾; 2) 参数σ_k²、ω_k、b_k、γ_k、t_c和β由其估计式依次更新; 3) 重复执行第2步, 直到满足预先设定的收敛标准为止。

1.4 统计检验

用两步选择方法^[17]对所有QTL进行似然比测验 (LRT)。

1.5 Monte Carlo模拟方案

用Zhang等^[28]的玉米群体来模拟4组试验, 以验证新方法的有效性。假定祖先亲本Founders是连锁平衡的。Non-founders基因型是由它们的亲本基因型产生的, 关联群体由Non-founders品种组成。假设在4条100 cM长的染色体上, 等距离设置44个标记, 每个标记有2个等比例等位基因。4个主效QTL (main-effect QTL, 简记为MQ)、7个上位性QTL (QTL-by-QTL interaction, 简记为QQ) 和2个QTL与环境互作 (QTL-by-environment interaction, 简记为QE) 的位置与效应参数设置见表 1。重复200次。字母Q标“*”表示该QTL含有主效, 即QQ^*和Q^*Q^*分别表示1个和2个位点有主效, Q*E表示主效QTL与环境互作。

在模拟试验1中, Founders数为100, Non-founders数为300, 表型分布设置为1:1:1:1:1、1:2:4:2:1、1:3:9:3:1和8:5:3:1:1共4种类型。在模拟试验2中, Founders数和Non-founders数同模拟试验1, 表型分类数设置为2(1:1)、4(1:3:3:1) 和6(1:3:6:6:3:1) 共3种类型。在模拟试验3中, 表型分布为1:2:4:2:1, Founders数为100, Non-founders数设置为200、300和500共3种类型。在模拟试验4中, 表型分布为1:2:4:2:1, Non-founders数为300, Founders数设置为10、50、75和100共4种情况。

QTL存在的显著性阈值为LOD>2.0。检测到的QTL样本数与重复总数200的比值为QTL检测的统计功效。假阳性效应个数与零效应总数的比值为假阳性率。

2 结果与分析 2.1 Monte Carlo模拟研究 2.1.1 模拟研究

4个模拟试验的研究结果见图 1。图 1-a是Founders数100、Non-founders数300、表型分类数5、表型分布1:2:4:2:1条件下, 不同遗传率 (0.01、0.03、0.05和0.10) 对MQ和QQ检测功效影响的模拟结果; 图 1-b是遗传率0.05, 其他参数同图 1-a下, MQ、QQ和QE三者功效比较的模拟结果。结果表明, MQ和QQ检测功效随遗传率增加而增加 (图 1-a); 遗传率小于5%时, 相同遗传率QE和QQ检测功效要高于MQ (图 1-b)。这种趋势在均匀分布和对称分布、分类数较少、遗传率较小、Founder数少和Non-founders数少情况下更明显^[29]。

新方法的假阳性率在不同参数设置情况下均较小 (图 2)。但是当样本容量太小 (如Non-founders数小于10) 时, 假阳性率会较高。这是因为样本容量越小, 数据所包含的信息量越少, 检测结果可靠度和准确度越低。新方法在检测MQ、QQ和QE时, QTL位置和效应以及阈值的估计值偏差都相对较小^[29]。

2.1.2 表型分布的影响

模拟试验1研究了表型性状分类数为5时, 不同表型分布对新方法的影响。模拟了均匀分布 (1:1:1:1:1)、对称分布 (1:3:9:3:1和1:2:4:2:1) 和偏态分布 (8:5:3:1:1)4种类型, 结果 (图 3-a和文献[29]) 表明:MQ、QQ和QE的遗传率较小 (如0.03) 时, 对称分布检测功效最高, 且对称分布变异越小 (如1:3:9:3:1) 检测功效越大; 其遗传率较大 (0.05或0.10) 时, 表型分布类型对检测功效的影响较小。

2.1.3 表型分类数的影响

模拟试验2是在表型分布对称条件下, 比较了表型分类数2(1:1)、4(1:3:3:1) 和6(1:3:6:6:3:1)3种情况对新方法的影响。结果 (图 3-b和文献[29]) 表明:MQ、QQ和QE的统计检测功效随分类数增加而增加。这是因为随着分类数增加, 表型观察值提供的生物信息量会增加, 离散表型更接近连续表型。

2.1.4 Non-founders数的影响

模拟试验3是在Founders数100时, 研究Non-founders数200、300和500对QTL检测结果的影响。结果 (图 3-c和文献[29]) 表明:Non-founders数越大, MQ、QQ的检测功效越大; 随着遗传率增加, MQ和QQ的检测功效增加; 遗传率为0.05时, QE的功效均较大, 受Non-founders数影响不显著。

2.1.5 Founders数的影响

模拟试验4比较了Founders数10、50、75和100时对QTL检测结果影响。结果 (图 3-d和文献[29]) 表明:随着Founders数增加, MQ和QQ的功效均增大, 但增加幅度逐渐减小; 而QE检测功效都较大 (遗传率为0.05), 受Founders数的影响不明显; QQ和QE检测功效大于MQ, 特别地, Founders数越小, 此关系越明显; 当Founders数较小 (如10和50) 时, MQ的存在会显著影响其相关QQ检测功效。

2.2 大豆耐盐碱性上位性QTL检测

大豆品种群体各SSR标记的等位基因数目不等, 很多位点等位基因数远大于2, 其效应数为等位基因数减1, 由此QQ效应数远大于样本容量, 导致遗传模型超饱和。本研究采用两步筛选方法以解决该问题。首先, 利用独立性χ²测验筛选标记以剔除遗传模型中明显不显著的效应, 以降低模型维数; 其次, 使用分层广义混合线性模型上位性检测方法分析大豆幼苗根长耐盐碱性数据。

在效应初步筛选时, 需要考虑筛选的显著水平。本研究设置显著水平分别为0.10、0.20、0.25和0.30, 以获得最适显著水平。结果表明, 显著水平大于0.25时的结果比较接近不筛选变量的结果 (表 1), 且假阳性率小于0.8%。因此, 选择0.25作为变量筛选的概率标准。

2.2.1 MQ、QE和QQ的检测

用新方法分析大豆耐盐碱指数分级数据, 检测到2个MQ (1个耐碱ATI和1个耐盐STI)、4个QQ (3个ATI和1个STI) 和3个QE (1个ATI和2个STI), 分别位于A1、B2、D1a、D2、F、I、M、N和O共9条染色体上 (表 2)。从表 2可见:耐碱性指数MQ、QQ和QE的累计表型方差解释率 (proportion of phenotypic variance explained, 简记为PVE) 分别为3.15%、26.37%和45.98%。QQ标记sat_274×satt348和QE标记sat_344的PVE均大于10%。耐盐性指数MQ、QQ和QE的累计PVE分别为3.33%、13.27%和14.47%。

与Feng等^[22]的结果比较, 有2个相似的MQ和3个相似的QE。前者为耐碱性关联标记sat_342(PVE=3.15%) 和耐盐性关联标记sat_274(PVE=3.33%); 后者为耐碱性关联标记sat_344(PVE=45.99%)、耐盐性关联标记satt270(PVE=7.38%) 和耐盐性关联标记satt022(PVE=7.09%)。Feng等^[22]采用主效QTL检测方法检测到碱性上位性sat_274×satt348的2个标记与MQ关联; ATI的QQ和STI的MQ同时检测到标记sat_274。标记satt580、satt149和satt022被EAM检测到^[24]; 标记sat_342、satt348、satt413和satt270被ECMLM^[30]检测到^[24]; 标记satt255、satt348、satt354被mrMLM^[31]检测到; 标记sat_274和satt354被BhGLM^[7]检测到; 根据soybase信息, 标记satt255与耐盐QTL (Salt tolerance 5-1) 关联^[32]。

2.2.2 优异等位基因挖掘和亲本组合预测

通过极大似然法, 估计品种群体关联标记等位基因效应, 筛选出效应最小即耐盐碱性最强的等位基因, 作为该位点的优异等位基因, 而携带该等位基因的品种为代表性载体材料, 结果见表 2。从表 2可见:与ATI关联的标记sat_342等位基因260 bp的效应最小 (-0.79), 为优异等位基因, 含该等位基因的‘遵义棕子豆’为其代表性载体材料。ATI的QQ (sat_274×satt348) 中, sat_274的412 bp × satt348的232 bp效应最小 (-6.20), 代表性载体材料为高耐的‘鲁豆1号’; STI的QE (satt270) 中, 等位基因252 bp与2009年环境的互作效应最小 (-0.56), 其代表性载体材料为‘临安八月白’, 2年的表型鉴定至少为中等耐性。

从包含所有被检测到优异等位基因的杂交组合中选出性状预测值最小的RILs作为最优亲本组合。结果表明, 任何两品种的组合均无法包含所有被检测到QTL的优异等位基因。为了将耐盐性或耐碱性优异等位基因全包含在内, 可考虑利用多个品种组合, 以实现亲本组合预测。例如‘周01042-10’ב白秋1号’ב鲁豆2号’包含所有被检测到QTL的耐碱性优异等位基因, 可作为改善耐碱性的最优亲本组合; ‘荣经四坪黄豆’ב大悟大黄豆’ב水城棕子豆’ב大黑脐’包含了所有被检测到QTL的耐盐性优异等位基因, 可作为改善耐盐性的最优亲本组合。

2.2.3 候选基因预测

分析了大豆耐盐碱性的潜在候选基因信息, 共发现12个与拟南芥耐盐碱性基因同源的大豆基因。它们分别与本文检测到的12个标记关联, 其中2个MQ、4个QQ和3个QE, 其物理距离在24.62~2 235.39 kb之间 (表 3)。发掘的3个与拟南芥耐盐碱基因AT1G69270共线性的大豆基因为Glyma01g31590、Glyma03g32260和Glyma13g04890, Glyma13g04890位于大豆F连锁群上标记satt149附近, 两者相距约175.71 kb。

3 讨论

本研究将品种群体多歧性状的多QTL分层广义混合线性模型关联分析方法^[22]推广到了上位性检测。新方法可以同时检测多歧性状的MQ、QQ和QE。在实际应用中, 面对超饱和线性模型实际数据资料, 可以结合卡方测验进行两步筛选, 以实现上位性检测。利用极大化正态逼近的似然和经验贝叶斯估计各关联标记效应, 并实现优异等位基因挖掘和优异亲本组合预测。本研究结果的可靠性已经在文中得以证实。

Monte Carlo模拟试验结果表明, 新方法在统计检测功效和参数估计精确性方面表现较好^[29], 且假阳性率比较低; 遗传率、表型分布、表型分类数、Non-founders数和Founders数对检测结果有一定影响, 这与Lü等^[17]和Feng等^[22]的结果相似。

新方法检测了257个大豆品种耐盐碱性的上位性, 除了标记sat_391, 其他关联标记均可被多QTL分层广义混合线性主效模型、EAM、ECMLM、mrMLM和BhGLM中的至少一种方法所证实。将本研究结果与该257份大豆连续数量性状的EAM和ECMLM定位结果^[24]对比, 共检测到7个相同标记; 2个MQ和3个QE被Feng等^[22]的主效多QTL方法检测到; 耐碱上位性的2个关联标记sat_274和satt348被Feng等^[22]的主效多QTL方法检测到为耐碱性MQ。

耐盐碱性优异等位基因代表性载体材料预测结果, 与所选材料的耐盐碱性试验表型鉴定结果一致。例如, 检测到耐碱性的3个QQ的载体材料:‘周01042-10’‘鸭蛋青’和‘鲁豆1号’至少在2009和2010两年的表型鉴定中为耐性; 检测到耐盐性MQ的载体材料‘荣经四坪黄豆’在两年表现鉴定中均为耐盐性, QQ的载体材料‘大悟大黄豆’至少在两年的表型鉴定中为耐盐性。

由于作物很多重要性状表现为多歧性状, 其表型数据是离散的。因此, 新方法的提出具有应用价值。目前, 关联分析主要是检测主效QTL, 上位性的检测方法还相对较少。因此, 新方法的提出对于推动植物关联分析上位性检测具有积极作用。随着标记数目的增加, 遗传模型中待估效应个数远大于样本容量, 模型的超饱和性越严重。因此, 品种群体多歧性状上位性QTL检测依然面临很大的挑战, 需要继续探索。