基于BSA-Seq和RNA-Seq方法鉴定大豆抗豆卷叶螟候选基因

曾维英 赖振光 孙祖东 杨守臻 陈怀珠 唐向民

基于BSA-Seq和RNA-Seq方法鉴定大豆抗豆卷叶螟候选基因

曾维英**赖振光**孙祖东*杨守臻 陈怀珠 唐向民

广西农业科学院经济作物研究所/ 农业农村部西南玉米大豆间套作区农业科学观测实验站, 广西南宁 530007

豆卷叶螟(Fabricius)是重要的大豆食叶性害虫, 挖掘大豆抗豆卷叶螟相关基因对大豆抗虫品种选育和遗传改良至关重要。本研究用大豆高抗豆卷叶螟材料赶泰-2-2和高感豆卷叶螟材料皖82-178进行杂交构建F2代分离群体, 从303个F2代单株中挑选出高抗豆卷叶螟和高感豆卷叶螟的单株各30株, 分别构建2个极端性状的DNA混合池用于全基因组重测序以分析控制豆卷叶螟相关的候选基因。结果表明, 4个样本中共有11,963,077个单核苷酸多态性(SNPs)标记, 根据SNP-index方法关联分析, 共有329个基因位于99%置信区间外, 这些基因主要集中在7号染色体5,601,065～5,865,237 bp区间(总长为0.26 Mb)、16号染色体2,975,110～6,336,096 bp区间(总长为3.36 Mb)、18号染色体44,366,115～54,297,600 bp区间(总长为9.93 Mb)等区域内。将BSA-Seq结果与转录组测序结果进行联合分析发现, 有12个基因相关联; 最后, 结合生物信息学分析、候选基因的表达模式和基因的同源注释, 锁定、转录因子、、丝氨酸/苏氨酸蛋白激酶、等12个基因为控制豆卷叶螟性状相关的候选基因。本研究结果不仅为解析大豆抗豆卷叶螟的分子机理奠定重要的基础, 也将为大豆抗虫基因的克隆奠定坚实的理论基础。

大豆; 豆卷叶螟; BSA-Seq; RNA-Seq; 候选基因

大豆是我国植物蛋白和食用油的主要来源, 虫害是影响大豆产量和品质的重要因素。豆卷叶螟(Fabricius)分布于东北、黄淮海及南方地区[1], 是南方和黄淮大豆产区主要食叶性害虫[2-3], 一年可发生4～5代, 有世代重叠现象, 在温度、湿度适宜时会大面积发生, 一般造成大豆减产15%～20%, 严重时可达30%以上。为害严重的年份, 叶片被取食后只剩叶脉和叶柄, 造成产量的严重损失[3-5]。

鉴于豆卷叶螟为害的严重性, 国内外在大豆抗豆卷叶螟方向展开了深入研究, 已鉴定出高抗豆卷叶螟材料赶泰-2-2、丰平黑豆等, 也鉴定出了高感材料皖82-178、山东大豆等[6-7]; 控制豆卷叶螟的抗性位点主要位于A2、C2、D1a、D1b、H、K和O等连锁群上[8-10]; 可溶性糖含量、茉莉酸含量、过氧化氢酶和多酚氧化酶活性只与豆卷叶螟诱导有关, 而超氧化物岐化酶活性、乙烯和脱落酸含量与虫害诱导和材料基因型都有关[11]; 大豆胰蛋白酶抑制剂、查尔酮异构酶4、脂氧合酶9等可能是大豆抗豆卷叶螟潜在的靶标蛋白(基因)[12-14], gma-miR156q、gma-miR166u、gma-miR166b、gma-miR319d、gma-miR394a-3p和gma-miR396e等miRNA可能参与大豆抗豆卷叶螟的抗性调控[15]。

目前, 大豆抗虫基因主要源于微生物的苏云金芽孢杆菌(, Bt)杀虫晶体蛋白基因(), 并转入大豆中获得了抗虫转基因大豆。郭东全等[16]将构建的和高效双价杀虫基因转入到吉林20与吉林27大豆品种中, 获得了高效抗虫且基本稳定遗传的转基因抗虫大豆; 陈秀华等[17]将基因转化到绥农14中, 得到抗虫效果较好的转化植株; 武小霞等[18]将基因转入黑农35中, 检测得到抗大豆食心虫转基因植株; 朱延明等[19]将具有豆荚特异性启动子Pmsg的29基因转入绥农28中, 获得了对鳞翅目和双翅目都具有抗性植株; 蓝岚等[20]将抗虫基因转入东农50中, 证明转基因植株对大豆食心虫具有明显抗性; 高嵩[21]将抗虫基因转入到吉农28中发现, 转基因植株对大豆食心虫具有明显抗性; 美国孟山都公司、陶氏公司将、等目的基因导入大豆中, 获得对鳞翅目类昆虫具有强烈杀虫活性的转基因抗虫大豆品种, 并在加拿大、日本、美国展开商业化种植[22-24]。但抗虫基因属于外源基因, 转抗虫基因大豆对非靶标生物造成了潜在的威胁。因此, 需鉴定出大豆自身含有的抗虫目的基因并进行遗传转化, 获得抗虫转基因材料。

随着高通量测序技术的兴起, 基于全基因组测序的BSA方法在不构建遗传图谱的情况下能高效快速地挖掘目标性状基因, BSA-Seq方法已被广泛应用在水稻、黄瓜、番茄、大豆等作物中, 并且成功地定位出水稻耐盐[25]、水稻稻瘟病[26]、黄瓜早花[27]、番茄果实重量[28]、大豆疫霉病[29]等基因。本研究利用高抗豆卷叶螟材料赶泰-2-2与高感材料皖82-178杂交构建F2代分离群体, 通过对亲本和2个极端表型子代混合池进行BSA-Seq全基因组重测序, 定位目标性状关联区域, 推测与抗豆卷叶螟有关的基因; 再结合前期的转录组测序结果, 关注抗虫性状候选区域的相关基因的表达水平, 进一步缩小候选基因范围, 鉴定出控制豆卷叶螟性状的候选基因, 旨在为更深入的阐明大豆抗豆卷叶螟的分子机制开拓新思路, 为大豆的抗虫分子育种奠定理论基础。

1 材料与方法

1.1 试验材料

2018年利用皖82-178 (高感豆卷叶螟)为母本和赶泰-2-2 (高抗豆卷叶螟)为父本进行杂交, 试验材料(F2代群体和亲本材料)于2020年3月种植在广西壮族自治区农业科学院院本部试验田防虫网室中(东经108.1146°、北纬22.5059°, 红壤土), 垄长1.5 m、行距0.5 m, 3行区播种, 每行10株, 在大豆整个生育期期间不喷施农药, 不施用肥料。待植株长至10片复叶时, 按每株3虫的密度分别接豆卷叶螟4龄幼虫3头(用镊子取3头虫接到第9片复叶的中间小叶上), 48 h后以卷包数和卷叶率为评价指标进行抗性鉴定(卷叶率0～10%为高抗、10%～30%为抗、30%～50%为中抗、50%～70%为感、>70%为高感)[7-8]。从303个F2分离群体中挑选30株极端抗豆卷叶螟单株和30株极端感豆卷叶螟单株构建子代高抗混合池和高感混合池, 用于抗豆卷叶螟的基因定位。田间采集相应植株幼嫩叶片。

1.2 试验方法

1.2.1 DNA提取 采用改良CTAB法[30]分别提取亲本材料(赶泰-2-2和皖82-178)和F2代群体极端个体的DNA, 将极端抗豆卷叶螟植株幼嫩叶片和极端感豆卷叶螟植株幼嫩叶片各30株的DNA等量混合, 制备总DNA 3 μg, 样品浓度≥20 ng µL-1, 构建高抗子代混合池、高感子代混合池。

1.2.2 文库构建及测序 由深圳华大基因科技服务有限公司完成样本建库和测序。将检测合格的基因组DNA样本用Covaris仪超声波随机打断成特定大小的片段, 打断后的样本用Agencourt AMPure XP-Medium kit进行片段选择, 使得样本条带集中在200～300 bp左右, 使用试剂盒Qubit dsDNA HS AssayKit, 500 assays检测纯化后DNA样本量; 对打断的DNA片段进行末端修复, 在3¢端连接A碱基, 将特定测序引物连接到DNA片段上; 纯化连接产物, 选择合适大小的DNA片段; 将DNA片段在cBot仪器上扩增制备文库; 对构建的文库进行质量检测; 将质量检测合格的文库上机测序。对亲本池和子代混合池分别开展20×和50×覆盖度的全基因组重测序。

1.2.3 生物信息学分析 对测序得到的Raw reads进行过滤, 将过滤后的Clean reads运用BWA软件与大豆参考基因组进行比对[31-32], 将比对结果进行排序、质量过滤和标记duplicate后进行变异类型检测。使用GATK软件查找样本相对于参考基因组之间的SNP, 用Breakdancer[33]软件来检测样本与参考基因组之间的SV位点, 用Control-FREEC[34]软件来查找样本相对于参考基因组之间的CNV区域。为得到更准确的变异信息, 对变异位点分别进行进一步的过滤, 使用变异位点注释软件Annovar[35]将变异位点在基因组上的位置区域注释出来。

1.2.4 SNP-index计算 子代群体与亲本之间的序列差异程度用SNP-index来表示, 计算经过质量控制后的每一个变异位点在高抗子代池和高感子代池中的SNP-index, 过滤掉2个子代极端混合池里面SNP-index都小于0.3的位点, 只保留双亲纯合且有差异、子代混池都为非缺少的SNP位点。同时, 将2个子代池的SNP-index相减得到Delta SNP-index。Delta SNP-index为1或者-1, 说明相对应的SNP来源于其中一个亲本, 这一位点所在的基因很有可能是与目的性状相关的基因。

1.2.5 目的基因的筛选与功能注释 使用滑窗法对滑窗内的Delta SNP-index进行计算, 取滑窗内所有SNP位点的Delta SNP-index的平均值, 滑窗大小为1.00 Mb, 步长为0.05 Mb, 对Delta SNP-index在各染色体上的分布进行作图。基于基因组Gff注释文件, 对所有候选位点进行来源注释, 挑选出位于基因区的SNP位点, 为缩小筛选范围, 选取99%置信区间作为筛选阈值, 单独把Delta SNP-index大于99%置信区间的SNP位点筛选出来进行注释。利用Blast2go_v2.5软件对所鉴定的基因进行GO注释分析, 计算得到的-value通过Bonferroni校正之后, 以Corrected-value≤0.05为阈值。利用Blast_v2.2.26软件在KEGG pathway数据库对所鉴定的基因进行Pathway富集分析, 计算得到的- value通过Bonferroni校正之后, Corrected-value≤0.05为阈值。

1.2.6 转录组分析 利用TRIzol Kit提取叶片总RNA, 采用Truseq RNA试剂盒进行mRNA纯化并构建cDNA文库, 对获得的cDNA文库进行PCR扩增富集, 2%琼脂糖凝胶电泳检测, 回收目的片段。然后使用Illumina Hiseq 2000测序仪进行对样本测序。以上试验均由深圳华大基因科技服务有限公司完成。测序后的数据经Base calling转化为Raw data或Raw reads, 随后对Raw reads进行质控。利用SOAPfuse软件[36]对原始测序数据进行去除处理得到Clean reads。再利用Tophat将Clean reads比对到大豆参考基因组, 允许有2个碱基的错配[37-38]。表达定量的结果以FPKM (Fragments per kilobase of exon model per million mapped fragments)为单位, 在得到差异检验的FDR值同时, 根据基因的表达量(FPKM值)计算该基因在不同样本间的差异表达倍数。以“FDR≤0.001和|log2Ratio (FC, Fold Change)|≥1”作为阈值来判断基因差异表达的显著性[39]。

2 结果与分析

2.1 测序质量评估

利用高通量测序技术对亲本及子代混合池共4个样本进行全基因组重测序, 共产生194.54 G原始测序数据, 经过SOAPnuke软件对数据进行接头过滤、低质量过滤和去N过滤后, 最终得到191.35 G高质量的Clean Reads数据。GC含量在35.46%～ 35.60%之间, 测序数据质量较高(Q20≥96.29%, Q30≥91.70%) (表1)。由此可知4个样本的数据量充足, 测序数据可靠, 测序质量合格, GC分布正常, 可进行下一步分析。

使用BWA软件(https://webblast.ipk-gatersleben. de/barley_ibsc/downloads/), 将4个样本中质控后的Clean reads比对到大豆参考基因组上, 然后使用Qualimap软件对比对结果进行统计。其中, 基因组大小是979,046,046 bp, 有效基因组大小为955,932,237 bp (参考序列中不含N), 参考基因组GC含量为33.95%, 4个样本的比对率介于99.76%和99.81%之间, 平均测序深度在30.78×和65.98×之间浮动(表2)。数据比对率高, 有利于下一步进行SNP检测。

将测序Reads比对到参考基因组以后, 统计参考基因组上不同染色体区域的覆盖度和测序深度分布情况。1×覆盖率百分比超过95.75%, 5×覆盖率百分比超过94.55%, 10×覆盖率百分比超过93.29%, 20×覆盖率百分比超过83.88% (表2)。通过Qualimap软件对比对结果进行统计, 参考基因组被均匀覆盖, 随机性良好, 有利于SNP过滤和筛选。

表1 原始测序数据统计表

表2 质控后测序数据与参考基因组比对结果统计分析

2.2 SNP检测和注释

SNP分析结果显示, 4个样本中共获得11,963,077个SNP位点, 其中高感子代混合池中最少, 皖82-178 (高感亲本)中最多。SNP位置和编码信息统计显示, 同类型材料间(亲本间, 子代池间)相同类型的SNP数和比例大体相当(表3)。纯合突变SNP比率明显高于杂合突变SNP比率, 纯合突变SNP 比率均高于99.73%。

2.3 候选基因筛选

利用SNP-index方法对检测到的SNP进行分析, 4个样本共鉴定出高质量的可信SNP位点11,963,077个。根据SNP-index方法关联阈值判断, 所有候选位点对应的基因数量为56,909个, 当拟合后Delta SNP-index置信度为99%时, 20条染色体上共有329个基因位于99%置信区间外, 主要集中在7号染色体5,601,065～5,865,237 bp区间(总长为0.26 Mb, 区间内基因数量为25个)、16号染色体2,975,110～6,336,096 bp区间(总长为3.36 Mb, 区间内基因数量为36个)、18号染色体44,366,115～ 54,297,600 bp区间(总长为9.93 Mb, 区间内基因数量为52个)等区域内, 区间内共包含113个基因(图1)。1号染色体和19号染色体上没有检测到基因, 2号、6号、8号、9号、10号、11号和12号染色体上检测到的基因数少于10个。

表3 多态性位点注释统计表

为进一步分析基因所在的亚细胞定位、分子功能及参与的生物学过程, 对所鉴定出的329个基因进行GO功能注释分析表明, 329个基因中有181个被注释到20个功能组中(图2), 包括9个生物学进程、7个细胞组分和4个分子功能, 在生物学进程中, 这些基因主要参与刺激反应、单一生物体过程、代谢过程、细胞过程等生物学过程; 在细胞组分中, 这些基因主要集中在细胞、细胞组分、膜、膜组分等组分; 在分子功能中, 这些基因主要参与结合、催化活性等功能。说明虫害胁迫下大豆对逆境的响应可能与细胞组织、细胞膜组织等的修复功能有关, 且一些催化活性的酶可能参与到虫害逆境胁迫中。

为鉴别329个基因参与富集的代谢途径, 利用KEGG通路数据库对候选基因进行Pathway富集分析表明, 329个基因中有30个基因主要富集在植物病原菌互作、代谢途径、次生代谢物的生物合成、RNA降解等通路中(图3)。表明当大豆受到豆卷叶螟危害以后, 植株体内的防御系统会立即响应刺激, 适当地增加体内的代谢活动, 产生防御物质, 如防御酶、蛋白质抑制酶等物质, 而且增强各种酶的活性来促进防御。

横坐标为染色体编号, 纵坐标为Delta SNP-index值。不同颜色的点表示不同染色体上筛选后的SNP, 红色曲线表示滑窗后的Delta SNP-index值, 蓝色直线表示99%置信区间值。

The X-axis is the chromosome number, Y-axis is the Delta SNP-index value. Different colored dots represent SNPs screened on different chromosomes, the red curve represents the Delta SNP-index value after the sliding window, the blue line represents the 99% confidence interval.

将BSA-Seq结果预测的基因和转录组测序结果[13]进行联合分析, 将位于99%置信区间以外SNP对应的基因与转录组数据进行比对, 获得转录组中的同源序列及基因的表达水平。比对发现共有12个基因在4个比较组中呈显著性差异表达, 这些基因分别位于2号、12号、16号、18号等9条染色体上(表4), 其中响应豆卷叶螟取食胁迫诱导的基因有8个, 这8个基因中共有7个为显著上调基因, 分别为环核苷酸门控离子通道4 ()、转录因子、赤霉素2-β加双氧酶8、氨基酸透性酶7 ()、Mdis1相互受体激酶2、丝氨酸/苏氨酸蛋白激酶、, 只有纤维素合成酶A催化亚基4 ()为显著下调基因; 类枯草杆菌蛋白酶、抗病蛋白RGA3和硫氰酸酶结构域蛋白酶10等3个基因在高感材料中的表达量显著高于高抗材料, 1个未知基因在在高感材料中的表达量显著低于高抗材料。最后, 结合基因的同源注释、基因差异表达分析和生物信息学分析等, 推测环核苷酸门控离子通道4 ()、转录因子、7、丝氨酸/苏氨酸蛋白激酶、等12个基因为抗豆卷叶螟有关的候选基因。

3 讨论

基于BSA-Seq的方法挖掘新基因在作物农艺性状相关的基因定位中的应用日益广泛, 该方法在构建遗传图谱的情况下高效快速地挖掘目的基因[40], 是一种简单快速、准确的基因定位方法, 已成为育种家进行快速QTL定位和挖掘各种不同作物目标性状基因的一种有效方法。如Takagi等[25]利用BSA-Seq方法找到控制水稻耐盐性的候选基因, 仅用2年时间培育出水稻耐盐性新品种; Lu等[27]利用BSA-Seq方法定位到84个控制黄瓜早花性状基因并确定基因为候选基因; Illa- Berenguer等[28]利用BSA-Seq方法定位到66个控制番茄果实重量性状基因; Guo等[41]利用BSA-Seq方法快速定位到29个辣椒抗黄瓜花叶病毒性状基因并挖掘出2个候选基因; Ma等[42]利用BSA-Seq方法定位到控制水稻低柱头外露突变候选基因; Zhao等[43]利用BSA-Seq方法鉴定到20个棉花细胞质雄性不育的育性恢复候选基因; 张尧锋等[44]利用BSA-Seq方法在甘蓝型油菜中定位到8个控制有限花序性状候选基因。BSA-Seq方法也被用于农作物的抗虫基因定位, Liang等[45]利用BSA-Seq方法定位到1个与黄瓜蚜虫抗性相关的区域, 该区域内包含40个基因, 其中16个基因可能与蚜虫抗性有关。利用全基因组重测序与BSA相结合同样能够高效和快速鉴定大豆中的目标性状基因[46], 本研究利用BSA-Seq分析和转录组测序分析对控制抗豆卷叶螟性状相关基因进行了初步定位, 筛选出12个可能与抗豆卷叶螟胁迫相关的候选基因, 再结合生物信息学分析、基因差异表达分析和基因同源性分析发现, 环核苷酸门控离子通道4 ()、转录因子、氨基酸透性酶7 ()、丝氨酸/苏氨酸蛋白激酶、、类枯草杆菌蛋白酶等6个基因已有报道与抗逆有关。而赤霉素2-b加双氧酶8、Mdis1相互受体激酶2、抗病蛋白RGA3、硫氰酸酶结构域蛋白酶10和1个未知基因可能是与抗大豆卷叶螟有关的新基因。

表4 候选基因信息表

HR代表抗豆卷叶螟材料赶泰-2-2和HS为高感豆卷叶螟材料皖82-178。数字0和48代表处理时间为0 h和48 h。a上调或者下调代表显著性上调或者显著性下调。

HR represents the highly resistant material Gantai-2-2, HS represents the highly susceptible material Wan 82-178. The numbers 0 and 48 represent the processing time of 0 hour and 48 hours.aThe ‘up’ or ‘down’ represents significantly different up or down.

环核苷酸门控离子通道()是一种非选择性的阳离子通道, 对一价和二价阳离子均有通透性, 是动植物细胞中非常重要的离子通道, 在植物生长发育、生物与非生物胁迫中都发挥着重要作用[47-48]。植物细胞中,是信号转导级联系统的组成部分, 是环核苷酸作用的一个主要靶标, 是连接环核苷酸和Ca2+的纽带。在病原菌入侵的早期, 它能够参与Ca2+内流的调控[49-50], 将细胞外信号通过阳离子流转变为胞内信号, 对细胞的生理活动进行调控[51]。当细胞受到外界刺激时, 定位于细胞膜上的受体蛋白识别相应刺激, 并激活胞质内腺苷酸环化酶, 生成的环核苷酸(cAMP)使通道打开, 胞外Ca2+内流[52]。转录因子是植物中最大的转录因子家族之一, 在调控植物的生长发育、生物及非生物胁迫等过程中具有十分重要的作用[53-55]。植物遭受昆虫胁迫后, 植物体内转录因子的表达水平会发生变化[56], 如拟南芥受到线虫危害时,基因被诱导[57]。氨基酸透性酶()是植物中的氨基酸转运蛋白, 在氨基酸转运过程中发挥着重要作用[58-59]。基因在植物抵御病原体侵害方面扮演着重要的角色, 如含基因和基因突变体植株能维持根结线虫生长率低于野生型[60]。丝氨酸/苏氨酸蛋白激酶是一类重要的信号分子, 当植物遭遇到如害虫取食、盐害、干旱等胁迫, 或是受到创伤以及其他细胞因子、激素刺激时, 丝氨酸/苏氨酸蛋白激酶迅速在丝氨酸、苏氨酸残基部位磷酸化而激活, 并进一步通过级联磷酸化作用活化下游的信号分子, 激活特定的信号传导途径, 最终将外界信号传递到细胞核, 激活或抑制特定基因的表达[61-63]。基因的表达对植物株型发育有重要的调控作用, 该基因的表达水平受植物激素和光照的调节, 最终导致植物株型的变化[64]。类枯草杆菌蛋白酶是一类在植物发育和信号级联中实现高度特异性功能的丝氨酸蛋白酶, 在植物和病原物互作过程中起到重要的免疫激发作用[65]。推测以上候选基因的功能和参与的代谢途径可能与大豆抗豆卷叶螟有关, 当大豆受豆卷叶螟诱导胁迫后, 一些候选基因能发生信号转导、相关基因表达及防御物质生成等一系列的反应来抵御虫害。大豆抗豆卷叶螟是一个复杂的生理过程, 目前尚缺乏遗传调控方面的研究, 本研究中通过BSA-Seq技术与RNA-Seq技术相结合, 挖掘出一些与抗豆卷叶螟相关的候选基因, 对于这些基因是否是调控豆卷叶螟的关键基因, 以及其调控机理将在后续的工作中进行分析和验证。

4 结论

通过BSA-Seq全基因组重测序发现共有329个基因位于99%置信区间以外, 主要集中在7号、16号和18号等染色体上。通过结合转录组测序结果分析、基因的同源注释、生物信息学分析等, 初步推测、转录因子、7、丝氨酸/苏氨酸蛋白激酶、等12个基因为大豆抗豆卷叶螟相关的候选基因, 这些基因可能在大豆抗豆卷叶螟过程中起着重要作用, 为大豆抗豆卷叶螟基因图位克隆奠定了基础。

[1] 中国农作物病虫图谱编写组. 中国农作物病虫图谱: 第五分册, 油料病虫(一). 北京: 中国农业出版社, 1982. pp 136–137. Editorial Committee of Plate of Chinese Diseases and Insects on Crop. Plate of Chinese Diseases and Insects on Crop, Fifth Fascicule, Diseases and Insects on Oil Crop (first). Beijing: China Agriculture Press, 1982. pp 136–137 (in Chinese).

[2] 崔章林, 盖钧镒, 吉东风, 任珍静. 大豆种质资源对食叶性害虫抗性的鉴定. 大豆科学, 1997, 16: 93–102. Cui Z L, Gai J Y, Ji D F, Ren Z J. A study on leaf-feeding insect species on soybeans in Nanjing area., 1997, 16: 93–102 (in Chinese with English abstract).

[3] 孙祖东, 杨守臻, 陈怀珠, 韦德卫. 南宁大豆食叶性害虫调查. 广西农业科学, 2001, 32: 104–106. Sun Z D, Yang S Z, Chen H Z, Wei D W. Investigation on leaf-feeding pests of soybean in Nanning., 2001, 32: 104–106 (in Chinese with English abstract).

[4] 孙祖东, 盖钧镒. 大豆对食叶性害虫抗性的研究. 中国农业科学, 1999, 32(增刊1): 81–88. Sun Z D, Gai J Y. Study on resistance of soybean to leaf-feeding insect., 1999, 32(S1): 81–88 (in Chinese with English abstract).

[5] 崔章林, 盖钧镒. 大豆抗食叶性害虫研究进展. 大豆科学, 1996, 15: 149–158. Cui Z L, Gai J Y. Advance of study on soybean leaf-feeding insects., 1996, 15: 149–158 (in Chinese with English abstract).

[6] 崔章林, 盖钧镒, 吉东风, 任珍静. 南京地区大豆食叶性害虫种类调查与分析. 大豆科学, 1997, 16: 12–20. Cui Z L, Gai J Y, Ji D F, Ren J Z. A study on leaf-feeding insect species on soybeans in Nanjing area., 1997, 16: 12–20 (in Chinese with English abstract).

[7] 孙祖东, 杨守臻, 陈怀珠, 李初英, 龙丽萍. 大豆对豆卷叶螟的抗性鉴定. 中国油料作物学报, 2005, 27(4): 69–71. Sun Z D, Yang S Z, Chen H Z, Li C Y, Long L P. Identification of soybean resistance to bean pyralid (Fabricius) and oviposition preference of bean pyralid on soybean varieties., 2005, 27(4): 69–71 (in Chinese with English abstract).

[8] Xing G N, Zhou B, Wang Y F, Zhao T J, Yu D Y, Chen S Y, Gai J Y. Genetic components and major QTL confer resistance to bean pyralid (Fabricius) under multiple environments in four RIL populations of soybean., 2012, 125: 859–875.

[9] 李广军, 程利国, 张国政, 何小红, 智海剑, 章元明. 大豆对豆卷叶螟抗性的主基因+多基因混合遗传. 大豆科学, 2008, 27: 33–36. Li G J, Cheng L G, Zhang G Z, He X H, Zhi H J, Zhang Y M. Mixed major-gene plus polygenes inheritance analysis for resistance in soybean to bean pyralid (Fabricius)., 2008, 27: 33–36 (in Chinese with English abstract).

[10] 李广军, 李河南, 程利国, 章元明. 大豆对豆卷叶螟抗性的QTL定位. 中国油料作物学报, 2009, 31: 365–369. Li G J, Li H N, Cheng L G, Zhang Y M. Mapping quantitative trait loci for resistance in soybean to bean pyralid (Fabricius)., 2009, 31: 365–369 (in Chinese with English abstract).

[11] 曾维英, 蔡昭艳, 张志鹏, 陈怀珠, 杨守臻, 唐向民, 赖振光, 孙祖东. 大豆抗豆卷叶螟的生理生化特性研究. 南方农业学报, 2016, 46: 2112–2116. Zeng W Y, Cai Z Y, Zhang Z P, Chen H Z, Yang S Z, Tang X M, Lai Z G, Sun Z D. Physiological and biochemical characteristics of(Fabricius)-resistant soybean., 2015, 46: 2112–2116 (in Chinese with English abstract).

[12] Zeng W Y, Sun Z D, Cai Z Y, Chen H Z, Lai Z G, Yang S Z, Tang X M. Proteomic analysis by iTRAQ-MRM of soybean resistance to., 2017, 18: 444.

[13] Zeng W Y, Sun Z D, Cai Z Y, Chen H Z, Lai Z G, Yang S Z, Tang X M. Comparative transcriptome analysis of soybean response to bean pyralid larvae., 2017, 18: 871.

[14] 曾维英, 孙祖东, 赖振光, 蔡昭艳, 陈怀珠, 杨守臻, 唐向民. 大豆抗豆卷叶螟的转录组和蛋白质组关联分析. 中国农业科学, 2018, 51: 1244–1260. Zeng W Y, Sun Z D, Lai Z G, Cai Z Y, Chen H Z, Yang S Z, Tang X M. Correlation analysis on transcriptomic and proteome of soybean resistance to bean pyralid ()., 2018, 51: 1244–1260 (in Chinese with English abstract).

[15] Zeng W Y, Sun Z D, Lai Z G, Yang S Z, Chen H Z, Yang X H, Tao J R, Tang X M. Determination of the miRNAs related to bean pyralid larvae resistance in soybean using small RNA and transcriptome sequencing., 2019, 20: 2966.

[16] 郭东全, 杨向东, 包绍君, 包绍君, 郭三堆, 康岭生, 尹爱萍, 钱雪艳, 赵桂兰. 转和双价抗虫基因大豆的获得与稳定表达. 中国农业科学, 2008, 41: 2957–2962. Guo D Q, Yang X D, Bao S J, Guo S D, Kang L S, Yin A P, Qian X Y, Zhao G L. Synchronous expression ofandgenes in soybean and analysis of their resistance to insect pests., 2008, 41: 2957–2962 (in Chinese with English abstract).

[17] 陈秀华, 柏锡, 潘欣, 翟红, 才华, 纪魏, 李勇, 朱延明. 转基因大豆的培育及抗虫性检测. 大豆科学, 2009, 28: 959–963. Chen X H, Bai X, Pan X, Zhai H, Cai H, Ji W, Li Y, Zhu Y M. Cultivation ofgene transformed soybean and insect resistant assay., 2009, 28: 959–963 (in Chinese with English abstract).

[18] 武小霞, 李静, 王志坤, 刘珊珊, 李海燕, 武天龙, 李文滨.基因转化大豆及抗虫性的初步评价. 上海交通大学学报(农业科学版), 2010, 28: 413–419. Wu X X, Li J, Wang Z K, Liu S S, Li H Y, Wu T L, Li W B. Transformation ofinto soybean and rough assessing its resistance to soybean pests.(Agric Sci), 2010, 28: 413–419 (in Chinese with English abstract).

[19] 朱延明, 郜庭, 张凤, 柏锡, 才华, 纪魏, 罗晓.29抗虫基因植物表达载体构建及对大豆的遗传转化. 东北农业大学学报, 2013, 44(1): 1–6. Zhu Y M, Gao T, Zhang F, Bo X, Cai H, Ji W, Luo X. Construction of plant expression vector of insect-resistant gene29and transformation intoL. Merr., 2013, 44(1): 1–6 (in Chinese with English abstract).

[20] 蓝岚, 吴帅, 申丽威, 王志坤, 孟凡立, 宋波, 拓云, 刘珊珊. 根癌农杆菌介导大豆转抗虫基因. 中国油料作物学报, 2013, 35: 29–35. Lan L, Wu S, Shen L W, Wang Z K, Meng F L, Song B, Tuo Y, Liu S S. Transgenic of soybean withgene mediated by., 2013, 35: 29–35 (in Chinese with English abstract).

[21] 高嵩. 抗虫基因在大豆中的遗传转化及抗虫性鉴定. 吉林农业大学硕士学位论文, 吉林长春, 2015. Gao S. Identification on Insect-resistant Gene ofTransformation intoL. Merr. MS Thesis of Jilin Agricultural University, Changchun, Jilin, China, 2015 (in Chinese with English abstract).

[22] Berman K H, Harrigan G G, Riordan S G. Compositions of seed, forage, and processed fractions from insect-protected soybean MON 87701 are equivalent to those of conventional soybean., 2009, 57: 11360–11369.

[23] Beazley K A, Burns W C, Cole II R H, Macrae T C, Miklos J A, Ruschke L G, Tian K R, Wei L P, Wu K S. Soybean transgenic event mon87751 and methods for detection and use thereof. 14303042. U.S. Patent Application 14/303, 042. 2014-06-12.

[24] Fast B J, Schafer A C, Johnson T Y, Potts B L, Herman R A. Insect-protected event DAS-81419-2 soybean (L.) grown in the United States and Brazil is compositionally equivalent to nontransgenic soybean., 2015, 63: 2063–2073.

[25] Takagi H, Tamiru M, Abe A, Yoshide K, Uemura A, Yaegashi H, Obara T, Oikawa K, Utsushi H, Kanzaki E, Mitsuoka C, Natsume S, Kosugi S, Kanzaki H, Matsumrua H, Urasaki N, Kamoun S, Terauchi R. MutMap accelerates breeding of a salt-tolerant rice cultivar., 2015, 33: 445–449.

[26] Takagi H, Abe A, Yoshide K, Kosugi S, Natsume S, Mitsuoka C, Uemura A, Utsushi H, Tamiru M, Innan H, Cano L M, Kamoun S, Terauchi R. QTL-seq: rapid mapping of quantitative trait loci in rice by whole genome resequencing of DNA from two bulked populations., 2013, 74: 174–183.

[27] Lu H F, Liu T, Joël K, Wang S H, Qi J J, Zhou Q, Sun J J, Zhang Z H, Weng Y Q, Huang S W. QTL-seq identifies an early flowering QTL located near flowering locus T in cucumber., 2014, 127: 1491–1499.

[28] Illa-Berenguer E, Houten J V, Huang Z J, Knaap E. Rapid and reliable identification of tomato fruit weight and locule number loci by QTL-seq., 2015, 128: 1329–1342.

[29] Zhong C, Sun S, Li Y, Duan C, Zhu Z. Next-generation sequencing to identify candidate genes and develop diagnostic markers for a novel Phytophthora resistance gene, RpsHC18, in soybean., 2018, 131: 525–538.

[30] Doyle J J, Doyle J L. Isolation of plant DNA from fresh tissue., 1990, 12: 13–15.

[31] Li H, Duibin R, Notes A. Fast and accurate short-read alignment with Burrows-Wheeler transform., 2009, 25: 1754–1760.

[32] Li H, Duibin R. Fast and accurate long-read alignment with Burrows-Wheeler transform., 2010, 26: 589–595.

[33] McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabrie S, Daly M, Depristo M A. The genome analysis toolkit: a MapReduce framework for analyzing next generation DNA sequencing data., 2010, 20: 1297–1303.

[34] Balentin A B, Tatiana P, Kevin B, Pierre C, Julie C, Gudrun S, Isabelle J L, Olivier D, Emmanuel B. Control-FREEC: a tool for assessing copy number and allelic content using next generation sequencing data., 2012, 28: 423–425.

[35] Wang K, Li M, Hakonarson H. ANNOVAR: Functional annotation of genetic variants from next-generation sequencing data., 2010, 38: e164.

[36] Jia W, Qiu K, He M, Song P, Zhou Q, Zhou F, Yu Y, Zhu D, Nickerson M L, Wan S, Liao X, Zhu X, Peng S, Li Y, Wang J, Guo G. SOAPfuse: an algorithm for identifying fusion transcripts from paired-end RNA-Seq data., 2013, 14: R12.

[37] Trapell C, Pacher L, Salzberg S L. TopHat: discovering splice junctions with RNA-Seq., 2009, 25: 1105–1111.

[38] Trapell C, Roberts A, Goff L, Pertea G, Kim D, Kelley D R, Pimental H, Salzberg S L, Rinn J L, Pachter L. Differential gene and transcript expression analysis of RNA-Seq experiments with TopHat and Cufflinks., 2012, 7: 562–578.

[39] Benjamini Y, Yektieli D. The control of the false discovery rate in multiple testing under dependency., 2001, 29: 1165–1188.

[40] 张之昊, 王俊, 刘章雄, 邱丽娟. 基于BSA-seq技术挖掘大豆中黄622的多小叶基因. 作物学报, 2020, 46: 1839–1849. Zhang Z H, Wang J, Liu Z X, Qiu L J. Mapping of an incomplete dominant gene controlling multifoliolate leaf by BSA-seq in soybean (L.)., 2020, 46: 1839–1849 (in Chinese with English abstract).

[41] Gao G J, Wang S B, Liu J B, Pan B G, Diao W P, Ge W, Gao C Z, Snyder J C. Rapid identification of QTLs underlying resistance to cucumber mosaic virus in pepper (s)., 2017, 130: 41–52.

[42] Ma X, Zheng Z, Lin F S, Ge T T, Sun H M. Genetic analysis and gene mapping of a low stigma exposed mutant gene by high- throughput sequencing., 2018, 13: e0186942.

[43] Zhao C P, Zhao G Y, Geng Z, Wang Z X, Wang K H, Liu S, Zhang H S, Guo B S, Geng J Y. Physical mapping and candidate gene prediction of fertility restorer gene of cytoplasmic male sterility in cotton., 2018, 19: 6.

[44] 张尧锋, 张冬青, 余华胜, 林宝刚, 华水金, 丁厚栋, 傅鹰. 基于极端混合池(BSA)全基因组重测序的甘蓝型油菜有限花序基因定位. 中国农业科学, 2018, 51: 3029–3039. Zhang Y F, Zhang D Q, Yu H S, Lin B G, Hua S J, Ding H D, Fu Y. Location and mapping of the determinate growth habit ofby bulked segregant analysis (BSA) using whole genome re-sequencing., 2018, 51: 3029–3039 (in Chinese with English abstract).

[45] Liang D N, Chen M Y, Qi X H, Xu Q, Zhou F C, Chen X H. QTL Mapping by SLAF-seq and expression analysis of candidate genes for aphid resistance in cucumber., 2016, 7: 1000.

[46] Song Q J, Jenkins J W, Hyten D L. Construction of high resolution genetic linkage maps to improve the soybean genome sequence assembly Glyma1.01., 2016, 17: 1.

[47] 王正朝, 黄瑞华, 潘玲梅, 李学斌, 石放雄. 环核苷酸门控离子通道的结构、功能及活性调节. 中国生物化学与分子生物学报, 2006, 22: 282–288. Wang Z C, Huang R H, Pan L M, Li X B, Shi F X. Molecular structures, physiological roles and regulatory mechanisms of cycli nucleotide-gated ion channels., 2006, 22: 282–288 (in Chinese with English abstract).

[48] 吴巨友, 薛亚男, 张绍铃. 植物环核苷酸门控离子通道基因的功能及其调控. 西北植物学报, 2020, 30: 1716–1720. Wu J Y, Xue Y N, Zhang S L. Function and modulation of plant cyclic nucleotide-gated channels., 2020, 30: 1716–1720 (in Chinese with English abstract).

[49] Dangl J L, Dietrich R A, Richberg M H. Death don’t have no mercy: cell death programs in plant-microbe interactions., 1996, 8: 1793–1807.

[50] Hetherington A M, Brownlee C. The generation of Ca2+signals in plants., 2004, 55: 401–427.

[51] Flynn G E, Johnson J P, Zagotta W N. Cyclic nucleotide-gated channels: shedding light on the opening of a channel pore., 2001, 2: 643–651.

[52] 刘海娇, 杜立群, 林金星, 李瑞丽. 植物环核苷酸门控离子通道及其功能研究进展. 植物学报, 2015, 50: 779–789. Liu H J, Du L Q, Lin J X, Li R L. Recent advances in cyclic nucleotide-gated ion channels with their functions in plants., 2015, 50: 779–789 (in Chinese with English abstract).

[53] Eulgem T, Rushton P J, Schemlzer E. Early nuclear events in plant defence signalling: rapid gene activation by WRKY transcription factors., 1999, 18: 4689–4699.

[54] Eulgem T, Rushton P J, Robatzek S, Somssich I E. The WRKY superfamily of plant transcription factors., 2000, 5: 199–206.

[55] Rushton P J, Somssich I E, Ringler P, Shen Q. WRKY transcription factors., 2010, 15: 247–258.

[56] Zhou X, Jiang Y, Yu D. WRKY22 transcription factor mediates dark-induced leaf senescence in., 2011, 31: 303–313.

[57] Grunewald W, Karimi M, Wieczorek K, Van D E, Cappelle E, Chnitzki E, Grundler F, Iize D, Beeckman T, Gheysen G. A role forin feeding site establishment of plant-parasitic nematodes., 2008, 148: 358–368.

[58] Tegeder M. Transporters for amino acids in plant cells: some functions and many unknowns., 2012, 15: 315–321.

[59] 李倩, 王罡, 张松皓, 杨丹, 王昱蓉, 季静, 安婷, 李辰, 马志刚, 史怀宇, 关春峰, 刘玉. 拟南芥基因的克隆与转化马铃薯的研究. 天津大学学报(自然科学与工程技术版), 2018, 51: 941–948. Li Q, Wang G, Zhang S H, Yang D, Wang Y R, Ji J, An T, Li C, Ma Z G, Shi H Y, Guan C F, Liu Y. Cloning and genetic transformation ofgene in potato.(Sci Technol), 2018, 51: 941–948 (in Chinese with English abstract).

[60] Marella H H, Nielsen E, Schachtman D P, Taylor C G. The amino acid permeasesandare involved in root-knot nematode parasitism of., 2013, 26: 44–54.

[61] Afzal A J, Wood A J, Lightfoot D A. Plant receptor-like serine threonine kinases: roles in signaling and plant defense., 2008, 21: 507–517.

[62] Hasegawa P M, Bressan R A, Zhu J K, Bohnert H J. Plant cellular and molecular responses to high salinity., 2000, 51: 463–499.

[63] Zipfel C. Pattern-recognition receptors in plant innate immunity., 2008, 20: 10–16.

[64] 史伟杰, 刘建伟, 张彦峰, 王晓峰. 拟南芥株型突变体的性状特征与基因鉴定分析. 西北植物学报, 2014, 34: 2153–2158. Shi W J, Liu J W, Zhang Y F, Wang X F. Characteristics and genetic analysis of plant type mutantin., 2014, 34: 2153–2158 (in Chinese with English abstract).

[65] 陈发晶, 谭蕊, 黄萌雨, 杜娇, 余洋, 杨宇衡, 毕朝位. 类枯草杆菌蛋白酶在植物和病原物互作中的研究进展. 分子植物育种, 2018, 16: 3146–3153. Chen F J, Tan R, Huang M Y, Du J, Yu Y, Yang Y H, Bi C W. Research progress of subtilases in the interaction of plant and pathogen., 2018, 16: 3146–3153 (in Chinese with English abstract).

Identification of the candidate genes of soybean resistance to bean pyralid (Fabricius) by BSA-Seq and RNA-Seq

ZENG Wei-Ying**, LAI Zhen-Guang**, SUN Zu-Dong*, YANG Shou-Zhen, CHEN Huai-Zhu, and TANG Xiang-Min

Institute of Economic Crops, Guangxi Academy of Agricultural Sciences / Southwest Experimental Station of Maize-Soybean Intercrop, Ministry of Agriculture and Rural Affairs, Nanning 530007, Guangxi, China

Bean pyralid is an important leaf-feeding insect in soybean. Identification of insect-tolerant genes from soybean has great significant to the crop insect-tolerant breeding and genetic improvement. In this study, an F2population with 303 individuals was constructed using insect-resistant line Gantai-2-2 and insect-sensitive line Wan 82-178. 30 F2insect-resistant individuals and 30 insect-sensitive individuals were selected respectively to construct two DNA pools which were used for the whole-genome re-sequencing. The results showed that there were a total of 11,963,077 single nucleotide polymorphism (SNPs) markers identified in two parental lines and two mixed pools. According to the association analysis of SNP-index method, a total of 329 genes were located outside the 99% confidence interval. These genes were mainly concentrated in the regions of 5,601,065–5,865,237 bp with a total of 0.26 Mb on chromosome 7, 2,975,110–6,336,096 bp with a total of 3.36 Mb on chromosome 16, and 44,366,115–54,297,600 bp with a total of 9.93 Mb on chromosome 18. Correlation analysis of BSA-Seq and transcriptome sequencing showed that 12 genes were correlated. Then, 12 candidate genes, including,transcription factor 16,7, serine/threonine protein kinase andwere identified by bioinformatics analysis, differential expression analysis, and homologous annotation. This study laid an important foundation for the analysis of the molecular mechanism of soybean resistance to bean pyralid and the cloning of anti-insect genes.

soybean; bean pyralid; BSA-Seq; RNA-Seq; candidate genes

10.3724/SP.J.1006.2021.04195

本研究由广西自然科学基金项目(2017GXNSFDA198037)和广西农业科学院科技发展基金项目(桂农科2020YM116, 2015YT58)资助。

This study was supported by the Natural Science Foundation of Guangxi Province, China (2017GXNSFDA198037) and the Science and Technology Development Fund of Guangxi Academy of Agricultural Sciences (Guinongke 2020YM116, 2015YT58).

孙祖东, E-mail: sunzudong639@163.com

**同等贡献(Contributed equally to this work)

曾维英, E-mail: zengweiying_1981@163.com; 赖振光, E-mail: 519229671@qq.com

2020-08-25;

2021-01-13;

2021-02-18.

URL: https://kns.cnki.net/kcms/detail/11.1809.S.20210216.1414.002.html