聶俊輝，王 平，張立新，王 通，郭建軍，曾 靜，袁 林

（1.江西省科學(xué)院微生物研究所，江西 南昌 330096；2.上饒市農(nóng)林水科學(xué)研究中心，江西 上饒 334000）

多粘類芽孢桿菌（Paenibacillus polymyxa）屬芽孢桿菌科（Bacillaceae）類芽孢桿菌屬（Paenibacillus），是一種產(chǎn)芽孢的革蘭氏陽性菌[1]。多粘類芽孢桿菌在自然界中分布廣泛，在玉米、高粱、小麥等多種植物體內(nèi)外以及根際土壤中均有被發(fā)現(xiàn)，是一類重要的根際有益菌[2]。多粘類芽孢桿菌對植物生長發(fā)育的有益作用主要表現(xiàn)在促進(jìn)植物產(chǎn)生激素、促進(jìn)植物吸收養(yǎng)分、促進(jìn)植物光合作用、增強(qiáng)植物固氮作用等方面[3]；同時(shí)，多粘類芽孢桿菌還能通過分泌抗菌殺蟲物質(zhì)如抗菌肽、細(xì)菌素、脂肽等增強(qiáng)植物抵御病蟲害的能力，是一種不可多得的生物防治原料，目前已廣泛應(yīng)用于農(nóng)業(yè)生產(chǎn)領(lǐng)域，具有較高的經(jīng)濟(jì)價(jià)值[4-5]。

近年來，高通量測序技術(shù)快速發(fā)展，使得細(xì)菌基因組的高效解析成為現(xiàn)實(shí)，為細(xì)菌的功能及其作用機(jī)制研究提供了重要手段[6-7]。趙興麗等[8]通過高通量測序技術(shù)發(fā)現(xiàn)土壤中有益微生物的含量在一定程度上可間接反映茶樹的抗病性，有益微生物菌群所占比率與發(fā)病程度呈負(fù)相關(guān)。陳德純等[9]對一株牦牛源產(chǎn)細(xì)菌素植物乳桿菌的基因組進(jìn)行測序，發(fā)現(xiàn)其基因組中包括多個(gè)細(xì)菌素基因以及細(xì)菌素轉(zhuǎn)運(yùn)系統(tǒng)的功能序列，為該菌株的應(yīng)用研究提供了重要的參考數(shù)據(jù)。高通量測序技術(shù)的普及有效地促進(jìn)了微生物的研究。

目前，有關(guān)多粘類芽孢桿菌的研究多集中在其生物特性及其對植物與生物防治的作用機(jī)理上，而對其基因序列的分析研究報(bào)道較少[10-12]。筆者從土壤中分離出一株多粘類芽孢桿菌LY1，利用基因組測序技術(shù)全面解析LY1 菌株的基因組序列，再通過生物信息學(xué)的方法比較不同菌株基因間的差異性和相似性，探索其基因結(jié)構(gòu)特征，明確其基因功能，從全基因組水平了解其系統(tǒng)發(fā)育進(jìn)化關(guān)系，以期為全面了解多粘類芽孢桿菌的遺傳背景，實(shí)現(xiàn)多粘類芽孢桿菌的遺傳改造提供基本線索。

1 材料與方法

1.1 供試菌株

多粘類芽孢桿菌LY1 從土壤中分離，由江西省科學(xué)院微生物研究所分離保存。

1.2 試驗(yàn)方法

1.2.1 測序及數(shù)據(jù)質(zhì)控取菌體委托上海生工生物工程股份有限公司使用Illumina Hiseq platform 進(jìn)行基因組測序。Illumina Hiseq?得到的原始圖像數(shù)據(jù)文件經(jīng)CASAVA 堿基識別 （Base Calling）分析轉(zhuǎn)化為原始測序序列（Sequenced Reads），對原始數(shù)據(jù)質(zhì)量值等信息進(jìn)行統(tǒng)計(jì)，并使用FastQC 對樣本的測序數(shù)據(jù)質(zhì)量進(jìn)行可視化評估。測序得到的原始數(shù)據(jù)，里面含有帶接頭的、低質(zhì)量的序列。為了保證信息分析質(zhì)量，必須對原始數(shù)據(jù)進(jìn)行過濾，得到Clean 數(shù)據(jù)。隨機(jī)從Clean 數(shù)據(jù)中抽取10 000 條序列與NCBI NT 數(shù)據(jù)庫進(jìn)行blastn 比對，取evalue ≤1E-10 并且相似度＞90%、coverage ＞80%的比對結(jié)果，計(jì)算其物種分布，同時(shí)進(jìn)行污染檢測。

1.2.2 數(shù)據(jù)拼裝使用SPAdes 對二代測序數(shù)據(jù)進(jìn)行拼裝。SPAdes 首先會對原始序列進(jìn)行序列錯(cuò)誤校正，然后通過多Kmer 值進(jìn)行組裝，最終綜合各Kmer 值組裝結(jié)果獲得最佳結(jié)果。再采用GapFiller 對拼接得到的contig 進(jìn)行Gap 修補(bǔ)，最后采用PrInSeS-G 進(jìn)行序列矯正，修正拼接過程中的剪輯錯(cuò)誤以及小片段的插入缺失。

1.2.3 基因預(yù)測與注釋采用Prokka 對組裝結(jié)果進(jìn)行基因元件預(yù)測。Prokka 是一系列基因元件預(yù)測工具的集合，調(diào)用Prodigal 預(yù)測編碼基因，Aragorn 預(yù)測tRNA，RNAmmer 預(yù)測rRNA，Infernal 預(yù)測miscRNA，預(yù)測出的各類基因元件匯總并完成初步注釋。采用Repeat Modeler 對組裝結(jié)果進(jìn)行重復(fù)序列Denovo 預(yù)測，再利用RepeatMasker 尋找基因組區(qū)段上各類型重復(fù)序列出現(xiàn)的位置和頻率。NCBI Blast+用于CDD、KOG、COG、NR、NT、PFAM、Swissprot、TrEMBL注釋?；赟wissprot 和TrEMBL 的蛋白注釋結(jié)果根據(jù)Uniprot 的注釋信息得到GO 注釋。使用KAAS，KEGG Automatic Annotation Server 進(jìn)行KEGG 注釋。

1.2.4 系統(tǒng)發(fā)育樹構(gòu)建將基因預(yù)測得到的16S rRNA序列與NCBI 的16S 數(shù)據(jù)庫進(jìn)行Blast 比對，設(shè)置參數(shù)identify ＞95。然后選取identify 最高的前30 條16S rRNA序列，利用muscle軟件進(jìn)行序列多重比對后，采用FastTree 軟件構(gòu)建系統(tǒng)發(fā)育樹。

2 結(jié)果與分析

2.1 測序數(shù)據(jù)質(zhì)控

二代測序和三代測序原始數(shù)據(jù)經(jīng)過去除帶接頭的、低質(zhì)量的序列，過濾后得到Clean 數(shù)據(jù)。Illumina Hiseq 原始測序數(shù)據(jù)過濾處理后，總共得到6 495 188個(gè)Clean Reads，總堿基數(shù)為951 905 276 bp，平均Read長度為146.56 bp；其中，Q10、Q20、Q30 的比例分別為100.00%（堿基數(shù)951 903 300 bp）、98.18%（堿基數(shù)934 538 164 bp）、93.79%（堿基數(shù)892 808 488 bp）；GC含量為45.63%（堿基數(shù)434 383 207 bp）。PacBio RSII原始測序數(shù)據(jù)過濾處理后，總共得到333 829 個(gè)Clean Reads，總堿基計(jì)數(shù)為925 726 342 bp，平均Reads 長度2 773.06 bp；其中，Reads ≥1 000 bp 、≥5 000 bp和≥10 000 bp 的分別占比47.34%、11.45%和5.30%，而Bases in Reads ≥1 000 bp、≥5 000 bp 和≥10 000 bp 的分別占比90.25%、61.52%和46.25%；GC 含量為44.40%。綜上所述，2 組數(shù)據(jù)質(zhì)量較高，可用于下一步基因組組裝。

2.2 基因預(yù)測與基因組基本特征

2.2.1 基因組的組成與信息2 組測序數(shù)據(jù)經(jīng)過序列矯正、Gap 修補(bǔ)、剪輯錯(cuò)誤修復(fù)、拼裝等程序后獲得菌株的基因組信息，如圖1 所示，多粘類芽孢桿菌LY1 的基因組為一條環(huán)狀閉合DNA，大小 5 765 474 bp，共預(yù)測到6 086 個(gè)蛋白質(zhì)編碼基因、162 個(gè) tRNA基因、42 個(gè)rRNA 基因、74 個(gè)small RNA ；最短基因長度為70 bp，最長基因長度為42 252 bp，基因平均長度為889.05 bp，長度≥500 bp 的基因有4 107 個(gè)，長度≥1 000 bp 的基因有1 916 個(gè)；重復(fù)區(qū)域計(jì)數(shù)為216 個(gè)，重復(fù)比例為1.45%；Low_complexity 為25 個(gè)，Simple_repeat 為117 個(gè)?；蜷L度和GC 含量的分布如圖2 所示。

圖1 多粘類芽孢桿菌LY1 的基因組圖

圖2 多粘類芽孢桿菌LY1 的基因長度（A）和GC 含量分布（B）

2.2.2 基因組注釋將預(yù)測的編碼基因與多個(gè)數(shù)據(jù)庫進(jìn)行比對，注釋到CDD（Conserved Domain Database）、KOG（euKaryotic Ortholog Groups）、NR（NCBI nonredundant protein sequences）、PFAM（Protein family）、Swissprot（A manually annotated and reviewed protein sequence database）、TrEMBL、GO（Gene Ontology）、KEGG（Kyoto Encyclopedia of Genes and Genomes）的基因比例分別為72.92%、64.27%、98.14%、69.57%、66.32%、97.52%、66.66%和37.92%，unigenes 總共5 813 個(gè)，同時(shí)注釋到NR、KEGG、Swissprot、COG的基因有2 133 個(gè)，注釋到所有數(shù)據(jù)庫的基因共有2 066 個(gè)，占unigenes 的35.54%。

2.2.3 NR 數(shù)據(jù)庫注釋由上面的數(shù)據(jù)可知，注釋到NR 數(shù)據(jù)庫的基因最為豐富，共有5 705 個(gè)，占比98.14%。NR 數(shù)據(jù)庫是一個(gè)非冗余的蛋白質(zhì)數(shù)據(jù)庫，內(nèi)容全面，通過與NR 數(shù)據(jù)庫的對比，可以查看物種轉(zhuǎn)錄本序列與相近物種的近似情況以及同源序列的功能信息[13]。由圖3 可知，鑒定到LY1 菌株有3 792 個(gè)基因與多粘類芽孢桿菌同源，與Paenibacillus polymyxaM1 菌株的相似性相對較高，在M1 的編碼基因中注釋到318 個(gè)基因，而在Paenibacillus polymyxaSQR-21、Paenibacillus polymyxaSC2 編 碼 基因中則分別注釋到76 和29 個(gè)基因。

圖3 多粘類芽孢桿菌LY1 基因的NR 數(shù)據(jù)庫比對

2.3 基因功能注釋

2.3.1 GO 功能注釋GO （Gene Ontology） 是一個(gè)國際標(biāo)準(zhǔn)化的基因功能分類體系，提供了一套動(dòng)態(tài)更新的標(biāo)準(zhǔn)詞匯表（controlled vocabulary）來全面描述生物體中基因和基因產(chǎn)物的屬性[14]。GO 總共有3 個(gè)ontology，分別描述基因的分子功能（MF，molecular function）、細(xì)胞組分（CC，cellular component）、參與的生物過程（BP，biological process）[15]。通過將基因進(jìn)行GO 注釋和分類，可判斷基因參與的主要功能。多粘類芽孢桿菌LY1 基因的3 大分類統(tǒng)計(jì)結(jié)果如圖4 所示。在生物過程中，分類到代謝過程（metabolic process）和細(xì)胞過程（cellular process）中的基因占比最多，分別為2 370（40.77%）和2 266 個(gè)（38.98%）；在細(xì)胞組分分類中，分類到細(xì)胞（cell）和細(xì)胞組成（cell part）中的基因占比最多，均為2 045 個(gè)（35.18%）；在分子功能分類中，分類到催化活性（catalytic activity）和結(jié)合（binding）中的基因占比最多，分別為2 345（40.34%）和2 080 個(gè)（35.78%）。

圖4 多粘類芽孢桿菌LY1 的GO 功能注釋

2.3.2 KEGG 注釋KEGG（Kyoto Encyclopedia of Genes and Genomes，京都基因與基因組百科全書）是基因組破譯方面的數(shù)據(jù)庫，在給出染色體中一套完整基因的情況下，它可以對蛋白質(zhì)交互網(wǎng)絡(luò)在各種細(xì)胞活動(dòng)起的作用作出預(yù)測，方便地尋找與行使某一類功能相關(guān)的所有注釋上的基因[16-17]。由圖5 可知，LY1 菌株的基因通過KEGG 分類主要分為5 大類，分別為細(xì)胞過程（Cellular Processes）、環(huán)境信息處理（Environmental Information Processing）、遺傳信息處理（Genetic Information Processing）、代謝（Metabolism）和有機(jī)系統(tǒng)（Organismal Systems），其中注釋到代謝的相關(guān)基因最多，達(dá)2 255 個(gè)，注釋到細(xì)胞過程、環(huán)境信息處理、遺傳信息處理和有機(jī)系統(tǒng)的基因數(shù)分別為99、613、334 和31 個(gè)。

圖5 多粘類芽孢桿菌LY1 的KEGG 功能注釋

2.3.3 COG 注釋蛋白質(zhì)直系同源簇?cái)?shù)據(jù)庫（COG，Cluster of Orthologous Groups of Prot-eins）是用于同源蛋白注釋的數(shù)據(jù)庫，是將細(xì)菌、藻類和真核生物等完整基因組的編碼蛋白根據(jù)系統(tǒng)進(jìn)化關(guān)系構(gòu)建而成，蛋白編碼序列與之比對后可以預(yù)測蛋白的功能[18-19]。多粘類芽孢桿菌LY1 的COG 功能注釋結(jié)果如圖6 所示，注釋到碳水化合物轉(zhuǎn)運(yùn)及代謝（Carbohydrate transport and metabolism）和轉(zhuǎn)錄（Transcription）分類的基因最多，分別為433（11.59%）和398 個(gè)（10.65%）。

圖6 多粘類芽孢桿菌LY1 的COG 功能注釋

2.4 毒力因子分析

把菌株LY1 的基因蛋白序列在VFDB 致病因子數(shù)據(jù)庫（Virulence Factors of Pathogenic Bacteria）中進(jìn)行比對，將基因與其對應(yīng)的毒力因子（Virulence Factors，VF）功能注釋信息相結(jié)合，發(fā)現(xiàn)注釋到SetB 組（full dataset，共30 053 個(gè)與毒力因子相關(guān)的基因）的蛋白序列共326 個(gè)，占比5.61%；注釋到SetA 組（core dataset，共2 585 個(gè)與毒力因子相關(guān)的基因）的蛋白序列共301 個(gè)，占比5.18%（表1）。

表1 多粘類芽孢桿菌LY1 的基因組基本特征

2.5 系統(tǒng)發(fā)育進(jìn)化分析

16S rRNA 基因是細(xì)菌上編碼rRNA 相對應(yīng)的DNA 序列，存在于所有細(xì)菌的基因組中，具有高度的保守性和特異性，是病原菌檢測和鑒定的一種有效參照[20-21]。如圖7 所示，LY1 菌株與Paenibacillus polymyxaM1、Paenibacillus polymyxaOSY-DF、Paenibacillus polymyxaSC2、Paenibacillus peoriae、Paenibacillus jamilae、Paenibacillus polymyxaE681 等菌株的距離較近。

圖7 菌株LY1 基于16S rRNA 的系統(tǒng)發(fā)育樹

3 小 結(jié)

目前，人們對多粘類芽孢桿菌的功能有了較為全面地了解，但其分子作用機(jī)制以及遺傳工程改造等方面的研究報(bào)道還不多[3]。該研究通過全基因組測序分析了多粘類芽孢桿菌LY1 的基因組序列，獲得其基因組基本特征、基因功能注釋及分類、系統(tǒng)發(fā)育進(jìn)化關(guān)系等關(guān)鍵信息，為其后續(xù)的功能挖掘、遺傳操作系統(tǒng)建立、基因工程改造等提供了研究基礎(chǔ)與前提條件。結(jié)果表明，多粘類芽孢桿菌LY1 的基因組長度為5 765 474 bp，共有6 086 個(gè)蛋白質(zhì)編碼基因、162 個(gè) tRNA 基因、42 個(gè)rRNA 基因，GC 含量約為45.23%；在GO 功能注釋結(jié)果中，在生物過程分類下，代謝過程和細(xì)胞過程中的基因占比最多，分別占比40.77%和38.98%；在細(xì)胞組分分類下，分類到細(xì)胞和細(xì)胞組成中的基因占比最多，為35.18%；在分子功能分類下，分類到催化活性和結(jié)合中的基因占比最多，分別為40.34%和35.78%；系統(tǒng)發(fā)育進(jìn)化分析顯示，菌株LY1 與Paenibacillus polymyxaM1、Paenibacillus polymyxaOSY-DF、Paenibacillus polymyxaSC2、Paenibacillus polymyxaE681 等菌株進(jìn)化距離較近。同時(shí)，根據(jù)多粘類芽孢桿菌LY1 的基因組信息及發(fā)育進(jìn)化信息，筆者認(rèn)為有望將其開發(fā)成新型生物防治制劑和農(nóng)業(yè)肥料添加劑等在農(nóng)業(yè)生產(chǎn)領(lǐng)域廣泛應(yīng)用。