一种基于全基因组关联分析和多环境预测模型的玉米分子育种方法
阅读说明：本技术 一种基于全基因组关联分析和多环境预测模型的玉米分子育种方法 (Corn molecule breeding method based on whole genome association analysis and multi-environment prediction model ) 是由 马娟 曹言勇 于 2020-12-17 设计创作，主要内容包括：本发明公开了一种基于全基因组关联分析和多环境预测模型的玉米分子育种方法,该方法包含表型和基因型缺失数据的处理方法,显著关联SNP表型变异解释率的计算,多环境预测中核函数、不同环境效应和基因型与环境互作效应模型构建,最佳多环境联合预测和最佳SNP预测体系确定以及根据育种值综合选择优良材料的筛选方法。本发明建立了针对多环境数据鉴定和筛选玉米优良材料的方法。相比常规育种,根据基因型选择多环境目标性状高的材料,可实现高效、定向和精确育种。(The invention discloses a corn molecule breeding method based on whole genome association analysis and a multi-environment prediction model, which comprises a method for processing phenotype and genotype deletion data, calculation of obvious association SNP phenotype variation interpretation rate, construction of a kernel function, different environmental effects and genotype and environmental interaction effect model in multi-environment prediction, determination of an optimal multi-environment combined prediction and optimal SNP prediction system and a screening method for comprehensively selecting excellent materials according to breeding values. The invention establishes a method for identifying and screening excellent corn materials aiming at multiple environmental data. Compared with conventional breeding, the method selects materials with high environmental target characters according to genotypes, and can realize efficient, directional and accurate breeding.)
技术领域
本发明属于植物分子育种领域，具体涉及一种基于全基因组关联分析和多环境预测模型的玉米分子育种方法。
背景技术
全基因组关联分析是基于连锁不平衡的定位方法，是全基因组水平解析数量性状遗传构成的重要方法。全基因组关联分析不需要专门构建作图群体，具有分辨率高，高通量的优势，有利于鉴定现有种质资源中的有利等位基因。
利用全基因组关联分析，科研工作者已经挖掘大量控制重要农艺性状的关键位点。但这些位点如何应用于田间育种，是目前亟待解决的一个重要课题。分子标记辅助选择只能利用效应值较大的共分离标记对数量性状进行标记辅助选择，对于一些微效的位点没有作用。然而，研究实践证明，重要的农艺性状多数是复杂数量性状，而且绝大多数都是受微效多基因控制。利用分子标记辅助选择对这些性状进行遗传改良收效甚微。
全基因组选择利用训练群体的基因型和表型数据建模，对只有基因型的育种群体进行表型预测和选择，是改良复杂农艺性状非常有效的一种方法。常见的模型包括最小二乘法、最佳线性无偏预测以及贝叶斯模型等。目前预测多数基于单一环境下的单一性状，忽略了不同环境间的关联。但实际育种中，评价一个材料或品种的优劣，往往需要多年多点试验。如果利用单环境模型处理，就忽略了不同环境间的效应和基因型与环境互作的影响，并不能准确评价品系或品种的表现。因此，我们提出一种基于全基因组关联分析和多环境选择模型的玉米分子育种方法，快速评价、筛选优良材料，促进育种发展。
发明内容
本发明的目的在于提供一种基于全基因组关联分析和多环境预测模型的玉米分子育种方法。
本发明的目的可以通过以下技术方案实现：
一种基于全基因组关联分析和多环境预测模型的玉米分子育种方法，该方法包括如下步骤：
(1)材料和田间试验设计
在多个不同环境种植多份玉米种质自交系，种植过程中采用随机区组试验设计，每个小区选取授粉较好的穗子，测量玉米穗子的行粒数和穗行数；
(2)表型数据统计分析
利用QTL IciMapping v4.0的AOV功能对不同环境计算穗子行粒数和穗行数的最佳线性无偏估计值即BLUE值，将其和多个环境用于全基因组关联分析和多环境联合预测。利用Excel数据分析工具的相关系数分析对不同环境的行粒数和穗行数进行相关性分析；
(3)基因型鉴定和分析
采用GBS(genotyping by sequencing)简化测序的方法对多份自交系进行基因型分型，测序仪为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的检测；以缺失率小于0.10，杂合率小于0.10，最小等位基因频率(MAF)大于0.05为筛选标准，筛选获得多个高质量SNP用于关联分析；
(4)全基因组关联分析
利用不同环境和BLUE环境的穗行数和行粒数进行全基因组关联分析筛选与行粒数表型和穗行数表型显著关联SNP，采用的方法为CMLM(compressed mixed linearmodel)、MLMM(multiple loci mixed model)和FarmCPU(fixed and random modelcirculating probability unification)，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方法计算，显著临界值设置为P＝1/58129＝1.72E-05；CMLM方法的表型变异解释率(PVE)由软件给出；MLMM和FarmCPU方法的PVE由线性回归模型计算，公式为：Y＝α+βX+ε，其中Y为表型，α为截距，β为斜率，X为标记编码(2，0，1)，ε为随机误差；回归模型计算的为MLMM和FarmCPU的PVE；其中，为观测值的估计值，为观测值的平均值；
(5)多环境全基因组选择模型和多环境联合预测的选择标准
由于基因型和表型数据均有缺失，在做全基因组预测时需要对其进行缺失填补，根据已知基因型的基因型频率，随机选择标记进行填补，对填补后的SNP基因型进行编码，高频纯合等位基因型编码为2，低频纯合等位基因型编码为0，杂合基因型编码为1，填补和编码均在R语言中实现；
表型缺失值利用R语言mice包的pmm(预测均值匹配)方法进行多重插补，多重填补法的填补矩阵数设为5，迭代次数为50次；
拟合不同环境间主效应模型即MM模型的基本模型如下：
y＝μ1+ZEβE+Zμμ+ε
y＝(y1，...，yj…，ys)′是观测值的向量，yj是第j个环境某个自交系的观测值向量。ZE是环境效应βE的设计矩阵，βE为固定效应。Zμ是不同环境间的遗传效应μ的设计矩阵。 是不同环境主效应的方差。即GBLUP核。X是标记矩阵，p是标记的个数。
拟合不同环境间主效应和基因型与环境互作效应模型即MDs模型的基本模型如下：
y＝μ1+ZEβE+Zμμ+μe+ε
μe为随机效应，其服从°是Haddamar内积。是基因型与环境互作的方差。其余同MM模型。其中，该模型也采用GBLUP核。
利用MM和MDs模型分别选择不同的环境组合进行多环境联合预测，根据穗行数和行粒数不同环境间的相关系数，对不同的环境组合进行多环境联合预测均设置高相关性环境和低相关性环境；
进行多环境联合预测时，所选择的不同环境组合的数量依据步骤(1)中所选择的不同环境的数量来确定，例如，步骤(1)中选择在四个不同环境种植多份玉米种质自交系和4个不同环境计算的BLUE环境，则分别选择两个环境、三个环境、四个环境和五环境进行联合预测。对两环境、三环境和四环境联合预测均设置高相关性环境和低相关性环境。
(6)SNP密度的选择标准
比较三种全基因组关联分析方法的检测功效，明确检测显著关联位点最多的模型。根据该模型下BLUE值计算的SNP与性状关联的显著性(P值)，从小到大排序，分别选择与行粒数和穗行数显著的前500个SNP、1000个SNP、5000个SNP、10000个SNP、20000个SNP、所有SNP(58129)和Bonferroni矫正选择的性状最显著SNP(5个行粒数显著SNP和21个穗行数显著关联SNP)进行多环境联合预测。
所有多环境联合预测均采用5-倍交叉验证方法，重复100次，利用100次重复计算的育种值与真实值间的相关系数均值作为评价预测准确性的标准；该模型利用R的BGGE包实现，迭代次数为15000次，预烧burnin为5000次，thin设置为1。
(7)优良材料的筛选标准
根据预测准确性，确定不同环境联合预测的最佳体系。在最佳预测体系下，利用BGGE中的cbind(fit$yHat，y)函数调出不同环境联合预测各自交系材料的育种值。利用Excel里的RANK.AVG计算每个材料育种值的秩，并计算每个材料的秩均值，选择秩均值最小的前几个材料作为重要种质。
优选的，步骤(1)中不同环境的数量为2～8个环境。
优选的，步骤(1)中所述随机区组试验设计的方法为：设3次重复，小区为2行区，行距60cm，株距25cm；每个小区，选取授粉较好的3个穗子。
优选的，步骤(4)中根据Bonferroni矫正方法，设置的显著性阈值为1.72E-05，在该阈值下，利用FarmCPU、CMLM和MLMM方法检测穗行数和行粒数显著关联的SNP。
本发明的有益效果：
本发明建立了一种基于全基因组关联分析和多环境预测模型的玉米分子育种方法，相比常规育种，该方法根据基因型选择目标性状高的材料，可实现高效、定向和精确育种。
附图说明
图1为MDs多环境模型预测行粒数不同标记密度的准确性。
图2为MM多环境模型预测行粒数不同标记密度的准确性。
图3为MDs多环境模型预测穗行数不同标记密度的准确性。
图4为不同多环境预测穗行数MM模型不同标记密度预测准确性。
图5为行粒数多环境最佳预测体系估计的育种值的秩、秩均值和育种值均值。
图6为穗行数多环境最佳预测体系估计的育种值的秩、秩均值和育种值均值。
具体实施方式
以下结合具体实施例对本发明做出详细的描述。根据以下的描述和这些实施例，本领域技术人员可以确定本发明的基本特征，并且在不偏离本发明精神和范围的情况下，可以对本发明做出各种改变和修改，以使其适用各种用途和条件。
实施例1
1材料与方法
1.1材料和田间试验设计
关联群体包括黄淮海骨干自交系、国内核心种质和美国GEM等种质309份。2017年种植在商丘虞城、新乡原阳和海南三亚。2019年仅种植在原阳。采用随机区组试验设计，3次重复。小区为2行区，行距60cm，株距25cm。每个小区，选取授粉较好的3个穗子，测量穗行数和行粒数。
1.2表型数据统计分析
利用QTL Icimapping v4.0的AOV功能对2017虞城、2017三亚、2017原阳和2019原阳计算穗行数和行粒数的最佳线性无偏估计值即BLUE值。BLUE值和4个环境均用于全基因组关联分析和多环境预测。利用Excel数据分析工具的相关系数分析对不同环境的行粒数和穗行数进行相关性分析。
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的检测。以缺失率小于0.10，杂合率小于0.10，最小等位基因频率(MAF)大于0.05为筛选标准，共获得58129个高质量SNP用于关联分析。
1.4全基因组关联分析
利用虞城2017、原阳2017、原阳2019、三亚2017和BLUE值的穗行数和行粒数进行全基因组关联分析，筛选与行粒数表型和穗行数表型显著关联SNP。全基因组关联分析采用的方法为CMLM(compressed mixed linear model)、MLMM(multiple loci mixed model)和FarmCPU(fixed and random model circulating probability unification)。3种方法均采用Q(群体结构)+K(亲缘关系)模型。
Q值利用Structure v2.3.4计算。首先设置群体内亚群数为1～8，length ofburn-in period设置为5000，蒙特卡罗重复个数设置为50000，每个亚群数迭代次数为3。根据ΔK，确定亚群数为2时的Q值用于关联分析。K值利用TASSEL v5.0的Centered_IBS方法计算。由于CMLM方法为单位点检测方法，需要进行多重检验来确定显著关联临界值。而FarmCPU和MLMM方法均为多位点检测方法，不需要进行多重检验。因此，选择中度Bonferroni矫正方法，即显著临界值设置为P＝1/58129＝1.72E-05。
CMLM方法的表型变异解释率(PVE)由软件给出。MLMM和FarmCPU方法的PVE由线性回归模型计算，公式为：Y＝α+βX+ε，其中Y为表型，α为截距，β为斜率，X为标记编码(高频纯合等位基因型编码为2，低频纯合等位基因型编码为0，杂合基因型编码为1)，ε为随机误差。线性回归模型计算的为MLMM和FarmCPU的PVE；其中，为观测值的估计值，为观测值的平均值；
1.5多环境全基因组选择模型和多环境联合预测的选择标准
由于基因型和表型数据均有缺失，在做全基因组预测时需要对其进行缺失填补。根据已知基因型的基因型频率，随机选择标记进行填补。对填补后的SNP基因型进行编码，高频纯合等位基因型编码为2，低频纯合等位基因型编码为0，杂合基因型编码为1。填补和编码均在R语言中实现。
表型缺失值利用R语言mice包的pmm(预测均值匹配)方法进行多重插补，多重填补法的填补矩阵数设为5，迭代次数为50次。
比较三种全基因组关联分析方法的检测功效，选择检测显著关联位点最多的模型。根据该模型下BLUE值计算的SNP与性状关联的显著性(P值)，从小到大排序，选择与行粒数和穗行数显著的前500个SNP、1000个SNP、5000个SNP、10000个SNP、20000个SNP、所有SNP(58129)和Bonferroni矫正选择的性状最显著SNP(5个行粒数显著SNP和21个穗行数显著关联SNP)进行两环境、三环境、四环境和五环境联合预测。采用的模型为拟合不同环境间主效应模型(MM)和不同环境间主效应和基因型与环境互作效应模型(MDs)，具体如下：
MM模型：y＝μ1+ZEβE+Zμμ+ε
y＝(y1，...，yj...，ys)′是观测值的向量，yj是第j个环境某个自交系的观测值向量。ZE是环境效应βE的设计矩阵，βE为固定效应。Zμ是不同环境间的遗传效应μ的设计矩阵。 是不同环境主效应的方差。即GBLUP核。X是标记矩阵，p是标记的个数。
MDs模型：y＝μ1+ZEβE+Zμμ+μe+ε
μe为随机效应，其服从°是Haddamar内积。是基因型与环境互作的方差。其余同MM模型。其中，该模型也采用GBLUP核。
利用MM和MDs模型分别选择两个环境、三个环境、四个环境和五个环境进行联合预测。根据穗行数和行粒数不同环境间的相关系数，对两环境、三环境和四环境联合预测均设置高相关环境和低相关环境。
1.6SNP密度的选择标准
比较三种全基因组关联分析方法的检测功效，明确检测显著关联位点最多的方法。根据该方法下BLUE值计算的SNP与性状关联的显著性(P值)，从小到大排序，分别选择与行粒数和穗行数显著的前500个SNP、1000个SNP、5000个SNP、10000个SNP、20000个SNP、所有SNP(58129)和Bonferroni矫正选择的性状最显著SNP(5个行粒数显著SNP和21个穗行数显著关联SNP)进行多环境联合预测。
采用5-倍交叉验证，将309份材料分为训练集和预测集，重复100次。根据预测集100次估计的育种值与观测值间的相关系数均值作为评价预测准确性的标准。该模型利用R的BGGE包实现，迭代次数为15000次，预烧burnin为5000次，thin设置为1。
1.7多环境联合预测综合筛选优良材料的标准
根据预测准确性，确定不同环境联合预测的最佳体系。在最佳预测体系下，利用BGGE中的cbind(fit$yHat,y)(fit为MM或MDs模型，yHat为估计的育种值，y为实际观测值)函数调出不同环境联合预测各自交系材料的育种值。利用Excel里的RANK.AVG计算每个材料育种值的秩，并计算每个材料的秩均值，选择秩均值最小的前几个材料作为重要种质。
2结果与分析
2.1不同环境间行粒数和穗行数的相关性
BLUE值是4个原始环境均值的最佳估计值，与4个环境的相关性均较高(表1和表2)。其中，行粒数BLUE与2017虞城的相关性最高(r＝0.79)，其次为2017三亚(r＝0.70)和2017原阳(r＝0.69)，2019原阳最低(r＝0.63)。穗行数BLUE与2019原阳的相关性最高(r＝0.73)，与2017虞城的相关性最低(r＝0.61)。两个性状4个原始环境间均表现出较低的相关性(r＝0.19-0.46)。
2.2全基因组关联分析结果
根据Bonferroni矫正方法，设置的显著性阈值为1.72E-05。在该阈值下，利用FarmCPU、CMLM和MLMM方法共检测到5个行粒数显著关联SNP(P<1.72E-05)(表3)。其中S1_173095105、S5_127421583、S2_35077012和S2_35076923均在3种模型中检测到。4个SNP解释行粒数表型变异的5.3％-9.0％。穗行数共检测到21个显著SNP(P<1.72E-05)，其中FarmCPU检测到18个显著SNP，MLMM检测到3个显著SNP(表4)。其中，S8_71716395、S9_10867079和S9_107695183解释穗行数表型变异率较高，分别为9.18％、8.65％和9.20％。
2.3两环境联合预测的最佳预测体系
由于FarmCPU模型检测到显著位点个数最多，因此以该模型下根据BLUE值计算的SNP与性状关联的显著性(P值)从小到大排序，选择与行粒数和穗行数显著的前500个SNP、1000个SNP、5000个SNP、10000个SNP、20000个SNP、所有SNP(58129)和Bonferroni矫正选择的性状最显著SNP(5个行粒数显著SNP和21个穗行数显著关联SNP)进行两环境联合预测。
根据行粒数相关分析的结果(表1)，BLUE与2017虞城的相关性最高，为0.79(表1)；2017虞城和2017原阳是原始环境相关性最高(0.46)的两个环境。因此，行粒数选择2017虞城+2017原阳(低相关环境)和2017虞城+BLUE(高相关环境)分别进行两环境联合预测。同理，穗行数选择的代表性两环境分别为2017原阳+2019原阳(低相关环境)和2019原阳+BLUE(高相关环境)。
除了穗行数21个SNP外，穗行数和行粒数的MM模型的预测准确性均高于MDs模型(图1-图4)。两个性状高相关的两环境预测准确性均高于低相关的两环境(图1-图4)。选择显著的500-20000个SNP预测准确性均高于利用所有SNP。相比其他标记密度，行粒数最显著的5个SNP和穗行数21个最显著的SNP预测准确性最低，分别为0.14-10.19和0.34-0.40(图1-图2)。在2017年虞城+BLUE环境的MM模型下，行粒数选择最显著的1000个SNP预测准确性最高(0.58)(图2)，为行粒数两环境的最佳预测体系。穗行数在2019原阳+BLUE环境的MM模型，选择最显著的10000个SNP预测准确性最高(0.62)，为穗行数两环境的最佳预测体系(图4)。
2.3三、四环境和五环境联合预测的最佳预测体系
不同环境选择和标记密度选择标准同两环境联合预测。同两环境，除了5个和21个最显著的SNP外，三环境、四环境和五环境MM模型预测准确性均高于MDs模型(图1-图4)。而且，同两环境联合预测，三环境和四环境中高相关环境联合预测的准确性均高于低相关环境预测。三环境预测模型中，在2017原阳+2017虞城+BLUE环境MM模型下，利用最显著的5000个SNP预测行粒数效果最好，预测准确性为0.60(图2)。穗行数利用2017原阳、2017虞城和BLUE环境联合预测，在MM模型的500个SNP预测效果最佳，其准确性为0.58(图4)。
四环境联合预测中，行粒数和穗行数均利用2017原阳+2017虞城+2019原阳+BLUE环境的MM模型效果最佳(图2和图4)。行粒数最佳的标记密度为5000个显著SNP，准确性为0.55(图2)。而穗行数最佳的标记密度为500个显著SNP，准确性为0.54(图4)。
五环境联合预测中，行粒数和穗行数MM模型最佳的SNP密度分别为5000和500，其预测的准确性分别为0.55和0.49(图2和图4)。
2.4多环境模型最佳预测体系选择的穗行数和行粒数较高的材料
根据多个不同环境联合预测的最佳体系，利用BGGE的cbind(fit$yHat,y)函数分别调出两环境、三环境、四环境和五环境最佳预测体系中309份自交系的育种值。根据育种值，利用Excel里的RANK.AVG计算每个材料育种值的秩，并计算每个材料育种值的秩均值，选择秩均值最小的前20个材料进行展示。前20个自交系穗行数和行粒数不同多环境预测育种值的秩、秩均值和育种值的均值如图5和图6。根据育种值的秩均值，行粒数选择的最高的5个材料是L10、L8、L20、L9和L248，其行粒数介于24.09-24.63(图5)。根据秩均值，穗行数选择的最高的5个材料是L85、L18、L101、L121和L96，其穗行数介于14.03-15.57(图6)。这些材料可以作为重要种质，开展进一步育种研究工作。
表1行粒数不同环境间相关系数
表2穗行数不同环境间相关系数
表3不同环境CMLM、MLMM和FarmCPU检测的行粒数显著关联SNP
表4不同环境MLMM和FarmCPU检测的穗行数显著关联SNP
- 上一篇：一种医用注射器针头装配设备
- 下一篇：一种基因定点敲入载体构建方法、系统及平台