玉米籽粒长度的全基因组关联分析

代力强，吴 律，董青松，闫 鸽，曲 静，王丕武

(吉林农业大学 农学院，吉林 长春130118)

粒长是玉米产量性状遗传改良的主要目标之一。选择相对较长的玉米籽粒，对提高玉米的单产水平具有重要意义。已有的研究表明，籽粒长度是高度遗传的数量性状，受环境的影响较小[1-3]。2005年樊庆琦等[4]采用多世代联合分析方法对普通玉米籽粒长度进行遗传分析，认为玉米粒长表现为主基因+多基因的混合遗传，只是在不同遗传背景的材料中控制籽粒长度的主基因数目表现出一定的差异。2011年Peng等[5]对2个F2:3玉米家系群体的研究结果显示，不同基因型玉米的粒长间存在着极显著差异，广义遗传力最高可达到0.89，表明该性状主要受到遗传因素的控制。2012年Macke[6]对12个美国玉米自交系及其杂交后代进行研究，认为玉米粒长性状的基因效应以加性为主并具有较高的遗传力，后代选择有效。所以在组配籽粒较深且产量高的玉米杂交种时，可选择粒长较长的玉米自交系。然而，控制玉米籽粒生长发育的因素较多且复杂，粒长发育形成的遗传机制尚不明确，还需进一步对其分子机理进行研究并找到相关的主效基因。

目前，科学家主要应用SSR和SNP等分子标记技术对玉米粒长性状进行QTL定位。2012年曹晓良等[7]以“农系531”和“SIL8”杂交所衍生的200个F2:3家系为作图群体，利用238对SSR标记共定位得到3个玉米粒长的QTL位点，最高可解释17.34%的表型贡献率。2014年Zhang等[8]利用SSR标记技术对243个永久F2代自交系进行遗传作图，经过两年两点的试验共检测出4个控制玉米粒长发育的QTL，这些位点分布在第3、第5和第10号染色体上。2016年Chen等[9]采用GBS简易基因组测序技术对173个重组自交系进行基因分型，在6个环境下共定位得到4个玉米粒长的QTL，解释的表型变异率在6.2%～8.0%。2016年Qin等[10]以Mo17和黄早四杂交所衍生的3个回交群体作为研究对象，利用SSR标记对之前得到的粒长QTL进行精细定位，结果显示，主效qKL1.07位点的遗传距离已缩小至1.6 Mb，为进一步挖掘相关的功能基因提供了依据。

尽管上述研究获得了大量控制玉米粒长性状的数量性状位点。但由于受到标记密度的影响，所检测到的QTL置信区间仍然较大，给后续功能基因的发掘带来了很大困难。所以，本研究采用第2代测序技术对80个核心玉米自交系进行基因组重测序，利用得到的高密度SNP标记对玉米籽粒长度的全基因组进行关联分析，筛选与粒长显著关联的SNP，深度挖掘玉米种质中控制粒长发育的优异等位基因，为分子标记辅助选择籽粒长、产量高的玉米材料及克隆相关的功能基因奠定基础。

1 材料与方法

1.1 试验材料与田间设计

本研究收集并筛选了80个核心玉米自交系作为关联群体，主要包含了在吉林省玉米育种中应用较多的自交系和近5年来在东北推广面积较大的玉米杂交种亲本。这些自交系的名称及系谱来源见表1。在2014年和2015年分别将其种植在吉林省长春市和梅河口市。采用完全随机区组设计，3行区，行长3 m，行距0.65 m，株距20 cm，设3次重复，两年的田间管理相同。

表1 80个核心玉米自交系的名称及系谱来源Table 1 Names and pedigree sources of 80 elite inbred lines

表1(续) Contiued table 1

1.2 表型鉴定

当玉米达到完熟期后每小区随机摘取10个果穗，风干后考种。按照以下公式计算玉米籽粒长度：粒长(cm)=(穗粗-轴粗)/2。其中穗粗和轴粗均为果穗中部的测量结果，每个果穗重复测量3次。各小区的粒长平均值用于后续的统计分析。

1.3 表型数据统计分析

使用SPSS 19.0(IBM Corp.,Armonk,NY,USA)软件分别对粒长数据进行描述性统计分析、方差分析和相关分析。按照Knapp等[11]提出的公式计算广义遗传力。

1.4 全基因组重测序及基因分型

采用康为世纪植物基因组DNA提取试剂盒提取80个玉米自交系叶片基因组DNA。用1%琼脂糖凝胶电泳、NanoDrop1000检测DNA质量，合格的DNA样品用于构建基因组文库，并利用Illumina Hiseq PE150进行测序。将获得的高质量测序数据通过BWA软件[12]比对到玉米自交系B73的参考基因组RefGen_v3(ftp://ftp.ensemblgenomes.org/pub/plants/release-29/fasta/zea_mays/dna/)上，比对结果经SAMTOOLS软件[13]去除重复，并利用贝叶斯模型进行群体SNP的检测。采用最小等位基因频率≥0.05、SNP检测缺失率<0.1等过滤条件进行筛选，最终获得1 490 007个高质量的SNP，用于后续的全基因组关联分析。

1.5 全基因组关联分析

应用软件PLINK[14]，通过LD测量参数r2(r2≥0.1)计算关联群体全基因组的平均LD衰减距离。运用FarmCPU方法[15]对高质量的SNP标记与籽粒长度数据间进行全基因组关联分析，其中关联群体的亲缘关系与5个代表群体结构的最大主成分，分别随机选取10% 的SNP标记，运用GAPIT软件[16]进行计算。采用Bonferroni法[17]对多重假设检验得到的P值进行校正，用以降低假阳性关联结果，最终确定当-lgP≥6时，判定SNP标记与玉米粒长性状关联显著。在显著性SNP位点的LD范围内(本试验所用关联群体平均LD衰减距离为5.2 kb)对控制粒长发育的候选基因进行扫描。候选基因的注释功能及相关信息参照玉米遗传学及基因组学数据库(http://www.maizegdb.org/)和美国国立生物技术信息中心网站(http://www.ncbi.nlm.nih.gov/)。

2 结果与分析

2.1 玉米粒长性状的统计分析

如表2所示，在2014年和2015年不同地点种植玉米自交系的籽粒长度变异范围在0.56～1.14 cm，标准差分别为0.11，0.15，0.12和0.14,说明玉米籽粒长度的遗传变异较大。通过方差分析得到的F值显示,不同玉米自交系的粒长间达到极显著差异水平(P<0.01);基因型与环境互作方差为0.023(P<0.01)。利用Knapp等[11]提出的公式计算玉米籽粒长度的广义遗传力为83.67%，与前人的研究结果[2,5,18]相近，说明该性状可以有效地通过表型进行选择。Pearson相关分析表明，不同环境下的玉米籽粒长度间达到极显著正相关(P<0.01)。

表2 80份玉米自交系籽粒长度的统计分析Table 2 Statistics of kernel lengths of 80 elite maize inbred lines in two years

注：**表示在P=0.01水平显著。

Note: **indicates significance atP=0.01.

2.2 玉米粒长关联群体的连锁不平衡

利用重测序得到的1 490 007个高质量的SNP标记对玉米粒长关联群体全基因组的LD进行估计。如图1所示，当r2=0.1时，关联群体全基因组的平均LD衰减距离约为5.2 kb，该结果小于已有研究得到的LD长度[19-20]，说明本研究所用的关联群体具有较快的LD衰减速率。

图1 玉米粒长关联群体全基因组范围内的连锁不平衡衰减Fig.1 Linkage disequilibrium decay across whole genome of maize kernel length association panels

2.3 玉米粒长的全基因组关联分析

本研究采用FarmCPU方法对80个玉米自交系的籽粒长度进行全基因组关联分析，结果显示除了2014年在长春未检测到与目标性状显著关联的标记外，其余3个环境下共检测到16个与玉米粒长显著关联(P<10-6)的SNP(图2)。这些标记分布在第1、第2、第3、第7和第10号染色体上，解释的表型变异率在1.68%～27.61% (表3)，其中位于染色体框10.03的SNP解释的表型变异率最低，而位于染色体框3.04物理位置为74171556的SNP解释的表型变异率最高，达到27.61%。同样位于染色体框3.04，物理位置为110936416的SNP标记在2015年的2个地点均被检测到，且该标记在长春环境下具有最小的P值，说明其与籽粒长度关联最显著。

表3 与玉米籽粒长度显著关联的SNP位点(P<10-6)Table 3 SNPs identified to be associated with maize kernel length (P<10-6)

注：已有研究中定位的粒长QTL与本研究检测到的显著性SNP具有一致性。

Note: There exists consistency between significant SNPs and QTL interval reported in previous studies.

2.4 玉米粒长候选基因分析

在5.2 kb的连锁不平衡范围内，对16个与玉米粒长显著关联(P<10-6)的SNP位点进行扫描分析，共得到3个候选基因(表3),其中GRMZM2G018607和GRMZM2G034882的基因座内各包含1个显著性SNP，它们均位于基因的外显子区域且为同义突变。基因GRMZM2G018607的编码产物是一个ABCF4转运蛋白(ABC transporter F family member 4)，该蛋白质包含ABCF_EF-3和ABC_tran_Xtn 2种结构域，可能参与调节应激相关蛋白的翻译合成过程[23]。基因GRMZM2G034882预测编码1种WIP2互作蛋白(WPP domain-interacting protein 2)，这种产物作为植物特有的外层核膜蛋白，在将RanGAP锚定在核膜上及维护细胞核形态的过程中发挥作用[24-25]。基因GRMZM2G322641编码1个含DUF2346结构域的蛋白质，目前该蛋白家族成员的生物学功能还未知。

3 讨论与结论

3.1 关联群体的连锁不平衡

关联群体的连锁不平衡水平不仅决定了全基因组关联分析的定位精度，而且对所选用的标记数量、密度以及所采取的试验方案等都会产生较大影响。2013年刘化龙等[26]研究认为，种质资源的连锁不平衡与其地理起源密切相关，群体遗传进化上的差异最终导致了LD的差异。2013年Romay等[19]利用681 257个SNP标记对2 815份来自美国玉米种质资源库的自交系进行了连锁不平衡分析，结果显示当r2=0.1时，整个群体的LD衰减距离在10 kb左右，其中热带玉米种质的LD接近于整个群体的连锁不平衡水平，而其他SS、NSS和ExPVP群体的LD衰减距离则达到50 kb以上。2012年Wang等[20]以144份我国玉米杂交种亲本作为研究群体，采用45 868个SNP对其全基因组的连锁不平衡进行分析，认为关联群体的平均LD衰减距离在200 kb左右，其中衰减最慢的染色体甚至达到了750 kb。该结果远大于本研究得到的5.2 kb的衰减距离，这可能是由于前者所选择的亲本自交系群体中基因或染色体片段受到的选择强度较大所致；另外，本试验采用的核心玉米自交系数量小于上述研究，供试种质的遗传构成较复杂，这也可能是导致群体LD衰减较快的另一个主要原因。

3.2 全基因组关联

自2011年Tian等[27]将全基因组关联分析用于玉米茎叶夹角的定位研究以来，该技术在玉米主要病害、开花期和抗逆性等其他复杂性状的遗传解析中也得到了广泛应用[28-30]。本研究采用FarmCPU方法对由80个核心玉米自交系组成的关联群体进行了全基因组关联分析，在4个环境下共检测到16个与玉米粒长显著关联(P<10-6)的SNP标记。这些位点中有7个位于前人已定位的粒长QTL置信区间内[7,21-22]。其中位于染色体框7.02(物理位置为96607424和111201213)、10.04(物理位置为124464498)及3.04(物理位置为74171556)的4个显著性SNP与Chen等[21]定位的籽粒长度QTL具有一致性；位于染色体框7.02，10.03(物理位置为37333210)和10.04的4个显著性SNP,与陈哲等[22]利用176个F2:3家系为作图群体定位的结果相吻合。2012年曹晓良等[7]利用1个F2:3分离群体在umc1917-umc1603处定位到粒长的QTL，本研究检测到的位于染色体框1.05的2个显著性SNP也在该区间的物理位置内。另外，2009年李永祥等[31]对玉米籽粒构型性状与产量性状进行相关性分析得出，玉米粒长与单穗产量的相关系数最高，达到了0.64，说明籽粒长度对玉米产量具有较大的贡献。所以，本研究对这16个显著性SNP与前人定位的产量性状QTL[32-36]进行比对，发现除了位于染色体框3.08及10.03的3个SNP外，其余的标记位点均位于已有产量性状的QTL置信区间内，说明控制这2个性状的基因可能存在连锁关系，或者在这些染色体热点区域有控制玉米产量及籽粒长度的多效基因存在。

3.3 候选基因分析

近年来，尽管研究者利用不同的分子标记技术定位了较多的控制玉米籽粒构型性状的数量性状位点[7-10]，但从这些位点中挖掘控制玉米粒长的功能基因还未见报道。本研究对16个与玉米籽粒长度显著关联(P<10-6)的SNP位点进行扫描分析，共得到3个候选基因。基因GRMZM2G018607编码一个ABCF4转运蛋白，该产物与ABC转运蛋白家族的其他成员有所不同，仅位于细胞质中，蛋白质无跨膜区域，因而不具有转运功能[37]。2011年Zeng等[23]对scord5和ila-3 2种拟南芥突变体进行研究得出，ABCF家族蛋白AtABCF3是细菌触发气孔关闭反应的重要组成部分，它可能在气孔关闭的早期反应通路中通过调节应激相关蛋白的翻译合成而发挥作用。基因GRMZM2G034882预测编码一种含有SMC_prok_A结构域的WIP2互作蛋白，它通过与特定的内层核膜蛋白进行互作，从而将其本身及RanGAP固定在核膜上，之后RanGAP催化RanGTPase对结合的GTP进行快速水解，精确调控核质间RanGTP和RanGDP的浓度梯度，有效促进细胞内的核质运输、纺锤体组装、核膜重建和异染色质的组装等[38-39]。2012年Zhou等[25]对3种wip拟南芥突变体的叶片表皮细胞进行DAPI荧光成像观察得出，与野生型植株相比，突变体叶片表皮细胞的细胞核长度显著降低，说明WIP2互作蛋白在维持细胞核特定形态的过程中也扮演着重要角色。基因GRMZM2G322641编码一个含有DUF2346结构域的未知功能蛋白。通过对其cDNA序列进行BLAST分析，结果显示在水稻、高粱和谷子等主要农作物的基因组中均发现GRMZM2G322641的直系同源基因，说明该基因的存在比较普遍，其生物学功能及参与的代谢途径还需要进一步验证研究。

候选基因的鉴定与所用的关联群体和多态性标记密度有直接关系。在不同的关联群体中等位基因频率有所差异，这种差异可能导致在特定基因组区域中多态性位点的不一致，从而使关联分析的结果不同[40]。另外，随着高通量测序技术的快速发展，应用不同测序技术获得的多态性信息含量差异较大，较高的标记密度在解析玉米复杂性状的遗传结构时具有更好的分辨率。所以，本研究采用第二代测序技术对由80个核心玉米自交系组成的关联群体进行重测序，经过全基因组关联分析共得到3个候选基因，其可能与玉米籽粒长度紧密相关，为进一步研究玉米粒长的遗传机理奠定了前期基础。