基于BSA-seq技術挖掘糙皮側耳抗螨候選基因

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

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

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

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

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

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

本研究以一株野生糙皮側耳菌株Pos045為試驗材料，構建擔孢子分離群體并開展抗螨性鑒定，分析其對腐食酪螨產(chǎn)生抗性的遺傳基礎，通過BSA-seq技術對潛在抗性基因進行初步定位，為開展糙皮側耳抗螨分子育種打下基礎。

1 材料與方法

1.1 試驗材料

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

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

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

1.2 試驗方法

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

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

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

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

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

2 結果與分析

2.1 糙皮側耳單核體群體抗性遺傳分析

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

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

2.2 測序數(shù)據(jù)量及數(shù)據(jù)質量

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

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

2.3 變異位點檢測

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

2.4 變異位點關聯(lián)分析

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

分析得到94 456個高質量可信InDel位點。ED值擬合關聯(lián)值分布如圖3所示。根據(jù)關聯(lián)閾值(0.69)判定得到2個候選區(qū)域，共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關聯(lián)抗螨候選基因的區(qū)間Fig.3 Identification of anti-mite candidate gene intervals using SNP and InDel association

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

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

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

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

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

4 結論

糙皮側耳抗螨性表現(xiàn)為數(shù)量性狀，基于BSA-Seq技術將候選基因定位在NW_023503160.1上總長度為1.75 Mb的相鄰2個候選區(qū)域，共包含605個基因信息，其中功能分析結果表明有26個基因值得優(yōu)先關注。這些候選基因為食用菌抗螨性基因挖掘和遺傳改良奠定了一定的基礎。