基于枯草芽孢桿菌多組學數(shù)據(jù)全基因組規(guī)模代謝模型與蛋白質降解模擬

袁 萍，江明鋒，官久強，安添午，張翔飛，羅曉林*

（1.西南民族大學青藏高原動物遺傳資源保護與利用，四川省教育部重點實驗室，成都 610041；2.四川省草原科學研究院，成都 611731）

枯草芽孢桿菌（Bacillus subtilis）是一種典型革蘭氏陽性細菌，分布廣泛，易分離培養(yǎng)，具有抑制病原菌生長、增強免疫應答、促進養(yǎng)殖動物生長、凈化水質等特點，在飼料、食品行業(yè)、基因工程和醫(yī)藥等領域應用廣泛[1]。關于枯草芽孢桿菌培養(yǎng)的微生物學、遺傳學和分子生物學方法非常成熟，已測定多株枯草芽孢桿菌基因組序列[2]?？莶菅挎邨U菌可利用低成本碳氮源大規(guī)模培養(yǎng)，因其可分泌大量胞外酶，如堿性蛋白酶（Extracellular alkaline protease，aprE）等，可將環(huán)境中不溶性蛋白質降解為多種寡肽和氨基酸，并從中衍生出其他含氮化合物供菌體生長[3]。基于此特性，有較強蛋白質降解能力的枯草芽孢桿菌被篩選用于降低水體蛋白質含量。

基因組規(guī)模代謝網(wǎng)絡模型（Genome-scale metabolic models，GEMs）通過基于約束的流量平衡算法，利用攝取和分泌通量有限數(shù)據(jù)預測代謝表型[4]。隨組學技術發(fā)展，GEMs提出與應用已近20年?；谄鋵蛐秃捅硇完P系綜合性描述能力，GEMs已廣泛應用于設計代謝工程策略和方面[5]。目前GEMs重構在模式微生物如大腸桿菌等發(fā)展和應用趨于成熟，非模式生物中應用待進一步研究[5-6]?？莶菅挎邨U菌代謝模型已有相關報道，但利用GEMs通過分子與蛋白水平預測枯草芽孢桿菌對環(huán)境中蛋白質降解能力研究較少。

隨水產(chǎn)養(yǎng)殖規(guī)模不斷擴大，養(yǎng)殖水體富營養(yǎng)化成為水產(chǎn)養(yǎng)殖主要問題。已有聯(lián)合利用枯草芽孢桿菌、光合細菌、硝化菌及反硝化菌等降解水體有機氮，改良水體富營養(yǎng)化相關研究。目前主要集中在調整不同細菌比例，根據(jù)有機氮或氨氮降解水平表型數(shù)據(jù)評估和篩選高效有機氮降解菌株，卻未見利用枯草芽孢桿菌基因組和代謝數(shù)據(jù)系統(tǒng)研究并基于計算推測其蛋白質降解能力的相關報道。本試驗首次利用全基因組規(guī)模代謝模型作為整合調控途徑與枯草芽孢桿菌碳氮循環(huán)代謝模擬不同碳氮條件下枯草芽孢桿菌降解蛋白質，使用菌株試驗數(shù)據(jù)和細胞內(nèi)通量數(shù)據(jù)評估模型準確性，為系統(tǒng)篩選并最大化枯草芽孢桿菌蛋白質降解能力提供參考。

1 材料與方法

1.1 試劑與培養(yǎng)基

本研究DNA操作所需酶和試劑盒均購自Takara Biotechnology（日本）和Parkson基因技術有限公司（中國上海）。其他試劑均購自國藥控股化學試劑有限公司（中國上海），具有最高分析級。

Luria-Bertani（LB）培養(yǎng)基和基本培養(yǎng)基用于細菌培養(yǎng)。LB培養(yǎng)基由NaCl 10 g·L-1，胰蛋白10 g·L-1，酵母提取物5 g·L-1和1 L H2O組成。根據(jù)需要使用1.5%瓊脂固化。

魚粉蛋白培養(yǎng)基用于蛋白降解菌分離篩選。魚粉蛋白培養(yǎng)基由含 NaCl 0.5 g·L-1，MgSO4·7H2O 0.5 g·L-1，NH4Cl 2 g·L-1，KH2PO40.2 g·L-1，K2HPO41 g·L-1，H2O 1 L和0.5 mL微量培養(yǎng)基（CoCl2·6H2O 0.1 g·L-1，MnCl2·4H2O 0.425 g·L-1，ZnCl20.05 g·L-1，NiCl2·6H2O 0.01 g·L-1，CuSO4·5H2O 0.015 g·L-1，Na2MoO4·2H2O 0.01 g·L-1，Na2SeO4·2H2O 0.01 g·L-1，pH 7.0）組成無機鹽培養(yǎng)基加入10 g·L-1魚粉蛋白作為唯一碳氮源。

改良葡萄糖蛋白胨培養(yǎng)基用于測定不同碳氮源條件下菌株生長速率和蛋白質降解速率，葡萄糖10 g·L-1，NaCl 5 g·L-1，不同濃度蛋白胨（40、30、20、10、5、3.3和2.5 g·L-1）。碳氮源比分別為1∶1，1∶2，1∶3，1∶4，2∶1，3∶1，4∶1。

1.2 菌株分離鑒定

本研究所用枯草芽孢桿菌菌株采自四川省成都二龍橋魚塘底部2～7 cm處泥層，用淤泥采樣器取50 g左右泥樣。采用文獻[7-8]方法從養(yǎng)殖水體中分離高效蛋白質降解菌株。將10 g底泥接種到500 mL魚粉蛋白培養(yǎng)基中，37℃180 r·min-1培養(yǎng)72 h，提取10 mL培養(yǎng)基接種于500 mL魚粉蛋白培養(yǎng)基中，37℃180 r·min-1培養(yǎng)24 h。將少量培養(yǎng)基涂布在固體LB培養(yǎng)基上，37℃培養(yǎng)24 h。挑取單克隆菌斑，接種于LB培養(yǎng)基擴大培養(yǎng)。使用通用引物27F（5'AGAGTTTGATCCTGGCTCAG 3'）和1492R（5'TACGGYTACCTTGTTACGACTT 3'）PCR擴增篩選得到菌株，獲得其16S rDNA。通過NCBI BLAST（https://blast.ncbi.nlm.nih.gov/）對該 16S rDNA序列線上比對，利用MEGA X構建其系統(tǒng)發(fā)育樹[9]。

1.3 菌株生長速率和菌株蛋白質降解速率鑒定

本研究通過測定培養(yǎng)基中菌體干重變化計算不同培養(yǎng)條件下菌株生長速率。首先，將接種菌株的培養(yǎng)基在37℃180 r·min-1分別培養(yǎng)0、12和24 h后搖勻抽取50 mL，離心取沉淀加純水50 mL混勻沖洗，然后離心取沉淀，重復3次。最后將沉淀烘干稱重。

使用凱氏定氮法測定培養(yǎng)基中蛋白含量。將接種菌株LB培養(yǎng)基或魚粉蛋白培養(yǎng)基在37℃180 r·min-1培養(yǎng)0和12 h后搖勻抽取50 mL，離心取上清液10 mL加入Viskase公司MD44透析袋（8 000～12 000 D）透析去除無機氮源、氨基酸和小分子寡肽。定量透析產(chǎn)物并使用凱氏定氮儀（上海赫冠儀器有限公司）測定蛋白質含量。

1.4 枯草芽孢桿菌HD-2菌株基因組DNA提取建庫測序及生物信息學分析

HD-2基因組DNA提取使用TaKaRa公司生產(chǎn)細菌基因組提取試劑盒。完成基因組DNA抽提后，利用1%瓊脂糖凝膠電泳檢測收集基因組DNA。使用Covaris M220將基因組DNA打斷為長度300～500 bp片段。使用TruSeqTMDNA Sample Prep Kit建庫試劑盒構建測序文庫并使用TruSeq PE Cluster Kit v3-cBot-HS試劑盒擴增。最終文庫通過Illumina Hiseq 2000測序平臺完成掃描測序。測序產(chǎn)生的DNA序列數(shù)據(jù)以FastQ格式記錄，對原始測序數(shù)據(jù)質控并去除測序錯誤和低質量測序片段。使用SOAPdenovo v2.04拼接軟件對測序片段組裝。使用Glimmer 3.02軟件預測細菌基因。

1.5 基因組規(guī)模代謝模型構建和生長模擬校驗

使用 RAST軟件（http://rast.nmpdr.org/）注釋HD-2基因組預測得到的基因序列并通過MODEL SEED軟件生成包含有完整Gene-Protein-Reaction（GPR）列表初始代謝草圖。根據(jù)KEGG和BioCyc代謝數(shù)據(jù)完善并設定反應方向。使用MEMOTE軟件評估模型完整性[10]，填補模型缺口，最終構建菌株HD-2基因組規(guī)模代謝模型。

本試驗使用Cytoscape2.8.2軟件作網(wǎng)絡可視化，并使用其Network Analyzer插件分析網(wǎng)絡拓撲結構，查找核心功能所在子網(wǎng)絡[11]。此后，使用COBRA2.0分析流平衡[12]。分別使用生長速率最大和單一產(chǎn)物量最大作為目標函數(shù)，函數(shù)如下：

maximize：biomass/target product

subject：S*v=0

0.01 ≤biomass≤0.02

通過模擬菌株在不同碳源和氮源條件下菌株最大生長速率，將模擬最大生長速率下碳氮比與試驗中結果比較，評估模型準確性。

1.6 枯草芽孢桿菌不同碳氮源降解蛋白質預測

利用上述構建的GEMs模型，通過改變碳源氮源輸入濃度，模擬aprE基因不同碳氮比營養(yǎng)條件下表達變化。將文獻報道不同碳氮營養(yǎng)條件下枯草芽孢桿菌生長數(shù)據(jù)[3]與試驗測量結果和模擬結果比較。將模擬結果中蛋白酶合成速率最大碳氮源比例等同于試驗中測定蛋白質降解速率最大時碳氮源比例，總蛋白酶合成量作為目標函數(shù)，最大葡萄糖流量設定為0.0050 g·DW（dry weight，干重）-1·h-1，即包含0.0020 g碳元素碳源設置為碳源輸入量，將包含0.0004 g氮元素 0.0025 g·DW-1·h-1總蛋白作為細胞外初始氮源輸入量（初始條件：碳氮比為5∶1）。此后，以0.0005 g·DW-1·h-1梯度逐步提高細胞外氮源輸入量至0.0125 g·DW-1·h-1作為細胞外最終氮源輸入量，模擬碳氮比1∶1營養(yǎng)條件。分別模擬蛋白質降解速率最大時碳氮比和蛋白酶合成量最大時碳氮比。調整模型中各反應對應流量變化范圍，使模擬結果與試驗結果數(shù)據(jù)變化趨勢一致。最后，利用調整后模型模擬分析連續(xù)梯度變化碳氮比條件下枯草芽孢桿菌蛋白質降解，確定枯草芽孢桿菌蛋白質降解最適碳氮源比例。

2 結果與分析

2.1 菌株選擇與鑒定

本研究分離鑒定得到1株不溶性蛋白質高效降解菌株HD-2，12 h培養(yǎng)基中不溶性蛋白降解率達89.76%。表明HD-2具有高不溶性蛋白質降解能力。以枯草芽孢桿菌HD-2基因組DNA為模板，利用細菌通用引物擴增16S rDNA，1.5 kb處有1條特異性條帶。純化16S rDNA的PCR擴增產(chǎn)物，并作序列測定。結果表明，枯草芽孢桿菌HD-2的16S rDNA長度為1 531 bp。在GenBank中BLAST比對分析，結果表明與菌株HD-2的16S rDNA序列相似性最高的是枯草芽孢桿菌BS49和BS34A菌株，達100%。利用其親緣關系標注系統(tǒng)發(fā)育樹見圖1，結果表明菌株HD-2屬于芽孢桿菌屬枯草芽孢桿菌。

2.2 枯草芽孢桿菌HD-2菌株基因組測序組裝和基因預測結果

枯草芽孢桿菌HD-2菌株基因組掃描測序共產(chǎn)生讀段（reads）15 412 018對，包含1 557 602 908個堿基，其中Q20以上reads占比92.77%。經(jīng)質控保留有效reads共15 049 172對，1 415 852 922個堿基。有效reads組裝得到206個長序列片段（scaffolds），包含5 753 251個堿基，測序reads覆蓋深度234x，詳細測序和組裝結果見表1。

枯草芽孢桿菌HD-2的因預測結果包含5 585個可能基因序列，共4 953 741 bp個堿基?；蚱骄L度886 bp，基因座平均分布密度為每千堿基分布0.97個基因。基因中平均GC含量為68.1%。全部基因序列儲存為ffn格式文件。基因組測序數(shù)據(jù)和組裝結果已提交至NCBI Genome數(shù)據(jù)庫，編號為SUB8287272。

圖1 利用HD-2菌株16S rDNA序列構建系統(tǒng)發(fā)育樹Fig.1 Phylogenetic analysis of the 16S rDNA of strain HD-2

表1 枯草芽孢桿菌HD-2菌株基因組組裝結果Table 1 Results of B.subtilis strain HD-2 genome assembly

2.3 枯草芽孢桿菌HD-2菌株基因組注釋

通過RAST在線快速注釋軟件對枯草芽孢桿菌HD-2菌株作基因組注釋。注釋結果中包含已知功能基因1 968個，分屬30個主要子系統(tǒng)，包括各類生物大分子合成和代謝、細胞結構合成、化合物轉運合成和生物防御運動等細胞過程。其中與蛋白酶合成直接相關氮代謝相關基因20個，蛋白質合成相關基因129個和蛋白質降解相關基因25個（見圖2）。

HD-2編碼外分泌蛋白酶是6個包含編碼多個開放閱讀框基因簇。這6種蛋白酶分別為胞外中性金屬蛋白酶（EC 3.4.24.28）、解封氨肽酶（EC 3.4.11.-）、YpdF氨肽酶、D-丙氨酰-D-丙氨酸羧肽酶（EC 3.4.16.4）、耐高溫羧肽酶1（EC 3.4.17.19）和氨肽酶PepA（EC 3.4.11.1）。

2.4 枯草芽孢桿菌HD-2基因組規(guī)模代謝模型iBac2019構建

本試驗構建HD-2菌株基因組尺度代謝模型iBac2019。該模型分為細胞外和細胞質兩個區(qū)域，共包含846個基因，1 394個反應和1 490種化合物，其中包含82個轉運反應模擬大分子跨膜運輸（見表2）。iBac2019還包含17個虛擬反應，模擬轉錄因子調控下游基因表達和被碳氮源供應量限制的表達水平。這17個虛擬反應包含存在于枯草芽孢桿菌中蛋白酶表達調控通路（見圖3），將轉錄調控過程轉化為虛擬生化反應，以便定量流平衡分析。

圖2 枯草芽孢桿菌HD-2注釋基因功能分類Fig.2 Gene annotation and classification in strain HD-2

表2 iBac2019模型基本參數(shù)Table 2 A summary of GEMs model iBac2019

圖3 枯草芽孢桿菌外泌蛋白酶表達調控通路Fig.3 Regulation pathways of exosprotease encoding genes in B.subtilis

2.5 iBac2019不同碳氮比下最大生長速率模擬

實驗室分別測定培養(yǎng)基碳氮比為1∶1、1∶2、1∶4、2∶1、4∶1條件下培養(yǎng)12 h菌體干重。根據(jù)試驗結果，所得不同碳氮條件下相應生長速率如下：碳氮比為1∶1時0.0124 g·DW-1·h-1；碳氮比為1∶2時0.0093g·DW-1·h-1；碳氮比為1∶4時0.0051g·DW-1·h-1；碳氮比為2∶1時0.0212 g·DW-1·h-1；碳氮比為4∶1時0.0083 g·DW-1·h-1（見表3）。使用基于約束的重構、分析和其軟件中內(nèi)置通量變異性分析，確定每個反應通量邊界。在此基礎上構建iBac2019模型Biomass方程，見附件2。

將目標條件設置為生長速率最大，模擬碳氮比為1∶3和3∶1條件下枯草芽孢桿菌最大生長速率分別為0.0072和0.0128 g·DW-1·h-1，與該條件下試驗測量值相差8.3%（0.0066 g·DW-1·h-1）和18.7%（0.0104 g·DW-1·h-1），模擬值與測量值變化趨勢吻合，表明模型仿真度較高。

2.6 通過目標產(chǎn)物最大化實驗模擬蛋白質降解最適碳氮比

通過MCODE插件尋找模型中堿性蛋白酶合成和轉運子網(wǎng)絡。子網(wǎng)絡拓撲分析表明蛋白酶合成通路和碳代謝通路關系緊密。通過與注釋信息比較，選擇模型中共有糖酵解、糖異生（30個反應），TCA循環(huán)（18個反應），丙氨酸、天冬氨酸代謝（15個反應）以及胞內(nèi)、胞外蛋白水解（22個反應）4個代謝通路共85個反應組成蛋白酶分泌及其相關代謝反應生物功能模塊。蛋白水解相關酶類共有胞內(nèi)蛋白水解酶16個和胞外蛋白水解酶6個（見表4）。

表3 不同碳氮源比例條件下枯草芽孢桿菌HD-2生長速率和蛋白質濃度變化Table 3 Growth rate of B.subtilis HD-2 and the changes in protein concentration under different carbon-nitrogen ratio

表4 iBac2019模型中與蛋白降解相關的蛋白（酶）Table 4 Proteins(enzymes)related to protein degradation in iBac2019 model

續(xù)表

將該模塊內(nèi)總蛋白酶合成量作為目標函數(shù)，將最大葡萄糖流量設定為0.0050 g·DW-1·h-1，即包含0.0020 g碳元素設置為碳源輸入量，將包含0.0004 g氮元素0.0025 g·DW-1·h-1總蛋白作為細胞外初始氮源輸入量，初始條件下碳氮比5∶1。在此條件上，以g·DW-1·h-1梯度逐步提高細胞外氮源輸入量至0.0125 g·DW-1·h-1總蛋白作為細胞外最終氮源輸入量，模擬碳氮比達1∶1。將胞外總蛋白上限設置為輸入蛋白量10倍，添加蛋白酶催化并轉運蛋白進入細胞質虛擬反應。模擬菌株分泌蛋白酶，降解蛋白為可利用肽鏈并增加氮源輸入流量。設置上述條件，利用流平衡分析模擬該模塊，并記錄蛋白酶合成量和菌株最大生長量（見圖4）。

圖4 不同氮源細胞流量下枯草芽孢桿菌代謝模型生長速率和蛋白酶合成速率模擬結果Fig.4 Simulation results of growth rate and protein synthesis in B.subtilis under different quantity of nitrogen source

結果表明，枯草芽孢桿菌在氮源進入細胞流量達到0.0650 g·DW-1·h-1時生長量達到頂峰，之后生長速率停止增加。即在碳氮比約為2∶1時，生長速率進一步增加開始受碳源供應量約束，單純增加氮源供應無法提高菌株生長速度。氮源進入細胞流量增加至0.0350 g·DW-1·h-1時達到最大單位蛋白酶合成速率，氮源進入細胞流量因氮源供應不足而限制菌株蛋白酶合成。單位蛋白酶合成速率達到最大時進一步供應氮源雖增加菌株生長速率，但抑制蛋白酶合成分泌。最終，在考慮菌株生長速度和單位蛋白酶合成速率兩個影響因素后，總蛋白酶合成量如圖4所示，氮源進入細胞流量達0.05000 g·DW-1·h-1時達到最大，此時碳氮比約為2.8∶1。進一步試驗驗證，在碳氮比達2.6∶1時，枯草芽孢桿菌HD-2菌株溶液蛋白質降解能力最高，達2.0000 mg·L-1·h-1。

3 討論

近年來，基因組規(guī)模代謝網(wǎng)絡模型廣泛應用于真核及原核生物。由于原核單細胞細菌有較為簡單的基因組結構及代謝通路，其代謝模型研究最為深入有效。目前，GEMs已在藥物生產(chǎn)、能源開發(fā)和健康監(jiān)測等多方面發(fā)揮重要作用。藥物生產(chǎn)方面，游動放線菌基因組規(guī)模代謝模型iYLW1028用于發(fā)酵模擬，顯著提高阿卡波糖含量[13]。能源開發(fā)方面，代謝網(wǎng)絡模型應用于提高藍細菌生物乙醇生產(chǎn)效率[14]。在健康監(jiān)測方面，腸道菌群在機體健康中發(fā)揮重要作用，腸道菌群全基因組規(guī)模代謝模型已廣泛應用于集體營養(yǎng)狀況檢測等[15]?？莶菅挎邨U菌作為一種模式細菌，其代謝網(wǎng)絡模型已構建并發(fā)展至今[16]。與傳統(tǒng)GEMs相比，本模型整合最新組學數(shù)據(jù)，顯著提高預測代謝通量分布和表型能力。

枯草芽孢桿菌降氮能力運用廣泛。據(jù)報道，枯草芽孢桿菌生物肥料可有效控制農(nóng)業(yè)氨氮排放，保持較高農(nóng)作物產(chǎn)量并減輕環(huán)境干擾[17]。此外，將枯草芽孢桿菌直接飼喂可降低犬糞便中氮發(fā)酵產(chǎn)物濃度和氣味[18]。在農(nóng)業(yè)及環(huán)境領域，篩選具有環(huán)境蛋白質降解能力枯草芽孢桿菌菌株意義重大。本研究在水產(chǎn)養(yǎng)殖環(huán)境中分離鑒定得到1株高效降解蛋白質菌株HD-2，利用其基因組掃描測序結果進一步探究其蛋白質降解潛能。

從機制上說，枯草芽孢桿菌可降解環(huán)境中不溶性蛋白質，通過分泌胞外酶（如aprE基因產(chǎn)物等）使蛋白質分解為多肽等，可進一步攝取產(chǎn)生自身所需含氮化合物[19]。同時可通過電子傳遞鏈將氨轉化為氮氣釋放到空氣中，起到對水環(huán)境等脫氮效果[20]。大量分泌的蛋白酶和其基因表達調控因子表明這些胞外蛋白酶對宿主枯草芽孢桿菌的重要性，但這也對枯草芽孢桿菌生長代謝造成負擔，因此需根據(jù)其細胞所處營養(yǎng)狀況和外界環(huán)境嚴格控制枯草芽孢桿菌胞外蛋白酶表達[4]。本研究利用系統(tǒng)計算方法，調整HD-2所處營養(yǎng)條件，預測其蛋白質降解能力，探究枯草芽孢桿菌HD-2最適蛋白質降解條件。與現(xiàn)階段通過大規(guī)模篩選試驗尋找蛋白質降解能力強菌株相比，本研究提出一種更高效便捷獲得環(huán)境氮降解菌株的方法。

4 結論

本研究成功分離得到枯草芽孢桿菌菌株HD-2，通過序列測定確定其親緣關系?；蚪M注釋分析得到20個與蛋白酶合成直接相關氮代謝相關基因、129個與蛋白質合成相關基因、25個與蛋白質降解相關基因以及6個分泌胞外蛋白酶。此外，利用枯草芽孢桿菌HD-2多組學數(shù)據(jù)重構基因組規(guī)模代謝網(wǎng)絡模型iBac2019。利用此模型，預測不同碳氮條件下蛋白水解相關胞內(nèi)蛋白水解酶16個，胞外蛋白水解酶6個，獲得蛋白質高效降解最優(yōu)生長狀態(tài)為碳氮比2.8∶1。通過試驗驗證預測結果準確性，發(fā)現(xiàn)碳氮比達2.6∶1時，枯草芽孢桿菌HD-2菌株溶液蛋白質降解能力最高。本研究成功將該模型應用于枯草芽孢桿菌對蛋白質降解能力預測，為從系統(tǒng)水平研究并篩選具蛋白質降解能力枯草芽孢桿菌提供理論基礎。