一种基于全基因组关联分析和全基因组选择选育玉米高产材料的方法

摘要

本发明提供一种是基于全基因组关联分析和全基因组选择建立快速筛选玉米优良自交系的分子育种方法。本发明通过3种全基因组关联分析方法即FarmCPU、MLMM和CMLM计算与玉米单穗产量、穗粒数和百粒重关联SNP的显著性。通过比较不同模型预测的准确性，确定最佳的预测模型。利用最佳模型对不同标记密度选择方案比较确定最佳的SNP信息。根据最佳预测模型和SNP信息，筛选单穗产量、穗粒数和百粒重育种值较高的优良材料。本发明可以快速鉴定候选材料的遗传变异，快速选育出聚合有利基因的优良材料，缩短育种进程。

说明书

技术领域

本发明属于植物分子育种领域，具体涉及一种基于全基因组关联分析和全基因组选择选育玉米高产材料的方法。

背景技术

玉米是重要的粮食、饲料和工业原料。单穗产量的提高是玉米育种的重要目标。单产是一个非常复杂的数量性状，主要受微效多基因控制。穗粒数和粒重是单穗产量的主要构成因子，也是重要的育种目标。

通过传统定位和全基因组关联分析，已经挖掘了大量控制单穗产量、穗粒数和百粒重的关键位点。然而这些位点尚停留在研究层面，并没有应用到育种中。由于单穗产量、穗粒数和百粒重主要受微效基因控制。利用分子标记辅助选择对这些性状进行遗传改良收效甚微。全基因组选择能够对所有遗传变异和遗传效应进行准确检测和估计，能够有效地利用微效位点。相比分子标记辅助选择，全基因组选择预测玉米散粉吐丝间隔、散粉期和吐丝期的准确性提高7％-43％。在产量和秸秆产量的遗传进度上，全基因组选择选育的玉米株系比分子标记轮回选择的高14％-50％。其选择株系的产量测交表现比传统系谱选择方法选择到的可以提高7.3％。因此，全基因组选择是利用单穗产量、穗粒数和百粒重等复杂数量性状的有效方法。

利用全基因组关联分析鉴定现有种质资源中的有利等位基因，挖掘关键位点，为玉米单穗产量遗传改良提供位点信息。在全基因组关联分析的基础上，选择与目标性状关联的位点，利用全基因组选择预测候选材料的遗传变异，能够快速选育出聚合有利基因的优良材料，创制高产新种质和新材料，从而缩短育种进程。

发明内容

本发明的目的在于提供一种基于全基因组关联分析和全基因组选择选育玉米高产材料的方法。

本发明的目的可以通过以下技术方案实现：

一种基于全基因组关联分析和全基因组选择选育玉米高产材料的方法，该方法包括以下步骤：

(1)材料和田间试验设计

在至少一个种植环境种植多份玉米种质自交系，种植过程中采用随机区组试验设计，每个小区选取授粉较好的穗子，测量玉米穗子的单穗产量、百粒重和穗粒数；

(2)表型数据统计分析

利用QTL IciMapping v4.0的AOV功能对各种植环境进行联合方差分析，并计算单穗产量、穗粒数和百粒重的广义遗传力和最佳线性无偏估计值即BLUE值；

(3)基因型鉴定和分析

采用GBS简化测序的方法对多份自交系进行基因型分型，利用BWA软件比对到参考基因组，采用SAMTOOLS软件进行群体SNP的检测，获得基因型数据的VCF；利用plink v1.9对基因型VCF筛选，筛选标准为缺失率小于0.10，杂合率小于0.10，最小等位基因频率(MAF)大于0.05；获得多个高质量SNP的用于关联分析；

(4)全基因组关联分析

利用各种植环境和BLUE环境的单穗产量、百粒重和穗粒数进行全基因组关联分析，根据Bonferroni矫正设置的显著性阈值(P<1.72E-05)，利用FarmCPU、MLMM和CMLM方法分别检测到单穗产量、穗粒数和百粒重显著关联的SNP；

详细过程为：利用各种植环境和BLUE环境的单穗产量、百粒重和穗粒数进行全基因组关联分析，筛选与单穗产量、百粒重和穗粒数显著关联的SNP；采用的方法为CMLM、MLMM和FarmCPU；3种方法均采用Q(群体结构)+K(亲缘关系)模型；Q值利用Structure v2.3.4计算，设置群体内亚群数为1～8，length of burn-in period设置为5000，蒙特卡罗重复个数设置为50000，每个亚群数迭代次数为3；根据ΔK，确定亚群数为2时的Q值用于关联分析，K值利用TASSEL v5.0的Centered_IBS方法计算，根据Bonferroni矫正方法，设置显著临界值设置为P＝1/58129＝1.72E-05；CMLM方法的PVE(表型变异解释率)由软件给出；MLMM和FarmCPU方法的PVE由线性回归模型计算，公式为：Y＝α+βX+ε，其中Y为观测值，α为截距，β为斜率，X为标记编码(2，0，1)，ε为随机误差；回归模型计算的

(5)最佳模型的确定

由于基因型和表型数据均有缺失，在做全基因组预测时需要对其进行缺失填补；根据已知基因型的基因型频率，随机选择标记进行填补；对填补后的SNP基因型进行编码，高频等位基因纯合型编码为2，低频等位基因型纯合编码为0，杂合基因型编码为1；填补和编码均在R语言中实现；

表型缺失值利用R语言mice包的pmm(预测均值匹配)方法进行多重插补，多重填补法的填补矩阵数设为5，迭代次数为50次；

全基因组选择模型采用BayesA、BayesB、BayesC、Bayes LASSO、GBLUP和RKHS；

通过比较BayesA、BayesB、BayesC、Bayes LASSO、GBLUP和RKHS模型预测的准确性，确定最佳的预测模型；

(6)标记密度选择标准和优良材料的筛选标准

利用最佳模型比较三种标记选择方案不同SNP密度的预测准确性，确定最佳的SNP信息，根据最佳预测模型和SNP信息，筛选单穗产量、穗粒数和百粒重育种值较高的优良材料。

三种标记选择方案：第一种标记选择方案是随机选择500个SNP、1000个SNP、5000个SNP、10000个SNP、20000个SNP和40000个SNP。第二种标记选择方案是根据与性状关联性即P值选择同等数量的SNP；该方案中，利用FarmCPU模型BLUE值计算的P值从小到大排序，选择最显著的500个SNP、1000个SNP、5000个SNP、10000个SNP、20000个SNP和40000个SNP。第一种方案也利用BLUE值作为真实观测值。因此，第一种和第二种方案预测准确性为BLUE与育种值间的相关系数均值。第三种标记选择方案是在3种全基因组关联分析方法下，将Bonferroni矫正阈值选择的与目标性状最显著的SNP为标记集，对各种植环境分别进行预测，评价性状显著关联SNP应用的价值。第三种方案预测准确性为各种植环境与育种值间的相关系数均值。利用1.5确定的最佳模型对3种SNP选择方案进行比较，确定预测效果最好的SNP信息。根据最佳预测体系估计的育种值，筛选单穗产量、穗粒数和百粒重育种值较高的优良材料。

作为一种详细技术方案，步骤(5)中所述的RKHS采用的是多核模型即核平均

优选的，步骤(4)中根据Bonferroni矫正方法，设置的显著性阈值为1.72E-05，在该阈值下，利用FarmCPU、CMLM和MLMM方法检测与单穗产量、百粒重和穗粒数显著关联的SNP。

优选的，步骤(6)中根据最佳模型和性状关联SNP信息，利用BGLR的cbind(fm$yHat,y)函数分别调出单穗产量5000个最显著SNP估计的育种值以及穗粒数和百粒重1000个最显著SNP估计的育种值结果，分别选择育种值最高的前10个材料并进行综合评价从中选择高产材料。

优选的，步骤(1)中所述种植环境的数量为1～6个环境。

优选的，步骤(1)中所述随机区组试验设计的方法为：设3次重复，小区为2行区，行距60cm，株距25cm；每个小区，选取授粉较好的3个穗子。

本发明通过3种全基因组关联分析方法即FarmCPU、MLMM和CMLM计算与玉米单穗产量、穗粒数和百粒重关联SNP的显著性。通过比较不同模型预测的准确性，确定最佳的预测模型。利用最佳模型对不同标记密度选择方案比较确定最佳的SNP信息。根据最佳预测模型和SNP信息，筛选单穗产量、穗粒数和百粒重育种值较高的优良材料。

本发明的有益效果：

本发明设计一种是基于全基因组关联分析和全基因组选择建立快速筛选玉米优良自交系的分子育种方法。本发明方法可以快速鉴定候选材料的遗传变异，快速选育出聚合有利基因的优良材料，缩短育种进程。

附图说明

图1为BLUE环境下单穗产量、穗粒数和百粒重不同模型的预测准确性。

图2为单穗产量、穗粒数和百粒重随机选择和性状-SNP显著性选择的不同标记密度的预测准确性。

图3为Bonferroni矫正(P<1.72E-05)选择的性状显著SNP利用RKHS方法预测不同环境单穗产量、穗粒数和百粒重的准确性。

具体实施方式

实施例1

1材料与方法

1.1材料和田间试验设计

关联群体包括黄淮海骨干自交系、国内核心种质和美国GEM等种质309份。2017年种植在郸城、虞城、原阳和三亚。2019年仅种植在原阳。采用随机区组试验设计，3次重复。小区为2行区，行距60cm，株距25cm。每个小区，选取授粉较好的3个穗子，测量单穗产量、百粒重和穗粒数。

1.2表型数据统计分析

利用QTL IciMapping v4.0的AOV对2017郸城、2017虞城、2017原阳、2019原阳和2017三亚进行联合方差分析，并计算单穗产量、穗粒数和百粒重的广义遗传力和最佳线性无偏估计值即BLUE值。

1.3基因型鉴定和分析

采用GBS(genotyping by sequencing)简化测序的方法对309份自交系进行基因型分型，测序仪为Illumina HiSeq PE150双端测序。利用BWA软件比对到参考基因组(ftp://ftp.ensemblgenomes.org/pub/plants/release-36/fasta/zea_mays/dna/Zea_mays.AGPv4.dna.toplevel.fa.gz)。采用SAMTOOLS软件进行群体SNP的检测，获得基因型数据的VCF。利用plink v1.9对基因型VCF筛选，筛选标准为缺失率小于0.10，杂合率小于0.10，最小等位基因频率(MAF)大于0.05。共获得58129个高质量SNP的用于关联分析。

1.4全基因组关联分析

利用郸城2017、虞城2017、原阳2017、原阳2019、三亚2017和BLUE环境的单穗粒重、百粒重和穗粒数进行全基因组关联分析。采用的方法为CMLM(compressed mixed linearmodel)、MLMM(multiple loci mixed model)和FarmCPU(fixed and random modelcirculating probability unification)。3种方法均采用Q(群体结构)+K(亲缘关系)模型。Q值利用Structure v2.3.4计算。利用plink将58129个SNP的VCF利用--recodestructure转化成Structure软件需要的格式。导入Structure软件后，首先设置群体内亚群数为1～8，length of burn-in period设置为5000，蒙特卡罗重复个数设置为50000，每个亚群数迭代次数为3。根据ΔK，确定亚群数为2时的Q值用于关联分析。K值利用TASSEL v5.0的Centered_IBS方法计算。根据Bonferroni矫正方法，设置性状与标记显著临界值为P＝1/58129＝1.72E-05。

CMLM方法的表型变异解释率由软件给出；MLMM和FarmCPU方法的PVE由线性回归模型计算，公式为：Y＝α+βX+ε，其中Y为观测值，α为截距，β为斜率，X为标记编码(高频纯合等位基因型编码为2，低频纯合等位基因型编码为0，杂合基因型编码为1)，ε为随机误差；回归模型计算的

1.5最佳模型的确定

由于基因型和表型数据均有缺失，在做全基因组预测时需要对其进行缺失填补。根据已知基因型的基因型频率，随机选择标记进行填补。对填补后的SNP基因型进行编码，高频纯合等位基因型编码为2，低频纯合等位基因型编码为0，杂合基因型编码为1。填补和编码均在R语言中实现。

表型缺失值利用R语言mice包的pmm(预测均值匹配)方法进行多重插补，多重填补法的填补矩阵数设为5，迭代次数为50次。

全基因组选择模型采用BayesA、BayesB、BayesC、Bayes LASSO、GBLUP和RKHS(再生核希尔伯特空间，reproducing kernel Hilbert space)。通过比较BayesA、BayesB、BayesC、Bayes LASSO、GBLUP和RKHS模型预测的准确性，确定最佳的预测模型。

其中RKHS采用的是多核模型即核平均

1.6标记密度选择标准和优良材料的筛选标准

设置三种选择标记的方案。第一是随机选择500个SNP、1000个SNP、5000个SNP、10000个SNP、20000个SNP和40000个SNP。第二是根据与性状关联性即P值选择同等数量的SNP。第二个方案中，利用FarmCPU模型BLUE值计算的P值从小到大排序，选择最显著的500个SNP、1000个SNP、5000个SNP、10000个SNP、20000个SNP和40000个SNP。第一种方案也利用BLUE值作为真实观测值。因此，第一种和第二种方案预测准确性为BLUE与育种值间的相关系数均值。第三种方案是在3种全基因组关联分析方法下，将Bonferroni矫正阈值选择的与目标性状最显著的SNP为标记集，对BLUE、2017郸城、2017虞城、2017三亚、2017原阳和2019原阳分别进行预测，评价性状显著关联SNP应用的价值。第三种方案预测准确性为不同环境与育种值间的相关系数均值。利用1.5确定的最佳模型对3种SNP选择方案进行比较，确定预测效果最好的SNP信息。根据最佳预测体系估计的育种值，筛选单穗产量、穗粒数和百粒重育种值较高的优良材料。

2结果与分析

2.1单穗产量、穗粒数和百粒重表型分析

多环境联合计算的单穗产量、穗粒数和百粒重的广义遗传力分别为0.60、0.61和0.70。多环境联合方差分析表明，309份材料间单穗产量、穗粒数和百粒重存在极显著遗传变异(表1)。不同环境间和基因型与环境互作也存在极显著差异。这说明3个性状受环境影响较大。

2.2全基因组关联分析结果

根据Bonferroni矫正设置的显著性阈值(P<1.72E-05)，利用FarmCPU、MLMM和CMLM方法分别检测到单穗产量、穗粒数和百粒重31个SNP、15个SNP和11个SNP(表2-表4)。单穗产量中，7个显著SNP在3个模型中均检测到，4个SNP在两个模型中检测到。其中，S3_62750920不仅在3个模型中检测到，而且在MLMM模型下的2个环境中也检测到。百粒重中，S4_215320252在3个模型中均检测到，S1_47210783在CMLM和MLMM方法中检测到。穗粒数8个SNP至少在2个模型中检测到。其中，S1_173095105、S3_120074687、S3_120074653和S1_143011251在3个模型中均检测到。单穗粒重和百粒重共有4个相同的位点，即S1_196234092、S1_143011251、S3_53872814和S3_120074687(表1和表3)。

2.3最佳全基因组预测模型的确定

对2017原阳、2017郸城、2017虞城、2017三亚和2019原阳的单穗产量、穗粒数和百粒重分别计算BLUE值，利用BLUE值进行预测。采用58129个SNP利用BayesA、BayesB、BayesC、Bayes LASSO、GBLUP和RKHS对单穗产量、穗粒数和百粒重进行5-倍交叉验证。对3个性状来说，RKHS模型的预测准确性最高(图1)。相比其他5个模型，RKHS模型对单穗产量预测准确性可以提高1.5倍，对穗粒数可以提高～1倍(图2)。相比其他5个模型，RKHS模型对百粒重的预测准确性可以提高0.02～0.05。

2.4优化的SNP密度

对3个性状来说，根据与性状关联显著性选择的SNP的预测效果均高于随机选择的相同数量的SNP，其中单穗产量可以提高0.75-1.32倍，穗粒数可以提高0.70-1.02倍，百粒重可以提高0.92-1.78倍(图2)。而且，与性状关联显著性选择的SNP的预测效果远高于Bonferroni阈值选择的SNP(图2和图3)。

对3个性状来说，利用最显著的500、1000和5000个SNP预测的准确性均高于利用最显著的10000、20000和40000个SNP(图2)。其中，单穗产量利用最显著的5000个SNP预测时，准确性最高为0.62，相比利用所有SNP预测的效果可以提高1倍。穗粒数和百粒重利用最显著的1000个SNP预测的准确性最高，分别为0.57和0.59。因此，对单穗产量、穗粒数和百粒重来说，最显著的5000个SNP，1000个SNP和1000个SNP为最佳的SNP密度。

2.5根据最佳预测体系选择育种值较高的优良材料

根据上述确定的最佳模型和最佳的SNP信息，利用BGLR的cbind(fm$yHat,y)函数(fm为RKHS模型，yHat为育种值，y为真实观测值即BLUE值)分别调出单穗产量5000个最显著SNP估计的育种值以及穗粒数和百粒重1000个最显著SNP估计的育种值结果，分别选择育种值最高的前10个材料。其中，L12和L80均在2个性状的最佳预测体系的前10名出现(表5-表7)。在单穗产量和百粒重的最佳预测体系中，L12的预测值分别73.71g和32.83g。在单穗产量和穗粒数的最佳预测体系中，L80的预测值分别为70.45g和293.25。因此，这两个材料可以作为候选高产材料，开展进一步育种研究工作。

表1单穗产量、穗粒数和百粒重多环境联合方差分析和遗传力

表2根据Bonferroni矫正设置的阈值(P<1.72E-05)检测到的单穗产量显著SNP

表3根据Bonferroni矫正设置的阈值(P<1.72E-05)检测到的穗粒数显著SNP

表4根据Bonferroni矫正设置的阈值(P<1.72E-05)检测到的百粒重显著SNP

表5根据最佳预测体系(RKHS模型、最显著的5000个SNP)预测的单穗产量最高的10个材料

表6根据最佳预测体系(RKHS模型、最显著的1000个SNP)预测的穗粒数最高的10个材料

表7根据最佳预测体系(RKHS模型、最显著的1000个SNP)预测的百粒重最高的10个材料