基于GWAS的西瓜种子性状候选QTL鉴定

高美玲，胡创然，袁成志，郭 宇，刘秀杰，刘继秀，高 越

（1.齐齐哈尔大学生命科学与农林学院，黑龙江齐齐哈尔 161006；2.黑龙江省抗性基因工程与寒地生物多样性保护重点实验室，黑龙江齐齐哈尔 161006；3.齐齐哈尔市园艺研究所，黑龙江齐齐哈尔 161006）

在大多数植物中，种子存在于生命周期的开始到结束，尤为重要[1-2]。作为重要的瓜类作物之一，西瓜在世界范围有着极广的栽培面积。西瓜分为籽用西瓜和瓤用西瓜，对于瓤用西瓜，种子小能提高食用口感；对于籽用西瓜，大种子内部营养物质含量高。所以对于两种用途的西瓜，种子大小都是其重要的育种目标，研究控制种子大小的机制有利于实现对其遗传操控[3-4]。

种子大小是数量性状，具有复杂性。目前已报道不少控制西瓜种子相关表型性状的基因，这类基因之间的相互作用形成了西瓜种子表型性状的遗传多样性[5]。J.Prothro等[6]利用F2分离群体和RIL群体，共定位到13个和种子大小性状相关的QTL，其中在连锁群LG2上定位到与种子长度、宽度和百粒重相关的QTL，表型贡献率在26.9%～73.6%范围之间，认为其是控制西瓜种子大小的主效QTL。与G.Meru和C.Mcgregor检测到的3个QTLs[7]定位在了相同的连锁群上。QTL定位的方法仅能对在分离群体的亲本材料间存在差异的基因效应进行分析，无法在全基因组范围广泛挖掘种子大小调控的基因[8]。而全基因组关联分析则是近年发展起来的作用于研究复杂性状的一种高效方法，具有高分辨率和高通量的优点，可同时对多个复杂性状进行关联，检测多个等位基因[9-11]。随着西瓜基因组测序的完成，利用GWAS研究，通过高通量数据研究西瓜一些复杂的数量性状变得更便捷。郭禄芹等[12]利用GWAS检测173份西瓜种质资源材料，两个环境中共同检测到与种子长度关联的QTL共4个，分别位于2、7和9号染色体；两个环境中共同检测到与种子宽度关联的1个QTL，位于11号染色体。

目前关于西瓜种子大小的GWAS研究尚少。本研究以实验室收集的144份纯合的西瓜种质资源为自然群体材料，通过全基因组重测序技术获得基因型，并调查自然群体表型数据，将基因型数据与表型数据结合进行全基因组关联分析，鉴定与西瓜种子大小显著关联的SNP，为进一步研究重要农艺性状奠定基础。

1 材料和方法

1.1 供体试材

试验材料为226份种子大小变异丰富的种质资源，来源于世界各地。经过纯度鉴定筛选后得到144份纯系材料，其中113份种质材料来自亚洲，11份来自非洲，10份来自北美洲，8份来自欧洲，2份来自南美洲。144份种质资源材料见表1。

表1 144份种质资源名称Table 1 The name of 144 germplasm resources

田间试验分3次，场地均为齐齐哈尔市园艺所。2017年4—8月进行第一次试验，种植了226份材料中的32份，株行距为55 cm×40 cm，吊蔓整枝，开花期自交授粉。第二雌花留瓜，并为了保证果实营养，每株只取一个瓜，严格进行田间管理。授粉后30 d获取种子，晾干后鉴定性状。2018年4—8月进行第二次田间试验，种植226份材料于基地大棚内，栽培和管理方式与之前相同。2019年4—8月进行第三次田间试验，种植鉴定为纯合的144份西瓜材料于基地大棚和露地，露地环境下，每个材料种植8株，株行距为55 cm×140 cm。大棚环境下，栽培和管理方式与之前相同。

1.2 自然群体表型数据调查

参照《西瓜种质资源描述规范和数据标准》[13]，对西瓜种子大小相关性状进行鉴定。鉴定标准如下：

种子长度（seed length，SL）：通过利用游标卡尺测量种子喙部边缘到底部边缘间的最大距离，单位为mm。

种子宽度（seed width，SW）：通过利用游标卡尺测量种面纵轴垂直方向两边缘间的最大距离，单位为mm。

种子百粒重（100-seed weight，100SWT）：利用万分之一天平（上海越平科学仪器有限公司，FA1004B）称量每20粒种子总质量，进行3次重复，最终折算成百粒重，单位为g。

平均值、变异系数、标准差使用Micosoft Excel 2016计算，相关性分析使用SPSS 20.0软件，多样性指数（Shannon-Wiener）计算公式为 H=ΣNi/N（ln Ni/N）[14]。

1.3 文库构建原理及流程

将基因组DNA用NEBNext DNA双链片段化酶随机打断成短的DNA片段，加入End Prep EnzymeMix和 End Repair Reaction Buffer（10×）进行平末端修复并在DNA片段两端连接dA尾，加入Blunt/TA Ligase Master Mix、NEBNext Adaptor for Illumina、Ligation Enhancer使测序接头连接。经过AMPure XP磁珠纯化，选择300～400 bp范围的片段进行PCR扩增，建好DNA文库。使用安捷伦2100或质检试剂盒DNA 1000 Assay Kit检测文库长度、浓度等。

1.4 测序及数据质量控制

将质检合格的文库使用Hiseq×10 PE150上机测序，得到原始图像数据后经base calling转化为序列数据，结果以FASTQ文件格式存储，对下机的clean reads再次过滤，得到high quality clean reads，用于后续的信息分析。以“97103 v2”作为参考基因组，运用对比软件BWA（0.7.15）的mem算法将过滤后的reads比对到参考基因组上。利用SAMtools将比对结果输出BAM格式文件，最后使用Picard（1.129）软件（Mark Duplicates）标记重复 reads，统计reads深度及覆盖度。使用软件GATK对样本的Variant检测，检测到变异（SNP和Indel）使用Variant Filtration过滤。

因为稀有等位（群体中频率很低的等位）、高缺失率、高杂合率等位点会引起群体分析和全基因组关联分析异常，使用自写的perl脚本，按照以下条件对SNP和Indel进行过滤：①非二等位位点去除；②第二等位基因频率（minor allele frequency，MAF）小于0.05的位点去除；③缺失率大于0.5的位点去除；④杂合比例大于0.8的位点去除。

1.5 全基因组关联分析

将经过全基因组筛选得到的984 827个SNP标记与种子大小相关性状结合进行关联分析，本试验使用TASSEL（v5.2.24）软件进行分析包括：一般线性模型（GLM模型）和混合线性模型（MLM模型）。一般线性模型可分为简单线性模型（LM或GLM模型）和广义线性模型[GLM（Q）模型]。混合线性模型包括加入K矩阵混合线性模型[简称MLM（K）模型]和加入Q矩阵和K矩阵的混合线性模型[MLM（QK）模型]。利用观测值p值，取p值-log10的对应值绘制Q-Q散点图（quantile-by-quantile plot）和曼哈顿图（Manhattan plot），通过比较各个模型的不同，选择最适的关联模型作为后续的试验分析的模型。

1.6 关联性状候选区段鉴定

基于群体的SNP标记，对自然群体进行LD分析，以r2作为主要参数，计算各条染色体及基因组SNP两两组合的LD强度，并绘制西瓜各条染色体及基因组水平的连锁不平衡衰减图。选择在曼哈顿图中有成簇SNP分布的位置为QTL区段，再根据连锁不平衡衰减距离分析关联性状的候选区段。

2 结果与分析

2.1 西瓜种子大小性状遗传分析

由表2可知，4次重复种子长度平均多样性指数为4.69，种子宽度为4.78，种子百粒重为4.80。4次重复种子百粒重平均变异系数最高，为75.19%；种子长度最低，为27.64%；种子宽度为28.56%。种子长度最大为17.64 mm，最小为5.01 mm，变异幅度为12.63；种子宽度最大为11.87 mm，最小为2.85 mm，变异幅度为9.02，种子长度和种子宽度变异系数相近；种子百粒重最大为29.65 g，最小为0.84 g，变异幅度为29.47。

表2 144份西瓜种质资源种子大小性状遗传分析Table 2 Genetic analysis of seed size traits in 144 watermelon germplasm resources

2.2 西瓜种子大小性状相关性分析

由表3可知，种子长度、种子宽度与种子百粒重两两之间呈极显著正相关。其中种子长度与种子宽度相关系数为0.967，种子长度与种子百粒重相关系数为0.947，种子宽度和种子百粒重相关系数为0.937。

表3 西瓜种子大小性状相关性分析Table 3 Analysis of the correlation of watermelon seed size

2.3 重测序数据的分析

对144份西瓜种质材料进行重测序，得到原始下机数据（raw data），经过过滤后的数据（clean data）共858.17 G，平均每个样品5.96 G。通过筛选去除含有adapter的reads；去除含N比例大于10%的reads，去除低质量reads（质量值 Q≤20的碱基数占整条read的50%以上），经过滤得到的有效数据（HQ clean data）共847.90 G，平均每个样品5.89 G。测序共获得标记5 800 304个，筛选后保留标记1134 013个，包括SNP位点984 827个，Indel位点149 186个（表4），SNP和Indel标记在全基因组范围内分布均匀（表 5）。

表4 标记过滤筛选结果统计表Table 4 The statistics results of markers after filtering

表5 标记在全基因组染色体上的分布Table 5 Ditribution of markers on chromosomes in the whole genome

SNP和Indel标记位于基因区和非基因区的位置情况见表6。

表6 标记位于基因区和非基因区的位置Table 6 Markers are located in genetic and non-genic regions

2.4 西瓜种子大小性状全基组关联分析

本研究基于4种统计模型进行关联分析（图1），将4种模型取-log10（p）后，进行Q-Q散点图绘制，如图所示，图中的MLM（QK）模型观测值较其他3个模型相比更近似于期望值，为更好地保证结果，减少假阳性，我们在全基因组关联分析中选择MLM（QK）模型进行后续试验。

图1 种子大小相关性状在GLM、GLM（Q）、MLM（K）和MLM（QK）模型的Q-Q图比较Figure 1 Comparison of Q-Q diagrams of seed size-related traits in GLM，GLM（Q），MLM（K）and MLM（QK）models

通过基因型与表型的关联分析（图2，表7），以0.05为筛选阈值，在此阈值基础上，MLM（QK）模型下种子百粒重性状共检测到显著性SNP 17个，分布在4、7和10号染色体上，p值的变化范围为7.36～8.54。用r2表示SNP对表型变异的解释率，表型变异的解释率范围16.97%～22.97%，均值为19.30%。

图2 144份西瓜种质材料在种子大小全基因组关联分析的曼哈顿图Figure 2 144 Manhattan plots of seed size genome-wide association analysis of watermelon germplasm materials

表7 关联分析检测到的17个相关性状的QTLTable 7 17 QTLs for related traits detected by association analysis

2.5 候选区段鉴定

当r2衰减至0.5时对应的物理距离为73.9 kb（约74 kb）（图3），因此该自然群体的衰减距离为74 kb，对显著性的SNP簇所在物理位置上下游加上衰减距离，形成的物理位置范围作为QTL物理区间，该区间也可以作为与性状相关联的候选区段。由此，得出与种子大小性状显著关联的区段结果（表8）。

表8 检测到与种子大小性状显著关联的QTL区段结果统计Table 8 Result statistics of QTL regions that are significantly associated with seed size traits detected

图3 LD衰减图Figure 3 Linkage disequilibrium decay

表8为检测到的与种子大小性状显著相关的QTL区段，2017年检测到区段 Chr4（18 555 982～18 703 982 bp）包含一个SNP位点，r2均值为17.85%；区段Chr7（1553 660～1719 419bp）包含5个SNP位点，r2均值为 19.22%；区段 Chr10（8 639 366～8832395bp）包含10个SNP位点，r2均值为19.12%。2018年大棚环境下，检测到区段 Chr7（1 553 030～1 701 030 bp）包含1个SNP位r点，r2均值为22.97%。由于7号染色体上的QTL区段被两次重复检测到，遗传稳定，且两次检测到的区段有大范围上的重叠，故合并为同一个区段 Chr7（1 553 030～1 719 419 bp）。

3 讨论

3.1 西瓜种子大小多样性与相关性分析

在本试验中，种子长度与种子宽度变异系数平均值分别为27.64%和28.56%，种子百粒重变异系数最高，达到了75.19%，种子长度、种子宽度和种子百粒重多样性指数平均值分别为4.69、4.78和4.80。王学征等[15]利用62份西瓜种质资源中，分析得出种子长度、种子宽度以及种子千粒重多样性指数4.20～4.32。郭禄芹利用173份西瓜材料，种子长度、种子宽度多样性指数为1.57～1.58[12]。可见本试验所用的144份种质资源种子大小表型变异更为丰富。各个性状在4次试验的变异系数的变化较小，说明材料种子大小相关性状表型遗传较稳定。

从本试验西瓜种子大小相关性分析结果来看，种子长度、种子宽度都与种子百粒重呈极显著正相关，在前人种子大小QTL定位中也得到相同结果，2014年刘传奇等[16]利用F2分离群体，将种子长度、宽度定位与种子百粒重定位到相同区域。这可能与基因一因多效遗传机制有关[17-18]。

3.2 控制西瓜种子大小性状基因定位分析

本试验利用144份西瓜种质资源种子长度、种子宽度与种子百粒重性状做全基因组关联分析，共检测到3个过阈值的QTL区段，均由与种子百粒重性状关联定位得到，分别位于4、7和10号染色体上。位于7号染色体上的QTL区段被两次重复检测到，遗传稳定。王学征等[15]利用62份西瓜种质资源，在1、4、6、7和10号染色体上定位到了种子大小相关QTL，以往的QTL定位研究多集中在1、6号染色体上，而对于4、7和10号染色体研究较少，种子大小QTL区段还未有人进行过精细定位[19]，染色体基因位点有待进一步确认。