利用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個候選基因, 這些基因主要參與植物激素生長調控及免疫系統的信號轉導．