马铃薯资源晚疫病抗性的全基因组关联分析

蒋 伟 潘哲超 包丽仙 周福仙 李燕山 隋启君,* 李先平,*

1 云南省农业科学院经济作物研究所, 云南昆明 650205; 2 云贵高原马铃薯与油菜科学观测实验站, 云南昆明 650200

晚疫病是马铃薯的主要病害, 每年可造成我国马铃薯减产约10%～15%, 经济损失高达20亿美元。西南地区低温高湿的气候环境, 晚疫病发生更严重,马铃薯减产可达 15%～40%[1-2]。提高品种晚疫病抗性是马铃薯育种重要目标之一, 而挖掘持久、高抗的晚疫病抗性基因对马铃薯抗病育种具有重要的意义。

20世纪50年代, 育种家将野生种Solanum demissum的抗性基因导入到栽培种中, 但抗性很快被新出现的晚疫病生理小种克服[3-4]。目前已经鉴定并克隆到一些对晚疫病表现出广谱抗性的R基因。通过构建BAC文库和长片段PCR扩增的方法, 从野生种S. bulbocastanum中克隆了第1个马铃薯晚疫病广谱抗性基因RB或Rpi-blb1[5-6]。通过构建S. bulbocastanum作图群体, 利用AFLP分子标记, 克隆到另1个广谱抗性基因Rpi-blb2[4]。利用野生种S. venturii的 5个作图群体, 结合分离群体分组分析法(bulked segregant analysis, BSA)找到与抗性基因连锁的AFLP分子标记, 通过构建BAC文库从而克隆了广谱抗性基因Rpi-vnt1[7]。对含有广谱抗性基因R8分子标记或QTL位点的群体构建BAC文库, 通过酶切PCR克隆到了该基因[8-9]。上述这些基因的鉴定主要是通过构建遗传群体、分子标记定位获得, 而通过自然群体全基因组关联分析定位晚疫病抗性基因的研究较少。

本研究选择的自然群体引自国际马铃薯中心,该中心从19世纪80年代起以淘汰马铃薯晚疫病垂直抗性而保留其水平抗性为目的, 筛选了一批群体材料。因此, 该群体中很可能保留着晚疫病广谱抗性基因。通过对四倍体马铃薯群体连续2年的田间晚疫病抗性评价, 结合群体简化基因组测序获得大量的SNP位点, 进行全基因组关联分析挖掘晚疫病抗性相关的遗传位点和候选基因, 旨在为晚疫病抗性品种选育和抗病机理研究提供一定的理论和材料基础。

1 材料与方法

1.1 材料

引自国际马铃薯中心的288份四倍体马铃薯晚疫病抗性群体材料(表1), 感病对照为Desiree (D214)。群体A中晚疫病抗性来源为改良的S. demissum衍生抗性, 包括来自安第斯栽培种S. tuberosumgroupsandigena、phureja和stenotomum, 以及野生种S.acaule和S. bulbocastanum抗性的材料 8份; 群体B3来源于群体A, 但筛选保留表现出晚疫病水平抗性的资源材料 82份; 群体 B1来自S. tuberosumgroupsandigena的材料15份; 群体LTVR具有马铃薯病毒病(PVY、PVX和PLRV)抗性、生长期短和适应温暖环境特性的材料111份; 群体B3-HT结合群体B3的晚疫病抗性、北美和欧洲品种的耐热性以及群体LTVR特性的材料34份; 群体B3-LTVR包括群体B3和群体LTVR的杂交后代22份; 群体BW具有细菌性枯萎病抗性的材料 12份; 群体PREBRED具有从野生种导入到四倍体B3或LTVR中的晚疫病抗性材料1份; VARIETY表示马铃薯品种 3份[10]。

1.2 田间试验和表型数据分析

2015—2016年, 群体材料种植于云南省曲靖市会泽县试验基地(26º06'N, 103º22'E), 海拔 2650 m。每个材料种植1行10株, 行距为0.7 m, 株距为0.3 m, 种植3次重复。采用Alpha试验设计(CIP提供软件), 整个生长季进行正常中耕、施肥管理等, 但不喷洒杀菌类农药。

调查田间晚疫病发病比率。发现初始病株, 每7 d调查1次, 共调查8次, 计算AUDPC和sAUDPC值。利用Microsoft Excel对所观测到的数据进行统计分析。

式中,n表示总的调查次数,xi表示第i次调查的晚疫病严重程度,ti表示第i次调查的时间。

1.3 DNA提取和测序方法

利用 CTAB法提取 DNA, 检测合格后利用Illumina HiSeq2500测序仪(Illumina, Inc; San Diego,CA, USA)测序。以马铃薯S. chacoenseM6的基因组测序信息为参考[11], 利用北京百迈客生物有限公司自主研发的简化基因组测序流程 SLAF-seq[12]进行电子酶切预测, 根据最佳酶切方案选择RsaI+HaeIII进行酶切, 酶切长度在314～414的序列被定义为SLAF标签, 获得多态性标签1,032,778个。以马铃薯S. chacoenseM6为参考基因组, 利用BWA软件[13]将 reads比对到参考基因组上, 利用 GATK[14]和SAMtools[15]软件筛选SNP。根据完整度大于0.8, 次要基因频率大于0.05过滤后, 得到353,753个SNP位点进行后续的分析。

1.4 晚疫病抗性全基因组关联分析

基于SNP位点, 利用TASSEL软件[16]的混合线性模型(compressed mixed linear model, CMLM), 通过公式y= Xα+Qβ+Kμ+e得到关联值。公式中, 通过Admixture软件[17]计算样品群体结构 Q矩阵,SPAGeDi软件[18]计算样品间亲缘关系K矩阵, X为基因型矩阵, y为表型, 最终获得每个SNP位点的关联值, 并计算它们的P值。以 AUDPC_2015、AUDPC_2016、sAUDPC_2015、sAUDPC_2016, 以及这 2年的 AUDPC和 sAUDPC平均值即AUDPC_Mean和 AUDPC_Mean为表型数据, 共分析6个性状。对这些性状利用GLM、MLM、CMLM、FASTLMM和EMMAX五种模型进行关联分析。利用 GGplot2软件[19]绘制 Quantile-Quantile散点图(Q-Qplot), 利用 QQman软件绘制曼哈顿图(Manhattan), 展示关联分析显著的SNP位点。

1.5 候选基因预测

对6个性状在5种模型下获得的显著水平是0.1和0.01的关联SNP位点进行注释。选取这些位点前后100 kb范围内的基因, 利用BLAST软件分别与NR、SwissProt、GO、COG和KEGG数据库进行比对, 获得关联区域内基因注释信息, 并根据注释信息寻找与晚疫病抗性相关的候选基因。

表1 自国际马铃薯中心引进的288份马铃薯群体材料Table 1 288 potato population resources from CIP

(续表 1)

(续表 1)

2 结果与分析

2.1 晚疫病抗性的表型分析

通过对288份马铃薯资源材料连续2年的晚疫病田间抗性评价发现, 晚疫病抗性评价值 AUDPC和 sAUDPC呈连续分布。2015年所有材料的AUDPC平均值为2244.1, 范围是121.3～4740.8, 变异系数为 61.73%, sAUDPC平均值为 3.9, 范围是0.2～8.3, 变异系数为 61.56%; 2016年所有材料的AUDPC平均值为1228.5, 范围是72.7～3371.3, 变异系数为80.47%, sAUDPC平均值为3.6, 范围是0.2～9.5,变异系数为80.21%。表明该群体表现出晚疫病抗性广泛的遗传差异, 且受多基因控制, 属于数量性状遗传。相关性分析中, 2015—2016年AUDPC值相关系数为0.8988, 高度线性相关; sAUDPC值相关系数为0.8974,高度线性相关(图1)。表明2年间群体材料田间晚疫病抗性表现一致, 抗性稳定。

2.2 马铃薯晚疫病抗性的全基因组关联分析

5种模型下对6个性状分析的QQ图结果表明,在所有模型下关联到显著相关的SNP位点都是可靠的(图2是2个性状AUDPC_mean和sAUDPC_mean在 5种模型下的 QQ图)。在 5种分析模型下, 共关联到与6种晚疫病抗性性状显著相关(显著水平为 0.1和 0.01)的 SNP位点 82个(表 2和图 3)。其中, 6个SNP位点分布在3号染色体(图3-A, B); 3个 SNP位点分布在 4号染色体(图 3-A, B); 1个SNP位点分布在5号染色体(图3-B, C); 分别有2个SNP位点分布在6号、7号染色体(图3-A, B, D);66个 SNP位点分布在 9号染色体(图 3-A～D); 各有1个SNP位点分布在11号和12号染色体(图3-A,B, D)。在82个SNP位点中, 20个SNP位点在所有模型分析中, 2种显著水平情况下, 与所有性状都相关。它们是31,375,673 (第9个)、31,707,832(第 23 个)、32,382,888 (第 25 个)、32,973,114 (第35 个)、33,118,463 (第 36 个)、33,145,927 (第 38个)、33,148,740 (第 39 个)、33,156,262 (第 40 个)、33,211,574 (第 43 个)、33,244,698 (第 48 个)、33,283,731 (第 49 个)、33,283,949 (第 50 个)、33,451,096 (第 53 个)、33,480,937 (第 55 个)、34,959,324 (第 62 个)、35,060,347 (第 66 个)、35,065,836 (第 67 个)、35,123,523 (第 72 个)、35,212,534 (第 77 个)、35,256,410 (第 80 个)。同一性状下, GLM模型分析定位到的SNP位点数目最多,而 MLM 模型分析定位到的 SNP位点数目最少。在GLM 模型下分析, 显著水平为 0.1时定位到与AUDPC_Mean相关的SNP位点数目(74个)最多(表2)。

2.3 候选基因预测

本研究在82个关联SNP位点上下游区域100 kb,共检测到922个候选基因。其中, 460个候选基因在所有分析(60种)中都存在, 其余基因在2～56种分析中出现。根据基因组注释, 筛选出54个已知或可能与晚疫病抗性相关的基因(表3)。其中, 23个基因为抗性基因包括晚疫病抗性基因R1同源基因、Sw-5同源基因(R8)和Rpi-vnt1以及编码多效性耐药蛋白基因; 5个基因编码MAPK蛋白和WRKY转录因子;1个基因参与茉莉酸途径; 3个基因与水杨酸途径相关; 6个基因是病程相关的基因; 3个基因参与苯基丙酸类合成途径; 其他与晚疫病抗性相关的基因如HMGR基因(2个)、细胞色素P450 (21个)。

表2 6个马铃薯晚疫病抗性性状在不同分析模型和显著水平下的显著关联标记Table 2 Significant correlation markers of six phenotype of potato resistance to late blight under different analytical models and significant levels

(续表 2)

(续表 2)

3 讨论

3.1 不同模型分析对全基因关联分析结果的影响

本研究采用了5种模型对基因组数据和表型数据进行关联分析。GLM 模型(general linear model)采用Q矩阵[16], MLM模型(mixed linear model)采用Q+K矩阵[16], CMLM模型(compressed mixed linear model)算法介于GLM和MLM之间, 寻找最优算法压缩计算时间[20], EMMAX模型(efficient mixed model association eXpedited)在MLM模型基础上校正关联群体的群体结构和遗传亲缘关系并加速这种运算过程[21-22], FASTLMM (factored spectrally transformed linear mixed model)同样修正了MLM模型的运算公式, 提高运算速度[23]。尽管不同的模型存在差异, 但是在所有模型下分析的表型性状得到的结果都比较可靠(图2)。5种模型下共关联到82个SNP位点, 而同一性状、同样的显著水平下, GLM模型分析关联到的SNP位点数量明显多于其他模型(表2)。与晚疫病抗性相关的候选基因中, 关联到48个 SNP位点, 其中 13个(第 6、14、16、31、32、56、60、65、68、69、70、75和81号) SNP位点只在GLM模型下关联到, 16个SNP位点在所有模型到关联到, 剩余的19个SNP位点在2～4个模型下获得。因此, 利用多种模型分析能关联到更多可能与晚疫病抗性相关的候选基因。

3.2 马铃薯晚疫病水平抗性群体中可能存在的抗性基因

栽培马铃薯晚疫病抗性主要来自安第斯栽培种S. tuberosumAndigenum group (包括andigena、phureja和stenotomun)[24-26], 并通过导入野生种如S.demissum[27]、S. venturii[7]、S. bulbocastanum[28]等的抗性基因不断改良马铃薯栽培种的晚疫病抗性。尽管国际马铃薯中心通过不断的轮回选择创制了一批晚疫病水平抗性群体[29], 但是这批材料中可能仍存在垂直抗性基因(R1-R11)[30]。本研究所用的试验材料即是来自国际马铃薯中心筛选的水平抗性群体,并且关联到12个R1类似基因, 3个番茄斑萎病抗性基因, 4个Rpi-vnt1抗性基因以及4个多效性耐药蛋白基因。

植物中抗性基因编码的最大的一类蛋白具有核苷酸结合位点和亮氨酸重复结构(nucleotide-binding site plus leucine-rich repeat, NB-LRR)[31], 而目前已知的马铃薯晚疫病抗性基因编码的都是NB-LRR蛋白[32]。R1基因定位在5号染色体上, 是植物抗性基因亮氨酸zipper/NBS/LRR中的一员[33]。根据前人的研究, 至少有4个抗性基因位于9号染色体, 分别是R8[8-9]、R9a[34]、Rpi-vnt1[7]和Rpi-moc1[35], 这 4 个基因均位于 9号染色体长臂末端, 而且遗传距离比较近[34,36]。本研究所关联到的12个R1类似基因位于9号染色体, 其长度为 35.34 Mb, 而关联到的 SNP的位置在 34.8 Mb左右, 因此, 这些抗性基因肯定不是R1, 而可能是已克隆的R8或其他抗性基因。Sw-5是番茄中一类广谱抗性的番茄斑萎病(ToMV)抗性基因[37], 位于 9号染色体上的抗性基因R8与Sw-5同源, 而Rpi-vnt1与番茄中另1个ToMV持久抗性基因 Tm-22同源[7,38]。本研究中另外关联到与RPi-vnt1基因类似的基因, 其位置在32.5 Mb附近,而关联到的ToMV抗性基因位于34.7 Mb左右, 根据R8与Rpi-vnt1的遗传位置[36],Rpi-vnt1在R8的上方, 因此关联到位于32.5 Mb的基因可能是Rpi-vnt1,而位于34.7 Mb的基因可能是R8。多效性耐药蛋白是一类 ABC转运蛋白, 可参与到拟南芥[39](PEN3)和小麦[40](Lr34)中非宿主抗性和对病原菌的持久抗性中。马铃薯中的多效性耐药蛋白也可参与到晚疫病菌的侵染过程中, 但是具体怎么调节这个过程尚不清楚[41]。本研究关联到的多效性耐药蛋白基因位于 31.7 Mb左右, 其对马铃薯晚疫病抗性贡献的大小, 有待进一步研究。

?

?

?

3.3 其他可能影响马铃薯晚疫病抗性的相关基因

本研究所用材料主要是表现出晚疫病水平抗性,即非R基因主导的持久抗性[29]。还存在其他途径影响晚疫病抗性[42-43], 包括参与信号转导和基因调控、茉莉酸途径、水杨酸途径、病程相关基因、苯基丙酸类合成途径、HMGR基因和细胞色素P450。马铃薯中MAPK激酶和WRKY转录因子可参与植物对晚疫病菌的防御过程如StMPK1[44]、StWRKY8[45]、StWRKY1[46]。本研究在 7号染色体上关联到 1个MAPK基因和WRKY基因, 在3号染色体上关联到3个WRKY基因, 它们对晚疫病水平抗性的贡献有待进一步研究。激素信号在植物病害防御中发挥着重要作用, 植物激素茉莉酸和水杨酸代谢途径中的因子不同程度地参与到马铃薯晚疫病水平抗性中[47-49]。本研究中关联到 1个与茉莉酸代谢途径可能相关和3个与水杨酸代谢途径可能相关的基因, 在其他马铃薯晚疫病抗性全基因组关联分析和 QTL定位研究中也定位到与茉莉酸和水杨酸代谢途径相关的基因[50-51]。有趣的是, 本研究中关联到 3个参与苯基丙酸类合成途径的查尔酮合成酶基因(CHS), 位于9号染色体。马铃薯中CHS基因至少有6个拷贝,研究较多的是位于5号染色体的CHS基因参与马铃薯块茎中花青素合成[52-54]。但是, 查尔酮合成酶能够参与植物病害防御反应[55-56], 表明位于 9号染色体的CHS基因可能是构成马铃薯晚疫病水平抗性因素之一。其他基因包括与病程相关的基因、HMGR基因和细胞色素 P450都有报道可以参与植物病害防御过程, 但是在马铃薯晚疫病水平抗性中所起的作用是不清楚的, 而本研究关联到的基因可以为马铃薯晚疫病水平抗性的形成机理提供参考[42,57-59]。

4 结论

基于5种模型的GWAS结果, 共检测到82个与晚疫病抗性显著相关的SNP位点, 其中48个位点关联到可能与晚疫病抗性相关的候选基因, 包括R1同源基因、Sw-5同源基因(R8)、Rpi-vnt1基因、与信号转导相关的MAPK基因、WRKY基因、参与植物激素茉莉酸和水杨酸代谢相关的基因、参与苯基丙酸类合成途径的查尔酮酶基因和其他病程相关的基因、HMGR基因以及细胞色素P450基因。