水稻(Oryza sativa)为全世界超过50%的人口提供了稳定的食物来源，水稻增产对保障我国乃至全球粮食安全做出卓越的贡献。20世纪60、70年代矮秆品种和杂交稻的培育和应用，使我国水稻产量实现了两次飞跃；但大量高产品种的培育和大面积推广，引发化肥、农药、水以及劳动力的投入大幅增长，导致水稻生产与资源环境的矛盾激增^[1]。针对我国水稻生产中病虫害严重、农药化肥使用过量、水资源短缺、产量持续徘徊不前等问题，张启发院士提出了培育“少打农药、少施化肥、节水抗旱、优质高产”的绿色超级稻育种目标^[2]。绿色超级稻的绿色性状，如抗病、抗旱、抗逆、氮磷高效利用等，大多是受多基因控制的复杂性状，受环境影响较大，常规选育的效果不太理想。随着高通量测序技术的发展，关联分析和基因组选择已成为对植物复杂数量性状进行遗传解析和改良的重要方法。

在关联分析和基因组选择模型中，标记的数量经常远远超过样本量，会产生自由度不足和多重共线性等问题。为了解决这些问题，前人发展出了许多统计方法^[3-4]。然而，目前的关联分析和基因组选择方法仍专注于对单一性状的分析，忽略了性状之间的遗传相关和环境相关。在绿色超级稻的育种工作中，需要培育同时具有高产优质、节水抗旱、抗病抗逆等绿色性状的品种，而不是追求单一目标性状的改良。为了满足绿色超级稻的育种需求，本文提出了关联分析和基因组选择的多性状方法。

1 关联分析的多性状方法

关联分析，是一种基于连锁不平衡，鉴定某一群体内目标性状与遗传标记或候选基因关系的分析方法。与传统的连锁作图相比，关联分析利用的是自然群体在长期进化中所积累的重组信息，解析精度较高，可实现对数量性状位点(quantitative trait locus, QTL)的精细定位，是解析数量性状的强有力的工具。自2001年以来，应用关联分析方法发掘植物数量性状基因就备受关注^[5]，大量的统计方法以及软件包应运而生。

目前应用最广泛的是基于单位点扫描的混合线性模型方法，最早由Yu等^[6]提出。该方法能够有效控制群体结构和多基因背景，从而降低关联分析的假阳性。为了处理较大样本并且提高运算效率，一系列基于混合线性模型的改进方法相继被提出，如EMMA^[7]、GEMMA^[8]等。为了获得更快的运算速度，一些近似算法也被提出，如EMMAX ^[9]、P3D ^[10]、CMLM ^[10]。虽然已被广泛应用，但单位点扫描方法仍存在一定的局限性，如要进行多重测验矫正，不能同时考虑所有位点信息等。因此，多位点的关联分析方法逐渐受到关注。在进行多位点关联分析时，样本量远大于标记数目(n)，许多方法被提出处理此类“大p小n”的问题，如逐步线性回归方法、岭回归分析、LASSO、主成分回归、偏最小二乘法(partial least squares, PLS)及贝叶斯方法等。为了更有效地控制群体结构，基于混合线性模型的多位点方法也逐渐被提出，如MLMM^[11]、LMM-Lasso^[12]和BSLMM^[13]等。

然而，上述的关联分析方法仍停留在单个目标性状上，而在绿色超级稻的育种工作中，需要培育出同时具备多个绿色性状的新品种。多性状关联分析既可以同时对多个目标性状进行联合分析，实现多个目标性状的协调发展，也可利用性状之间的遗传相关和环境相关信息，提高关联分析的功效。要想获取绿色性状基因，必须准确鉴定出与目标性状关联的基因变异位点，并精确估计出该位点对目标性状的遗传效应和贡献率。因此，本文拟介绍一种基于偏最小二乘的多性状联合关联分析方法，并利用模拟实验验证该方法的可行性。

1.1 分析方法

本文采用以下策略进行基于偏最小二乘的多性状联合关联分析，步骤如下。(1)采用非线性迭代偏最小二乘法求解多变量模型的回归系数。(2)计算变量投影重要性(variable importance in projection, VIP)，分别计算M个自变量的VIP指标并从大到小排序，即：${\rm{VI}}{{\rm{P}}_{\left(1 \right)}} \ge {\rm{VI}}{{\rm{P}}_{\left(2 \right)}} \ge {\rm{VI}}{{\rm{P}}_{\left(M \right)}}$。(3)根据回归系数选择对应每个表型性状的重要变量$S_k^v$，然后根据VIP指标同样选择重要变量子集${S_{{\rm{VIP}}}}$，均选择前$\frac{n}{2}-1$(n为偶)或$\frac{{n-1}}{2}$(n为奇)个效应变量，从而综合对各个性状都相对重要的自变量集合。整合两类重要效应变量集合，生成对应每个表型性状的重要效应变量的预选择${S_k} = S_k^V \cup {S_{VIP}}$。其中S_k中的变量个数小于n–1，经过变量压缩，线性模型由过饱和变为普通的线性模型。(4)采用双向筛选的逐步回归进行重要效应变量的选择，并用BIC信息准则作为增加或删除变量的标准。最终根据最优模型中的效应变量，以普通线性回归模型计算效应大小及相应的p值，具体算法参见[14]。

1.2 模拟研究

为了探究该方法的可行性，本团队在不同的条件下比较了该方法和两个单性状的多位点方法LASSO^[15]和基于偏最小二乘的多位点关联分析(PLS-based MLAS)^[16]的检验功效。本文以两个相关数量性状的联合分析为例进行。假设基因组中存在10 000个位点，每个位点均有2个等位基因，其中10个位点具有遗传效应，分别为Q1~Q10。Q1~Q4仅控制性状1的表达，Q7~Q10仅控制性状2的表达，Q5和Q6设为一因多效基因；为了更加符合实际情况，对不同位点设置不同水平的多态信息含量(PIC)。10个位点的效应值及PIC设置如表 1所示。同时，考察样本容量n和候选基因遗传力h₂两个因素，其中样本容量设定为3个水平，分别为100、300和500；遗传力设定3个水平，分别为30%、50%和70%。性状1和性状2的平均值均设为10，两个性状间的剩余误差相关系数r_e设定为0.5。

结果表明，样本容量和遗传力均对检测功效具有明显影响(图 1)。当样本容量一定时，随着遗传力的增加，检测功效逐步提高；与此类似，当遗传力一定时，随着样本容量的增加，检测功效也相应提高。从表 2可以看出，PIC值对其检测功效具有较大的影响：位点的PIC值越大，其被检测到的概率越高；PIC值较小的位点在样本较少和遗传力较低时被检测到的概率较低，甚至检测不到。由表 2还可以看出，效应的大小对检测功效有较大影响：在PIC保持恒定的情况下，QTL的效应值的绝对值愈大，其被检测到的概率愈高，然而效应真值的正负对检测功效影响不大。样本容量和遗传力对效应估计的准确度和精确度均具有一定影响，随着样本容量和遗传力的增加，位点效应的估计值越趋近于真值。

采用多性状联合分析、LASSO以及PLS-based MLAS方法获得的统计功效均列于表 2，从中可以看出，多性状联合分析时，一因多效位点的统计功效与单性状分析相比有了大幅度的提高。此外，大多数情况下，即使一个位点仅控制一个性状，联合分析仍具有一定优势。但是对于PIC值极低的位点，联合分析似乎没有优势。

2 基因组选择的多性状方法

基因组选择(genomic selection, GS)能够估计全基因组上所有标记的遗传效应，进而实现对品种更加可靠的选择^[17-18]。特别是在水稻等作物的杂种育种中，杂交种的基因型可以由亲本基因型进行推断，GS的优势更加突出^[19-20]。然而传统的GS方法专注于单一环境下单个性状的预测，忽视了性状之间的遗传相关和不同环境之间的关联，而相互关联的性状或者不同环境下的作物品种之间可能拥有共同的遗传或生物学基础^[21]，从而使多性状联合分析拥有更高的预测能力^[22]。在绿色超级稻的育种实践中，一些性状可能难以度量或者观测的成本过于昂贵(小区产量、抗逆性状、肥料的吸收利用和根系相关性状等)，多性状联合分析提供了一种新的策略，即利用更易鉴定的相关表型预测难以获取的重要性状^[23]。

此外，水稻等作物的产量遗传力较低，所以基因组预测的效果不佳。如Xu等^[24]利用GBLUP模型预测了水稻杂交种的4个性状，其中产量的预测能力最低，只有0.13。另一方面，育种家真正关心的是水稻等作物多个农艺性状的综合表现，千粒重等农艺性状的遗传力较高，相应的基因组预测精度也较高。绿色超级稻的育种目标有多个，要求“少打农药、少施化肥、节水抗旱、优质高产”等，实际的GS工作中可以利用一些较高遗传力的性状构建选择指数，实现对品种的综合预测。而一般的多性状预测方法往往是基于多变量模型的建立，未能与选择指数理论^[25]相结合，从而无法利用与目标性状相关的多个其他性状构建选择指数，对作物进行更加全面的选择。

选择指数是农作物多目标育种选择的常用方法，其效率高于对所有性状的逐一选择，能够被用来同时改良多个数量性状。GS的快速发展为选择指数带来了新的前景。Dekkers ^[26]利用全基因组高密度标记估计的育种值构建选择指数；Jesus Ceron-Rojas等^[27]利用GS方法和选择指数对多个性状同时进行选择；Schulthess等^[28]使用黑麦中的两个性状建立选择指数，并将其看做单一性状用GS法进行预测；Lyra等^[29]将玉米杂交种在不同氮胁迫下的性状组合以构建选择指数，然后用GS方法进行预测，结果表明其方法是有效的。但是，已有的选择指数研究大多专注于对动植物多个目标性状的加权选择，在将选择指数和GS结合时，指数本身的预测精度往往成为关注的对象，而忽视了利用选择指数聚集辅助性状信息以预测目标性状的能力。在对复杂目标性状进行预测的过程中，与目标性状相关的辅助性状能够提供大量有用的信息，将选择指数和GS方法相结合，有助于充分利用这些信息，帮助育种家对目标性状进行更加准确的预测，从而有效开展选种和配种工作。

因此，本研究开展了多性状联合预测的GBLUP模型研究，并利用与目标性状相关的多个辅助性状建立选择指数，实现对水稻性状的综合选择，提高预测精度。本研究使用的水稻数据集来自武汉大学^[20]，以115份水稻纯系品种为父本，5份不育系作母本，采用NCII交配设计配制575份杂交种。考察了8个性状，包括单株产量(grain yield per plant, GY)、千粒重(thousand-grain weight, TGW)、有效穗(productive panicle number per plant, PN)、株高(plant height, PH)、一次枝梗(primary branch number, PB)、二次枝梗(secondary branch number, SB)、主穗实粒数(grain number per panicle, GN)和穗长(panicle length, PL)。本研究模型中使用到的水稻表型数据是两种环境下两次重复的平均值。同时，对120份水稻亲本进行高通量测序，共获得3 299 150个SNP标记。

2.1 基因组选择方法的比较

自从Meuwissen等^[30]首次提出GS的概念后，大量的方法涌现，包括各种参数、半参数和非参数的方法。参数方法主要有GBLUP、LASSO、PLS以及贝叶斯方法，非参数方法有随机森林、SVM和RKHS等。比较各种GS方法的异同，了解每种方法所适用的条件，以更好地服务于育种工作，是GS研究中的重要内容。已有一些研究对这些方法进行了比较。de los Campos等^[31]对参数方法的预测准确性进行了比较，发现对于大多数性状GBLUP都有较好的预测准确性，BayesB对于大效应QTL控制的性状预测效果较好。Riedelsheimer等^[32]用5种不同预测方法对玉米3个农艺性状和3个代谢性状进行了预测，得出的结论是这些方法在预测准确性上差异不大。Heslot等^[33]利用了10种不同的参数和非参数方法对从不同物种中测量的18个性状进了预测准确性比较，发现RKHS在不同的性状和物种中都有较好的预测表现。Howard等^[34]利用模拟数据进行了参数方法和非参数方法的预测准确性比较，发现参数方法预测受加性效应控制的性状时效果比非参数方法要稍好一些。然而在实际水稻数据中，这些基因组预测方法的比较尚少有报道。

本课题组基于上述杂交水稻数据集，对6种常用的预测方法进行比较^[35]，包括4种参数方法LASSO、GBLUP、BayesB和PLS以及2种非参数方法RKHS和SVM。通过对6种预测方法的预测准确性进行方差分析发现，不同方法的预测准确性之间存在极显著差异，而且预测准确性在方法、性状之间的互作也存在极显著差异。进一步进行多重比较发现，6种方法中，GBLUP方法的预测准确性最高，而PLS和SVM的预测准确性最低，其余方法的准确性介于两者之间。此外，不同的方法适用于预测不用的性状，如GBLUP方法对于主穗实粒数、株高和千粒重的预测准确性最高，LASSO对于一次枝梗和二次枝梗的预测准确性最高，SVM是穗长的最佳预测方法。尽管没有一种方法能适用于所有性状，但GBLUP方法在所有性状的预测中表现最稳健。因此，后续的多性状的GS研究将基于GBLUP模型开展。

2.2 基因组选择模型的多变量方法

传统的GS方法主要用于单一环境下单个性状的预测。虽然迄今为止已经发表的大部分GS研究结果使用的都是单变量模型，但是一些研究建议使用多变量模型^{[28, 36]}。本课题组将单性状加-显模型扩展为多性状加-显模型，并且把非遗传的残差项分解为完全相关的部分和完全独立的部分，得到了多性状加-显模型MV-AD以及包含共同环境效应的多性状加-显模型MV-ADE。此外，本研究利用辅助变量构造的关系矩阵开发了一种高效的多变量模型MV-ADV^[20]。各种模型的预测能力使用交叉验证的方法衡量。

与单性状模型相比，多性状模型能够利用到群体更多的性状信息。以GY和TGW分别为目标性状，使用MV-AD、MV-ADE和MV-ADV对水稻杂交种进行的两性状联合预测结果见表 3。对于GY，MV-ADE和MV-ADV的预测能力明显优于MV-AD和UV-AD。虽然MV-ADV相较MV-ADE的优势绝对值较小，但是成对比较表明MV-ADV的优势是显著的。然而无论以哪个性状为辅，都无法提高对TGW的预测能力。表明多性状模型能够提高对GY等低遗传力性状的预测能力，对于TGW等高遗传力的性状，多性状模型的预测效果并无明显改进，此时只需应用单性状模型进行预测即可。

本研究不仅考察了两性状联合预测的情况，还将所有的8个性状放在一起进行了联合预测。图 2对比了单性状UV-AD、两性状MV-ADV和八性状MV-ADV的预测能力。其中八性状联合的预测能力显著高于单性状预测和两性状预测，特别是GY和PN的预测能力提高最大，TGW和PH的预测能力提高较少，再次说明多性状联合分析的策略更适用于低遗传力的性状。两性状的平均预测能力比单性状预测要高10.2%，而八性状的平均预测能力比单性状预测要高50.6%。很明显，八性状的预测能力相对于两性状和单性状预测获得了更为显著的增加，说明结合更多的辅助性状信息，可以较大程度提高目标性状的预测能力。

2.3 综合性状的基因组选择方法

本课题组以杂交种的SNP数据为基础，针对基于选择指数的多性状GS方法开展了大量模拟研究，随机模拟了具有100个QTL的多个性状。其中遗传力为0.3的表型性状T3为待预测的目标性状，遗传力为0.7和0.4的多个性状为辅助性状，用以构建选择指数。本研究同时以模拟表型和实际表型数据为例说明多性状选择指数GS方法的应用，并对其效果进行评估。

本研究针对基于选择指数的多性状GS方法，选定某一目标性状(如GY)，对与之相关的s个辅助性状构建了选择指数SI。选择指数能够综合多个辅助性状与目标性状的相关遗传信息，具有较强的预测价值，能够实现对水稻多个性状的综合选择。此外，如果预测的目标是产量等单一性状，可以将选择指数的预测结果与传统GS的预测结果相结合，以获得更高的精度。

为了衡量利用选择指数进行预测的效果，本研究定义预测值与遗传效应值之间的决定系数为预测精度，与表型值之间的决定系数为预测能力。研究中使用5倍交叉验证的方法衡量各种模型的预测效果。在模拟研究中，每次以T3为目标性状，可选取不同组合的其他性状为辅助性状，然后建立选择指数。本研究对比11个不同的辅助性状组合，仅利用训练集数据构造选择指数SI_train。指数预测能力反映了选择指数本身所具备的预测能力，在模拟的情形下，11个组合的平均预测能力为0.7035，远高于T3预测能力(0.2470)，说明利用选择指数进行预测能实现更加精确的综合选择。而且，辅助性状组合与目标性状遗传相关程度越高，指数预测能力越高。此外，辅助性状的遗传力对指数预测能力也存在较大影响。因为辅助性状的遗传力越高，其中包含的遗传信息比重越大，以其为基础构建的选择指数也必然获得更高的遗传力，从而拥有较高的指数预测能力。水稻实际表型数据的研究中，指数预测能力的平均值为0.5479，也远远高于GY的预测能力(0.1344)。

当目标性状缺失时，可以用构建的选择指数直接预测目标性状。若选取不同的辅助性状组合构建选择指数，并直接利用指数预测T3，在大多数情况下，预测的精度无法超越目标性状的预测精度，但是可以十分接近这一水平，这在目标性状缺失的情况下不失为一种有效的方法。但是，如果已经拥有了目标性状的训练群体，就没有必要直接用选择指数进行预测。选择指数能够综合多个辅助性状与目标性状的相关遗传信息，可以将选择指数的预测结果与传统GS的预测结果相结合，提高对目标性状的预测精度和预测能力。

对本研究选取的多个辅助性状组合，其预测结果见表 4。虽然大部分组合的指数直接预测精度低于T3预测精度，但是，所有组合的指数辅助预测精度都优于T3预测精度，提高的百分比为0.6%~ 8.3%。同时，辅助性状与T3的遗传相关越高，指数辅助预测精度越高。指数辅助预测能力也有类似的特点，只是由于误差项的干扰而波动较大。

除了模拟的情形，本课题组还以水稻杂交种的GY为目标性状，TGW、PN、PH、PB、SB、GN和PL为辅助性状，考察了选择指数辅助预测的情况。因为实际性状的育种值是未知的，所以无法得到真实的预测精度，只能用指数辅助预测能力来衡量指数的优劣。为了找到最佳的辅助性状组合，本研究采用逐步剔除的方法，从全部辅助性状组合出发，每次去除一个辅助性状，然后将去除前和去除后的所有组合进行对比，找到指数辅助预测能力最大的组合。结果表明，各组合所构建选择指数的辅助预测能力均高于GY的预测能力，说明指数辅助预测是有效的。其中，TGW、PN、PH、PB和SB等5个性状构建的选择指数辅助预测能力为0.1461，是所有组合中最大的，比GY的预测能力高8.7%。

3 结语与展望

针对绿色超级稻育种实践中多个绿色性状的遗传改良问题，本文提出了关联分析和基因组选择的多性状方法。

首先，基于多元偏最小二乘和两阶段变量选择策略，提出了用于优异等位基因挖掘和遗传效应估计的多性状关联分析方法，并通过大量的模拟研究验证了该方法的可行性。研究表明，随着样本容量、候选基因遗传力以及PIC的增加，检测功效以及效应估计的准确度和精确度均明显提高，由于本文提出的方法充分利用了性状之间的遗传相关与环境相关信息，并且对所有的基因位点同时分析，因此在进行关联分析时具有明显的优势。与单性状分析相比，不论候选基因是同时控制多个性状的表达，还是仅控制其中一个性状的表达，绝大多数情况下具有更高的检测功效和更准确的效应估计。对于一因多效基因，多性状联合分析的优势更为明显，这与Xiao等^[37]对主基因的分离分析、Jiang等^[38]对多性状联合QTL定位的研究结论一致。然而对于PIC值很小的位点，例如Q1，基于偏最小二乘的多性状联合分析的检测功效较单性状分析差，这可能是偏最小二乘对成分的提取损失了部分遗传变异方差信息的缘故。本文提出的模型目前仅适用于群体结构不明显的自然群体，下一步的工作设想是对模型作进一步改进，考虑消除群体结构对分析结果的影响，将群体结构纳入协变量，使之适用于任何的自然群体，进一步提高模型的适用范围。此外，本研究的模型仅涉及到候选基因的主效应，该分析方法可以进一步扩展到主效应加互作效应的模型，模型的构建可以参考Zeng等^[39]的方法进行。

其次，本文提出的多性状基因组选择方法，特别有利于对一些低遗传力或者难以观测的重要性状进行选择^[40]。利用这一策略能够提高预测的精度，特别是性状之间存在较强的关联时^[41]。在绿色超级稻育种中，很多绿色性状(如氮磷高效利用、抗逆性状和小区产量)较难测定，此时利用目标群体的其他辅助性状进行预测具有十分重要的现实意义。但是如果所预测的群体未得到种植，即没有任何性状得到鉴定时，可以使用基于选择指数的GS方法进行预测。另一方面，传统的水稻育种关注产量水平的提高，对绿色性状的重视不够。目前绿色超级稻育种要求多个绿色性状同时得到发展，这就给全基因组选择带来新的挑战。基于选择指数的GS技术适应了这一需要，它能实现对水稻多个性状的同时选择。不论是模拟表型数据还是实际表型数据，选择指数的预测能力都大大高于传统GS的预测能力，充分说明选择指数具有更加精确的综合选择能力，对于绿色超级稻等项目的多目标育种具有重要的意义。此外，基于选择指数的GS方法能够利用与目标性状相关的多个辅助性状建立选择指数，实现对某个目标性状的辅助预测，提高对目标性状的预测精度和预测能力。本研究中选择指数的建立都是基于不同性状间的线性关系，而实际中这种关系可能是非线性的，在今后的研究中考虑性状间更加复杂的非线性联系，可进一步提高选择指数预测的效果。对于水稻杂交种实际数据选择指数的研究，虽然表型性状局限于GY、TGW等性状，分析结果对水稻其他性状的应用可能有一些限制，但是本研究的方法和结果为水稻绿色性状选择指数的构建，以及利用选择指数预测目标性状，提供了重要的技术参考。