木麻黄青枯病抗性与EST-SSR标记的关联分析*

许秀玉 张 勇 甘先华 张华新 张应中 王明怀 毛亦杨

（1 广东省森林培育与保护利用重点实验室/广东省林业科学研究院，广东 广州 510520；2 中国林业科学研究院 国家林业局盐碱地研究中心， 北京 100091； 3 中国林业科学研究院 热带林业研究所， 广东 广州 510520； 4 广东省森林资源保育中心， 广东 广州 510173）

植物青枯病（Ralstonia solanacearum）抗性性状受多基因控制，属于数量性状，采用常规育种方法选育抗病品系难度大，进展慢[1]。木麻黄科（Casuarinaceae）植物自然分布区广，不同树种、种源及无性系间遗传变异非常丰富，在青枯病抗性上也存在着极显著分化[2]。近年来，随着越来越多木本植物[3-6]完成全基因组测序及大量分子标记的开发，数量性状的研究，特别是利用关联分析挖掘数量性状相关功能基因已经越来越成为国际植物基因组学研究的热点。林木繁衍以自然杂交为主，群体遗传多样性水平高，连锁不平衡水平低，特别适用于关联分析[7]。目前，许多林木树种开展了数量性状与分子标记的关联分析，但大多集中在生长性状与木材性状[8]，利用关联分析挖掘林木抗病基因位点的研究极少见报道，仅有火炬松（Pinus taeda）-抗镰刀球菌（Fusarium circinatum）[9]、抗脂溃疡病（Fusarium circinatum）[10]，石榴（Punica granatum）-抗白叶枯病（Xanthomonas axonopodis）[11]，巨尾桉（Eucalyptus grandis × E. urophylla）-抗锈病（Puccinia psidii）[12]等少数几个树种病害关联分析。木麻黄不仅尚未开发足量的分子标记，在分子辅助选育、遗传图谱构建、抗病基因定位等分子生物学领域的研究严重滞后，木麻黄青枯病抗性相关分子标记研究国内外也极少见报道。本研究利用木麻黄223 个EST-SSR 标记，对我国现有木麻黄种质资源青枯病抗性进行关联分析，旨在挖掘与木麻黄青枯病抗性密切相关的标记位点或具有特定功能的基因位点，为分子标记辅助育种提供有益参考，加速抗病种质资源的选育进程。

1 材料与方法

1.1 试验材料

供试53 份木麻黄种质来源于广东、福建、海南，是广东、福建、海南各省根据不同育种目的选育并保存下来的无性系，由广东省林业科学研究院培育。

供试材料青枯病抗性鉴定方法、病情指数、抗性分级等数据来源于许秀玉等[13]研究结果，为了统计分析方便，本文将病情指数换算成抗病指数，抗病指数=100-病情指数。供试木麻黄种质材料信息见表1。由表1 可见，供试材料抗病指数变异范围大，介于27.2～100.0 之间，基本呈连续正态分布，适合用于抗病关联分析。

1.2 DNA 提取与标记分型

采取木麻黄新鲜嫩枝，采用改良的 CTAB 法[14]提取基因组DNA。检测浓度与纯度后稀释至100 ng/μL，作为PCR 扩增模板。

选取杂交、30、41、45、W6、G1、平5 和X1等8 个亲缘关系较远、抗病指数大于95 的无性系构建抗病基因池；选取2、34、37、59、82、G88、K18 和C7 等8 个亲缘关系较远、抗病指数小于50的无性系构建感病基因池。利用223 个木麻黄ESTSSR 标记[15]对抗、感池进行标记分型与筛选。初步筛选标记方法[16]：比较两个基因池的标记-等位基因频率，相差80%以上的标记入选，筛选出用于关联分析的标记。

反应体系、分型过程及参试223 个标记具体信息参见Xu 等[15]研究结果。

1.3 数据分析

1.3.1 群体结构及关联分析 标记分型与筛选采用GenAlEx 6.4.1 检测[17]，筛选出在抗感基因池中能差异表达且不偏离哈-温平衡的标记用于关联分析。利用Structure 2.3.4[18]软件分析木麻黄供试材料群体结构，设定亚群数（K）的取值范围1～16，将MCMC（Markov Chain Monte Carlo）的Length of Burn-in Period（不作数迭代）设为 100 000 次，每个K 值重复运行10 次，其余参数采用默认设置。根据 lnP(D)计算ΔK，以ΔK 最大原则选取最适合的亚群数（K），并计算各参试材料的基因型属于第K 亚群的概率（Q 值）。利用SPAGe-Di1-5a 软件[19]获取供试材料个体间亲缘关系系数K 矩阵。基于Q+K 模型，选择TASSEL 3.0 软件[19]的MLM 模型，确定P＜0.05 时的关联位点，计算标记对表型变异的解释率（R2）。

表1 供试无性系及抗病指数Tab.1 The clones and the disease index

1.3.2 等位变异表型效应分析 采用无效等位基因（null allele）法[20]计算各增效/减效等位变异表型效应值（Ai），确定关联标记的增效/减效等位变异，将某一标记增效/减效等位变异的平均效应作为该标记的增效/减效应值（AAE），计算表型效应比（AN）。Ai=∑xij/ ni-∑Nk/ nk；AAE= ∑ac/ nc；AN=AAE / ( ∑Nk/ nk) × 100%。（xij：第j 个具有i 等位片段个体的表型值；ni：具有i 等位片段的个体数；Nk：第k 个具有无效等位基因个体的表型值；nk：具有无效等位基因的个体数；ac：关联位点内第c 个增效/减效等位变异表型效应值；nc：关联位点内增效/减效等位变 异数。）

2 结果与分析

2.1 木麻黄EST-SSR 标记分型与筛选

检测发现，223 个木麻黄EST-SSR 标记中，共有49 个标记在抗、感基因池中等位基因频率差异达80%以上，且CASeSSR244 与CASeSSR064两个标记偏离哈-温平衡（表2），因此选择其余47 个标记作为关联分析标记，利用这47 个标记对参试材料进行群体结构分析与亲缘关系系数计算。

2.2 群体结构分析

根据模型后验概率lnP(D)计算ΔK，通过ΔK 来确定亚群数（K），以ΔK 值最大时的K 值即为最佳亚群数。由图1A 可知，当K=2，ΔK 值最大，即参试群体明显存在2 个亚群，表明将53份木麻黄种质资源分为2 个组群时，具有最大的后验概率ln P(D)，组群间具有最大的遗传差异，组群内植株具有相对一致的遗传背景。Q 值计算结果发现：所有参试材料在各自所属组群中的Q值全部大于0.5，这表明所有参试材料都可归入2个组群中某一个，群体结构简单，遗传组分相对单一，有利于减少关联分析的假阳性，提高关联分析效果。

表2 49 个标记哈-温平衡卡方检验结果Tab.2 Results of 49 EST-SSRs Chi-square tests for Hardy-Weinberg equilibrium

构建的群体结构图如图1B，亚群1 包含34个无性系，亚群2 包含19 个无性系，来源不同的参试材料在不同的亚群中均有分布，这表明亚群的划分与地理来源之间没有必然联系，这与前期木麻黄种质资源亲缘关系分析结果相一致[21]。

2.3 关联分析

由表3 关联分析结果可见：有10 个标记与青枯病抗性显著相关（P ＜ 0.05），相关标记对表型变异的解释率较高，介于12.6%～60.1%之间。其中有5 个标记具有增效效应，标记CASeSSR359、CASeSSR052 与CASeSSR056 增效效应值最大，分别为19.1、18.1、11.7，其相对的表型变异解释率分别为19.1%、18.1%、43.7%，推断这些标记与某个贡献率较大的抗病位点连锁的可能性较大，可利用这些标记对木麻黄种质资源的青枯病抗性潜力做出评价。除了CASeSSR359、CASeSSR052标记外的其余8 个标记均具有减效效应，CASeSSR022、CASeSSR242、CASeSSR176 减效效应值最大，分别为32.2，22.0，21.4，相应的表型变异解释率分别为35.9%，20.7%，12.6%，初步推断这些标记与感病位点连锁的可能性较大，可利用这些标记淘汰易感植物材料。

2.4 位点和等位片段的表型效应

由表4 可见，具有增效表型效应等位变异的 关 联 位 点 有CASeSSR052、CASeSSR056、CASeSSR241、CASeSSR263 及CASeSSR359，其表型效应范围分别为：2.9～37.1、6.1～17.3、2.0～5.3、2.1～13.8、 11.5～27.3，其优异增效等位变异分别为：CASeSSR052-255 bp（PE=37.1）、CASeSSR052-237 bp（PE=37.1）、CASeSSR359-370 bp（PE=24.0）、CASeSSR359-386bp（PE=27.3）。具有减效表型效应的关联位点有CASeSSR022、CASeSSR056、CASeSSR176、CASeSSR241、CASeSSR242、CASeSSR253、CASeSSR258、CASeSSR263，其中CASeSSR022-196 bp、 CASeSSR022-190 bp、CASeSSR022-188 bp、CASeSSR253-202 bp、CASeSSR242-227 bp、CASeSSR253-186 bp、CASeSSR176-204 bp 的等位变异减效表型效应均高于20.0。育种中，可优先利用表型效应值大的优异等位变异，利用优异增效等位变异进行优质亲本及抗病植物材料的筛选，为后续标记辅助育种及杂交育种提供有益参考。

图1 木麻黄53 个无性系群体结构Fig.1 Population structure of Casuarinaceae

表3 关联标记的平均增效效应、减效效应Tab.3 EST-SSR associated with R. solanacearum resistance could be detected with MLM model in Casuarinaceae and average positive (negative) allele effect

表4 木麻黄无性系青枯病抗病关联位点及其等位变异Tab.4 Associated loci and phenotypic effect for R. solanacearum resistance

3 结论与讨论

木麻黄分子生物学研究基础薄弱，早期利用农杆菌开展了木麻黄基因工程研究[22]，利用FISSR、ISSR、SSR、SCAR 等标记开展了木麻黄遗传多样性[23]、群体结构[23-24]、群体亲缘关系[21]、木麻黄物种鉴定[25]等少量分子育种研究，2018 年开 发 了223 个 木 麻 黄EST-SSR 标 记[15]，2019 年首次公布了短枝木麻黄（Casuarina equisetifolia）的基因组数据[3]，本研究在此基础上首次开展了木麻黄青枯病抗病性状与分子标记的关联分析，共发掘出CASeSSR241、CASeSSR263、CASeSSR056 等10 个与木麻黄青枯病抗病性状显著关联的EST-SSR 位点及CASeSSR052-255 bp、CASeSSR052-237 bp、CASeSSR359-370 bp、CASeSSR359-386 bp 等优异增效等位变异，为今后预测控制性状的主效QTL、新基因发掘、标记辅助育种、杂交育种等研究提供参考。其次，CASeSSR253-202 bp、CASeSSR022-190 bp、CASeSSR022-196 bp 等多个等位片段对木麻黄青枯病抗性具有较高的减效效应，利用这些等位变异可以淘汰木麻黄感病材料。

本研究从223 个木麻黄EST-SSR 标记中筛选得到10 个与青枯病抗性相关的标记，关联标记数量偏少，这可能由于与木麻黄青枯病抗性相关的基因原本就不多，也可能由于所用的EST-SSR 标记数量不够多，今后可扩大标记数量与种类以获得更多的关联位点。本研究筛选出的10 个关联标记表型变异解释率介于12.6%～60.1%，高于大多数植物抗病虫害关联标记，如辣椒（Capsicum annuum）疫病（Phytophthora capsici Leonian）抗性关联分析中表型变异解释率2.98%～16.05%[26]，烟草（Nicotiana tabacum）青枯病抗性关联分析中表型变异解释率5%～14%[19]，陆地棉（Gossypium hirsutum）黄萎病（Verticillium dahliae）抗性关联分析表型变异解释率3%～13%[27]等，这可能由于本次试验群体数量少造成的，也可能由于这些标记位于基因编码区或基因表达调控区，是青枯病抗性性状重要位点。

木麻黄为外来树种，本研究参试群体为我国海南、福建、广东各省根据抗风、速生、抗旱等不同育种目的选育出来的，为我国现存的木麻黄重要种质资源，因此群体数量少。利用这个群体进行关联分析难免存在着LD 偏差及表型取样误差[26]，为了提高关联分析的准确性，今后要进一步加强木麻黄育种群体的收集、积累、选育及杂交品系的创制，不断扩大用于木麻黄关联分析的群体数量，对获得的标记位点与优异等位变异进一步验证。其次，本研究中木麻黄青枯病抗性表型数据为室内接种试验所得，而林木对青枯病抗性的强弱除了受基因型影响外，还受自身生长势、侵染方式及环境条件等因素的影响，今后还应开展多地点田间试验，获得不同地点、不同生境条件下，不同木麻黄资源对青枯菌的抗病性能，研究环境与基因型的互作，以得到更多更准确的抗性相关位点与优异等位变异。