基于BSA-seq技术挖掘糙皮侧耳抗螨候选基因

李辉平， 骆 昕， 侯子强， 蒋 宁， 林金盛， 侯立娟， 徐 平， 马 林， 曲绍轩,

(1.江苏省农业科学院蔬菜研究所，江苏省高效园艺作物遗传改良重点实验室，江苏南京210014；2.江苏大学生命科学学院，江苏镇江212013)

糙皮侧耳(Pleurotusostreatus)，又称平菇、蚝菇，其栽培方式简单粗放、抗病虫害能力强、产量高，栽培技术需求低，受到全球食用菌生产者的欢迎，是世界重要栽培食用菌之一[1]。又因其富含膳食纤维、蛋白质、矿物质、维生素和多种活性物质，越来越受到消费者的关注和欢迎[2-4]。2020年全国平菇产量达到6.829 6×106t，是中国第三大栽培食用菌[5]。

腐食酪螨(Tyrophagusputrescentiae)是世界性仓储害虫之一，取食各种干制食品和生物原材料，对30多种食用菌造成危害，是食用菌栽培过程中需要防控的主要螨虫。该螨以刺吸式口器取食菌丝体和子实体，同时在培养基质上传播病原菌[6]。发菌期暴发螨虫可造成大规模减产甚至绝收，出菇期暴发螨虫则可减产30%～40%[7-8]。

目前，食用菌害螨主要依赖化学防治。但由于螨虫个体微小、生长周期短、变异快，已对多种药剂产生抗药性，越来越难防治[9]。培育和利用抗螨食用菌品种，可有效地解决田间螨害问题，降低防控成本，推进食用菌绿色安全生产。本研究团队前期研究发现，腐食酪螨在中国不同栽培糙皮侧耳品种上繁殖力、发育周期、产卵量、若螨死亡率等性状存在显著差异，表明国内不同平菇菌株对腐食酪螨存在不同的抗性水平[10]。然而，关于糙皮侧耳对螨虫抗性的来源和遗传机制尚不清楚，国内外也尚无抗螨性食用菌品种的相关报道。

群体分离分析法(Bulked segregant analysis，BSA)通过表型差异构建分离群体快速筛选和定位目标性状相关基因[11]。该方法在动植物功能基因定位研究中应用十分广泛[12-13]，在食用菌功能基因定位研究中也逐步得到了应用[14]。BSA-seq技术基于高通量测序技术，具有实验周期短、定位准确，开发的分子标记分布均匀、密度高等优点[15-16]。

本研究以一株野生糙皮侧耳菌株Pos045为试验材料，构建担孢子分离群体并开展抗螨性鉴定，分析其对腐食酪螨产生抗性的遗传基础，通过BSA-seq技术对潜在抗性基因进行初步定位，为开展糙皮侧耳抗螨分子育种打下基础。

1 材料与方法

1.1 试验材料

1.1.1 供试菌株 糙皮侧耳Pos045由江苏省农业科学院提供。

1.1.2 供试螨虫 腐食酪螨由江苏省农业科学院提供，室内继代饲养[10]。

1.1.3 培养基 PDA固体培养基，购自台湾MDBio公司，按说明书配制灭菌备用。棉籽壳出菇培养基，按质量百分比棉籽壳55%，玉米芯30%，麸皮15%，混匀后加水搅拌至含水量65%，用330 ml培养瓶装至瓶肩，121 ℃高压蒸汽灭菌90 min备用。平菇抗螨性鉴定培养基，棉籽壳200 g(煮水30 min过滤)，葡萄糖20 g，琼脂15 g，加水定容到1 L后121 ℃蒸汽高压灭菌30 min备用。

1.2 试验方法

1.2.1 孢子单核体群体构建 糙皮侧耳Pos045经PDA培养基活化后接种至棉籽壳出菇培养基，25 ℃下避光培养至菌丝体长满培养瓶，降温至18 ℃进行催蕾出菇。待菌盖开始舒展时，在超净台中悬于灭菌的锥形瓶中静置24 h，用200 μl无菌水冲洗瓶底孢子即为孢子悬浮液，4 ℃短期保存待用。将孢子悬浮液稀释至500 CFU/mL以下，取100 μl涂布90 mm平板， 25 ℃下避光培养5 d后，挑取单个菌落进行继代培养，显微镜检查有无锁状联合。挑取不少于400个单核菌株建立Pos045单核体群体。

1.2.2 抗螨性鉴定 以本实验室建立的平皿菌丝体接种螨虫取食孔洞计数法进行抗螨性鉴定。菌丝活化后接种至抗螨性鉴定培养基，菌丝长满平板后接种30 mg约500头螨虫，透气膜封闭25 ℃下继续培养，5个重复。接种螨虫48 h后，借鉴玉米抗玉米螟鉴定技术规范(NY/T 1248.5-2006)，观测螨虫取食菌丝后形成的孔洞数量和直径大小，统筹划分为害级别1～9级。计算为害程度的平均值，为害程度的平均值=∑(为害级别×该级别菌株数)/调查总数，并据此划分抗性等级，共分5级：1.0～2.9为高抗；3.0～4.9为抗；5.0～6.9为中抗；7.0～8.9为感；9.0为高感。

1.2.3 极端混池构建并测序 群体其他表型数据尽量一致，根据抗性数据，选择高抗/高感极端株各20株构建极端表型混池。单株分别使用MiniBEST Plant DNA纯化试剂盒(宝生物公司产品)提取DNA，等量混合构建DNA混池。利用Illumina HiSeq平台(北京百迈克生物科技有限公司产品)完成建库测序，单池测序数据2 G以上。

1.3 数据分析

基于NCBI公布的糙皮侧耳参考基因组(ASM1446616v1)信息[17]，通过BWA进行数据比对[18]，Picard去重复[19]。过滤掉变异质量值低于30、QD值低于2.0、MQ值低于40、FS值高于60的位点[20-21]。通过SnpEff完成SNP和InDel的注释[22]。对过滤掉多重和无差异后的位点通过欧式距离(Euclidean distance，ED)方法进行关联分析，然后采用DISTANCE方法对ED值进行拟合，统计并根据关联阈值判定候选区域[23]。对SNP和InDel 2个关联区域取交集缩小候选区间，应用BLAST进行NR、Swiss-Prot、GO、KEGG等多个数据库的基因注释[24]。

2 结果与分析

2.1 糙皮侧耳单核体群体抗性遗传分析

以抗螨性的糙皮侧耳Pos045为亲本，共挑取479个Pos045担孢子单核菌株建立分离群体，其对螨虫抗性呈现近似正态分布，表明糙皮侧耳抗螨性符合数量性状遗传(图1)。其中鉴定为高感的菌株48株，高抗54株，分别从中选取菌落大小、菌丝疏密等其他外观性状相似的菌株各20株，构建高感/高抗混池。

图1 糙皮侧耳Pos045单核体群体抗性分布Fig.1 Resistance distribution of Pleurotus ostreatus Pos045 monokaryotic population

2.2 测序数据量及数据质量

高感/高抗混池通过测序、过滤后，得到5 G高质量数据，平均Q30值达到94.70%，平均G+C碱基含量46.86%。read平均比对效率达到81.31%，总体测序深度超过46×(总碱基数与物种基因组大小的比值)，基因组覆盖率为97.20%，文库片段大小呈单峰正态分布(图2)。所得测序数据质量较好，可靠性较高。

图2 测序数据基因组覆盖深度分布及插入片段长度分布Fig.2 Genome coverage depth distribution and insert size distribution

2.3 变异位点检测

与参考基因组比对后，高抗池共检测到612 863个SNP，高感池检测到542 869个SNP，两个混池之间去掉相同位点后共获得105 323个SNP。根据检测SNP位点所在的基因位置信息，SNP位点分布在13个区域，其中位于同义突变、基因上游、基因下游、非同义突变、内含子的SNP位点最多，其余的则分布于可变剪切位点、终止子提前、终止子丢失、基因区间、3′非翻译区、5′非翻译区等区域中，其中13 325个SNP引起非同义突变。高抗池共检测到102 721个InDel，其中11 602个在基因编码区。高感池检测到91 711个InDel，有10 372个在基因编码区。两个混池之间去掉相同位点后共获得22 025个InDel。根据基因注释位置信息，InDel位点分布在基因上/下游、内含子、移码突变等15个功能区或类型，其中引起移码突变的InDel位点1 571个。

2.4 变异位点关联分析

过滤后共获得556 886个高质量SNP位点。ED值拟合关联值分布如图3所示。根据关联阈值(0.70)进行判定，得到2个可信候选区域，合计总长度为1.87 Mb，包含658个基因注释信息，其中有370个基因存在非同义突变位点。

分析得到94 456个高质量可信InDel位点。ED值拟合关联值分布如图3所示。根据关联阈值(0.69)判定得到2个候选区域，共1.75 Mb，含611个基因，其中89个基因具有移码突变位点。

a:NW_023503156.1;b:NW_023503157.1;c:NW_023503158.1;d:NW_023503159.1;e:NW_023503160.1;f:NW_023503161.1;g:NW_023503162.1;h:NW_023503163.1;i:NW_023503164.1;j:NW_023503165.1;k:NW_023503166.1;l:NW_023503167.1;m:NW_023503168.1。图3 SNP和InDel关联抗螨候选基因的区间Fig.3 Identification of anti-mite candidate gene intervals using SNP and InDel association

对SNP和InDel两种方法关联抗螨候选基因区域取交集，得到2个与抗螨候选基因相关的候选区域，都位于参考基因组(ASM1446616v1)中scaffold NW_023503160.1上。关联区域总长度为1.75 Mb，共包含基因605个，其中非同义突变基因353个，移码突变基因89个。

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

对候选区间内的605个基因进行多个数据库注释，NR注释到605个、NT注释到605个、trEMBL注释到602个、SwissProt注释到269个、GO注释到409个、KEGG注释到302个、COG注释到198个。KEGG和GO富集部分结果见(图4、图5)。在KEGG通路富集中，富集到ABC转运蛋白、MAPK信号通路、吲哚生物碱合成、吲哚二萜生物碱合成、萜类骨架合成等相关基因26个(表1)。这些基因都有可能与平菇对螨虫抗性相关，需要进一步验证。

图4 注释基因在KEGG通路中富集结果Fig.4 Enrichment results of annotated genes in KEGG pathway

图5 注释基因在GO通路中富集结果Fig.5 Enrichment results of annotated genes in GO pathway

3 讨论

BSA分析方法利用构建极端表型群体的策略来进行经济且快速的连锁标记筛选或QTL定位。随着近十几年来高通量测序技术的不断发展，测序成本逐步降低，加速了BSA-seq分析方法在不同物种中的应用[25-26]，包括酵母[27-28]、农作物[29-32]、园艺作物[33-34]、果树[35-36]、脊椎动物[37-39]和昆虫[40-41]等。研究中，我们针对Pos045的单核体群体进行抗螨性BSA分析，充分利用了担子菌孢子单核体易获取、易保存的优势，避免了双核体和多核体杂合基因型对表型数据的干扰，有助于提高定位准确度和精准度。无需亲本测序的分析方法，可以省却亲本基因纯合培育的时间，加快抗性机制研究进程。

表1 候选区间内相关基因信息

本研究对野生糙皮侧耳Pos045的单核体群体进行抗螨性鉴定，发现螨虫对群体内菌株取食造成的危害程度不一，呈正态分布，说明糙皮侧耳抗螨性是数量遗传，由多个基因控制。进一步，我们将糙皮侧耳的抗螨性基因定位到一条染色体(NW_023503160.1)上，共得到2个与抗螨性相关的相邻候选区域。已报道的植物抗虫候选基因在染色体上距离靠近并成簇状分布，如已定位的水稻抗褐飞虱基因中，22个抗性基因在染色体上成簇存在[42]。大量的植物抗虫性机理研究结果表明，害虫取食激活了多种信号途径，如植物激素(水杨酸和茉莉酸等)的调控、MAPK通路的激活、细胞质钙离子浓度升高、活性氧含量升高等，继而调节防御相关基因的活性，代谢出防御物质，从而对害虫产生直接抗性或间接抗性[43]。本研究团队前期在食用菌和螨虫互作过程研究中同样发现，腐食酪螨取食可以诱导食用菌产生大量萜烯类物质，这些挥发性倍半萜在螨虫识别宿主真菌方面发挥了重要作用[44]。而机械损伤和螨虫取食食用菌菌丝体，对其转录组和蛋白质组的影响存在显著差异。这些结果都喻示着，食用菌在螨虫长期取食胁迫下，可能产生了基于次级代谢产物抗螨的机制。本研究中，我们在初定位候选区间内经功能注释以及GO和KEGG通路富集分析，发现26个基因参与了信号传导、防御过程和次级代谢相关通路，推测这些候选基因可能与抗螨性相关，但具体的功能仍需进一步分析和验证。

4 结论

糙皮侧耳抗螨性表现为数量性状，基于BSA-Seq技术将候选基因定位在NW_023503160.1上总长度为1.75 Mb的相邻2个候选区域，共包含605个基因信息，其中功能分析结果表明有26个基因值得优先关注。这些候选基因为食用菌抗螨性基因挖掘和遗传改良奠定了一定的基础。