陆地棉抗黄萎病性状与SSR标记的关联分析

王龙，李黎贝，马启峰，耿洪伟，陈全家，曲延英，范术丽

(1.新疆农业大学农学院，乌鲁木齐 830052；2.棉花生物学国家重点实验室/中国农业科学院棉花研究所，河南安阳 455000)

0 引 言

【研究意义】棉花黄萎病是以大丽轮枝菌浸染为主的维管束病害，对棉花的产量和纤维品质等农艺性状有很大影响。1914年在美国发现棉花黄萎病，随后在澳大利亚等主要产棉国也陆续发现了棉花黄萎病[1]。1935年，棉花黄萎病由美国传入中国，20世纪90年代，棉花黄萎病在我国主要产棉区逐年蔓延，其中西北内陆植棉区、黄河流域植棉区、长江流域植棉区等地的范围不断扩大。棉花黄萎病成为阻碍棉花实现优质、高产、稳产等育种目标的主要因素，已成为提高棉花的产量、纤维品质的主要障碍[2]。【前人研究进展】有研究认为，棉花黄萎病抗性由单基因控制的质量性状[3-5]。多数学者认为，棉花黄萎病抗性由多个微效基因共同控制的数量性状[6-9]。产生两种结论的原因，可能是研究人员所选不同供试材料而造成的，并且棉花抗黄萎病性状受环境影响较大。棉花黄萎病抗性QTL(Quantitative trait locus)定位在国外的相关研究较少，国内相关研究较多；SSR分子标记已成为研究黄萎病的主要技术，研究群体由先前的陆陆群体转变为海陆群体，因为后者的多态性较好。作图群体也由初级的F2、BC1群体向更适合精细定位的RIL(Recombinant inbred lines)、NIL(Near isogenic lines)等群体过度。高玉千等[10]利用海陆组合的F2群体、通过RAPD(Random amplified polymorphic DNA)和SSR(Simple sequence repeat)标记找到2个主效QTLs和一个微效QTL；王红梅等[11]利用F2/F2∶3群体定位到6个与黄萎病抗性相关的QTLs；昝伟等[12]研究结果表明：在F2群体中，通过SSR分子标记检测到与棉花黄萎病抗性相关的2个主效QTLs，并且定位到第一连锁群上；李磊[13]研究结果表明：在构建的两个F2群体上进行定位，找到6个与黄萎病抗性相关的QTLs，其中包括3个抗病QTLs，3个感病QTLs。目前，关联分析已经成为一项技术体系比较完善、相对比较成熟的研究功能基因定位的一种方法。关联分析要比常规QTL定位技术优势更明显，(1)以自然群体为材料，无需构建专门的作图群体，而QTL作图至少需要花费2年以上的时间来完成群体的构建；(2)可实现对其作图群体(自然群体)一个基因座上所有等位基因的考察，而绝大部分QTL作图所利用的群体是两亲本杂交、重组自交后代，每一基因座一般只能涉及两个等位基因；(3)定位精度高，可以达到单基因水平。关联分析的方法在玉米、大麦、大豆、小麦等多个作物的一些功能基因定位研究中得到广泛应用[14-24]。近年来，通过关联分析定位与棉花农艺性状相关功能基因的报道已有一些。Abdurakhmonov[18-19]利用285份陆地棉材料构建自然群体，并结合使用202个SSR标记，研究了该群体纤维品质性状的关联分析，研究结果表明：与纤维长度显著关联的标记有23个，与纤维比强度显著关联的标记有19个。贺道华等[20]用92份栽培棉构建自然群体，对其进行扫描，并结合均匀分布于26条染色体上的132个 SSR标记，获得21个与纤维品质性状相关的SSR位点。梁冰等[21]构建68份陆地棉材料为自然群体，以160对SSR引物对该群体进行基因型检测，发现了与棉花纤维品质等主要农艺性状相关的SSR位点。Kantartzi等[22]以56份亚洲棉品种构建自然群体，通过98个SSR标记对该研究群体的纤维品质性状进行关联分析，研究结果表明：与纤维长度显著关联的标记3个，与纤维比强度显著关联的标记有2个。Wang等[23]以55个海岛棉品系构建自然群体，使用258个SRAP标记、170个SSR标记与该研究群体的产量、纤维品质等性状进行关联分析， 研究结果表明，与产量性状相关联的标记有26个，与纤维品质性状相关联的标记有46个。【本研究切入点】近年来，关于研究棉花黄萎病抗性方面的关联分析报道较少。Zhao等[24]构建自然群体以158个陆地棉材料为研究对象，利用与棉花抗黄萎病性状相关的212个SSR标记进行关联分析，研究结果表明，与黄萎病抗性相关联的标记位点有42个，其中新发现的有32个标记位点。郭志军等[25]构建自然群体以178份材料为研究对象，利用74对SSR引物与该群体进行关联分析，检测到位于NAU1514- NAU2754区间内与棉花抗黄萎病性状相关的1个QTL。【拟解决的关键问题】以186份陆地棉材料组成的自然群体为研究对象，对其抗黄萎病性状进行4个播期的表型鉴定，利用筛选出的140个SSR标记进行基因分型，分析群体结构和品种间亲缘关系，利用TASSEL4.0软件的广义线性模型(General line model，GLM)(P<0.001)和混合线性模型(Mixed linear model，MLM)(P<0.01)软件，对棉花农艺性状和标记进行关联分析，发掘棉花农艺性状显著关联的SSR标记，为棉花抗黄萎病分子标记辅助选择育种提供依据。

1 材料与方法

1.1 材 料

试验以186份陆地棉材料构建自然群体。表1

表1 186份陆地棉材料来源
Table 1 The origin regions of 186 upland cotton accessions

中国农业科学院棉花所试验基地位于安阳市白璧镇，将供试材料种植在试验基地，试验设计为2年(2015和2016年)、2个不同播期(春播和夏播)一共4个环境的种植模式。186份供试材料播种期分别标记为 M(2015年4月26日)、 W(2015年5月23日) 、 Y(2016年5月1日)和 Z(2016年5月25日)，一共4个环境。采用随机区组排列，设置3次重复，行长6 m，行距0.7 m，株距0.3 m，常规大田农事管理。

在棉花黄萎病发病高峰期进行黄萎病表型病级鉴定，分别为2015年8月1日、2016年8月2日进行，试验地棉花黄萎病的表型病级鉴定采用5级制的划分方法[26]。表2

表2 棉花叶片黄萎病级划分标准
Table 2 The classification standard of cotton leaf verticillium wilt

1.2 方 法

1.2.1 SSR标记

以棉花幼苗八叶期叶片为试验材料，提取DNA参考宋国立等[27]改进的CTAB法。选取覆盖陆地棉26条染色体的140个SSR标记位点，对186份供试材料进行基因型检测，在该群体中检测到355个等位变异。所用SSR引物由上海生工生物工程股份有限公司合成。

PCR反应体系为10 μL： 包括模板DNA 30 ng，10×PCR Buffer 1.0 μL，10 mmol/L dNTP 0.2 μL，10 mmol/L SSR引物2 μL，5 U/μLTaqDNA聚合酶0.1μL，加ddH2O补足至10 μL。试验所用的Taq聚合酶和dNTP均由北京全式金生物技术有限公司提供。PCR反应程序：95℃预变性5 min；30个循环(94℃变性30 s，58℃退火温度45 s，延伸45 s)；72℃延伸10 min；10℃保存。使用浓度为8的聚丙烯酰胺凝胶电泳对PCR扩增产物进行分离，待跑完电泳以后，银染方法依据张军等[28]的方法进行并照相保存。

1.3 数据处理

带型统计采用0、1统计法，在同一迁移率位置上，有带记为1，无带记为0，缺失记为‘？’。根据DNA Marker结合Quantity One软件的方法，以此来判断多态性条带的大小。

利用Power Maker V3.25软件计算引物的多态性信息含量(Polymorphism information Content，PIC)[29]。

参照Henderson[30]提出的BLUP(Best Linear Unbiased Prediction)法，通过R语言中的“lme4”包进行四个环境下的育种值计算。

对186份陆地棉材料构建的自然群体基于数学模型进行的类群划分，通过使用 STRUCTURE2.3.4分析软件[31]，进一步对该群体的群体结构进行估算，并计算得出的186份研究群体中每个个体所相应的Q值。使用该软件程序设定K为1至10，将MCMC(Markov Chain Monte Carlo)开始时的不作数迭代(Length of burn- in period)设为10 000次，再将不作数迭代后的MCMC设为100 000次，参照文自翔等[32]的设置对程序进行处理，并进行10次重复计算[24]。通过SAS运算出最大似然值，根据运行结果确定合适的K值，若似然值LnPD随K值的增大而持续增大，则可参照Evanno等[33]的方法根据ΔK确定合适的K值。

通过使用TASSEL4.0软件模型，利用广义线性模型(General linear model，Q模型)是以个体Q值作为协变量，而混合线性模型(Mixed liner model，Q+K模型)是以Q值亲缘关系(Kinship)作为协变量，运用两种模型对4个环境下抗黄萎病表型数据与SSR位点进行回归分析。利用SPAGeDi-1.3[34]软件对亲缘关系进行计算。通过运用Bresghellohe和Sorrells等[35]提出的“无效等位变异”方法，分析了棉花抗黄萎病表型效应值情况。

2 结果与分析

2.1 田间鉴定

供试材料种植在中国农业科学院棉花所试验基地，在棉花黄萎病发病高峰期进行黄萎病表型鉴定，分别为2015年8月1日、2016年8月2日进行，黄萎病的表型病级鉴定采用马存[26]提出的5级制划分方法。图1

图1 田间棉花叶片病级划分标准
Fig.1 The standard phenotype of cotton leaf in classifying verticillium wilt resistance in the field

2.2 陆地棉材料的群体结构

研究186份供材料的群体结构，通过使用Structure version2.3.4软件对其进行分析，结果显示：LnP(D)值随K值的增大而增大，发现没有一个拐点，因此，无法确定该研究群体结构的K值；随后采用通过ΔK法来确定群体结构的K值[33]。当K=2时，平均似然值的斜率改变，将186份研究材料分成2个亚群。一个亚群有96份材料，占总数的51.60% ；另一亚群有90份材料，占总数的48.40%。

利用Structure软件，计算出当K=2时，可将186份研究材料划分为2个亚群，得到对应的Q值，为关联分析协变量的分析提供依据。图2

图2 186份陆地棉材料的群体结构
Fig.2 The population structure 186 cultivated upland cotton

2.3 单个环境下关联

通过使用140对SSR引物做分子标记，在该群体中检测到355个等位变异，平均每对引物的等位变异为2.54，其中平均值为0.76，引物多态性信息含量变化范围：0.50～0.99；将186份供试材料两年、两个播期共四个环境的棉花黄萎病指数，结合140个SSR标记，通过GLM(P<0.001)模型和MLM(P<0.01)模型同时进行关联分析，在检测到的355个位点中，同时能在两个模型下检测到且与黄萎病抗性相关的位点有22个。有3个标记位点与M播期相关联，有6个标记位点与W播期相关联，有3个标记位点与Y播期相关联，有11个标记位点与Z播期相关联。其中，NAU998和CGR5258两个位点不仅与Y播期、Z播期的黄萎病病指相关联，还分别于M播期、W播期的黄萎病病指相关联。NAU998和CGR5258两个位点表型变异解释率分别为5.53%和12.07%。表3

表3 陆地棉抗黄萎病性状关联标记的贡献率及P值

2.4 多点环境下关联

将两年、两个播期下共四个环境的表型数据构建BLUP(Best linear unbiased prediction)模型，同时进行该模型与单个环境下的GLM模型以及MLM模型下的关联分析。研究表明，在GLM模型下关联到14个位点，在MLM模型下关联到3个位点，最终发现两个位点同时存在于2个环境以上，分别为NAU998和CGR5258，该结果与单个环境下关联分析的结果相吻合，说明这两个标记可能在棉花黄萎病抗性的分子鉴定方面有较好的效果，对棉花分子标记辅助育种有一定的应用价值。图3

图3 关联分析流程
Fig.3 The flow chart of association analysis

2.5 抗黄萎病性状表型效应值

与棉花抗黄萎病性状相关的两个位点，分别为NAU998和CGR5258。依据文自翔[32]提出的“无效等位变异”判断方法，得出SSR位点对表型效应的影响。位点NAU998在四个环境下的表型效应分别为0.03/0.08/-0.005/0.033，增效效应更为显著，说明这个位点的存在有利于提高棉花黄萎病的抗病性；位点CGR5258在四个环境下的表型效应都表现为-0.001/-0.001/-0.001/-0.001，均为减效效应，这个位点的存在会使棉花黄萎病的感病性提高。在育种过程中，可选择具有NAU998位点同时缺失CGR5258位点的中间材料，以提高选育材料的黄萎病抗性。表4

表4 抗黄萎病性状的表型效应
Table 4 Phenotypic effects of resistance of Verticillum wilt

3 讨 论

目前，有关棉花关联分析的研究大多数于集中于主要农艺性状等方面，然而关于棉花抗黄萎病性状的报道仍然较少[21,25]。陆地棉黄萎病抗性的遗传方式为数量性状，由多个微效基因共同控制，且致病菌有变化快的特点，而且棉田的发病易受湿度和温度等条件的影响。研究利用两年、两个播期的四环境数据，发掘与陆地棉抗黄萎病性状相关联的SSR位点，为培育棉花抗黄萎病新品种提供可靠的分子标记。在检测到与抗黄萎病性状相关联的22个位点中，两个环境中同时被关联到的位点NAU998和CGR5258，并且在前人研究结果中尚未出现。因此，可通过构建分离群体对这2个位点进行验证，并应用于改良棉花抗黄萎病育种实践。

研究在进行关联分析过程中运用了GLM和MLM两种模型，GLM通过群体结构Q值作协变量；而使用MLM模型在进行关联分析时，不但使用了Q值作协变量，还引用了亲缘关系K；因此，在一定程度上能降低关联结果的假阳性，使结果更为准确，无论是单个环境下的关联结果还是多个环境下的关联结果，MLM模型所关联到的位点始终少于GLM模型。研究在进行关联分析时使用了Henderson提出了BLUP法，即最佳线性无偏预测法，这种方法最初运用于动物育种当中，而20世纪80年代中期以后，BLUP方法开始应用于牧草、农作物及园林植物的遗传改良研究，在关联分析研究中用来构建模型处理多年多点的数据，可进一步消除多年多点对数据关联结果的影响，且已在大豆，玉米等[36-38]的关联分析研究中得到应用。

结合前人对棉花黄萎病抗性QTL定位的研究结果，利用e-PCR程序[39]，与前人研究结果相比，NAU998和CGR5258不能比对到Zhang等[40]构建的连锁群中。但是在最近的研究中，分别与抗病性QTLs位点qDI-C3-3的近距离标记MUSS59相近，与QTLs位点qDI-C15-3的近距离标记NAU3337相近[41]，且分别位于NAU3775与BNL0663和NAU3684与BNL3971之间[42]。出现这种差异的原因可能是由于定位的材料不同而导致的，NAU998和CGR5258标记的发现，将对棉花抗黄萎病基因的克隆和功能验证奠定坚实的基础[43]。

4 结 论

将四个环境下的表型数据进行BLUP之后进行关联分析，无论是GLM模型还是MLM模型所关联到的位点都少于单个环境下所关联到的位点，BLUP的精准性，研究得到的NAU998和CGR5258两个位点具有较好的可靠性。

近年来，有关棉花关联分析主要以SSR标记居多，传统的研究提供基因组多态性信息相对较少，未来的发展更趋向于通过基于简化基因组测序技术和重测序技术进行全基因组关联分析[44-46]。研究通过使用140对SSR做分子标记，检测到355个等位变异，进行关联分析，标记数量需要进一步加大，研究结果只是粗略意义上的全基因组关联分析， 因此，还需增加标记的密度，可以结合SNP标记，开展更加详实的棉花抗黄萎病GWAS研究；第二，要对已发掘的SSR位点进一步功能验证，通过构建分离群体进行标记与性状的相关性分析，可进行QTL分析，以此来提高相应位点的准确性。通过对棉花抗黄萎病性状的关联分析，发掘了一些与黄萎病抗性相关联的SSR位点，为陆地棉抗黄萎病材料进行全基因组扫描提供了更多的分子标记。