解淀粉芽孢桿菌群體遺傳特征分析

王依婷，李素珍，毛志濤，廖小平，許慧敏，邱藝華，閆雪崧，楊鳳妍，李貞景，郭慶彬*，劉歡歡*

（1.天津科技大學 食品科學與工程學院，天津 300457；2.思念食品有限公司，河南 鄭州 450000；3.中國科學院 天津工業生物技術研究所，天津 300308；4.濱州市沾化區海洋和漁業發展服務中心，山東 濱州 256600）

解淀粉芽孢桿菌（Bacillus amyloliquefaciens）是一類革蘭氏陽性菌，最初被定義為可產α-淀粉酶的枯草芽孢桿菌，廣泛用于生物防治、傳統發酵食品等領域[1-3]。B.amyloliquefaciens 具有較廣的抑菌譜，能通過自身代謝產生抑菌蛋白、聚酮化合物等來達到抑制病原菌生長繁殖的作用。劉冬等[4]從大豆根際土壤分離得到的菌株JDF3，可抑制大豆疫霉孢子囊的產生，對大豆疫病的防治效果達到70.70%。Rasiya 等[5]從印度紅樹內生細菌中分離到的菌株RKEA3，其產生的抗菌脂肽能夠抑制Staphylococcus aureus、Streptococcus pyogenes和Cryptococcus neoformis。同時，B.amyloliquefaciens 也可用于醬油[6]、白酒[7]、豆豉[8]等傳統食品的發酵。此外，B.amyloliquefaciens 也被發現對黃曲霉毒素B1 誘導的小鼠盲腸炎癥具有保護作用[9]，表現出一定的益生活性。

由于B.amyloliquefaciens 在農業與食品領域中的廣泛應用，越來越多的菌株已完成了基因組測序。截止目前，NCBI 數據庫已收錄超過150 個基因組序列，為更好地理解該菌株的生理生化特征以及工業化提供了支撐，同時也為解析該物種群體遺傳的多樣性和系統進化以及功能基因的多態性、基因組變異與表型關聯、正向基因篩選等提供了數據基礎。近些年發展起來的泛基因組學技術（pan-genome）[10]，通過群體內大量的基因組來完整描述一個物種全部的遺傳信息[11]，側重于分析特定菌株基因的存在/缺失對種群功能的獲得/缺失變異（presence/absence variation，PAV）之間的關聯。其中，該種群所有樣本共有的基因稱為核心基因，一般與物種生物學功能和主要表型特征相關；僅在部分樣本中存在的基因稱為附屬基因，一般與物種對特定環境的適應性或特有的生物學特征相關，反映了部分物種的特性。泛基因組是全面解析物種多樣性和研究表型差異背后分子機制的強有力工具[12]。

本文從泛基因組學角度出發，通過Prokka 基因組注釋[13]、Roary 泛基因組構建[14]、ProCAST 蛋白功能注釋、antiSMASH 次級代謝產物基因簇注釋[15]以及ABRicate 抗生素耐受性基因注釋等分析流程，構建B.amyloliquefaciens 的泛基因組并進行深入分析，挖掘B.amyloliquefaciens 物種的普遍性遺傳特征及碳水化合物的利用、植物促生、生物防治等方面的遺傳特異性，以期為深度利用B.amyloliquefaciens 這種微生物資源作參考。

1 材料與方法

1.1 材料

從NCBI 數據庫（https：//www.ncbi.nlm.nih.gov/datasets/genomes）下載組裝水平標識為“Complete”的66株B.amyloliquefaciens 基因組序列（fasta 格式，截止2022 年8 月17 日）。

1.2 方法

1.2.1 平均核苷酸一致性分析

平均核苷酸相似度（average nucleotide identity，ANI）是在核苷酸水平比較兩個基因組親緣關系的指標。FastANI 可將基因組序列分割為短序列片段[16]，使用Mashmap 算法計算同源映射并估計ANI，在本研究中，參數fragLen 設定為500 bp。之后對所有基因組的ANI 數值矩陣在R 語言中進行層級聚類，用iTOL v6[17]（https：//itol.embl.de/）繪制聚類樹圖。

1.2.2 泛基因組構建及分析

采用Prokka 注釋流程[13]及默認參數，對所有基因組進行基因結構預測和功能注釋，其中獲得的gff 格式文件作為泛基因組的輸入文件。Roary 通過blastp 對gff文件中的序列進行直系同源基因分析[14]，獲得B.amyloliquefaciens 泛基因組（pan-genome）和核心基因組（core genome）。采用Origin 軟件對泛基因組/核心基因組編碼基因數目與基因組數目之間的關系進行擬合[18]，擬合公式如下。

y=A×xB+C

式中：y 為泛基因組編碼基因數目；x 為新增基因組數目；A、B、C 為系數。對于核心基因組的變化，設新增核心基因組數目為x，泛基因組大小為y，A、B、C 為系數，擬合公式如下。

y=A×eBx+C

使用mafft 對核心基因組進行序列比對，并采用FastTree 及GTR 算法對比對文件構建基于核心基因組序列的系統發育樹[19]。同時，對Roary 生成的泛基因組PAV 矩陣，進行層級聚類并構建樹圖。

1.2.3 核心基因組功能注釋及富集分析

利用ProCAST 平臺（https：//procast.biodesign.ac.cn）對Roary 分析獲得的核心基因組蛋白序列進行功能注釋。該管路是由中科院天津工業生物技術研究所開發的綜合蛋白質注釋的無服務器網絡工具，涵蓋蛋白家族、域、分類、結構、基因本體和途徑信息等26 個數據庫，并提供多個工具用于結果的可視化分析。

1.2.4 質粒、抗生素耐藥性等基因的鑒定

采用ABRicate 流程（Seemann T，Abricate，https：//github.com/tseemann/abricate），整合CARD[20]、BacMet[21]、PHAGE、PlasmidFinder[22]、NCBI AMRFinderPlus[23]數據庫，對66 個菌株基因組的所有蛋白編碼序列（coding sequence，CDS）進行功能注釋。其中，CARD 與NCBI AMRFinderPlus 為抗生素耐藥性信息數據庫；BacMet為抗菌殺菌劑和金屬抗性基因數據庫；PHAGE 為噬菌體基因數據庫；PlasmidFinder 包含質粒復制子序列，可用于原始和組裝測序數據的復制子序列分析。

1.2.5 次級代謝產物基因簇的鑒定

采用antiSMASH 本地化軟件流程[15]，對66 個菌株基因組注釋產生的gbk 文件進行次級代謝產物基因簇挖掘，參數采用--cb-general、--cb-subclusters、--cbknownclusters、--asf 及--pfam2go。

1.2.6 碳水化合物水解酶基因的鑒定

采用本地化的CAZy 數據庫對66 個菌株基因組的所有蛋白編碼序列（CDS）進行功能注釋[24]。

1.2.7 解磷途徑及吲哚乙酸（indoleacetic acid，IAA）合成途徑相關基因的鑒定

根據京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）中IAA 合成途徑中的蛋白質酶編號[25]，以及解磷途徑涉及的酶編號[26]，在UniProt 數據庫中檢索已校驗的fasta 蛋白序列，并通過blast 建庫，對66 個基因組蛋白序列進行比對，identity閾值設定為40%，鑒定潛在的IAA 合成和解磷途徑。

2 結果與分析

2.1 平均核苷酸一致性分析結果

ANI 被定義為兩個微生物基因組同源片段之間平均的堿基相似度，是從全基因組水平評估物種間親緣關系，尤其是近緣物種的重要工具，通常以ANI>95%作為劃分物種邊界的閾值[16]。本研究通過使用fastANI對66 個B.amyloliquefaciens 基因組進行同源性評估，結果如圖1 所示。

圖1 基于ANI 分析的B.amyloliquefaciens 層級聚類樹Fig.1 Hierarchical clustering based on ANI values from B.amyloliquefaciens genomes

如圖1 所示，所有基因組可聚類為2 個主要分支Clade I 與II，其中，在分支內兩兩基因組之間的ANI數值均高于95%，而分支間的兩個基因組ANI 數值均高于90%而低于95%[16]。因此，Clade I 和Clade II 之間已經出現較為明顯的物種邊界，未來可將這2 個分支進行獨立命名。

2.2 泛基因組功能分析

利用Prokka 對所有基因組進行基因結構預測，結果見表1。

表1 B.amyloliquefaciens 基本信息Table 1 Basic information of B.amyloliquefaciens

由表1 可知，B.amyloliquefaciens 的基因組介于3.70～4.44 Mb，最小為RD7-7，最大為XJ5，GC 含量為45.50%～46.50%，編碼基因數目為3 716～4 445 個，平均每個基因組編碼3 997 個基因。結合Roary 泛基因組分析流程，對上述基因組的核心基因組和泛基因組進行鑒定，結果見圖2。

圖2 核心基因數量和總基因數量隨基因組數目增長的變化趨勢Fig.2 Variation trend of core gene number and total gene number with the increase of genome number

由圖2 可知，B.amyloliquefaciens 泛基因組基因總數為13 051 個，核心基因組中包含基因2 165 個，占基因組平均編碼基因數目的54.16%。

由于物種受到生存環境與生態位差異、有效群體規模、多態性水平差異等因素的影響，泛基因組呈現開放型和閉合型兩種模式。開放型泛基因組特點是物種內隨著測序菌株的增加，基因個數也隨之增加，較難到達平臺期；閉合型泛基因組特點是物種內隨著有限的菌株個體增加，基因個數迅速達到平臺期[27-28]。

在對泛基因組編碼基因數目與基因組數目之間的擬合公式中，y=3 925x0.29-343.8（R2=0.998）。其中，0.29為非零正數，泛基因組中基因總數會隨著新的基因組的加入而呈指數性增加，因此推斷B.amyloliquefaciens的泛基因組呈現一定的開放性。

核心基因組編碼基因數目與基因組數目之間的擬合公式為y=1 388e-0.1571x+2 235（R2=0.980）。根據擬合方程，隨著B.amyloliquefaciens 基因組測序數量的增加，核心基因組編碼基因數目可能穩定在2 235 個。上述分析表明，B.amyloliquefaciens 易從外界獲得新的基因，但同時核心基因組具有較強的穩定性。

泛基因組的缺失與存在變異（PAV）進行聚類分析，泛基因組中核心基因序列的系統發育分析結果如圖3 和圖4 所示。

圖3 B.amyloliquefaciens 基因存在/缺失變異（PAV）矩陣的層次聚類Fig.3 Hierarchical clustering based on gene presence/absence variation（PAV）matrix of B.amyloliquefaciens

圖4 基于B.amyloliquefaciens 核心基因組序列的系統發育樹Fig.4 Phylogenetic tree based on core genome sequence alignment of B.amyloliquefaciens

由圖3 與圖4 可知，此結果與ANI 結果一致，表明66 株菌在全基因組ANI 相似性、PAV 聚類以及核心基因組序列相似性中均可被明確劃為兩個分支，Clade II中的15 株菌與Clade I 中的51 株菌不屬于一個種，將來可作為獨立物種進行研究。

為進一步解析該物種所具有的核心基因功能特征，即菌株的普遍性遺傳特征，本研究采用ProCAST中集成的21 個數據庫對核心基因組蛋白序列進行了注釋。其中NR、SwissPort 和PFAM 數據庫的注釋基因占核心基因數的比例均超過90%。

基于COG 數據庫的功能注釋結果如圖5 所示。

圖5 核心基因組的COG 注釋結果Fig.5 Main annotations for COG annotation

如圖5 可知，R（11.14%）、E（8.56%）、J（7.30%）、K（7.03%）、P（5.71%）、C（5.65%）、M（5.10%）、G（4.83%）等占據了主要的功能模塊。核心基因組的PATHWAY注釋結果見圖6。

圖6 核心基因組的PATHWAY 注釋結果Fig.6 Main annotations for PATHWAY annotation

在PATHWAY 注釋中（包含KEGG pathway、Reactome、MetaCyc 等數據庫），Aminoacyl-tRNA biosynthesis、Mitochondrial translation elongation 以及Purine metablism 等模塊聚集了大量的功能基因，顯示出B.amyloliquefaciens 在轉錄、翻譯、氨基酸代謝、離子運輸、能量代謝等基礎代謝方面功能保守。

CAZy 是一個整合了碳水化合物酶活性的數據庫。B.amyloliquefaciens 的核心基因組注釋結果見圖7。

圖7 核心基因組的CAZy 注釋結果Fig.7 Main annotations for CAZy annotation

如圖7 所示，糖苷轉移酶類（glycoside transferases，GT）數量為50 個，糖苷水解酶（glycoside hydrolases，GH）數量為49，占CAZy 酶總數的37.04%和36.30%。GH 和GT 是參與糖的合成和代謝的主要酶類，其中GH13 家族具有α-淀粉酶活性，占CAZy 酶總數的5.19%；GH3 和GH51 兩個家族具有β-葡聚糖酶活性，占CAZy 酶總數的3.70%；GH3 和GH43 兩個家族具有木聚糖酶活性，共占CAZy 酶總數3.70%。表明B.amyloliquefaciens 具有較好的水解淀粉、β-葡聚糖和木聚糖潛力。這為B.amyloliquefaciens 的糖苷水解酶和糖苷轉移酶類開發利用提供了遺傳基礎。已有研究表明，B.amyloliquefaciens 可以產生β-1，3-葡聚糖酶等病程相關蛋白，通過作用于真菌細胞壁的不同成分來破壞真菌細胞壁，以達到抵抗病原真菌侵染的目的[29-30]。核心基因組的基因本體論注釋結果見圖8。

圖8 核心基因組的GO 注釋結果Fig.8 Main annotations for GO annotation

如圖8 所示，在基因本體論（Gene Ontology，GO）注釋中，B.amyloliquefaciens 大量核心基因聚集在oxidoreductase activity（7.14%）、structural constituent of ribosome（5.42%）、integral component of membrane（2.97%）等涉及氧化還原、遺傳信息傳遞和生物膜合成等生物學過程模塊；在細胞組成模塊中，主要涉及catalytic activity（9.25%）、protein binding（1.98%）、cytoplasm（0.66%）等；而在分子功能模塊中，大量功能基因聚集在DNA binding（10.58%）、catalytic activity（9.25%）、transmembrane transporter activity（7.87%）等方面，表明B.amyloliquefaciens 在催化分解、蛋白質的合成以及轉運方面較為保守。上述核心基因組的功能注釋證明B.amyloliquefaciens 在基礎代謝相關的基因中顯示出一定的遺傳和功能保守性，為研究該菌的普遍代謝特征提供了基礎。

此外，利用TMHMM 數據庫對蛋白質跨膜螺旋結構域進行預測，同時采用SIGNALP 數據庫對氨基酸序列中的信號肽進行注釋。將含有信號肽且無跨膜結構域的蛋白質識別為分泌蛋白。結果顯示，在B.amyloliquefaciens 的核心基因組中，共識別出99 個分泌蛋白，占核心基因總數的4.57%，并且在分泌蛋白中發現了溶桿菌素的蛋白質序列，表明B.amyloliquefaciens普遍能夠分泌溶桿菌素這一抗菌肽，賦予了該物種廣譜抗菌的特性。

2.3 所有B.amyloliquefaciens 菌株基因組功能模塊的PAV 分析

因B.amyloliquefaciens 在生物防治以及食品等領域應用廣泛，因此對所有66 株菌的CAZy、解磷途徑、IAA 合成途徑、次級代謝產物基因簇以及抗生素耐藥性基因進行全面比較分析，為比較B.amyloliquefaciens種群內生物功能提供理論基礎。

將66 株B.amyloliquefaciens 與其他芽孢桿菌（例如：B.subtilis 168、B.pumilus Ha06YP001 等）對比，結果見圖9。

圖9 CAZyme 編碼基因的分布Fig.9 CAZyme coding gene distribution

如圖9 所示，B.amyloliquefaciens GH3、GH43 以及GH51 基因數較多。其中，GH3 和GH51 家族為β-葡聚糖酶，能催化水解谷物細胞壁中的β-葡聚糖[29]；GH43 是一類降解木聚糖的酶系，可降解自然界中大量存在的木聚糖類半纖維素。此外，B.amyloliquefaciens ARP23 等30 株菌的GH13 家族編碼基因數普遍較多，而GH13 家族為支鏈淀粉酶，能夠催化淀粉等多糖化合物中的α-1，6-糖苷鍵水解，形成直鏈淀粉。以上結果表明66 株B.amyloliquefaciens 具有豐富的β-葡聚糖和木聚糖降解酶的基因資源。依賴色氨酸的IAA 生物合成途徑見圖10。解磷途徑和IAA 合成途徑相關酶的基因簇分布見圖11。

圖10 依賴色氨酸的IAA 生物合成途徑Fig.10 Tryptophan dependent IAA biosynthesis pathway

圖11 解磷途徑和IAA 合成途徑相關酶的基因簇分布Fig.11 Gene distribution of enzymes related to phosphorus hydrolysis and IAA synthesis pathway

如圖11 所示，吲哚丙酮酸（indole pyruvic acid，IPyA）途徑為66 株菌所共有，而其它途徑均不完整或缺失。如圖11 所示，15 株菌（Clade II）在醛脫氫酶（NAD+）的基因存在與缺失上與其余51 株菌（Clade I）呈“互補”狀態。在IPyA 介導的IAA 合成途徑中，色氨酸通過氧化轉氨作用形成吲哚丙酮酸，隨后在苯丙酮酸脫羧酶作用下脫羧形成吲哚乙醛，再經過醛脫氫酶（NAD+）脫氫變成IAA，而IAA 對于植物生長具有重要影響。

此外，通過blast 分析比對，所有菌株的基因組中共含有5 種編碼堿性磷酸酯酶（4 個Pho D 和1 個Pho A）基因和1 種植酸酶（phy）編碼基因。由圖11 可知，這些基因在66 株菌中均被發現。微生物的解磷過程十分復雜，目前被廣泛接受的是有機磷（organophosphorus，OPB）酶解和無機磷（inorganic phosphorus，IPB）酸解兩個方式[31]。堿性磷酸酶是解磷菌參與OPB 礦化的重要酶類，分為phoA、phoD、phoX 3 種類型[32]，其中在海洋和土壤微生物中phoD 的含量更高[33]。磷酸酶通過去磷酸化或斷裂有機磷物質中的磷酯鍵釋放有效磷，植酸酶則通過釋放植酸磷來增加有效磷的含量[34]。

2.4 抗性基因的挖掘與次級代謝產物基因簇

為揭示不同菌株在次級代謝產物合成以及抗生素耐藥性等方面的特征，本研究利用本地化的anti SMASH 及ABRicate 注釋流程對每個B.amyloliquefaciens 的基因組進行了功能挖掘，結果如圖12 所示。

圖12 B.amyloliquefaciens 次級代謝產物和抗性、噬菌體以及質粒基因簇分布Fig.12 Distribution of B.amyloliquefaciens secondary metabolites resistance，phage and plasmid gene clusters

如圖12 所示，在所有B.amyloliquefaciens 基因組中，共鑒定出5 種抗藥性基因和3 種銅離子抗性基因，菌株覆蓋率高達98.48%，只有B.amyloliquefaciens B408未鑒定出抗性基因。其中YFMP、YFMO、FABL、CUER數量最多，出現頻次各為66 次；其次為rphC、satAbsub、clbA、rphB、tet（L）、CSOR、COPZ，共254 個。YFMP、YFMO、CUER、CSOR、COPZ 蛋白表現出對銅的抗性；FABL 蛋白表現出對酚類化合物的抗性；rphB 和rphC蛋白對利福霉素具有耐受性；clbA 和satA-bsub 對鏈霉素具有耐受性，且clbA 對林可酰胺以及大環內酯類抗生素也具有耐受性，tet（L）對四環素以及多西環素具有耐藥性。

噬菌體和質粒的存在可以增強細菌對環境的適應性。通過PlasmidFinder 數據庫注釋發現，B.amyloliquefaciens HM618、TL106、LL3、IT -45、S499、LS1 -002-014s 和MBE1283 中含有質粒，覆蓋率達10.61%；PHAGE 數據庫對66 株B.amyloliquefaciens 進行噬菌體基因注釋，共檢測出16 種噬菌體基因，覆蓋率達100%，涉及參與細菌素生產或免疫、小酸溶性孢子蛋白和硫氧還蛋白的合成等功能，為B.amyloliquefaciens 的抗菌、遺傳穩定、環境耐受等特性提供遺傳物質基礎。

antiSMASH 注釋結果表明，66 株B.amyloliquefaciens 共預測到815 個基因簇，每個基因組內約包含10～15 個次級代謝基因簇。其中，與已知基因簇同源性大于60%的有441 個，同源性為100%的有354 個，主要涉及桿菌烯（bacillaene）、桿菌肽（bacillibactin）和溶桿菌素（bacilysin）等的生產（圖12），這可能是B.amyloliquefaciens 作為生防菌株的主要遺傳基礎。

上述分析表明，作為生防菌株，B.amyloliquefaciens有望通過次級代謝產物遏制病原菌來實現抑菌效果，并通過抗藥性/噬菌體/質粒增強自我保護與環境適應性；B.amyloliquefaciens 對銅離子具有較好的抗性，表明在土壤受到銅污染時，具有自我保護的能力。江春玉等[35]從重金屬污染土壤中分離到Cu2+抗性B.amyloliquefaciens，可顯著增加培養液中游離Cu2+的含量，從而提高植物吸收Cu2+速率，增強Cu2+污染土壤治理的效率，但相關機制尚不明確。

3 結論

本研究通過對B.amyloliquefaciens 進行泛基因組構建以及對核心基因組的功能注釋，從基因組層面評估了其作為生物防治菌株的遺傳潛力。ANI 以及聚類分析表明66 株B.amyloliquefaciens 分為相對獨立的兩個分支；COG、GO、PATHWAY 等對核心基因組基礎注釋表明B.amyloliquefaciens 具有較為保守的能量生產、蛋白質合成、基因表達、化合物轉運等功能，這為B.amyloliquefaciens 的代謝穩定性提供了遺傳基礎。

在功能基因注釋上，CAZy 注釋結果表明B.amyloliquefaciens 具有較強的水解多糖類物質的能力。PATHWAY 和SwissProt 注釋結果表明B.amyloliquefaciens 可通過IPyA 途徑合成IAA 以及通過分泌植酸和堿性磷酸酶降解環境中的磷，促進植物的生長。

在66 株B.amyloliquefaciens 的次級代謝產物合成基因簇分析中預測到大量的桿菌烯和桿菌肽合成基因簇。同時抗生素耐藥性分析表明，多數B.amyloliquefaciens 基因組中含有大量編碼利福霉素、鏈霉素、林可酰胺以及大環內酯類抗生素等耐藥性基因以及銅抗性基因，上述基因元件進一步增強了該菌在遏制病原菌、聯合用藥、自我保護和土壤改良方面的優勢和潛力。總之，本文通過構建B.amyloliquefaciens 泛基因組，系統揭示了該物種內普遍性及特異性遺傳特征，為將來深度利用B.amyloliquefaciens 這種微生物資源奠定了基礎。