基于元分析的大豆胞囊线虫病3号生理小种抗性基因挖掘

韩英鹏，田利峥，姜海鹏，包冬芳，王 俊，赵 雪

（东北农业大学大豆研究所，大豆生物学教育部重点实验室，农业农村部东北大豆生物学与遗传育种重点实验室，哈尔滨150030）

大豆胞囊线虫病（Soybean cyst nematode，SCN,Heterodera glycines Ichinohe）是大豆生产重要病害，主要分布于中国、美国、巴西、日本等国家。大豆胞囊线虫病害分布广泛，可造成20～30%产量损失，病害严重时高达80%，甚至绝产，严重危害大豆生产[1]。

大豆胞囊线虫病生理小种分化明显，Riggs 等利用大豆4个鉴别寄主和1个对照感病寄主调查胞囊线虫反应，将其划分为16 个生理小种[2]，根据7个鉴别寄主区划分为HG 类型[3]。目前，我国大豆主产区发现1～7 号、9 号、11 号、13 号、14 号等11个SCN生理小种，其中3号生理小种在我国分布最广，主要分布在东北三省和内蒙古；4号生理小种浸染能力更强，主要分布在黄淮和西北地区；6号生理小种主要发现于黑龙江省内[4]。抗病品种（系）选育是防治SCN 最有效方法，挖掘SCN 抗性基因是分子辅助育种基础，而QTL 定位是挖掘抗性基因最有效方法之一。

Weismann 等首次将抗性基因Rhg4 定位在pBLT24 与pBLT56 区间[5]。Concibido 定位多个SCN抗性QTL，且整合1994～2004年发表文献报道的60多个SCN 抗性QTL，分布在除2、9、10 等以外17个染色体[6]。目前，国内外文献已报道大量不同SCN生理小种抗性QTL，但因遗传抗性材料和接种方法不同，且不同地理环境和遗传背景对QTL 稳定性影响较大，导致QTL 基因组区域与实际所在基因组区域不完全相同。

Goffinet 建立可整合不同遗传背景和不同材料一致性QTL 元分析方法，基于最大邻接法整合不同来源QTL，缩小置信区间，提高QTL 定位结果有效性和准确性[7]。通过meta分析在QTL区段发现相关性状候选基因在大豆和其他作物中被广泛利用。李长育等利用BioMercatorV2.1软件通过meta分析获得2个一致性QTL位点，得到与大豆根瘤相关的8 个候选基因[8]。在玉米中，通过meta 分析方法，在847 个有关玉米籽粒蛋白质含量候选基因中，发现13个影响该性状的关键基因[9]。

大豆SCN抗性研究已定位和鉴定160多个QTL（http://www.soybase.or），其中，一些QTL未经重复验证，不同试验条件下，鉴定结果不一致。目前，已克隆和鉴定2 个主要SCN 抗性位点rhg1 和Rhg4，并揭示新的大豆抗病机制。目前已有栽培品种对SCN 抗 性 主 要 来 自PI88788、PI437654 和Peking，SCN 在不断适应过程中，易克服大豆寄主单一抗性，因此需鉴定新的SCN抗性位点和基因。

本研究使用图谱整合软件BioMercatorV4.2 将不同来源SCN Race 3 抗性QTL 整合到遗传图谱GmComposite 2003，并通过元分析方法缩小QTL置信区间，进一步提高QTL 准确性并获得一致性QTL，选取区间内所有基因并作基因注释，进一步筛选与SCN Race 3 相关抗性基因，分析其抗病机理，最终为分子辅助育种尤其选育抗SCN 品种（系）提供种质资源和理论依据。

1 材料与方法

1.1 大豆抗胞囊线虫QTL信息收集与整理

收集已发表的SCN Race 3 抗性QTL 信息，包括QTL 位置、置信区间、LOD 值、遗传贡献率、作图群体类型、群体规模和临近标记等信息。其中，我国品种（系）有中品03-5373 和L-10，中品03-5373 对SCN Race 3 抗性来自Peking、PI437654和灰皮支黑豆等3 个抗源，L-10 对9 个胞囊线虫小种均有抗性；美国品种（系）有PI437655、PI438489B、PI567516C、 PI89772、 Bell、 Msoy8001 和PI437654，均具有多种胞囊线虫小种抗性，其中PI437655、 PI438489B、 PI567516C、 PI89772 和PI437654等均为野生型大豆。

1.2 SCN Race 3抗性QTL映射与meta-QTL分析

利用BioMercatorV4.2 软件完成各作图群体SCN Race 3 抗性QTL 映射及meta-QTL 分析。基于每个QTL 图谱名称（Map name）、QTL 所在染色体（Chromosome）、群 体 规 模（Size）、群 体 类 型（Type）、QTL 名称（QTL name）、QTL 位置（Position）、遗传贡献率（R2）、LOD 值和各原始图谱与公共图谱上相关标记，构建一致性遗传图谱（Consensus map），随后利用ConsMap 和QTLProj 功能将每个原始图谱中QTL 映射到GmComposite2003 上。利用Veyrieras元分析作meta-QTL分析。

meta-QTL 分析后得到5 个模型，每个模型均按照最大似然函数比方法，利用高斯定理推导出在染色体上最可能排列位置和置信区间。根据AIC（Akaike-type-criteria-values）最小值原则，当在每条染色体上选取meta-QTL 数目时，首先考虑AIC 最小值对应的meta-QTL 数目模型，参考一致性图谱上原始QTL 分布及数目，选取合适数目，获取一致性置信区间。

1.3 SCN Race3 抗性meta-QTL 区段候选基因发掘及分析

利用已获得的QTL 一致性置信区间，根据基因组预测注释信息，挖掘SCN Race 3 抗性相关基因。查找已克隆的与大豆胞囊线虫病抗性相关基因，下载其表达序列，分析其结构域特征。将本研究所获置信区间内基因作为候选基因，同时根据其包含的不同功能结构域，将抗性基因分类，将所得抗性基因与大豆EST 数据库比对、注释。

1.4 SCN Race 3抗性相关基因蛋白结构域分析

在Phytozome网站（https://phytozome.jgi.doe.gov/）利用Glycine max Williams82 基因组数据库，获得基因参考序列，通过网站Functional Annotation 工具分析蛋白结构域，并使用绘图软件IBS绘图。

2 结果与分析

2.1 SCN Race 3抗性QTL信息收集与整理

分析收集的SCN Race 3 抗性QTL 信息可见，共涉及7 个作图群体，获得33 个SCN Race 3 相关QTL 信息（见表1），分布于11 个连锁群，其中G（chr18）分布较多，为5 个；F（chr13）、I（chr20）较少，均为1个。

表1 大豆胞囊线虫病3号小种抗性QTL相关信息Table 1 Information related to QTL resistance of soybean cyst nematode race 3 published

2.2 SCN Race 3抗性QTL映射与meta-QTL分析

共有21 个报道成功的QTL 映射到大豆参考图谱，利用BioMercatorV4.2 软件中的QTL meta-analyses程序分析计算，每次运算后以AIC最小值为原则，根据QTL模型得到3个一致性QTL，分别位于M（chr7）、E（chr15）和G（chr18）（见表2、图1）等3个连锁群中，其中置信区间指数（CI）为95%。这些一致性QTL 可为分子标记辅助育种、QTL 精细定位、挖掘抗性候选提供重要位置信息。

2.3 SCN Race 3 抗性meta-QTL 区段候选基因发掘及分析

挖掘利用BioMercatorV4.2 获得的一致性置信区间内候选基因，根据基因组测序注释信息，共发现18 个SCN Race 3 抗性相关基因（见表3），17个蛋白激酶（PK）型基因和1 个亮氨酸重复序列类受体蛋白激酶（LRR）型基因。

表2 候选QTL位点Table 2 Candidate QTL loci

表3 抗病基因类型统计Table 3 Types statistics of resistance genes

2.4 SCN Race 3相关抗性基因蛋白结构域分析

通过Functional Annotation工具分析基因结构域，发现17个PK型基因分为丝/苏氨酸蛋白激酶、半胱氨酸蛋白激酶、组/精/赖氨酸蛋白激酶3类（见图2）。

丝/苏氨酸蛋白激酶包括Glyma.18G074200、Glyma.18G096500、Glyma.18G026700、Glyma.18G0 26800、Glyma.18G026900、Glyma.18G062500、Glyma.07G003200、 Glyma.18G060900、 Glyma.18G042200、Glyma.18G038500 和Glyma.07G021100， 其 中Glyma.18G074200和Glyma.18G096500具有钙依赖蛋白激酶活性；Glyma.07G003200 和Glyma.18G060900具有分裂原活化蛋白激酶活性，该结构域参与大豆PTI和ETI两种防御机制，从而影响大豆对胞囊线虫抗性。半胱氨酸蛋白激酶包括Glyma.18G0 46100、Glyma.18G046200、 Glyma.18G046300、 Glyma.8G0 46600 和Glyma.18G046400，该结构域通过诱导大豆PTI防御机制从而影响大豆对胞囊线虫抗性；组/精/赖氨酸蛋白激酶包括Glyma.18G054100；Glyma.15G155600 为LRR 型基因，具有亮氨酸重复序列类受体蛋白激酶活性，该结构域通过参与大豆ETI防御机制从而影响大豆对胞囊线虫的抗性。

3 讨 论

利用精细定位方法可从已定位QTL 中准确确定候选基因，但耗时耗力。利用元分析方法获得QTL一个较小置信区间，将目的基因锁定在较小置信区间内，再根据抗病基因保守性挖掘目标基因，最后利用EST 注释信息进一步确认所获基因，该方法省时省力且费用极低。常玮等利用BioMercator2.1 软件，通过已知SCN 抗性QTL 和参考图谱GmComposite2003获得一致性位点，根据大豆20个连锁标记（SSR、SNP 和RFLP）序列和大豆全基因组数据构建物理图谱，与已获得一致性QTL 位点，根据hmmsearch 搜索结果获得77 个抗性基因[20]。在77个候选基因中，多数基因已证明与SCN抗性有关。因此，通过元分析计算确定一致性QTL位点从而准确筛选候选基因，为大豆分子辅助育种尤其是选育抗SCN 品种（系）提供一种新的思路和途径。

Li等研究中品03-5373（抗性）与中黄13（感性）及其后代所构成群体，根据其基因型和表型，经全基因组关联分析（GWAS），发现基因Glyma.18G026900具有SCN 抗性连锁的单核苷酸多态性（SNP）位点（C/T）[21]。因此，可推断其同源基因Glyma.18G026700、Glyma.18G026800 和Glyma.18G062500 等也可能具有SCN抗性连锁的单核苷酸多态性（SNP）位点。但其余15个有关SCN抗性基因尚无报道。

Glyma.18G026900及其同源基因具有丝氨酸/苏氨酸蛋白激酶活性，属于PK 酶（蛋白激酶）的一种。植物蛋白激酶在逆境胁迫信号传导中发挥关键作用，维持植物体正常生命活动，根据底物特异性，植物蛋白激酶分为丝/苏氨酸蛋白激酶、组/精/赖氨酸蛋白激酶、酪氨酸蛋白激酶、半胱氨酸蛋白激酶、天冬酰胺基/谷氨酰胺基蛋白激酶等5类[22]。其中丝/苏氨酸蛋白激酶包括分裂原活化蛋白激酶（MAPK）和钙依赖蛋白激酶（CDPK）：当大豆受体蛋白参与植物体PTI 防御机制时，MAPK 在其下游的防御机制中发挥重要作用，植物体特异抗病蛋白，用于识别效应蛋白，即防御机制ETI，大豆通过MAPK在大豆胞囊线虫防御机制ETI中发挥重要作用[23]；CDPK 主要调节植物体中盐信号传导，其在大豆胞囊线虫中抗性机理还需进一步研究。半胱氨酸类受体激酶（Cysteine-rich receptorlike-kinases，CRKs）属于蛋白激酶，在植物体受病原菌胁迫时，CRK 诱导分泌甲基SA 诱导植物体PTI防御机制，使植物局部过敏性坏死，保护植物体[24]。其余蛋白激酶在植物非生物胁迫中发挥作用，但在生物胁迫中，尤其在SCN抗性中尚无报道。

LRR 型蛋白因其结构域中富含亮氨酸重复序列而得名。参与植物体防御反应的LRR 型蛋白可分为类受体蛋白激酶、抗病基因编码蛋白质、多聚半乳糖醛酸酶抑制蛋白和伸展蛋白等4类[25]。类受体蛋白激酶主要参与植物体PTI防御过程；植物抗病基因编码蛋白质主要参与植物体ETI 过程[26]；多聚半乳糖醛酸酶抑制蛋白（Polygalac-turonase-inhibiting proteins，PGIPs）通过保护植物细胞壁抵御病原菌入侵；伸展蛋白含有组成植物体细胞壁的结构蛋白，是防止病原物侵入和扩散的屏障。Kanazin 等在16 号染色体上定位到一个具有LRRTM 结构域的基因，已证实其与大豆胞囊线虫抗性相关[27]。

本研究在整合前人QTL 信息基础上，采用图谱整合软件BioMercatorV4.2 挖掘SCN Race 3 抗性基因，将收集和整理的33 个QTL 整合在遗传图谱GmComposite2003上，并通过元分析方法得到一致性QTL，获得3 个抗性候选区段，发现18 个与SCN Race 3 抗性相关基因。本研究选取的SCN Race 3抗性QTL信息及所构建遗传图谱描述较为详细，因此可使用meta 分析方法发现与SCN Race 3抗性相关基因，对发现其他SCN 生理小种抗性基因具有借鉴意义。该方法为分子辅助育种提供新的思路和途径，也为后期筛选候选基因，分析其抗病机理，最终为选育抗SCN 品种（系）提供种质资源和理论基础。

4 结 论

采集美国、中国和巴西等发表的33 个与SCN Race 3抗性相关QTL及其相关信息，利用QTL-meta分析得到3 个一致性QTL，从中发现18 个与SCN Race 3 抗性相关基因，其中17 个PK 型基因，1 个LRR 型基因。17 个PK 型基因位于连锁群M（chr7）和连锁群G（chr18），1个LRR型基因位于连锁群E（chr15）。通过蛋白结构域分析，在17 个PK 基因中，1个基因具有组/精/赖氨酸蛋白激酶活性、5个基因具有半胱氨酸蛋白激酶活性、11 个基因具有丝/苏氨酸蛋白激酶活性。