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

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

多粘類芽孢桿菌（Paenibacillus polymyxa）屬芽孢桿菌科（Bacillaceae）類芽孢桿菌屬（Paenibacillus），是一種產(chǎn)芽孢的革蘭氏陽性菌[1]。多粘類芽孢桿菌在自然界中分布廣泛，在玉米、高粱、小麥等多種植物體內(nèi)外以及根際土壤中均有被發(fā)現(xiàn)，是一類重要的根際有益菌[2]。多粘類芽孢桿菌對(duì)植物生長發(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]。

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

目前，有關(guān)多粘類芽孢桿菌的研究多集中在其生物特性及其對(duì)植物與生物防治的作用機(jī)理上，而對(duì)其基因序列的分析研究報(bào)道較少[10-12]。筆者從土壤中分離出一株多粘類芽孢桿菌LY1，利用基因組測(cè)序技術(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 測(cè)序及數(shù)據(jù)質(zhì)控取菌體委托上海生工生物工程股份有限公司使用Illumina Hiseq platform 進(jìn)行基因組測(cè)序。Illumina Hiseq?得到的原始圖像數(shù)據(jù)文件經(jīng)CASAVA 堿基識(shí)別 （Base Calling）分析轉(zhuǎn)化為原始測(cè)序序列（Sequenced Reads），對(duì)原始數(shù)據(jù)質(zhì)量值等信息進(jìn)行統(tǒng)計(jì)，并使用FastQC 對(duì)樣本的測(cè)序數(shù)據(jù)質(zhì)量進(jìn)行可視化評(píng)估。測(cè)序得到的原始數(shù)據(jù)，里面含有帶接頭的、低質(zhì)量的序列。為了保證信息分析質(zhì)量，必須對(duì)原始數(shù)據(jù)進(jìn)行過濾，得到Clean 數(shù)據(jù)。隨機(jī)從Clean 數(shù)據(jù)中抽取10 000 條序列與NCBI NT 數(shù)據(jù)庫進(jìn)行blastn 比對(duì)，取evalue ≤1E-10 并且相似度＞90%、coverage ＞80%的比對(duì)結(jié)果，計(jì)算其物種分布，同時(shí)進(jìn)行污染檢測(cè)。

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

1.2.3 基因預(yù)測(cè)與注釋采用Prokka 對(duì)組裝結(jié)果進(jìn)行基因元件預(yù)測(cè)。Prokka 是一系列基因元件預(yù)測(cè)工具的集合，調(diào)用Prodigal 預(yù)測(cè)編碼基因，Aragorn 預(yù)測(cè)tRNA，RNAmmer 預(yù)測(cè)rRNA，Infernal 預(yù)測(cè)miscRNA，預(yù)測(cè)出的各類基因元件匯總并完成初步注釋。采用Repeat Modeler 對(duì)組裝結(jié)果進(jìn)行重復(fù)序列Denovo 預(yù)測(cè)，再利用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ù)測(cè)得到的16S rRNA序列與NCBI 的16S 數(shù)據(jù)庫進(jìn)行Blast 比對(duì)，設(shè)置參數(shù)identify ＞95。然后選取identify 最高的前30 條16S rRNA序列，利用muscle軟件進(jìn)行序列多重比對(duì)后，采用FastTree 軟件構(gòu)建系統(tǒng)發(fā)育樹。

2 結(jié)果與分析

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

二代測(cè)序和三代測(cè)序原始數(shù)據(jù)經(jīng)過去除帶接頭的、低質(zhì)量的序列，過濾后得到Clean 數(shù)據(jù)。Illumina Hiseq 原始測(cè)序數(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原始測(cè)序數(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ù)測(cè)與基因組基本特征

2.2.1 基因組的組成與信息2 組測(cè)序數(shù)據(jù)經(jīng)過序列矯正、Gap 修補(bǔ)、剪輯錯(cuò)誤修復(fù)、拼裝等程序后獲得菌株的基因組信息，如圖1 所示，多粘類芽孢桿菌LY1 的基因組為一條環(huán)狀閉合DNA，大小 5 765 474 bp，共預(yù)測(cè)到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è)。基因長度和GC 含量的分布如圖2 所示。

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

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

2.2.2 基因組注釋將預(yù)測(cè)的編碼基因與多個(gè)數(shù)據(jù)庫進(jìn)行比對(duì)，注釋到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ù)庫的對(duì)比，可以查看物種轉(zhuǎn)錄本序列與相近物種的近似情況以及同源序列的功能信息[13]。由圖3 可知，鑒定到LY1 菌株有3 792 個(gè)基因與多粘類芽孢桿菌同源，與Paenibacillus polymyxaM1 菌株的相似性相對(duì)較高，在M1 的編碼基因中注釋到318 個(gè)基因，而在Paenibacillus polymyxaSQR-21、Paenibacillus polymyxaSC2 編 碼 基因中則分別注釋到76 和29 個(gè)基因。

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

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ù)庫，在給出染色體中一套完整基因的情況下，它可以對(duì)蛋白質(zhì)交互網(wǎng)絡(luò)在各種細(xì)胞活動(dòng)起的作用作出預(yù)測(cè)，方便地尋找與行使某一類功能相關(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)建而成，蛋白編碼序列與之比對(duì)后可以預(yù)測(cè)蛋白的功能[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)行比對(duì)，將基因與其對(duì)應(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 相對(duì)應(yīng)的DNA 序列，存在于所有細(xì)菌的基因組中，具有高度的保守性和特異性，是病原菌檢測(cè)和鑒定的一種有效參照[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é)

目前，人們對(duì)多粘類芽孢桿菌的功能有了較為全面地了解，但其分子作用機(jī)制以及遺傳工程改造等方面的研究報(bào)道還不多[3]。該研究通過全基因組測(cè)序分析了多粘類芽孢桿菌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)用。