张素君，唐丽媛，李兴河，王海涛，刘存敬，张香云，张建宏

(河北省农林科学院 棉花研究所，农业部黄淮海半干旱区棉花生物学与遗传育种重点实验室，国家棉花改良中心河北分中心，河北 石家庄 050051)

陆地棉(GossypiumhirsutumL.)是世界上最重要的天然纤维作物，也是我国重要的经济作物之一[1]。黄萎病是一种严重的土传病害，寄主范围广泛[2]。目前，国内缺乏高抗黄萎病的陆地棉品种，且没有有效治愈黄萎病的杀菌剂，使得黄萎病成为危害棉花的最大病害之一[3-4]。2009-2010年，超过我国棉花种植面积的50%受到黄萎病的影响(美国国家棉花委员会-疾病数据库)[5]。培育和推广抗病品种是防治棉花黄萎病的有效方法[3]。

由于当前陆地棉品种遗传背景狭窄，且抗病性和产量、纤维品质之间存在一定程度的遗传负相关性[6]，实现三者的同步改良是一项艰巨的考验。分子数量遗传学发展为实现遗传改良提供了有效的工具，通过挖掘抗病性相关的QTLs(Quantitative trait loci)应用于分子辅助育种可加速选育抗病品种的进程。

为了挖掘抗病性相关的QTLs，不同研究者利用2个不同黄萎病抗性的棉花材料构建分离群体，通过连锁作图定位到众多与抗黄萎病性状相关的QTLs[7-13]。同时，随着遗传学的发展，基于连锁不平衡的关联分析在植物基因组学研究领域取得了成功的应用[14]。与连锁作图相比，全基因组关联分析(Genome-wide Association Study,GWAS)具有精确度高、不需专门构建群体、同时分析同一位点的多个等位基因等优点[15-17]。

近年来，关联分析方法已开始应用于棉花黄萎病抗性的研究。郭志军等[18]采用 74 个 Simple sequence repeat(SSR)标记对 172 份陆地棉栽培种的基因组变异进行扫描，通过关联分析发现 12个与黄萎病抗性显著相关(P<0.05)的位点；李黎贝[19]利用140 个SSR标记在 186 份陆地棉材料中共检测出 355 个等位变异，挖掘与棉花黄萎病抗性状相关联的位点22个(P<0.01)。Zhao等[20]搜集了世界各地的158份陆地棉资源，在全基因组内均匀筛选212个多态性较好的SSR标记，在成株期和苗期对陆地棉的黄萎病抗性进行了鉴定，利用关联分析挖掘了42个与黄萎病抗性关联的位点，其中有30个为首次挖掘，此外，本研究在16号染色体上发现了与黄萎病抗性相关的QTL簇。目前，通过关联分析挖掘陆地棉黄萎病抗性相关位点取得了一定的进展，但是与棉花其他性状比如纤维品质、产量相比[1,6,21-24]，相关研究尚不充足。

本研究通过分析214份陆地棉资源的遗传多样性，利用黄萎病病圃对其多年发病情况进行调查，利用均匀分布于棉花基因组上的237个SSR标记进行全基因组扫描，发掘与棉花黄萎病抗性显著关联的标记位点，旨在为棉花抗黄萎病性状的分子检测及标记辅助选择育种提供参考信息。

1 材料和方法

1.1 试验材料

214份陆地棉资源均由河北省农林科学院棉花研究所分子育种室提供，包括我国陆地棉栽培种(107个)、搜集于国内的资源(55份)、国外的资源(7份)、本实验室自选的育种材料(45个)。这些材料来源广泛，表型性状比如株高、纤维品质、熟性、衣分、产量、叶色、株形等方面差异显著。2013-2015年将上述材料种植于河北省农林科学院棉花研究所小安舍试验站田间病圃(河北石家庄)，该田间病圃是多年植棉又经过人工改造的黄萎病病田，经多年田间观察统计，具有与黄萎病病圃相当的致病效果。试验采用3行区，行长7.5 m，行距0.8 m。定苗时每行留27株左右。采用局部随机排列，设3个重复。组内每隔10个品种设感病对照品种(冀棉11)一次。

1.2 田间黄萎病抗病性鉴定

黄萎病鉴定采用田间成株期鉴定，参考国家棉花品种区试中抗黄萎病的鉴定方法[25]。根据黄萎病发病高峰期分别于每年8月底至9月初进行田间鉴定，每重复取中间行20株调查黄萎病发病情况，感病对照为冀棉11，采用3个重复相对病情指数(Relative disease index，RDI)平均值表示各个材料的抗病性强弱[25]。

1.3 棉花总DNA的提取与检测

CTAB法提取基因组DNA[26]，采用1%琼脂糖检测DNA完整性和杂质，采用分光光度计Nanodrop 2000检测，稀释至终浓度25～40 ng/μL。

1.4 PCR体系与电泳检测

本试验从6 000余对SSR引物中选取237个扩增效果稳定、多态性好且分布在棉花26条染色体上的SSR标记用于群体扩增，引物在染色体上的位置参考Liang等[27]构建的连锁图谱。所有引物信息均可从CottonGen(http://www.com.cottongen.org/)和CMD(http://www.cottonmarker.org/)查取。PCR(Polymerase chain reaction)反应体系10 μL，模板DNA 1.2 μL，天根2×TaqRCR Master Mix 5 μL，上下游引物各0.5 μL，ddH2O 3 μL。PCR反应程序为：95 ℃预变性4 min；94 ℃变性45 s，55 ℃退火45 s，72 ℃延伸1 min，30个循环；72 ℃ 10 min，8%的垂直板非变性聚丙烯酰胺凝胶电泳分离，参考张军等[28]方法银染。

1.5 数据统计与分析

人工读带，同一位置上清晰且重复性好的记为“1”，无带记为“0”，带弱或模糊不清的记为缺省值；利用Powermarker 3.25[29]计算引物的多态性位点、多态性信息含量和多样性指数，并利用其Phelogeny功能，计算Nei氏遗传距离，构建群体材料的Neighbor-Joining(NJ)聚类图；利用Structure 2.2 软件[30]估测该自然群体遗传结构，计算获得Q值。K值为1～10，将MCMC(Markov chain monte carlo)开始时的不作数迭代(Length of burn-in period)设为10 000次，采用Evanno等[31]提出的ΔK值方法确定合适K值；利用Tassel 2.1软件的一般线性模型(General linear model，GLM)进行关联分析，并计算各标记对表型变异的解释率[32]。

2 结果与分析

2.1 陆地棉资源的黄萎病抗性统计

对214份陆地棉资源进行抗病性鉴定表明，陆地棉在黄萎病抗性上表现出广泛的变异。在田间病圃鉴定条件下，2013年的RDI变化为5.6～91.6，平均为30.6，2014年的RDI变化为1.5～98.3，平均为26.0，2015年的RDI变化为0～90.0，平均为33.3(图1)。方差分析表明：相同年份3个重复内的相对病情指数(RDI)不存在显著差异(P>0.5)，而不同材料间RDI存在极显著差异(P<0.01)。不同年份之间进行比较发现，2014年黄萎病整体发病情况与其他2年相比较轻，年份效应显著。3年均表现抗病(RDI<20)的材料有18份，其中，本实验室自选的育种材料6份(抗8、冀棉18选系、冀选系258、冀选系120、Z56-1系、JMX-20)，栽培品种11份(冀棉958、中棉所69、冀151、冀122、新陆中44、德利农5号、仁和39号、中棉所40、冀优768、冀棉616、冀228)，资源材料1份(马里4号)。

2.2 SSR分子标记多态性信息

利用PowerMarker v3.25计算软件进行遗传多样性分析，237个SSR引物共检测到280个多态性位点，涉及695个等位变异，变异范围 2～6 个，平均等位变异数为2.479 0个(表1)。JESPR153、BNL169和BNL3442检测到的等位变异数最多为6个，多态性信息含量(PIC)变幅为0.004 9(NAU3816)～0.610 8(JESPR153)，平均0.213 3；基因多样性指数变异为0.004 9(NAU3816)～0.665 5(JESPR153)，平均为0.255 8。SSR引物在所选棉花品种中检测的等位变异数目和基因多样性的跨度较大，但平均值较低，说明所选择的棉花材料在基因组水平上变异比较丰富，但整体上所选陆地棉的遗传基础比较狭窄。

图1 田间病圃鉴定条件下供试材料病指变化频率分布Fig.1 Histogram of relative disease index of tested lines identified in field disease nursery

表1 SSR 标记多态性信息统计Tab.1 Statistics of SSR marker polymorphism information

2.3 聚类分析

利用PowerMarker v3.25软件中的Phelogeny功能，根据280个具有多态性位点分子标记数据，计算Nei氏遗传距离；并基于该遗传距离，构建群体材料的Neighbor-Joining(NJ)聚类图。采用Nei氏距离对该自然群体进行聚类分析(图2)，可以分为3个亚群，分别包含8，79，127份材料。此分类与材料类型无对应关系。8份材料中有5份资源材料、3个品种，其中，7份材料在群体结构分析中划分为第1亚群；79份材料中有33份育种中间材料、31个品种和15份资源材料，在群体结构划分中均被划分为第1亚群；127份材料中有73个栽培品种、42份资源材料和12份育种中间材料，在群体结构划分中有70份划分为第1亚群，57份材料划分为第2亚群。聚类结果与材料的地域来源无显著对应关系，说明不同育种单位之种质流动性较大，打破了种质的地域界限。

图2 214份陆地棉资源基于SSR标记的聚类Fig.2 Clustering map of 214 upland cotton based on SSR markers

2.4 群体结构分析

利用 Structure 2.2软件中的混合模型聚类分析法对214份陆地棉的群体结构进行划分，在K=1～10时，随着K值的增大，模型的后验概率(Lnp Posterior probability，Ln P(D))也持续增大，因而，无法依据最大似然值的原则确定合适的K值。参照 Evanno 等[31]的方法计算ΔK值，ΔK在K=2时出现顶点，推测该自然群体分为2个亚群(图3),2个亚群分别包含156，58份材料。此结果与利用NJ聚类方法分析结果略有差别，但共同表明供试群体结构简单，来源较为单一，在一定程度上反映了陆地棉品种遗传基础狭窄的问题。

图3 214份陆地棉材料的群体结构Fig.3 The population structure of 214 upland cotton germplasm

2.5 陆地棉SSR位点间的连锁不平衡分析

利用Tassel 2.1进行连锁不平衡分布(LD)分析，共有204个标记覆盖了整个棉花基因组2 918.7 cM，标记间的平均距离为14.38 cM，平均每条染色体上有7.85个位点，另有33个标记尚未定位到染色体上。237个标记在供试群体中检测到了280个多态性位点，这些位点间的连锁不平衡分布如图4所示。在280个位点的39 060种位点组合中，连锁不平衡位点广泛存在于共线性组合和非共线性组合中，概率统计极显著支持(P<0.01)的成对位点占总组合数的3.73%(图4)，所选陆地棉基因组内的连锁不平衡水平较低。

2.6 关联分析

通过关联分析，在不同年份中发掘与黄萎病抗性显著关联(P<0.01)的SSR位点27个(表2)，其中有2个位点(BNL3442和BNL1064)在3年均被重复检测到，表型变异解释率最高分别为12.10%和8.02%。3个位点(Gh433-2、NAU3607和JESPR153)在2013,2014年同时检测到与黄萎病抗性显著关联；6个位点(BNL3902、NAU5428、BNL1646-1、BNL1646-2、CGR5996和HAU1774-1)在2014,2015年同时检测到与黄萎病抗性显著关联，其中，NAU5428的2年表型变异解释率较高，分别为6.93%和5.57%；其余16个位点仅在一个年份检测到与黄萎病抗性相关联。这些在不同年份中稳定存在的SSR标记位点有可能与抗病基因紧密连锁，有望用于棉花黄萎病抗性材料筛选和抗病基因挖掘中。

图4 SSR位点间的连锁不平衡Fig.4 The linkage disequilibrium between SSR loci

表2 与黄萎病抗性相关的SSR位点Tab.2 SSR loci associated with Verticillium wilt-resistance

本研究通过对抗黄萎病性状进行关联分析的同时，挖掘出一些与棉花抗黄萎病稳定关联的SSR位点(BNL3442和BNL1064)，可以为抗黄萎病候选位点做进一步验证，并挖掘优异等位变异，对棉花抗黄萎病分子育种具有重要的理论和实践意义。同时，本研究通过3年陆地棉资源的黄萎病抗性统计，鉴定了一批抗病性突出的育种材料和资源，为棉花抗病育种提供了参考。

3 结论与讨论

3.1 黄萎病鉴定方法

目前，对于棉花黄萎病抗性的鉴定，主要有田间病圃鉴定、营养钵苗期鉴定和黄萎病毒素法鉴定[33]。其中，田间病圃鉴定是棉花种质抗黄萎病鉴定常用且能真实反映鉴定材料抗性强弱的方法，但由于大田鉴定存在不同材料成熟期早晚和接种量的差异，以及试验小区的干扰等，存在不同年份、不同棉区筛选的抗黄萎病材料的抗病性差异较大的情况。通过设置合适的感病对照品种，利用相对病情指数来描述棉花的抗病性是减少环境变异的有效措施[4,34]。本研究中，使用冀棉11同一批种子作为感病对照，最大程度保证了标准品种年度间抗病性的一致。本研究结果表明，尽管努力减少试验的环境误差，仍然存在部分材料在不同年份间鉴定结果差异较大的情况，可能是由于陆地棉黄萎病抗性属数量性状遗传[35]，容易受到环境变化的影响。通过上述鉴定方法筛选出3年均表现为抗黄萎病的材料18份，这些材料中可能含有多环境稳定的抗黄萎病优势等位基因，需要在分子水平进行进一步挖掘。这些材料经其他鉴定方法重复检验后，有望作为抗黄萎病材料用于棉花育种工作中。

3.2 群体结构分析

群体结构是影响关联分析的一个重要因素，研究一个群体的遗传结构组成时，会因为亚群的混合导致整个群体的连锁不平衡水平增强，假阳性结果增加[36]，如果群体结构比较简单，假阳性关联结果出现的可能性就会减少[37]。因此，在进行关联分析前对自然群体进行遗传多样性和群体结构分析很有必要。本研究中，聚类结果显示，所选群体被分为3个亚群，群体结构分析表明，供试材料被分为2个亚群，2个聚类结果差异并不大，共同说明了所选陆地棉整体上遗传背景狭窄。这也解释了本研究中SSR引物在所选棉花品种中检测的等位变异数目和基因多样性的跨度较大，但平均值较低的问题。由于群体结构较为简单，本研究仅利用GLM模型进行了关联分析。

3.3 黄萎病相关分子标记的获得及利用

本研究通过关联分析获得了与黄萎病性状显著关联的SSR标记位点27个，与前人黄萎病抗性相关关联分析结果进行比较[18-20]，Zhao等[20]利用212个SSR标记与陆地棉成株期和苗期的黄萎病抗性进行了关联分析，鉴定出NAU5428与黄萎病极显著相关(P<0.000 1)，解释表型变异率达10.93%，本研究中，此标记在2年田间鉴定中与黄萎病显著相关，解释表型变异率分别为6.93%和5.57%，此结果表明了本研究结果与前人研究结果的一致性。

随着黄萎病抗性相关的QTLs的不断挖掘，已经有一些分子标记陆续应用于棉花抗病基因分子标记辅助选择育种上[4,38-41]。王芙蓉等[38]研究发现，标记 NAU751 和BNL1395 的抗性基因型聚合后可显著提高后代的抗性水平；孔祥瑞等[39]研究发现，BNL1721、BNL2733和 BNL3452这 3个标记的黄萎病抗性选择效应能够在不同世代间稳定遗传；祁伟彦等[40]利用和黄萎病抗性紧密连锁的 SSR 标记对棉花黄萎病抗性进行了辅助筛选，提出将人工病圃筛选和分子标记辅助育种相结合是选育棉花抗黄萎病材料可行、高效的育种方法。赵云雷等[4]研究表明，抗黄萎病优异等位基因的累加具有明显提高抗病性的作用，材料抗病性的强弱在一定程度上可以通过优异等位基因数目和效应值之和来反映，这就为分子方法鉴定抗病性强弱提供了依据，同时为黄萎病关联位点利用提供了重要参考。

通过关联分析挖掘了多个与陆地棉黄萎病抗性关联的分子标记位点，为棉花抗黄萎病性状的分子检测及标记辅助选择育种提供信息。今后应进一步挖掘与棉花抗黄萎病性状显著关联的优异等位基因位点，逐步实现抗病基因型的直接选择，有效缩短抗病性的鉴定周期。