蘇云金芽胞桿菌G03菌株的比較基因組分析

姚萌 王奎 束長龍 杜立新 李海濤 丁明月 劉榮梅

（1. 東北農(nóng)業(yè)大學生命科學學院，哈爾濱 150030；2. 中國農(nóng)業(yè)科學院植物保護研究所植物病蟲害生物學國家重點實驗室，北京 100193；3. 河北農(nóng)林科學院植物保護研究所，保定 071000）

作為目前研究最多、應用最廣的生防微生物，蘇云金芽胞桿菌（Bacillus thuringiensis，Bt）廣泛分布于土壤、水體、塵埃、病死昆蟲、植被等基質中，其最大的特點是在形成芽胞的同時生成伴胞晶體，這種伴胞晶體由不同的殺蟲晶體蛋白組成（Insecticidal crystal proteins，ICPs），包括 Cry和Cyt蛋白［1］，對多種昆蟲具有特異的高效殺滅活性，如鱗翅目、鞘翅目、雙翅目、半翅目，及線蟲、螨類［2-7］等，并且對非靶標天敵安全、對人畜無害、不污染環(huán)境，因此Bt廣泛用于防治農(nóng)、林、醫(yī)、倉貯等領域的害蟲［8］。

目前國際上Bt商品化產(chǎn)品約有數(shù)百種，其中用于防治鱗翅目害蟲的Bt產(chǎn)品主要源自庫斯塔克亞種（B. thuringiensissubsp.kurstaki，血清型H3a，3b，3c，簡稱Btk）菌株，如HD1菌株（Bt殺蟲劑Dipel的生產(chǎn)菌株）。該菌株已經(jīng)在全世界應用幾十年，并且是目前用于評估鱗翅目害蟲特異Bt殺蟲劑毒力效價的標準參考菌株［9］。在國家有關項目的大力支持下，我國研究機構也分離出一系列對鱗翅目高毒力的Bt菌株。例如，河北農(nóng)林科學院植物保護研究所分離到的鲇澤亞種（B. thuringiensissubsp.aizawai，血清型7，簡稱Bta）菌株G03，同樣對多種鱗翅目害蟲有效，中國農(nóng)業(yè)科學院植物保護研究所在此基礎上構建出了同時對鱗翅目和鞘翅目（葉甲科）害蟲均具有高毒性的工程菌株G033A［10-11］，并且目前已經(jīng)獲得了農(nóng)藥登記證（PD20171726），商標為“禁衛(wèi)軍”。

為了進一步分析G03等鱗翅目害蟲高效菌株的特點，本研究對G03菌株進行了基因組測序，并與已公布序列的Bt生產(chǎn)菌株HD1以及Btk模式菌株HD73進行比較基因組的研究，相關研究結果對G03應用以及今后Bt菌株遺傳改良均具有重要的指導意義。此外，為了分析不同殺蟲活性特性Bt菌株之間的進化差異，本研究選擇了兩株蠐螬有效的菌株（Bt185、HBF18），以及對鱗翅目害蟲有效但是殺蟲活性相對較差的菌株HD12進行了基因組進化分析。

1 材料與方法

1.1 供試菌株

G03菌株是實驗室保存菌株。從NCBI數(shù)據(jù)庫（https：//www.ncbi.nlm.nih.gov/）下載菌株 HD1、HD73，Bt185、HBF18以及HD12基因組數(shù)據(jù)用于后續(xù)分析。

1.2 方法

1.2.1 基因組制備與測序 按照張彥蕊等方法提?。?2］G03菌株的基因組DNA并檢測基因組DNA的質量，合格后送交公司測序。利用TruSeq DNA無PCR文庫制備試劑盒構建測序文庫，并交由華大基因公司利用Illumina HiSeq 2500 測序平臺測序。

1.2.2 測序數(shù)據(jù)的組裝、注釋 對測序獲得的原始Reads進行質量控制，刪除低質量的Reads，利用IDBA_UD軟件包［13］對剩余高質量的pair-end clean reads數(shù)據(jù)進行基因組序列組裝。并利用NCBI 原核生物基因組注釋流程PGAP（Prokaryotic Genome Annotation Pipeline，https：//www.ncbi.nlm.nih.gov/genome/annotation_prok/）對組裝基因組數(shù)據(jù)進行基因預測與注釋。

利用泛基因組分析PGAP軟件包［14］進行基因家族聚類。利用維恩圖在線繪制網(wǎng)站（http：//bioinfogp.cnb.csic.es/tools/venny/）進行菌株間共享基因及特有基因的統(tǒng)計分析。利用Blast軟件包［15］、TBtools軟件［16］進行Go功能富集。并利用Blast軟件包［15］進行殺蟲基因及功能基因的注釋。利用MEGAX軟件包［17］和CVTree軟件包［18］分別進行蛋白質系統(tǒng)發(fā)生分析和全基因組系統(tǒng)發(fā)生分析。

2 結果

2.1 基因組測序、組裝與注釋

利用IDBA_UD軟件包對G03菌株進行基因組組裝。結果顯示，基因組大小為6.38 M，組裝contig數(shù)量為395，GC含量為34.71%，Contig N50為45.68 kb，Contig L50為41。相關組裝結果已經(jīng)提交到NCBI GenBank數(shù)據(jù)庫，基因組登錄號為NHNQ00000000.1。進一步利用NCBI 原核生物基因組注釋流程PGAP對基因組進行了注釋，結果顯示G03菌株包含6916個基因，其中蛋白質編碼基因（Coding sequence，CDS）6 848 個。

2.2 CVTree系統(tǒng)進化分析

CVTree（Composition vector tree，簡稱 CVTree）是郝柏林研究組于2003年提出的一種全新的用于推測物種之間親緣關系與分類的研究方法［18］。這種方法基于全基因組，不需要進行序列比對，根據(jù)全基因組序列中不同核苷酸堿基或氨基酸殘基短串的出現(xiàn)頻率來構建系統(tǒng)發(fā)育樹。迄今為止，CVTree已經(jīng)廣泛用于分析多種不同物種及數(shù)據(jù)，包括古生菌［19-20］、原核生物［18，21］、真菌［22］、病毒［23］、葉綠體序列［24］、tRNA 序列［25］、癌癥序列［26］等。

我們利用CVTree對包括G03、HD1、HD73菌株在內的7株Bt菌株的全基因組構建組分矢量并作出Neighbor-joining系統(tǒng)發(fā)生樹（氨基酸短串K值取6）。結果（圖1）發(fā)現(xiàn)，Btk菌株HD1和HD73聚在同一個分支，菌株G03與Btk菌株HD1和HD73所在分支相聚較近，屬于同一大分支。而同樣對鱗翅目害蟲有效的菌株HD12則相距較遠。此外，對鞘翅目害蟲有效的菌株Bt185與HBF18與G03所在分支也有較大的距離。

圖1 菌株CVTree系統(tǒng)進化分析樹圖

2.3 比較基因組分析

我們對CVTree分析中親緣關系較近的Bta菌株G03以及Btk菌株HD1、HD73這3株Bt菌株進行了比較基因組的分析，首先利用PGAP軟件包對其進行基因家族聚類，并利用維恩圖在線繪制網(wǎng)站（http：//bioinfogp.cnb.csic.es/tools/venny/）對這 3株菌株的共享基因及特有基因家族進行繪圖統(tǒng)計分析，如圖2所示。PGAP分析共鑒定出了6 860個基因家族，其中核心基因家族共有 4 720個（全部菌株均含有的基因家族），占全部基因家族的68.8%。而G03、HD1和HD73菌株分別有540、433和327個特有基因（僅在一株菌株中存在的基因）。G03與HD1有5122基因家族是共有的，這大于HD1與HD73共有基因家族（5101）。

2.4 菌株功能分類比較分析

我們利用Blast程序包分別將G03和HD1菌株的全部基因、共有基因以及特有基因序列與uniprot數(shù)據(jù)庫（https：//www.uniprot.org/）進行序列比對，并利用TBtools軟件進行Go功能注釋以及l(fā)evel 2水平上的功能解析。圖3統(tǒng)計了菌株G03基因組Go功能注釋在level 2水平上的統(tǒng)計情況，紅色是菌株G03基因組特有基因的統(tǒng)計結果，藍色部分是G03和HD1菌株都有的基因的統(tǒng)計結果。

圖2 3株Bt菌株共享及特有基因的維恩圖分析

結果分析發(fā)現(xiàn)生物學途徑（biological process）相關功能部分，G03特有基因在cellular process和metabolic process兩個方面數(shù)量最多；在細胞學組件（cellular component）相關功能部分，G03特有基因在cell和cell part兩個方面數(shù)量最多；在分子功能（molecular function）方面，G03特有基因在binding和catalytic activity兩個方面數(shù)量最多。其中生物學途徑與分子功能相關的Go富集反映了G03在生物學水平與分子水平上的特點。

cellular process是指在細胞水平上發(fā)生的過程相關功能富集。例如，細胞對環(huán)境的應答過程，在該方面的功能有較多Go富集顯示菌株G03在環(huán)境適應等細胞過程方面有較多的特點。metabolic process是指生物體轉化化學物質的化學反應和途徑相關功能富集，包括合成代謝和分解代謝，這方面較多Go富集說明G03菌株在物質分解利用或化合物合成方面有較特別的能力。

Binding和catalytic activity分別指結合與催化方面的功能富集，這兩個功能是緊密相關的，G03特有基因在這兩個方面有較多的富集，顯示其在分子識別催化方面有特別的功能。

2.5 G03菌株殺蟲基因分析

我們對Bta及Btk菌株包含的殺蟲基因種類以及數(shù)量做了比較分析，HD73菌株僅有cry1Ac一個殺蟲基因，菌株G03與HD1均含有cry1Aa、cry1Ac、cry1Ia、cry2Ab、vip3Aa這5個高毒力殺蟲基因［3，27-29］。已有研究表明 Cry1A、Cry1I、Cry2A、Vip3A這幾類殺蟲蛋白具有不同的受體，并且殺蟲蛋白之間還有協(xié)同增效的作用，這可能是G03與HD1具有高活性的原因。此外，與HD1相比，G03菌株還含有cry1Ca、cry1Da、cry9Ea基因。

圖3 G03菌株Go功能注釋

2.6 G03、HD1菌株特有基因分析

利用基因家族分析的方法，對菌株G03與 HD1進行比較分析，結果顯示G03有597個特有基因家族，而HD1有814個特有基因家族。其中，G03菌株特有的基因家族中有232個可以注釋到功能，而HD1菌株特有的基因家族中有404個可以注釋到功能。進一步，對這些注釋到功能的基因家族進行分析，發(fā)現(xiàn)噬菌體與轉座酶相關的基因家族最為豐富。HD1菌株特有的基因家族中有34個基因家族與噬菌體相關，而G03菌株特有的基因家族中只有 10個與噬菌體有關。

此外，特有基因家族中，HD1菌株轉座酶基因也比較多，有26個基因家族屬于轉座酶，而G03菌株特有的基因家族中只有8個。轉座酶是轉座子執(zhí)行轉座功能的酶，通常由轉座子編碼。本研究進一步比較了兩個菌株含有轉座酶的數(shù)量與多樣性。結果顯示，G03菌株僅僅含有33個轉座酶或片段，而HD1多達216個轉座酶或片段。這些轉座酶按照序列一致性（按照序列一致性50%設定閾值）分為58個家族，其中16個基因家族是兩個菌株共有的。這些基因家族中，最大的基因家族含有45個轉座酶基因，其中有4個來源于菌株G03，41個來源于菌株HD1。該家族轉座酶含有TNP1 DDE 結構域，是 IS4、IS421、IS5377、IS427、IS402、IS1355 以及IS5等多種插入序列（Insert sequence，IS）的轉座酶。進一步用UPGMA對這些含有TNP1 DDE 結構域的轉座酶進行系統(tǒng)進化分析（去除了4個轉座酶片段），結果（圖 4）顯示這些轉座酶分為3個類群，其中類群1是兩個菌株共有的；類群2 有6個基因，是HD1特有的；類群3 只有1個基因，是G03特有的。

噬菌體與插入序列、轉座子是微生物獲得外源基因的重要途徑，同時也是微生物基因組變異的重要因素。上述分析結果表明G03含有較少的噬菌體與插入序列、轉座子相關的基因，顯示菌株G03與HD1相比，應該具有較低的變異概率。

3 討論

目前市場上常用的Bt生產(chǎn)菌株種類相對較少，主要以Btk菌株HD1為主，相關研究報道較多。Bta菌株G03是我國科學家自主分離并且已經(jīng)用于害蟲防控的菌株，已通過基因工程改造拓展了其殺蟲譜，在鱗翅目和鞘翅目（葉甲科）害蟲防治中發(fā)揮重要作用，然而該菌株遺傳學、基因組相關研究鮮見報道。本研究利用Illumina測序技術進行了菌株G03基因組測序，并進行了基因功能注釋。在此基礎上與NCBI數(shù)據(jù)庫中已經(jīng)公開的Btk生產(chǎn)菌株HD1、Btk模式菌株HD73、兩株蠐螬有效的菌株（Bt185、HBF18）以及含有多種殺蟲基因類型的菌株HD12進行了比較基因組研究，旨在從基因組角度分析G03等鱗翅目害蟲高效菌株的特點、比較與Btk菌株的異同，為進一步對Bt菌株進行改良、推進我國自主分離研發(fā)的Bt菌株及其殺蟲劑產(chǎn)品的應用奠定基礎。

圖4 含有TNP1 DDE結構域轉座酶進化關系分析。箭頭標出的是G03轉座酶

本研究首次利用基因組測序技術，從基因組水平上全面的分析了G03殺蟲相關的基因。與PCR殺蟲基因鑒定技術相比，基因組測序技術避免了PCR鑒定技術由于引物設計不全面、PCR擴增條件不合適等原因，導致基因鑒定的誤差。通過基因組分析發(fā)現(xiàn)，菌株G03含有8個殺蟲基因（cry1Aa、cry1Ac、cry1Ca、cry1Da、cry1Ia、cry2Ab、cry9Ea、vip3Aa），其中cry1Ca、cry1Da是菌株HD1不具備的。目前，表達Cry1A、Cry2A和Vip3A等蛋白的轉基因的抗蟲作物已經(jīng)在多個國家商業(yè)化種植（International Service for the Acquisition of Agri-biotech Applications，www.isaaa.org/），使用表達與上述蛋白無交互抗性的Bt殺蟲劑在害蟲抗性治理方面具有重要的意義。利用小菜蛾抗Cry1Ac種群研究表明，Cry1Ca與Cry1Ac蛋白沒有交互抗性［30］，因此，菌株G03比HD1在害蟲的Bt殺蟲劑抗性風險控制方面更具優(yōu)勢。

Bt菌株中與殺蟲活性相關的免疫抑制因子A（Immune Inhibitor A，InhA）、 雙 效 菌 素 ZwA（Zwittermicin A，ZwA）、?；呓z氨酸內酯水解酶（AHLs Hydrolase，aiiA）以及幾丁質酶（Chitinase）對Bt菌株殺蟲活性及環(huán)境適應性都有重要意義。InhA是Bt分泌的一種金屬蛋白酶，能夠降解昆蟲產(chǎn)生的抗菌肽從而逃避宿主的免疫系統(tǒng)［31-32］。aiiA可以水解細菌的群體感應相關的信號分子?；呓z氨酸內酯（Acylated Homoserine Lactones，AHLs），對多種細菌有抑制作用，因此有助于提高Bt在昆蟲腸道內競爭優(yōu)勢［33-34］。Chitinase是一種可溶性的胞外蛋白質類殺蟲活性物質，可以幫助Bt菌株降解昆蟲腸道圍食膜中的幾丁質成分，使其能通過穿孔的腸道進入血腔引起昆蟲敗血癥，對Bt殺蟲蛋白具有增效作用［35-36］。前期研究發(fā)現(xiàn)，菌株G03與HD1都含inhA、aiiA、chitinase基因以及合成ZwA的基因簇，并且aiiA、chitinase基因進行系統(tǒng)發(fā)生分析發(fā)現(xiàn)鱗翅目害蟲高效的菌株aiiA、chitinase基因序列比較接近，說明G03菌株在感染寄主、適應環(huán)境方面與HD1類似［37］。

此外，本研究還發(fā)現(xiàn)與HD1相比，G03含有較少的噬菌體與轉座子基因，說明G03基因組應該具有更好的遺傳穩(wěn)定性。作為生產(chǎn)菌株，遺傳穩(wěn)定性對產(chǎn)品發(fā)酵生產(chǎn)尤為重要，基因組比較分析結果顯示G03菌株在該方面更有優(yōu)勢。

綜上所述，本研究對生產(chǎn)菌株G03進行的基因組測序、注釋與比較基因組分析，不僅可以從殺蟲基因角度揭示生產(chǎn)菌株高活性的原因，還可以為進一步菌株改良提供基因組數(shù)據(jù)支撐。