利用BSA-Seq筛选甘蓝型油菜菌核病抗性候选基因

夏吉春， 彭于洲， 王东， 李桐舟， 张凯， 徐新福

西南大学 农学与生物科技学院，重庆 400715

由核盘菌引发的油菜菌核病是油菜的第一大病害, 近年来随着机械化和双低油菜的推广, 该病害呈加重趋势． 油菜菌核病的致病机理较为复杂, 尚未有此方面系统全面的报道[1-3]．

现有研究表明, 多种水解酶及草酸等毒素在致病中起到了重要作用． 核盘菌通过分泌释放纤维素酶、 果胶酶、 水解酶等裂解植物细胞壁[4], 从而获得孢子萌发及菌丝生长、 扩散所需的营养成分[5]． 早在19世纪就有学者发现草酸在核盘菌侵染中起重要作用[6]． 在后续研究中发现草酸会酸化、 腐败寄主植物组织, 主要作用机理: ① 降低pH值, 促进核盘菌分泌水解酶, 同时抑制寄主中起抗病作用的酶[7]． ② 结合细胞中的二价金属离子, 破坏细胞结构与功能, 如损坏细胞壁与核糖体、 阻塞导管、 阻碍叶绿素合成、 改变膜的通透性、 影响细胞信号转导等, 致使细胞正常生理功能受到影响[8-9]． ③ 刺激保卫细胞对K+的吸收, 改变膜两侧的渗透压, 使保卫细胞无法正常工作, 致使气孔打开, 叶片脱水, 为核盘菌侵入创造有利条件[10]． ④调控寄主活性氧(ROS)代谢． 侵入前期, 抑制活性氧生成, 阻止植物自身的防御进程; 侵染中后期, 诱导氧爆反应, 诱发细胞程序性死亡(PCD), 死亡细胞可作为营养源供病原菌生长、 繁殖[11]．

油菜通过多种途径响应来应对核盘菌的侵染． Wang等[12-13]发现油菜依次激活水杨酸(SA)途径、 茉莉酸/乙烯(JA/ET)途径来产生对核盘菌的抗性, 通过表达丝裂原活化蛋白激酶 BnMPK4, 提高了油菜植株对菌核菌的抗性, 说明 MAPK 途径参与油菜抗菌核病反应． 在次生代谢物方面, 宋志荣等[14]发现高硫苷的油菜菌核病抗性更高． 另外, 植保素、 木质素、 酚类化合物也参与了油菜的抗病响应[15]．

当前, 对油菜菌核病抗病基因的发掘主要从中心法则的3个方面开展研究． 在DNA水平上, 通常采用数量性状位点(QTL)定位进行研究． Zhao等[16]利用107个分子标记在甘蓝型油菜128-F2: 3家系群体中检测到3个与苗期叶片抗性相关的QTL和3个与成株期茎秆抗性相关的QTL, 马田田[17]在F2群体中使用复合区间作图法检测到3个叶片抗性的QTL和2个茎秆抗性的QTL． 在这些研究中, 苗期与成株期抗病QTL未检测出重合区域, 说明油菜菌核病的抗性具有特异性, 在不同时期或不同器官中可能存在不同的抗性位点． 值得一提的是, 梅家琴[18]在甘蓝抗感分离的F2无性系群体中检测到了一个叶片抗性与茎秆抗性重叠的QTL, 这与之前的研究出现了差异． Wu等[19]2013年鉴定出一个与茎秆相关的重要QTL位点, 具有较高的遗传效应, 并选定BnaC.IGMT5.a为候选抗病基因．

在RNA水平上, 多数研究者利用cDNA芯片技术研究核盘菌侵染后诱导表达的基因, 或不同抗性材料在接种核盘菌后的基因表达差异． 除去指导合成各种抗性相关的蛋白质外, 相关研究主要集中于转录因子、 信号转导及次生代谢物合成基因等[20]． 张卡[21]发现BnaA03.WRKY28和BnWRKY33两个转录因子均能结合BnWRKY33启动子, 前者为抗性负调控因子, 后者的大量表达可增强抗病反应． 许李明[22]的研究表明,BnEIN3是正向调控抗性的乙烯信号转导关键调节因子．

在蛋白质水平上, 早期学者们通过检测接种核盘菌前后植株体内防御酶的变化规律来研究抗性相关性, 王汉中等[23]发现几丁质酶、 苯丙氨酸解氨酶(PAL)、 多酚氧化物酶(PPO)、 过氧化物酶(POD)等与菌核病抗性有关, 之后则采用蛋白组学的方法研究抗病基因． Wen等[24]采用蛋白组双向电泳技术, 根据蛋白质表达谱差异得到了甘蓝型油菜20个抗性相关蛋白． 此外, 对油菜菌核病病程蛋白(PR)的研究也具有广阔的发展前景, 石美娟等[25]发现甘蓝型油菜类甜蛋白基因BnTLP1正调控菌核病抗性, 且JA/ET途径可能早期参与其抗性响应．

本研究是在极端抗性材料构建的F2分离群体基础上, 利用高通量测序方法筛选鉴定可能的甘蓝型油菜菌核病抗性基因, 为后续基因克隆、 功能分析及油菜菌核病在分子机理层面的剖析奠定基础．

1 材料与方法

1.1 试验材料

本试验所有材料均来源于重庆市油菜工程技术研究中心． 高抗父本21Y490发病率在30%左右, 相对于ZY821为中抗品种, 病情严重度为1级(病斑长1～2 cm); 高感母本21Y689发病率高达100%, 相对于ZY821为高感品种, 病情严重度在3级以上(病斑长3 cm以上)． 将高抗父本21Y490与高感母本21Y689杂交, 杂种F1代自交后获得F2分离群体, 于2021年在重庆市武隆区仙女山种植．

将核盘菌在PAD培养基(20%马铃薯、 2%葡萄糖及1.5% 琼脂)平板上复苏后, 取平板边缘菌丝转接, 对第二轮转接后的菌丝进行暗培养3 d, 待菌丝刚布满平板时用灭菌后的打孔器(内径0.6 cm)沿平板边缘打孔, 获得的带菌丝琼脂块作为接种体, 用于接种[26]．

1.2 试验方法

1.2.1 茎秆离体接种鉴定

于终花期离地面10 cm处截取约30 cm长的茎段, 去除多余枝叶并将两端敷上浸过水的棉花球, 用保鲜膜封严保湿． 通过无菌打孔器(内径 0.6 cm)在茎秆距两端约 10 cm 处破坏表皮, 产生与接种体大小一致的伤口, 将接种体带菌面紧贴创伤处, 覆盖保鲜膜, 接种后置于培养箱内72 h测量菌斑长度[26]．

1.2.2 构建混池进行重测序

根据侵染菌斑长度从F2分离群体中挑选30株高抗植株和30株高感植株使用十六烷基三甲基溴化铵(CTAB)法分别提取DNA, 构建DNA混合池, 两混池、 两亲本池均送北京百迈克生物科技有限公司利用Illumina HiSeq平台建库并进行全基因组重测序．

1.2.3 混合分组分析法(BSA)分析

由百迈克生物科技公司完成BSA分析, 主要步骤如下:

(1) 测序数据质控

根据碱基质量、 类型等数据确认测序质量, 并删除原始测序序列中低质量Reads(带接头、 N比例大于10%、 质量值Q≤10), 过滤得到Clean Reads．

(2) 参考基因组比对

利用bwa软件将Clean Reads与甘蓝型油菜参考基因组(ZS11油菜参考基因组)比对[27]． 通过比较确定了Clean Reads在参考基因组上的位置, 并统计测序深度和基因组覆盖度．

(3) SNP和InDel的检测与注释:

SNP和InDel检测使用GATK软件工具包, 根据前面比对定位的结果, 预处理后使用GATK进行SNP和Small InDel检测, 过滤后得到最终结果[28]． 具体流程参考GATK官网: https: //www. broadinstitute.org/gatk/guide/best-practices.php． 根据比对结果, 将4个混池间所有差异的变异位点汇总, 统计出样品间的SNP和InDel位点．

SNP和InDel注释使用SnpEff软件, 根据变异位点在参考基因组上的位置及信息, 可以得到变异位点在基因组发生的区域, 以及变异产生的影响[29]．

(4) 关联分析

对过滤后的SNP和InDel, 分别用欧式距离(Euclidean Distance, ED)算法和SNP-index方法分析关联区域． 根据公式计算各位点的ED值后, 取ED值的5次方为关联值以消除背景噪音, 采用DISTANCE方法对ED值进行拟合[30]． SNP-index算法同样采用DISTANCE方法对ΔSNP-index进行拟合, 然后根据关联阈值, 选择阈值以上的区域作为与性状相关的区域．

(1)

式(1)中:Amut为A碱基在突变混池中的频率,Awt为A碱基在野生型混池中的频率;Cmut为C碱基在突变混池中的频率,Cwt为C碱基在野生型混池中的频率;Gmut为G碱基在突变混池中的频率,Gwt为G碱基在野生型混池中的频率;Tmut为T碱基在突变混池中的频率,Twt为T碱基在野生型混池中的频率．

SNP-index(aa)=Maa/(Maa+Paa)

(2)

式(2)中:Maa表示aa池来源于母本的深度;Paa表示aa池来源于父本的深度;

SNP-index(ab)=Mab/(Mab+Pab)

(3)

式(3)中:Mab表示ab池来源于母本的深度;Pab表示ab池来源于父本的深度;

ΔSNP-index=SNP-index(aa)-SNP-index(ab)

(4)

(5) 候选区域筛选与功能注释

取两种算法的交集得到候选区域, 参考NR[31],GO[32],KEGG[33],COG[34]等数据库, 通过BLAST[35]软件对候选区间内的编码基因进行注释．

1.2.4 菌核病抗性基因预测

本研究主要通过基因注释文件, 根据蛋白质结构域和同源基因等信息, 结合油菜抗性机制进行抗性基因预测．

2 结果与分析

2.1 甘蓝型油菜抗菌核病表型鉴定

于2021年9月1日在重庆市武隆区仙女山收取F2分离群体材料, 接种核盘菌培养后每个材料量得两个菌斑长度, 取平均值． 由统计结果可得, F2分离群体菌核病抗性表型符合正态分布规律, 也符合微效多基因控制的数量性状遗传规律(图1)． 其中, 最大的菌斑长9.5 cm, 最小为1 cm． 选取菌斑平均值最大和最小各30株材料构建混池． 图2为极端材料表型对照图片．

图1 菌核病表型数据统计

图2 极端材料对比

2.2 BSA分析

2.2.1 测序数据质控

测序得到的原始测序序列(Sequenced Reads), 里面含有带接头的、 低质量的Reads, 为了保证信息分析质量, 对数据进行过滤, 得到过滤后数据, 用于后续信息分析(表1)．

表1 样品数据评估

2.2.2 与参考基因组比对统计

参考基因组为Brassica_napus: Bra_napus_v2.0(基因组大小为979 Mb, GC%含量为36.47%)． 4个样品同参考基因组的比对效率在98.82%以上, 将Reads定位至参考基因组, 统计碱基覆盖情况, 亲本池覆盖深度为13 X,16 X, 抗病池、 感病池覆盖深度为31 X,34 X, 平均覆盖深度约为23.50 X, 基因组至少覆盖1 X的比率约为95.12%． 以上结果表明, 测序正常且随机性好(表2)．

表2 比对结果统计

2.2.3 SNP和InDel的检测与注释

亲本间获得基因间区SNP 531 100个, 基因内SNP 627个, 非同义突变SNP 101 340个, 总计获得1 879 742个SNP． 混池间获得基因间区SNP 160 720个, 基因内SNP 149个, 非同义突变SNP 16 850个, 总计获得440 478个SNP(表3)．

表3 SNP注释统计

亲本间获得移码突变Small InDel 6 768个, 基因间区Small InDel 95 392个, 基因内Small InDel 601个, 总计获得502 652个Small InDel． 混池间获得移码突变Small InDel 1 958个, 基因间区Small InDel 32 569个, 基因内Small InDel 127个, 总计129 915个Small InDel(表4)．

表4 InDel注释统计

2.2.4 关联分析

SNP的ED法分析, 关联阈值为0.06, 共得到18个区域, 位于A04,A06,C02,C03,C05,C07,C09共7条染色体上, 总长度为 28.64 Mb, 其中最长区段大小为9.95 Mb, 包含615个基因, 最短区段大小为0.07 Mb, 包含3个基因, 具有非同义突变位点的基因550个． InDel的ED法分析, 关联阈值为0.09, 共得到13个区域, 位于C02,C03,C05,C06,C09共5条染色体上, 总长度为24.46 Mb, 其中最长区段大小为9.05 Mb, 包含616个基因, 最短区段大小为0.1 Mb, 包含1个基因, 具有移码突变位点的基因87个(图3)．

NC_027757.2-NC_027775.2依次为甘蓝型油菜A01-A10、 C01-C09染色体．图3 ED关联值在染色体上的分布

对于SNP(InDel)-index的计算结果, 当置信度为0.99时, 并没有关联到区域． 为了最大限度地利用这些数据, 本文将阈值调整为0.24, 共得到14个区域, 全部位于C09染色体上, 总长度为8.40 Mb, 最长区段大小为2.33 Mb, 包含123个基因, 最短区段大小为0.03 Mb, 包含1个基因, 具有非同义突变位点的基因124个． InDel共得到8个区域, 也同样全部位于C09染色体上, 总长度为 8.58 Mb, 最长区段大小为2.30 Mb, 包含122个基因, 最短区段大小为0.02 Mb, 包含1个基因, 具有移码突变位点的基因25个(图4)．

NC_027757.2-NC_027775.2依次为甘蓝型油菜A01-A10,C01-C09染色体．图4 SNP(InDel)-index关联值在染色体上的分布

2.2.5 候选区域筛选与功能注释

取以上4个部分的交集得到8个候选关联区域, 全部位于甘蓝型油菜C09染色体上, 主要分布在该染色体上17.8～21.2 Mb及31.6～39.1 Mb两个区间内(表5)．

表5 候选关联区域统计

参考NR,GO,KEGG,COG等数据库, 注释结果表明候选基因主要集中在代谢及遗传信息进程, 而富集分析表明在二萜生物合成、 硫代谢、 萜类骨架生物合成等通路中富集程度较高(图5)．

图5 候选区域内基因的通路分布图

2.3 抗病基因预测

在获得的C09染色体上8个关联区域中重点筛选111个非同义突变或移码突变基因, 根据基因注释文件, 结合油菜菌核病抗性机制, 预测8个基因为油菜菌核病抗性候选基因(表6), 分别是:BnaC05g50540D,BnaC09g33900D,BnaA04g09130D,BnaA01g29170D,BnaC09g15640D,BnaC09g35310D,BnaA03g09220D,BnaCnng15460D．

表6 候选基因功能注释、 同源分析及参与的调控途径

3 讨论

随着现代分子生物学技术的快速发展, 以二代测序技术为背景衍生出的BSA-seq技术凭借其高效的优势成为目前基因定位的热门方法． 本研究利用BSA-seq技术定位油菜菌核病茎秆抗性基因, 发现关联区域全部位于C09染色体上, 而此前定位的油菜菌核病QTL大多出现在A基因组上, 且与Mei等[36]的研究相互印证, 证实了油菜C基因组对菌核病抗性研究的重要价值．

本研究在候选区间共筛选出8个候选基因, 其中BnaC05g50540D基因包含MATH结构, 有研究表明该结构域与植物同真菌的共生效应、 抵抗病毒有关, 对信号转导也会产生影响[37], 其在拟南芥中的同源基因USP12也参与了茉莉酸反应及脱落酸信号转导． 已有研究表明, 茉莉酸(jasmonic acid, JA)、 水杨酸(salicylic acid, SA)及活性氧(Reactive Oxygen Species, ROS)是植物体内抵杀病原菌的重要信号途径, 其中JA为植物创伤反应调节子, 能根据植物与病原菌相互作用的特异类型促使症状朝正或负的方向发展, 以控制病害程度[38-39]．BnaC09g15640D的同源基因GH3.12也与茉莉酸调节有关, 且其GH3基因家族被证实在水稻抗白叶枯病中起正向调节作用[40]．BnaC09g33900D是拟南芥MAPKKK17的同源基因,MAPK途径参与了油菜菌核病的抗性响应, 超表达BnMAPK4能促进茉莉酸抗病信号转导途径中重要基因PDF1.2的表达[13]．BnaA04g09130D和BnaA01g29170D都含有植物NLR免疫受体具有的ZF-BED结构域, 该结构域在水稻白叶枯病和小麦条锈病对抗性方面起重要作用, 而在油菜中有关ZF-BED结构域还未有报道． 预测其可能在油菜对抗核盘菌入侵过程中发挥作用, 并为油菜抗病育种提供新的方向[41]．BnaC09g35310D的同源基因GA1[42]及BnaA03g09220D的同源基因GASA10都与赤霉素调节有关, 说明赤霉素可能在油菜菌核病抗性方面发挥重要作用, 值得更深入地研究．BnaCnng15460D含有FBD结构域, 属于F-box基因家族, 而该家族C端结构域能和许多蛋白特异性结合, 导致其功能多样性． 左蓉等[43]已对油菜F-box-LLR基因家族进行了研究, 表明该基因家族与植物抗病响应有紧密的联系, 且极大可能具有菌核病抗性功能． 值得注意的是,BnaA04g09130D同样属于F-BOX蛋白, 这可能暗示ZF-BED结构域和FBD/LLR结构域存在相互作用关系, 且可能在植物抗病系统中起重要作用．

油菜对于菌核病的抗性是一个复杂的过程, 由多基因协同调控． 在对候选基因的挖掘中发现大多数基因都参与了茉莉酸信号转导, 表明油菜在面对核盘菌侵染时, 通过调节茉莉酸通路来获得植物对病原菌的抗性, 从而形成自我保护． 同时, 菌核病抗性基因还和其他植物激素及信号转导相关, 也可能同植物自身免疫体系有着紧密的联系, 往后在油菜菌核病的研究中可以更加关注植物激素和免疫体系两个方面． 后续对候选基因功能验证等研究值得期待, 可能会对油菜抗菌核病育种起到极大的推动作用．

4 结论

将甘蓝型油菜高抗父本21Y490与高感母本21Y689杂交, 对F2群体室内接种核盘菌并进行表型鉴定, 筛选极端表型． 利用BSA-seq技术关联到两个显著区域, 即C09染色体的17.8～21.2 Mb和31.6～39.1 Mb区域, 共399个基因, 包含111个非同义突变或移码突变基因, 从中预测了8个候选基因, 这些基因主要参与植物激素生长调控及免疫系统的信号转导．