一种高效预测切花菊株高性状SNP位点及其方法和应用

一种高效预测切花菊株高性状snp位点及其方法和应用
技术领域
1.本发明属于生物技术领域，涉及一种高效预测切花菊株高性状snp位点及其方法和应用。


背景技术：

2.菊花(chrysanthemum morifolium ramat.)是菊科(asteraceae)菊属的多年生宿根花卉，是我国传统名花和四大切花之一，具有很高的观赏价值和经济价值(luo c et al.ssr analysis of genetic relationship and classification in chrysanthemum germplasm collection.horticultural plant journal,2018.4(2):73-82)。株高是决定菊花株型的关键性状之一，也是切花菊品质分级的重要标准。优质切花品种株高通常大于80cm(马婉茹等.多头切花菊品种茎、枝特性评价体系构建与品种评价.中国农业科学,2019.52(14):2515-2524)，但是株高过高又会导致易倒伏、生产耗能高等问题。因此，理想株高一直是切花菊育种的重要目标。节间数、节间长度和茎粗作为株高的主要构成因子，也是切花菊理想株型的重要指标。菊花株高是受多基因控制的数量性状，具有较高的遗传力(zhang f et al.mapping single-locus and epistatic quantitative trait loci for plant architectural traits in chrysanthemum.molecular breeding,2011.30(2):1027-1036；su j et al.combining ability,heterosis,genetic distance and their intercorrelations for waterlogging tolerance traits in chrysanthemum.euphytica,2017.213(2):e42)。此外，菊花株高与节间数、节间长度、茎粗、花径、花期等性状均存在显著的表型相关(韩勇等.标准切花菊花径与其他重要数量性状的相关与通径分析.南京农业大学学报,2012.35(1):33-37.)。
3.全基因组关联分析(genome wide association study，gwas)是一种通过检测全基因组遗传标记与表型变异之间关联程度定位与目标性状相关的遗传位点的有效方法。目前，利用gwas分析已鉴定到与菊花株型(li p et al.genetic diversity,population structure and association analysis in cut chrysanthemum(chrysanthemum morifolium ramat.).molecular genetics and genomics,2016.291(3):1117-1125)、花型(chong x et al.identification of favorable snp alleles and candidate genes responsible for inflorescence-related traits via gwas in chrysanthemum.plant molecular biology,2019.99(4-5):407-420)、抗蚜性(fu x et al.genetic variation and association mapping of aphid(macrosiphoniella sanbourni)resistance in chrysanthemum(chrysanthemum morifolium ramat.).euphytica,2018.214(2):e21)、耐旱性(li p et al.association analysis of drought tolerance in cut chrysanthemum(chrysanthemummorifolium ramat.)at seedling stage.3biotech,2018.8(5):e226)和耐涝性(su j,et al.genome-wide association study identifies favorable snp alleles and candidate genes for waterlogging tolerance in chrysanthemums.horticulture research,2019.6:e21)等性状相关的有利变异位点，但是
由于菊花的重要园艺性状大多数是数量性状，同时受多个微效基因控制，且鉴定到的位点大多属于环境特异性位点，因此在实际育种过程中尚未得到应用。不同于全基因组关联分析，全基因组选择(genomic selection，gs)不需要预先确定具有显著效应的分子标记，而是通过对基因组中所有标记遗传效应的评估预测个体的估计育种值，尤其对微效多基因控制的复杂性状具有较好的预测效果(juliana p et al.improving grain yield,stress resilience and quality of bread wheat using large-scale genomics.nature genetics,2019.51(10):1530-1539)。随着高通量测序技术的快速发展，越来越多的研究表明将gwas鉴定到的关联位点纳入gs模型可在一定程度上提高gs的预测能力，该方法在动植物育种方面取得了一定成果，但在菊花株高相关性状中尚未报道。


技术实现要素：

4.发明目的：针对现有技术存在的缺陷和不足，本发明提供了一种高效预测切花菊株高性状snp位点，该位点对株高、节间数、节间长度和茎粗的预测准确度分别达到0.97、0.98、0.97和0.96，可用于切花菊株高相关性状的预测，缩短了育种周期，加速育种进程。
5.本发明提供了一种高效预测切花菊株高性状snp位点的方法。
6.本发明还提供了一种高效预测切花菊株高性状snp位点及其方法在预测切菊花株高中的应用。
7.技术方案：为了实现上述目的，本发明提供了一种高效预测切花菊株高性状snp位点，所述切花菊性状snp位点及其特异性核苷酸如下所示：
8.[0009][0010]
本发明所述高效预测切花菊株高性状snp位点的方法，包括以下步骤：
[0011]
步骤1：选取来源不同且无直接亲缘关系的切花菊品种，于生长末期统计该群体在不同环境下的株高相关性状；
[0012]
步骤2：对步骤1统计的切花菊性状进行描述性统计及方差分析，并计算出各性状的广义遗传力；
[0013]
步骤3：采用gbs(genotyping-by-sequencing)对步骤1中的切花菊品种进行基因型分型，利用bwa软件将clean reads比对到菊花参考基因组，采用samtools软件进行snp calling，利用vcftools软件对检测到的snp过滤，获得的高质量snp用于全基因组关联分析；
[0014]
步骤4：基于步骤1获得的表型数据和步骤3获得的基因型数据进行全基因组关联分析，利用emmax法进行全基因组关联分析，依据显著性阈值筛选与各个性状显著关联的snp；
[0015]
步骤5：对切花菊性状进行全基因组选择分析，通过比较不同snp标记集的预测准确度，确定最佳标记集，筛选符合育种值需求材料。
[0016]
进一步地，所述切花菊株高相关性状为切花菊株高、节间数、节间长度或茎粗中任意一种。
[0017]
进一步地，所述步骤1中统计切菊花性状方法为选取不同环境下每个品种长势一致的植株进行表型测量。
[0018]
优选地，所述步骤1中的环境数量为2个，环境一为南京湖熟e1，环境二为南京锁石e2，于2020年10月中旬对每个品种选取长势一致的5个植株进行表型测量。
[0019]
进一步地，所述步骤2中的切菊花表型数据分析主要利用r语言完成，方差分析模型为y
ij
＝μ+gi+yj+gy
ij
+e
ij
。
[0020]
作为优选，所述方差分析模型y
ij
＝μ+gi+yj+gy
ij
+e
ij
中y
ij
为表型观测值，μ为多个环境下的均值，gi为基因型效应，yj为环境效应，gy
ij
为基因型与环境互作效应，e
ij
为随机误差。依据方差分析的结果计算各性状的广义遗传力(h2)，其表达式为h2＝σ2g/(σ2g+σ2ge/n+σ2e/rn)，σ2g为基因型方差，σ2ge为基因型与环境互作方差，σ2e为残差方差，n为环境个数，r为单个环境的重复数。
[0021]
进一步地，所述步骤3中snp利用vcftools软件按照
“‑‑
mindp 6
‑‑
min-alleles 2
‑‑
max-alleles 2
‑‑
max-missing 0.85
‑‑
maf 0.05”的标准过滤。
[0022]
进一步地，所述步骤4中emmax方法将主成分分析(pca)和亲缘关系k矩阵作为协变量计入模型，利用emmax软件对不同环境下的切菊花性状进行全基因组关联分析，计算出每个snp位点的显著度即p值，将p《1
×
10-5
作为显著关联的阈值筛选出潜在的候选snps。
[0023]
进一步地，所述步骤5中具体包括以下步骤：
[0024]
(1)由于个别品种的表型数据有缺失，在全基因组选择分析时需要对其进行缺失填补；若表型数据有缺失，利用r软件hmisc包中impute()函数进行自动插补，若基因型数据有缺失，利用rrblup包中a.mat()函数自动填补，并将snp的基因型格式编码为{-1,0,1}基因型格式；
[0025]
(2)设置不同标记方案对不同环境下的目标性状进行预测；
[0026]
(3)采用五倍交叉验证法，利用rrblup模型进行全基因组选择，基因组估计育种值(gebv)与实际观测值的相关系数作为评价预测准确度的指标，相关系数最大的标记集为最佳snp标记集。
[0027]
作为优选，所述rrblup模型的表达式为：y＝wgu+ε，y是表型向量，u～n(0,iσ
2u
)是标记效应向量，w表示基因型和表型间的关联矩阵，g表示遗传矩阵，ε是随机误差，该模型的解为：其中z＝wg,λ＝σ
2e
/σ
2u
指的是残差方差与标记效应方差的比值。
[0028]
进一步地，所述不同标记方案具体为：依据所述emmax模型计算的p值对snps从小到大排序，选择显著阈值下的snp以及与目标性状显著相关的前100～2000个snps作为方案一，或从上述步骤3中获得的高质量snp中随机选择与方案一同等数量的snps作为方案二；利用这两种方案分别对不同环境下的目标性状进行预测，确定最佳的snp标记集。
[0029]
本发明所述高效预测切花菊株高性状snp位点在预测切花菊株高中的应用。
[0030]
本发明所述高效预测切花菊株高性状的方法在预测切花菊株高中的应用。
[0031]
进一步地，所述高效预测切花菊株高相关性状snp位点或高效预测切花菊株高相关性状snp位点方法在预测切花菊株高中的应用，包括以下步骤：
[0032]
(1)根据筛选出的最佳snp标记集，利用公式y＝μ+xg+e估计200个不同切花菊品种株高育种值gebv(基因组估计育种值)；
[0033]
(2)当切花菊品种株高的育种值大于80cm时，该品种可作为优质切花菊品种。
[0034]
本发明基于现有种质资源群体，利用gwas(全基因组选择关联分析)筛选与株高相关性状显著关联的snp标记，为菊花株型的遗传改良提供位点信息。此外，在gwas的基础上，利用与目标性状相关的分子标记，通过gs(全基因组选择)对个体的育种值进行估计，实现切花菊株高性状的早期预测和选择，从而缩短繁育周期，加速育种进程。
[0035]
本发明的有益效果：与现有技术相比，本发明具有以下优点：
[0036]
(1)本发明使用数量多、分布广泛、能够稳定遗传且易于高通量基因分型的snp标
记进行切花菊株高相关性状的全基因组关联分析，通过全基因组关联分析鉴定到与株高、节间数、节间长度、茎粗四个性状相关的snp标记，为菊花株高相关性状功能标记的开发和遗传改良提供了位点信息。
[0037]
(2)本发明在全基因组关联分析的基础上，结合鉴定到的snp标记进行全基因组选择，在不同snp标记集下筛选出对切花菊株高相关性状进行有效预测的最佳snp标记，这些最佳标记对株高、节间数、节间长度和茎粗的预测准确度分别为0.97、0.98、0.97和0.96，初步认为可以应用于切花菊株高相关性状的早期预测，便于快速筛选出具有理想株高的切花菊优异材料，缩短育种周期，对于培育理想株型切花菊新品种具有重要理论和实践意义。
附图说明
[0038]
图1为切花菊在环境e1和环境e2中的表型分布直方图，其中a图为株高、b图为节间数、c图为节间长度、d图为茎粗；
[0039]
图2为切花菊株高在随机选择和性状相关snp下的全基因组预测准确度，其中图a为e1环境，图b为e2环境。
[0040]
图3为切花菊节间数在随机选择和性状相关snp下的全基因组预测准确度，其中图a为e1环境，图b为e2环境。
[0041]
图4为切花菊节间长度在随机选择和性状相关snp下的全基因组预测准确度，其中图a为e1环境，图b为e2环境。
[0042]
图5为切花菊茎粗在随机选择和性状相关snp下的全基因组预测准确度，其中图a为e1环境，图b为e2环境。
[0043]
图6为基于最佳snp标记集和两个环境表型均值的株高、节间数、节间长度和茎粗四个性状的全基因组预测准确度。
[0044]
图7为一种基于全基因组选择高效预测切花菊株高相关性状的方法的流程图。
具体实施方式
[0045]
下面结合实施例对本发明做进一步的详细说明。
[0046]
实施例1
[0047]
(一)试验材料与田间设计
[0048]
选取200份不同来源且无直接亲缘关系的切花菊品种为供试材料(表1)，所有材料来自于南京农业大学中国菊花种质资源保存中心，本领域技术人员可从“中国菊花种质资源保存中心”中获得上述种质。于2020年7月分别种植在南京农业大学湖熟基地(e1)和锁石基地(e2)。小区行距和株距均为0.10m，水肥药及其他田间管理技术措施按菊花优质生产技术操作。在营养生长末期(十月下旬)，选取各品种长势一致的5个植株，测量株高、节间数、节间长度和茎粗。测量方法：株高(ph/cm):茎基部到花冠顶部的高度；节间数(in/个):茎基部到花冠顶部所有的节间数量；节间长度(il/cm)：植株中上部节间长度的平均值；茎粗(sd/mm)：植株中上部茎的直径。
[0049]
表1 200份供试切花菊品种
[0050]
[0051][0052]
(二)表型数据统计分析
[0053]
在e1和e2环境下，200份切花菊品种株高、节间数、节间长度和茎粗四个性状的表型分布如图1所示，利用r语言的
‘
pastecs’包对步骤(一)e1和e2环境下的株高、节间数、节间长度和茎粗四个性状进行描述性统计，利用
‘
lm4’包进行方差分析，方差分析模型为y
ij
＝μ+gi+yj+gy
ij
+e
ij
，其中y
ij
为表型观测值，μ为多个环境下的均值，gi为基因型效应，yj为环境效应，gy
ij
为基因型与环境互作效应，e
ij
为随机误差，依据方差分析的结果计算各性状广义遗传力(h2)，其表达式为h2＝σ2g/(σ2g+σ2ge/n+σ2e/rn)，σ2g为基因型方差，σ2ge为基因型与环境互作方差，σ2e为残差方差，n为环境个数，r为单个环境的重复数。
[0054]
结果表明200份切花菊材料的株高、节间数、节间长度和茎粗四个性状存在广泛的变异，变异系数为16.49～30.18％(表2)。株高、节间数、节间长度和茎粗的广义遗传力分别为79.61％，86.06％，84.42％和79.21％(表2)，说明切花菊株高、节间数、节间长度和茎粗四个性状主要受遗传因素的影响。方差分析结果表明除基因型效应外，环境效应和基因型与环境互作效应对株高相关性状也存在显著的影响(表3)。
[0055]
表2株高、节间数、节间长度和茎粗四个性状的描述性统计
alleles 2
‑‑
max-missing0.85
‑‑
maf 0.05”的标准进行过滤，最终获得330,710个高质量snps用于进一步的全基因组关联分析。
[0062]
(四)全基因组关联分析
[0063]
基于步骤(一)获得的表型数据和步骤(三)获得的基因型数据进行全基因组关联分析，利用emmax方法分别对e1和e2下的株高、节间数、节间长度和茎粗进行全基因组关联分析。该方法将主成分分析(pca)和亲缘关系k矩阵作为协变量计入模型。采用emmax软件的“emmax-kin”子程序计算该群体之间的亲缘关系k矩阵。利用gcta软件进行主成分分析，选取前10个主成分作为协方差矩阵用于关联分析。利用emmax软件计算出每个snp位点的显著度即p值，将p《1
×
10-5
作为显著关联的阈值筛选出潜在的候选snp(表4)。
[0064]
如表4所示，共检测到12个与株高相关的snp位点，其中有3个snp来自e1，剩余的9个snp来自e2。除了chr16__103728540and chr21__9532188两个位点外，其他的位点均表现出增效效应。共检测到12个与节间数相关的snp位点，其中只有2个snp来自e2。共检测到10个与节间长度相关的snp位点，其中5个snp来自e1，5个snp来自e2，表型效应值的范围为0.2-0.45。共检测到8个与茎粗相关的snp位点，其中只有chr8__245595683表现出减效效应。然而，同一个性状在两个环境下没有检测到相同的位点，说明环境效应和基因型与环境互作效应在个体的遗传变异中起着重要作用。此外，株高和节间数共有2个相同的关联位点，即chr12__270808434和chr12__270808478，且均表现出增效效应，说明会这些性状可能由多效性基因或紧密连锁的基因位点控制。
[0065]
表4与株高、节间长度、节间数和茎粗四个性状显著相关snp位点
[0066]
[0067][0068]
注：a加粗snp位点为重复检测到的等位位点；b为表型效应值；**和*分别表示在0.01和0.05水平差异显著，ns表示差异不显著。
[0069]
(五)全基因组选择
[0070]
由于表型数据和基因型数据均有缺失，表型数据有缺失，利用r软件hmisc包中impute()函数进行自动插补，基因型数据有缺失，利用rrblup包中a.mat()函数自动填补，并将snp的基因型格式编码为{-1,0,1}基因型格式。设置两种标记方案，依据emmax模型计算的p值对步骤(三)中330,710个snps从小到大排序，选择显著阈值下的snps以及与目标性状(株高、节间数、节间长度和茎粗)显著相关的前100、300、500、1000、5000、10000、15000和20000个snps作为方案一；从步骤(三)获得的330,710个snps中随机选择与方案一同等数量的snps作为方案二；利用这两种方案分别对e1和e2环境下的目标性状进行预测，确定最佳的snp标记集。以e1和e2环境下的表型均值为真实值，以两个环境下最佳的snp标记的并集为标记集，利用全基因组选择对株高、节间数、节间长度和茎粗四个性状的预测准确度进行验证。全基因组选择利用r语言的“rrblup”包完成，该模型的表达式为：y＝wgu+ε，y是表型
向量，u～n(0,iσ
2u
)是标记效应向量，w表示基因型和表型间的关联矩阵，g表示遗传矩阵，ε是随机误差，该模型的解为：是随机误差，该模型的解为：其中z＝wg,λ＝σ
2e
/σ
2u
指的是残差方差与标记效应方差的比值。采用五倍交叉验证法重复500次以提高模型的泛化能力，利用rrblup模型进行全基因组选择，基因组估计育种值(gebv)与实际观测值的相关系数作为评价预测准确度的指标，相关系数最大的标记集为最佳snp标记集。
[0071]
全基因组选择结果表明，无论是在e1还是e2环境下，性状相关snp的预测准确度远远高于同等数量下随机选择snp的预测准确度。对于株高、节间数、节间长度、茎粗四个性状来说，随着snp标记个数的增加至1000时，预测的准确度呈现显著增加的趋势，但当snp标记个数从1000增加至20000时，预测准确度的变化趋势呈现近乎水平状态(表5，图2-图5)。对于株高来说，e1和e2环境下，最显著相关的1000个snps的预测准确度分别为0.97和0.96；对于节间数来说，e1和e2环境下，最显著相关的1000个snps的预测准确均为0.97；对于节间长度来说，e1和e2环境下，最显著相关的1000个snps的预测准确度均为0.96；对于茎粗来说，e1和e2环境下，最显著相关的1000个snps的预测准确分别为0.95和0.96。可见，最显著相关的1000个snps以相对较少的标记即可获得相当高的预测准确度，将其确定为最佳snp标记集。基于两个环境下最佳的snp标记的并集和两个环境的性状表型均值进行全基因组选择分析，株高、节间数、节间长度和茎粗四个性状的预测准确度分别为0.97，0.98，0.97，0.96(图6)，进一步说明了本发明开发的snp标记用于预测切花菊株高、节间数、节间长度和茎粗的可行性。
[0072]
表5株高相关性状在不同分子标记集下的全基因组预测准确度
[0073]
[0074][0075]
注：性状相关_sig表示显著阈值下的snps(株高：sig＝12；节间数：sig＝12；节间长度：sig＝10；茎粗：sig＝8)；随机选择_ran表示随机选择的与显著阈值同等数量的snp(株高：ran＝12；节间数：ran＝12；节间长度：ran＝10；茎粗：ran＝8)。
[0076]
实施例2根据筛选的最佳snp标记集选择优良育种材料
[0077]
根据实施例1中选取的最佳snp标记集(表5)和全基因选择预测的流程图(图7)，利用公式y＝μ+xg+e估计200个不同切花菊品种株高育种值gebv(基因组估计育种值)，其中y为预测表型值，μ为训练群体表型均值，x为测试群体表型均值，g为标记效应矩阵，e为随机效应矩阵。将育种值(预测表型值)大于80cm的品种作为切花菊优选品种。
[0078]
表6基于gs最佳snp标记集筛选出的株高符合切花菊生产标准的优质品种
[0079][0080]
如表6所示，共有40个品种株高的育种值大于80cm，可作为优质切花菊的候选材料，用于后续株型育种工作。

技术特征：
1.一种高效预测切花菊株高性状snp位点，其特征在于，所述切花菊株高性状snp位点及其特异性核苷酸如下所示：
2.一种权利要求1所述高效预测切花菊株高性状snp位点的方法，其特征在于，包括以下步骤：步骤1：选取来源不同且无直接亲缘关系的切花菊品种，于生长末期在不同环境下统计该群体的株高相关性状；步骤2：对步骤1统计的切花菊性状进行描述性统计及方差分析，并计算出各性状的广义遗传力；步骤3：采用gbs(genotyping-by-sequencing)对步骤1中的切花菊品种进行基因型分型，利用bwa软件将clean reads比对到菊花参考基因组，采用samtools软件进行snp calling，利用vcftools软件对检测到的snp过滤，获得的高质量snp用于全基因组关联分析；步骤4：基于步骤1获得的表型数据和步骤3获得的基因型数据进行全基因组关联分析，依据显著性阈值筛选与各个性状显著关联的snp；步骤5：对切花菊株高相关性状进行全基因组选择分析，通过比较不同snp标记集的预测准确度，确定最佳标记集，筛选符合育种值需求材料。3.根据权利要求2所述的高效预测切花菊株高相关性状snp位点的方法，其特征在于，所述步骤1切花菊株高相关性状为切花菊株高、节间数、节间长度或茎粗中任意一种。4.根据权利要求2所述的高效预测切花菊株高相关性状snp位点的方法，其特征在于，所述步骤2中的切菊花表型数据分析方差分析模型为y
ij
＝μ+g
i
+yj+gy
ij
+e
ij
。5.根据权利要求2所述的高效预测切花菊株高相关性状snp位点的方法，其特征在于，所述步骤3中snp利用vcftools软件按照
“‑‑
mindp 6
‑‑
min-alleles 2
‑‑
max-alleles 2
‑‑
max-missing 0.85
‑‑
maf 0.05”的标准过滤。6.根据权利要求2所述的高效预测切花菊株高相关性状snp位点的方法，其特征在于，所述步骤4对不同环境下的切菊花性状进行全基因组关联分析，计算出每个snp位点的显著度即p值，将p<1
×
10-5
作为显著关联的阈值筛选出潜在的候选snps。7.根据权利要求1所述的高效预测切花菊株高相关性状snp位点的方法，其特征在于，所述步骤5中具体包括以下步骤：(1)若表型数据有缺失，利用r软件hmisc包中impute()函数进行自动插补，若基因型数据有缺失，利用rrblup包中a.mat()函数自动填补，并将snp的基因型格式编码为{-1,0,1}基因型格式；(2)设置不同标记方案对不同环境下的目标性状进行预测；(3)采用五倍交叉验证法，利用rrblup模型进行全基因组选择，基因组估计育种值(gebv)与实际观测值的相关系数作为评价预测准确度的指标，相关系数最大的标记集为最佳snp标记集。
8.一种权利要求1所述的高效预测切花菊株高性状snp位点在预测切花菊株高中的应用。9.一种权利要求2～7任一项所述的高效预测切花菊株高性状snp位点的方法在预测切花菊株高中的应用。10.根据权利要求8所述的高效预测切花菊株高性状snp位点在预测切花菊株高中的应用，其特征在于，包括以下步骤：(1)根据筛选出的最佳snp标记集，利用公式y＝μ+xg+e估计200个不同切花菊品种株高育种值gebv；(2)当切花菊品种株高的育种值大于80cm时，该品种可作为优质切花菊品种。

技术总结
本发明公开了一种高效预测切花菊株高性状SNP位点，还公开了高效预测切花菊株高性状SNP位点方法及其应用。本发明使用数量多、分布广泛、能够稳定遗传且易于高通量基因分型的SNP标记进行切花菊株高相关性状全基因组关联分析，结果表明，与株高性状相关的最佳标记对株高、节间数、节间长度和茎粗的预测准确度分别为0.97、0.98、0.97和0.96，可以应用于切花菊株高、节间长度、节间数和茎粗的早期预测，缩短育种周期，为菊花株高相关性状功能标记的开发和遗传改良提供了位点信息，对于选育不同株型菊花优良材料具有重要理论和实践意义。菊花优良材料具有重要理论和实践意义。


技术研发人员：张飞 张雪峰 苏江硕 管志勇 房伟民 陈发棣
受保护的技术使用者：南京农业大学
技术研发日：2023.06.20
技术公布日：2023/9/20