基于BSA-seq技術(shù)挖掘糙皮側(cè)耳抗螨候選基因

李輝平， 駱 昕， 侯子強(qiáng)， 蔣 寧， 林金盛， 侯立娟， 徐 平， 馬 林， 曲紹軒,

(1.江蘇省農(nóng)業(yè)科學(xué)院蔬菜研究所，江蘇省高效園藝作物遺傳改良重點(diǎn)實(shí)驗(yàn)室，江蘇南京210014；2.江蘇大學(xué)生命科學(xué)學(xué)院，江蘇鎮(zhèn)江212013)

糙皮側(cè)耳(Pleurotusostreatus)，又稱平菇、蠔菇，其栽培方式簡(jiǎn)單粗放、抗病蟲害能力強(qiáng)、產(chǎn)量高，栽培技術(shù)需求低，受到全球食用菌生產(chǎn)者的歡迎，是世界重要栽培食用菌之一[1]。又因其富含膳食纖維、蛋白質(zhì)、礦物質(zhì)、維生素和多種活性物質(zhì)，越來越受到消費(fèi)者的關(guān)注和歡迎[2-4]。2020年全國(guó)平菇產(chǎn)量達(dá)到6.829 6×106t，是中國(guó)第三大栽培食用菌[5]。

腐食酪螨(Tyrophagusputrescentiae)是世界性倉儲(chǔ)害蟲之一，取食各種干制食品和生物原材料，對(duì)30多種食用菌造成危害，是食用菌栽培過程中需要防控的主要螨蟲。該螨以刺吸式口器取食菌絲體和子實(shí)體，同時(shí)在培養(yǎng)基質(zhì)上傳播病原菌[6]。發(fā)菌期暴發(fā)螨蟲可造成大規(guī)模減產(chǎn)甚至絕收，出菇期暴發(fā)螨蟲則可減產(chǎn)30%～40%[7-8]。

目前，食用菌害螨主要依賴化學(xué)防治。但由于螨蟲個(gè)體微小、生長(zhǎng)周期短、變異快，已對(duì)多種藥劑產(chǎn)生抗藥性，越來越難防治[9]。培育和利用抗螨食用菌品種，可有效地解決田間螨害問題，降低防控成本，推進(jìn)食用菌綠色安全生產(chǎn)。本研究團(tuán)隊(duì)前期研究發(fā)現(xiàn)，腐食酪螨在中國(guó)不同栽培糙皮側(cè)耳品種上繁殖力、發(fā)育周期、產(chǎn)卵量、若螨死亡率等性狀存在顯著差異，表明國(guó)內(nèi)不同平菇菌株對(duì)腐食酪螨存在不同的抗性水平[10]。然而，關(guān)于糙皮側(cè)耳對(duì)螨蟲抗性的來源和遺傳機(jī)制尚不清楚，國(guó)內(nèi)外也尚無抗螨性食用菌品種的相關(guān)報(bào)道。

群體分離分析法(Bulked segregant analysis，BSA)通過表型差異構(gòu)建分離群體快速篩選和定位目標(biāo)性狀相關(guān)基因[11]。該方法在動(dòng)植物功能基因定位研究中應(yīng)用十分廣泛[12-13]，在食用菌功能基因定位研究中也逐步得到了應(yīng)用[14]。BSA-seq技術(shù)基于高通量測(cè)序技術(shù)，具有實(shí)驗(yàn)周期短、定位準(zhǔn)確，開發(fā)的分子標(biāo)記分布均勻、密度高等優(yōu)點(diǎn)[15-16]。

本研究以一株野生糙皮側(cè)耳菌株P(guān)os045為試驗(yàn)材料，構(gòu)建擔(dān)孢子分離群體并開展抗螨性鑒定，分析其對(duì)腐食酪螨產(chǎn)生抗性的遺傳基礎(chǔ)，通過BSA-seq技術(shù)對(duì)潛在抗性基因進(jìn)行初步定位，為開展糙皮側(cè)耳抗螨分子育種打下基礎(chǔ)。

1 材料與方法

1.1 試驗(yàn)材料

1.1.1 供試菌株 糙皮側(cè)耳Pos045由江蘇省農(nóng)業(yè)科學(xué)院提供。

1.1.2 供試螨蟲 腐食酪螨由江蘇省農(nóng)業(yè)科學(xué)院提供，室內(nèi)繼代飼養(yǎng)[10]。

1.1.3 培養(yǎng)基 PDA固體培養(yǎng)基，購(gòu)自臺(tái)灣MDBio公司，按說明書配制滅菌備用。棉籽殼出菇培養(yǎng)基，按質(zhì)量百分比棉籽殼55%，玉米芯30%，麩皮15%，混勻后加水?dāng)嚢柚梁?5%，用330 ml培養(yǎng)瓶裝至瓶肩，121 ℃高壓蒸汽滅菌90 min備用。平菇抗螨性鑒定培養(yǎng)基，棉籽殼200 g(煮水30 min過濾)，葡萄糖20 g，瓊脂15 g，加水定容到1 L后121 ℃蒸汽高壓滅菌30 min備用。

1.2 試驗(yàn)方法

1.2.1 孢子單核體群體構(gòu)建 糙皮側(cè)耳Pos045經(jīng)PDA培養(yǎng)基活化后接種至棉籽殼出菇培養(yǎng)基，25 ℃下避光培養(yǎng)至菌絲體長(zhǎng)滿培養(yǎng)瓶，降溫至18 ℃進(jìn)行催蕾出菇。待菌蓋開始舒展時(shí)，在超凈臺(tái)中懸于滅菌的錐形瓶中靜置24 h，用200 μl無菌水沖洗瓶底孢子即為孢子懸浮液，4 ℃短期保存待用。將孢子懸浮液稀釋至500 CFU/mL以下，取100 μl涂布90 mm平板， 25 ℃下避光培養(yǎng)5 d后，挑取單個(gè)菌落進(jìn)行繼代培養(yǎng)，顯微鏡檢查有無鎖狀聯(lián)合。挑取不少于400個(gè)單核菌株建立Pos045單核體群體。

1.2.2 抗螨性鑒定 以本實(shí)驗(yàn)室建立的平皿菌絲體接種螨蟲取食孔洞計(jì)數(shù)法進(jìn)行抗螨性鑒定。菌絲活化后接種至抗螨性鑒定培養(yǎng)基，菌絲長(zhǎng)滿平板后接種30 mg約500頭螨蟲，透氣膜封閉25 ℃下繼續(xù)培養(yǎng)，5個(gè)重復(fù)。接種螨蟲48 h后，借鑒玉米抗玉米螟鑒定技術(shù)規(guī)范(NY/T 1248.5-2006)，觀測(cè)螨蟲取食菌絲后形成的孔洞數(shù)量和直徑大小，統(tǒng)籌劃分為害級(jí)別1～9級(jí)。計(jì)算為害程度的平均值，為害程度的平均值=∑(為害級(jí)別×該級(jí)別菌株數(shù))/調(diào)查總數(shù)，并據(jù)此劃分抗性等級(jí)，共分5級(jí)：1.0～2.9為高抗；3.0～4.9為抗；5.0～6.9為中抗；7.0～8.9為感；9.0為高感。

1.2.3 極端混池構(gòu)建并測(cè)序 群體其他表型數(shù)據(jù)盡量一致，根據(jù)抗性數(shù)據(jù)，選擇高抗/高感極端株各20株構(gòu)建極端表型混池。單株分別使用MiniBEST Plant DNA純化試劑盒(寶生物公司產(chǎn)品)提取DNA，等量混合構(gòu)建DNA混池。利用Illumina HiSeq平臺(tái)(北京百邁克生物科技有限公司產(chǎn)品)完成建庫測(cè)序，單池測(cè)序數(shù)據(jù)2 G以上。

1.3 數(shù)據(jù)分析

基于NCBI公布的糙皮側(cè)耳參考基因組(ASM1446616v1)信息[17]，通過BWA進(jìn)行數(shù)據(jù)比對(duì)[18]，Picard去重復(fù)[19]。過濾掉變異質(zhì)量值低于30、QD值低于2.0、MQ值低于40、FS值高于60的位點(diǎn)[20-21]。通過SnpEff完成SNP和InDel的注釋[22]。對(duì)過濾掉多重和無差異后的位點(diǎn)通過歐式距離(Euclidean distance，ED)方法進(jìn)行關(guān)聯(lián)分析，然后采用DISTANCE方法對(duì)ED值進(jìn)行擬合，統(tǒng)計(jì)并根據(jù)關(guān)聯(lián)閾值判定候選區(qū)域[23]。對(duì)SNP和InDel 2個(gè)關(guān)聯(lián)區(qū)域取交集縮小候選區(qū)間，應(yīng)用BLAST進(jìn)行NR、Swiss-Prot、GO、KEGG等多個(gè)數(shù)據(jù)庫的基因注釋[24]。

2 結(jié)果與分析

2.1 糙皮側(cè)耳單核體群體抗性遺傳分析

以抗螨性的糙皮側(cè)耳Pos045為親本，共挑取479個(gè)Pos045擔(dān)孢子單核菌株建立分離群體，其對(duì)螨蟲抗性呈現(xiàn)近似正態(tài)分布，表明糙皮側(cè)耳抗螨性符合數(shù)量性狀遺傳(圖1)。其中鑒定為高感的菌株48株，高抗54株，分別從中選取菌落大小、菌絲疏密等其他外觀性狀相似的菌株各20株，構(gòu)建高感/高抗混池。

圖1 糙皮側(cè)耳Pos045單核體群體抗性分布Fig.1 Resistance distribution of Pleurotus ostreatus Pos045 monokaryotic population

2.2 測(cè)序數(shù)據(jù)量及數(shù)據(jù)質(zhì)量

高感/高抗混池通過測(cè)序、過濾后，得到5 G高質(zhì)量數(shù)據(jù)，平均Q30值達(dá)到94.70%，平均G+C堿基含量46.86%。read平均比對(duì)效率達(dá)到81.31%，總體測(cè)序深度超過46×(總堿基數(shù)與物種基因組大小的比值)，基因組覆蓋率為97.20%，文庫片段大小呈單峰正態(tài)分布(圖2)。所得測(cè)序數(shù)據(jù)質(zhì)量較好，可靠性較高。

圖2 測(cè)序數(shù)據(jù)基因組覆蓋深度分布及插入片段長(zhǎng)度分布Fig.2 Genome coverage depth distribution and insert size distribution

2.3 變異位點(diǎn)檢測(cè)

與參考基因組比對(duì)后，高抗池共檢測(cè)到612 863個(gè)SNP，高感池檢測(cè)到542 869個(gè)SNP，兩個(gè)混池之間去掉相同位點(diǎn)后共獲得105 323個(gè)SNP。根據(jù)檢測(cè)SNP位點(diǎn)所在的基因位置信息，SNP位點(diǎn)分布在13個(gè)區(qū)域，其中位于同義突變、基因上游、基因下游、非同義突變、內(nèi)含子的SNP位點(diǎn)最多，其余的則分布于可變剪切位點(diǎn)、終止子提前、終止子丟失、基因區(qū)間、3′非翻譯區(qū)、5′非翻譯區(qū)等區(qū)域中，其中13 325個(gè)SNP引起非同義突變。高抗池共檢測(cè)到102 721個(gè)InDel，其中11 602個(gè)在基因編碼區(qū)。高感池檢測(cè)到91 711個(gè)InDel，有10 372個(gè)在基因編碼區(qū)。兩個(gè)混池之間去掉相同位點(diǎn)后共獲得22 025個(gè)InDel。根據(jù)基因注釋位置信息，InDel位點(diǎn)分布在基因上/下游、內(nèi)含子、移碼突變等15個(gè)功能區(qū)或類型，其中引起移碼突變的InDel位點(diǎn)1 571個(gè)。

2.4 變異位點(diǎn)關(guān)聯(lián)分析

過濾后共獲得556 886個(gè)高質(zhì)量SNP位點(diǎn)。ED值擬合關(guān)聯(lián)值分布如圖3所示。根據(jù)關(guān)聯(lián)閾值(0.70)進(jìn)行判定，得到2個(gè)可信候選區(qū)域，合計(jì)總長(zhǎng)度為1.87 Mb，包含658個(gè)基因注釋信息，其中有370個(gè)基因存在非同義突變位點(diǎn)。

分析得到94 456個(gè)高質(zhì)量可信InDel位點(diǎn)。ED值擬合關(guān)聯(lián)值分布如圖3所示。根據(jù)關(guān)聯(lián)閾值(0.69)判定得到2個(gè)候選區(qū)域，共1.75 Mb，含611個(gè)基因，其中89個(gè)基因具有移碼突變位點(diǎn)。

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關(guān)聯(lián)抗螨候選基因的區(qū)間Fig.3 Identification of anti-mite candidate gene intervals using SNP and InDel association

對(duì)SNP和InDel兩種方法關(guān)聯(lián)抗螨候選基因區(qū)域取交集，得到2個(gè)與抗螨候選基因相關(guān)的候選區(qū)域，都位于參考基因組(ASM1446616v1)中scaffold NW_023503160.1上。關(guān)聯(lián)區(qū)域總長(zhǎng)度為1.75 Mb，共包含基因605個(gè)，其中非同義突變基因353個(gè)，移碼突變基因89個(gè)。

2.5 候選區(qū)域篩選與功能注釋

對(duì)候選區(qū)間內(nèi)的605個(gè)基因進(jìn)行多個(gè)數(shù)據(jù)庫注釋，NR注釋到605個(gè)、NT注釋到605個(gè)、trEMBL注釋到602個(gè)、SwissProt注釋到269個(gè)、GO注釋到409個(gè)、KEGG注釋到302個(gè)、COG注釋到198個(gè)。KEGG和GO富集部分結(jié)果見(圖4、圖5)。在KEGG通路富集中，富集到ABC轉(zhuǎn)運(yùn)蛋白、MAPK信號(hào)通路、吲哚生物堿合成、吲哚二萜生物堿合成、萜類骨架合成等相關(guān)基因26個(gè)(表1)。這些基因都有可能與平菇對(duì)螨蟲抗性相關(guān)，需要進(jìn)一步驗(yàn)證。

圖4 注釋基因在KEGG通路中富集結(jié)果Fig.4 Enrichment results of annotated genes in KEGG pathway

圖5 注釋基因在GO通路中富集結(jié)果Fig.5 Enrichment results of annotated genes in GO pathway

3 討論

BSA分析方法利用構(gòu)建極端表型群體的策略來進(jìn)行經(jīng)濟(jì)且快速的連鎖標(biāo)記篩選或QTL定位。隨著近十幾年來高通量測(cè)序技術(shù)的不斷發(fā)展，測(cè)序成本逐步降低，加速了BSA-seq分析方法在不同物種中的應(yīng)用[25-26]，包括酵母[27-28]、農(nóng)作物[29-32]、園藝作物[33-34]、果樹[35-36]、脊椎動(dòng)物[37-39]和昆蟲[40-41]等。研究中，我們針對(duì)Pos045的單核體群體進(jìn)行抗螨性BSA分析，充分利用了擔(dān)子菌孢子單核體易獲取、易保存的優(yōu)勢(shì)，避免了雙核體和多核體雜合基因型對(duì)表型數(shù)據(jù)的干擾，有助于提高定位準(zhǔn)確度和精準(zhǔn)度。無需親本測(cè)序的分析方法，可以省卻親本基因純合培育的時(shí)間，加快抗性機(jī)制研究進(jìn)程。

表1 候選區(qū)間內(nèi)相關(guān)基因信息

本研究對(duì)野生糙皮側(cè)耳Pos045的單核體群體進(jìn)行抗螨性鑒定，發(fā)現(xiàn)螨蟲對(duì)群體內(nèi)菌株取食造成的危害程度不一，呈正態(tài)分布，說明糙皮側(cè)耳抗螨性是數(shù)量遺傳，由多個(gè)基因控制。進(jìn)一步，我們將糙皮側(cè)耳的抗螨性基因定位到一條染色體(NW_023503160.1)上，共得到2個(gè)與抗螨性相關(guān)的相鄰候選區(qū)域。已報(bào)道的植物抗蟲候選基因在染色體上距離靠近并成簇狀分布，如已定位的水稻抗褐飛虱基因中，22個(gè)抗性基因在染色體上成簇存在[42]。大量的植物抗蟲性機(jī)理研究結(jié)果表明，害蟲取食激活了多種信號(hào)途徑，如植物激素(水楊酸和茉莉酸等)的調(diào)控、MAPK通路的激活、細(xì)胞質(zhì)鈣離子濃度升高、活性氧含量升高等，繼而調(diào)節(jié)防御相關(guān)基因的活性，代謝出防御物質(zhì)，從而對(duì)害蟲產(chǎn)生直接抗性或間接抗性[43]。本研究團(tuán)隊(duì)前期在食用菌和螨蟲互作過程研究中同樣發(fā)現(xiàn)，腐食酪螨取食可以誘導(dǎo)食用菌產(chǎn)生大量萜烯類物質(zhì)，這些揮發(fā)性倍半萜在螨蟲識(shí)別宿主真菌方面發(fā)揮了重要作用[44]。而機(jī)械損傷和螨蟲取食食用菌菌絲體，對(duì)其轉(zhuǎn)錄組和蛋白質(zhì)組的影響存在顯著差異。這些結(jié)果都喻示著，食用菌在螨蟲長(zhǎng)期取食脅迫下，可能產(chǎn)生了基于次級(jí)代謝產(chǎn)物抗螨的機(jī)制。本研究中，我們?cè)诔醵ㄎ缓蜻x區(qū)間內(nèi)經(jīng)功能注釋以及GO和KEGG通路富集分析，發(fā)現(xiàn)26個(gè)基因參與了信號(hào)傳導(dǎo)、防御過程和次級(jí)代謝相關(guān)通路，推測(cè)這些候選基因可能與抗螨性相關(guān)，但具體的功能仍需進(jìn)一步分析和驗(yàn)證。

4 結(jié)論

糙皮側(cè)耳抗螨性表現(xiàn)為數(shù)量性狀，基于BSA-Seq技術(shù)將候選基因定位在NW_023503160.1上總長(zhǎng)度為1.75 Mb的相鄰2個(gè)候選區(qū)域，共包含605個(gè)基因信息，其中功能分析結(jié)果表明有26個(gè)基因值得優(yōu)先關(guān)注。這些候選基因?yàn)槭秤镁跪曰蛲诰蚝瓦z傳改良奠定了一定的基礎(chǔ)。