基于宏基因組測序技術解析高溫快速發酵豆豉菌群結構與功能注釋

冼芳瑩，趙文鵬，王思宇，李 浩，楊慧林，王筱蘭

（江西師范大學生命科學學院，江西南昌 330022）

豆豉是以黑豆或黃豆為原料，經微生物發酵制成的一種風味獨特的調味品。曲霉、毛霉、細菌等是發酵豆豉的主要菌種[1]。曲霉型豆豉的生產主要包括原料處理、制曲、（后）發酵三個階段。制曲主要利用米曲霉等分泌的多種酶對原料初步分解，發酵則主要利用細菌分泌的多種酶對洗曲后的原料再次分解。研究表明，豆豉發酵可產生大量的功能營養物質，如游離氨基酸、大豆異黃酮、大豆低聚糖、纖溶酶、降壓肽等[2-4]，提升了大豆的營養價值及藥用價值。另外，微生物在發酵階段廣泛參與醇發酵、酯化、美拉德和脫氨等反應，形成獨特的豉香風味[5]。值得注意的是，傳統豆豉室溫發酵周期以30～60 d居多，且發酵過程中涉及的背景菌群組成復雜，生產成本高，產品批次穩定性差[1]。解析豆豉發酵的關鍵功能菌并明確其在發酵中扮演的角色，可為豆豉多菌混合發酵提供新思路，也將為縮短發酵周期保障風味與營養提供技術參考。

利用可培養技術結合變性梯度凝膠電泳（DGGE）、一代測序技術分析發酵食品的主要功能菌群是傳統的研究思路[6-7]，但上述傳統方法通量低，靈敏度受到局限，如DGGE方法，不同實驗條件（核酸提取、PCR擴增程序及引物的選擇）可能導致帶型譜圖的復現力欠佳，從而給序列信息的探針設計和菌群結構分析帶來負面影響[8]。宏基因組高通量測序可在免培養的情況下對樣品中全部微小生物的基因組DNA進行測定，從而能快速確定微生物菌群結構、進行功能基因分析與代謝特征分析等[9-11]。本課題組前期已采用高通量測序技術揭示曲霉型豆豉在室溫發酵過程中的菌群多樣性[12]，以及菌群動態跟蹤[13]。本研究涉及的高溫快速發酵曲霉型豆豉與傳統室溫發酵豆豉不同的是，其發酵溫度大于50 ℃，且大大縮短發酵周期，而對于其高溫快速發酵過程中其菌群結構分布、功能預測還尚未有報道。

本研究以經高溫發酵的曲霉型豆豉為研究對象，采用宏基因組測序技術，得到菌群物種結構信息，重點探尋豆豉高溫發酵過程菌群功能變化，輔以優勢菌屬的功能貢獻度分析，探討主要菌屬在功能代謝中扮演的角色。為豆豉的多菌混合發酵工藝的改進提供依據，也可為后續運用多組學技術探討豆豉微生物與功能風味內在機制的研究工作奠定一定基礎。

1 材料與方法

1.1 材料和儀器

OMEGA Soil DNA Kit Omega公司；瓊脂糖Biosharp公司。

電泳儀 北京市六一儀器廠；凝膠成像儀Biorad公司；HiSeq4000 Illumina公司；臺式高速離心機 上海生工生物工程有限公司；NanoDrop ND-2000超微量分光光度計 ThermoScientific公司。

1.2 樣品采集

曲霉型豆豉采自南昌稻香園調味食品有限公司，選取發酵1、3、5和7 d的高溫發酵豆豉（制曲第3 d，以55～60 ℃發酵至第8 d）樣品，裝入已滅菌的自封袋，共4袋。每袋樣品均是從豆豉發酵罐不同位置（四個對稱點及中心）自上而下各取出約20 g后，充分混勻后形成的混合樣，每袋約100 g。采集完畢后及時將所采豆豉樣品冰盒保存，并立即送返實驗室，于-80 ℃超低溫備用。

1.3 豆豉宏基因提取

稱取所采豆豉樣品10.0 g，于無菌研缽中研磨充分后，取1.0 g，根據E.Z.N.A. Soil DNA Kit 試劑盒實驗說明書步驟，提取樣品中微生物宏基因組DNA。對所得宏基因組DNA樣本進行質檢，包括完整性、濃度以及純度（1.0%瓊脂糖凝膠電泳檢測，NanoDrop超微量分光光度計）。收集OD260/OD280值介于1.8～2.0的合格DNA樣品（電泳條帶清晰、無明顯拖尾），置于-20 ℃備用。

1.4 PE文庫構建及基因組測序

將宏基因組樣品片段化（Covaris M220儀器），至約300 bp。為構建PE文庫，使用TeuSeqTMDNA Sample Prep Kit試劑盒進行試驗，并通過變性獲得單鏈DNA片段。后進行橋式PCR（HiSeq 3000/4000 PE Cluster Kit），以獲取Cluster簇，并將所獲簇線性化為單鏈。借助Illumina HiSeq 4000平臺、以HiSeq 3000/4000 SBS Kits完成宏基因組測序。逐次記錄并統計熒光信號結果，獲得模板DNA的片段序列數據。

1.5 序列質控與組裝

剪切序列的3’和5’接頭序列（軟件Seqprep）后，去除經質量剪切后長度＜50 bp、平均質量值＜20且含N堿基的片段（Sickle軟件），保留高質雙端測序數據。將所測數據與大豆（Glycine max）基因組序列比對，去除相似度較高的宿主序列片段（BWA軟件）。基于質控后的序列數據和de Brujin組裝策略，以Megahit軟件從小到大迭代片段的k-mer值，質控后將測序深度低和測序深度高的片段逐步進行拼接，保留＞300 bp的重疊群（contigs），選擇最優組裝結果[11]。將未映射的序列片段再次拼接組裝，對＜1000 bp的片段再次用Newbler軟件拼接，得到＞1000 bp的contigs后，與之前各樣品contigs進行合并，為最終拼裝結果。

1.6 基因預測與非冗余基因集構建

針對所獲contigs，以MetaGene軟件預測其開放閱讀框（Open Reading Frame，ORF），對核酸長度≥100 bp的基因進行翻譯，統計各樣品的基因預測結果。以CD-HIT軟件對預測的基因序列進行聚類分析（90% coverage、95% identity），以每類中最長基因為代表序列，進行非冗余基因集構建。使用SOAPaligner軟件將各樣品中的優質序列與非冗余基因集（95% identity）進行比較，并統計各樣品中基因的豐度。

1.7 物種與功能注釋及物種貢獻度分析

本研究基于I-Sanger生物信息云平臺（www.isanger.com）進行非冗余基因集的數據庫比對。借助Diamond軟件[14]，進行非冗余基因集BLASTP比對（e-value≤10-5），包含GenBank的非冗余蛋白質序列（Non-Redundant Protein Sequence Database，NR）和京都基因百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）數據庫。由NR物種注釋結果，獲得物種豐度信息；基于KEGG比對結果，利用KOBAS 2.0軟件功能注釋后獲得KO、Pathway、EC等類別的基因豐度信息后，根據KEGG level 3注釋結果，進一步采用ipath v3.0工具標記代謝通路總豐度前50的pathway，以變異系數（Coefficient of Variation，CV）對豆豉各發酵階段其代謝通路的差異進行描述性統計。采用碳水化合物活性酶（Carbohydrate active enzyme database，CAZy）數據庫的hmmscan，比對非冗余基因集（CAZyme 6.0數據庫、evalue≤10-5），獲得該基因集所對應的碳水化合物活性酶注釋的豐度數據。采用定制的R軟件包分析主要菌群（屬水平）對每個功能途徑的貢獻。由于四個樣本的測序數據的大小和測序深度在一定程度上存在不同，故采用RPKM（Reads Per Kilobase Per Million Mapped）計算物種與功能注釋、物種貢獻度分析中的基因豐度或相對豐度，以避免直接使用基因數目（Reads Numbers）造成的統計性誤差。

1.8 數據處理

本研究涉及的4樣品的宏基因組原始序列信息已提交至NCBI網站上的SRA數據庫。項目檢索號：SRP193072，4樣品檢索號：SRR8931515～SRR 8931518。

2 結果與分析

2.1 宏基因組測序與質控結果

利用Illumina測序平臺得到宏基因原始序列，在經過質控、去宿主處理后，4組樣品得到reads數各約為5400、5600、4300、4700萬條（表1），總長各約為7.51、7.76、5.99、6.46 Gbp。對優化后序列拼接組裝后，均獲得10萬條以上的contigs。對總contigs進行ORF預測，共獲得1324268個ORF，均長610.36 bp，其中201至600 bp長度的基因比例占比達到52.75%?；贠RF的預測基因進一步聚類、構建，共得到982672個unigenes，均長657.51 bp。N50表示將contigs由大到小疊加至contigs總長的50%時對應的contig長度，通常其值越大，組裝效果越好[15]；本次組裝結果顯示，各組N50均大于500 bp，最大值達7907 bp，表明本次組裝良好，可用于進行后續分析。

表1 宏基因組測序的質控及組裝信息Table 1 Overview of metagenomics sequencing quality control and assembly information

2.2 發酵豆豉宏基因組NR物種注釋

NR是非冗余蛋白質序列數據庫，將非冗余基因集與其進行比對從而獲得物種分類學信息。NR注釋結果表明，在域水平上，細菌的相對豐度高達98.9%～99.8%，真菌與古細菌僅0.03%～0.3%，表明在高溫快速發酵過程中細菌極有可能是豆豉風味及其他特征的主要承擔者[16]。以下集中從門、屬及種三水平分析豆豉高溫發酵階段菌群結構變化趨勢。

在門水平上，發酵中優勢門為厚壁菌門（Firmicutes）、擬桿菌門（Bacteroidetes）、放線菌門（Actinobacteria）及變形菌門（Proteobacteria）（圖1a），這與張菊所報導的門水平菌群結構相似[17]。但在發酵過程中以上菌門的相對豐度均存在較大波動，如厚壁菌門的相對豐度在發酵前期（Day 1）高達87.1%，此后連續下降，在發酵中期（Day 3～5）達到最低27.3%，但在發酵后期（Day 7）迅速回升至61.8%。擬桿菌門的趨勢恰好相反，發酵前期相對豐度不足0.1%，但在發酵中后期躍升為兩大優勢門之一，其相對豐度于中期Day 5達到最高值64.6%后在后期回落，表明兩個門可能存在競爭關系。

屬水平上，擬桿菌屬（Bacteroides）、魏斯氏菌屬（Weissella）、葡萄球菌屬（Staphylococcus）、乳桿菌屬（Lactobacillus）、真桿菌屬（Eubacterium）、柔嫩梭菌屬（Faecalibacterium）、普雷沃氏菌屬（Prevotella）為優勢屬（圖1b），這與過往的報道非常相似[18]。魏斯氏菌屬、乳桿菌屬、片球菌屬（Pediococcus）及腸球菌屬（Enterococcus）4種乳酸菌僅在發酵前期擁有優勢地位，相對豐度占比總和高達69.91%，之前的研究表明[19]，乳酸菌對發酵食品的防腐與良好風味形成具有重要作用。擬桿菌屬是中后期唯一全程優勢屬，其豐度呈先升后降的趨勢，于Day 5達到最高相對豐度（52.5%）。發酵前期的特有屬數目達324屬，遠大于其他三階段（圖1c），而發酵中后期共有屬達713個。

圖1 豆豉各發酵階段菌群的分布Fig.1 Distribution of flora during various fermentation stages in Douchi

種水平注釋結果表明，發酵前期種數目達5900種，中后期回落至約3800種，表明豆豉發酵過程可能存在微生物量減少的現象；明顯優勢種為泰國魏斯氏菌（Weissella thailandensis）、戊糖片球菌（Pediococcus pentosaceus）、普拉梭菌（Faecalibacterium prausnitzii）、同形擬桿菌（Bacteroides uniformis）、腸道普氏菌（Prevotella copri）（圖1d）。其中，泰國魏斯氏菌在前期為第一優勢菌，該菌參與細菌素和酸的代謝，不僅有助于風味物質形成，而且能抑制低耐酸性病原體的生長[12]。與屬水平類似的是，發酵前期6大優勢種即占據了51.8%的相對豐度；中后期優勢種數目更多且分布更均一。需要指出的是，由于種水平上菌種遺傳同源性高、差異小，諸多種水平注釋并不能與現有菌種庫匹配，因此完全無法分類的種亦占一定比例。

2.3 豆豉發酵階段宏基因組KEGG功能注釋

2.3.1 發酵階段KEGG通路的分布 KEGG是將基因組與功能信息聯系起來，進行系統分析基因功能的數據庫。經KEGG比對，共有85.94%的unigenes得到注釋，其中得到KEGG Pathway注釋的unigenes超過27萬條（27.54%）。KEGG通路（level 1）注釋顯示，發酵中后期的絕對功能豐度均遠低于發酵前期（RPKM=699879），這與NR注釋的種水平分類結果呼應，證明發酵中菌群總數具有減少的趨勢[20]（圖2），從而導致代謝總量減少；各階段的基因豐度特征均為代謝相關基因占比最多且均超過50%，說明豆豉發酵過程中代謝活性較高。與發酵前期相比，中后期屬于代謝通路的基因比例均不同程度上升，而屬于環境信息處理的比例卻明顯下降。對代謝通路相關基因展開二級KEGG注釋，結果（45類）表明，在整個發酵過程中，相對豐度主要集中在碳水化合物代謝，全局和概覽、氨基酸代謝、核苷酸代謝、復制與修復等通路。各功能類別變化趨勢為，發酵中后期菌群的氨基酸代謝、全局和概覽、能量代謝等功能的相對豐度較發酵前期來說明顯上升，而與增殖相關的功能類別如核苷酸代謝、復制與修復、膜轉運等明顯下降。氨基酸代謝與豆豉風味的形成密切相關[21]，因此發酵中后期可能生成大量風味物質；相對前期，高溫發酵的中后期增殖代謝相關基因表達量減少，這可能是中后期的發酵環境對菌群繁殖不利所致。

圖2 不同后酵階段豆豉樣品的KEGG 代謝通路（level 1、level 2）分布Fig.2 Distribution of KEGG Pathway predicted in different fermentation stage Douchi samples at level 1 and 2

2.3.2 發酵階段主要代謝通路的分布 基于KEGG通路（level 3）基因注釋結果，利用ipath工具對發酵中總豐度前50條的代謝通路進行標注，以了解豆豉發酵的主要代謝通路分布規律。結果顯示（圖3），就整個發酵過程來看，波動偏大（CV＞15%）的代謝通路（紫色、紅色、藍色）共18條，其中各階段波動差異均偏大（紫色）僅3條；其余15條代謝通路具有階段偏向性，即這些通路在中后期波動偏小，但在前期至中期間波動明顯。其中有5條（紅色）在發酵前期相對豐度明顯更高，如丙酮酸代謝（ko00620）、丁酸代謝（ko00650）、甘油酯代謝（ko00561），為豆豉各類風味物質形成奠定基礎[19]。據前文所述，乳酸菌為前期主要優勢菌，因此，乳酸菌可能與上述5條代謝通路在前期的高豐度表現有關。在中后期相對豐度更高的10條代謝通路（藍色），反而是與氨基酸或能量代謝有關，包括Phe、Tyr和Trp的生物合成（ko00400）、泛酸與輔酶A（CoA）的生物合成（ko00770）、卟啉與葉綠素代謝（ko00860）、Arg的生物合成（ko00220）等。這種具有階段偏向性的表現意味著菌群的整體代謝特征在發酵前期到中期階段發生較大轉變，而在中期至后期代謝組成相對穩定，結合NR物種注釋結果，這可能與菌群結構在前中期存在較大波動、而中后期相對穩定有關。其中芳香族氨基酸（Phe、Tyr和Trp）可通過轉氨、脫氫、脫羧、還原等反應產生酸、醇、醛、酯等風味物質[22]，表明在發酵中后期可能生成較多風味物質。

圖3 不同發酵階段KEGG主要代謝通路的全局（ipath）差異分析Fig.3 Global metabolic pathway （iPath） difference analysis for top 50 KEGG metabolic pathways during various fermentation stages

2.3.3 厭氧發酵過程關鍵功能酶基因豐度的分布發酵過程中微生物所分泌的酶催化一系列代謝反應，產生各種風味物質，從而形成豆豉獨特的風味[23]。為更好地了解高溫快速發酵豆豉的發酵代謝特征，基于KEGG酶基因豐度注釋結果，對發酵中某些關鍵酶基因豐度進行了比較。結果顯示，發酵中后期，纖維素及淀粉酶（EC 3.2.1.1、EC 3.2.1.4）基因豐度較發酵前期明顯更高，但與部分寡糖分解酶（如EC 3.2.1.21）橫向對比，其豐度明顯較低，這表明菌群可能更傾向于分解寡糖（圖4），所生成單糖可供菌體用于后續代謝，參與風味物質的生成。在發酵中期，Day 5的纖維素酶基因豐度較Day 3反而降低，此前研究表明[24]，經過適應環境后菌群可能提高纖維素酶表達量以將纖維素分解成寡糖，隨著環境中寡糖量上升，纖維素酶表達減少，而寡糖分解酶表達上升。從厭氧發酵主要的有機酸產物（乳酸、乙酸、琥珀酸、丁酸等）來看，主要有機酸合成的關鍵酶基因豐度變化趨勢不一。中后期乳酸的同型或異型發酵的關鍵酶基因豐度較前期下降70%以上，產乙酸、琥珀酸關鍵酶基因的變化趨勢與乳酸類似；而產丁酸的關鍵酶基因豐度雖很低，但變化趨勢與乳酸相反。此外，乙偶姻作為吡嗪類物質的前體，與其相關的功能酶（EC 4.1.1.5、EC 1.1.1.303/304）基因在前期數目較多，自Day 3起急劇回落并維持近零水平。結合NR注釋結果，上述生成有機酸關鍵酶基因低豐度可能與發酵中后期魏斯氏菌屬、乳桿菌屬、葡萄球菌屬等5種菌屬失去優勢地位存在一定的相關性。

圖4 豆豉厭氧發酵不同階段重要功能酶的基因豐度的分布Fig.4 Distribution of gene abundance related to some vital functional enzymes at various stages in Douchi anaerobic fermentation summarily

2.4 豆豉宏基因組的CAZy功能注釋

KEGG注釋結果提示，作為與碳水化合物代謝密切相關的CAZy家族，其中各族數目極有可能與樣品中所注釋基因的豐度存在一定聯系，故從CAZy數據庫對unigenes展開注釋。CAZy注釋結果顯示（圖5a），糖苷水解酶（Glycoside Hydrolases，GHs）家族被注釋到的基因數目最多，糖苷轉移酶（Glycosyl Transferases，GTs）家族次之排列。與發酵前期相比，GH家族在發酵中后期的豐度占比增加，GTs變化趨勢則較為平穩。研究表明[25]，GHs負責水解糖苷鍵，生成單糖、寡糖或糖復合物，在寡糖、芳香基糖苷的合成、氨基酸和多肽的糖基化方面具有重要作用。而GTs負責催化糖分子與其它分子（如蛋白、寡糖、脂等）連接[26-28]。豐富的GHs和GTs為單糖、寡糖的形成、轉移及后續代謝提供了基礎。

針對CAZy各亞家族分布特征展開統計，結果發現，GT兩大亞家族GT2、GT4豐度最高，且二者在發酵中穩定維持在較高水平（圖5b）。部分GH亞家族（如GH2、GH3、GH92等）表現出在前期豐度較低、中后期反而較高的現象，它們主要參與纖維素、寡糖的分解，負責甘露糖苷、葡萄糖苷、乳糖苷等糖苷殘基的剪切[29]，表明在發酵中后期碳水化合物代謝中可能存在著較高水平的糖類分解代謝，生成豐富的單糖資源，為豆豉風味的形成提供可利用的還原糖和風味前體物質[30]。

圖5 豆豉發酵樣品的碳水化合物活性酶的功能分類（科或綱水平）Fig.5 Carbohydrate-active enzyme (CAZy) functional classification of Douchi fermentation sample at Family or Class level

2.5 豆豉微生物的功能貢獻度分析

通過計算微生物各屬水平在各代謝通路的相對豐度，探索豆豉發酵菌群與KEGG代謝通路之間的功能貢獻度。結果顯示，在碳代謝、氨基酸的生物合成、氨基糖和核苷酸糖代謝等5大代謝通路中，貢獻度前10屬皆分屬厚壁菌、擬桿菌兩大優勢門（圖6）。結合門水平的NR注釋結果，推測菌屬的相對豐度是其對代謝通路功能貢獻度重要因素。發酵中后期的優勢菌屬擬桿菌屬在Day 5對各大代謝通路貢獻度均超50%，對氨基糖和核苷酸糖代謝的貢獻度更是高達55.54%。

圖6 豆豉發酵過程核心菌群與部分KEGG類別的貢獻度分析Fig.6 Analysis of contribution of dominant flora and patial KEGG categories flora during various fermentation stages in Douchi

另外，就相對豐度前15的菌屬對各CAZy家族的功能貢獻度來看，發酵前期的魏斯氏、葡萄球菌等4個優勢屬即囊括CAZy六家族70%以上的功能貢獻度，其中乳酸菌魏斯氏屬與乳桿菌屬的貢獻總和明顯較高，表明發酵前期可能存在一定量的乳酸積累[19]；其中魏斯氏菌屬與多糖裂解酶類（Polysaccharide Lyases，PLs）相關基因的貢獻度達37.60%，表明此時魏斯氏屬可能與多糖的裂解存在較大的相關性；乳桿菌屬對PLs的功能貢獻度的分布與魏斯氏菌屬基本相反（圖7）。在發酵中后期，擬桿菌屬對PL家族的相對貢獻度總和最高，說明在后兩階段多糖的裂解主要與擬桿菌相關。

圖7 豆豉發酵過程核心菌群與碳水化合物活性酶各家族的貢獻度分析Fig.7 Analysis of contribution of dominant flora and all CAZy families during various fermentation stages in Douchi

3 討論

3.1 物種豐度分布

當下，宏基因組數據庫比對已然成為一項普遍而有力的基因功能分析技術，幫助科研人員從系統探究了解基因所對應功能，以進一步揭示菌群構建的微生態所具有的生物學意義[10]。與以往分別進行細菌或真菌的16S/18S測序不同，本研究基于同一統計水平獲得了豆豉高溫發酵豆豉細菌及真菌的相對豐度。結果中顯示真菌豐度遠不及以往報道所述，基于該豆豉工藝及高溫發酵特征，推測洗曲過程存在降低真菌豐度的可能性[31]；同傳統型豆豉發酵工藝相比，快速發酵工藝采用高溫發酵，以此推測該結果或與真菌及細菌間高溫耐受及生長速率的差異有關。以往一些16S測序結果表明[12,31]，豆豉發酵中主要含厚壁菌門、變形菌門、放線菌門的細菌，這與本研究中發酵前期門的分布較為類似。但發酵中期的擬桿菌門，在過去往往未檢出或僅某階段少量出現，這或許與發酵工藝（發酵溫度、發酵濕度、發酵容器）因素有關。具體到屬來看，前期主要菌屬（乳桿菌屬、魏斯氏菌屬、片球菌屬等），均為室溫主產乳酸的細菌，與以往豆豉的菌屬分布也較大程度類似[18,31]。而中后期階段，以擬桿菌屬、普雷沃氏菌屬為代表的優勢菌屬在其他發酵食品卻鮮有報道，有趣的是，這些屬在諸多腸胃菌群研究中卻是?？蚚32-33]。將發酵與腸道環境特征對比后發現[34],偏酸（約pH6.0）、缺氧、營養基質豐富等均為兩類環境的共性。但驚訝的是，這些菌屬在55～60 ℃的高溫發酵環境仍可存活并增殖，或許與發酵前期的優勢菌因高溫（＞50 ℃）消失后預留較大的營養空間及高濕環境的散熱保護有關。而該高溫發酵豆豉在發酵中后期優勢屬豐度的波動過程，提示各菌屬可能在相繼對高溫環境有所適應后，緊接逐步進入生存空間爭奪的狀態。

種水平的注釋表明，豆豉發酵中種數目達3500～5900種，而采用可培養法及DGGE法篩選到的菌種僅有數十種[23,35]。這樣的差距既表明宏基因測序法的優越，也從種水平揭示了豆豉發酵菌群的豐富性。另外，豆豉發酵中十幾個甚至幾個主要種以較少的種數即占據了大多（50%～80%）的生存空間，這些優勢種（泰國魏斯氏菌、嗜酸乳桿菌、同型擬桿菌等）的豐度特點，或許可為今后豆豉多菌混合純種發酵提供新思路；而某些（條件）致病菌種（如Staphylococcussciuri）的檢出，側面說明發展科學、可控的豆豉發酵工藝將是必然趨勢。

3.2 功能豐度分布

KEGG代謝通路分析揭示了在發酵中后期菌群的氨基酸、能量代謝、輔酶與維生素代謝等與生長發育相關的通路較前期都得到一定提升，暗示菌群可能增強了其代謝強度以維持自身生存，而核酸代謝、復制與修復等與菌群繁殖相關的功能分布呈下降的趨勢則說明，或存在高溫環境迫使菌群采取偏保守的繁殖策略的可能。

從ipath注釋結果可以看出，中后期的許多氨基酸代謝通路（如ko00400，ko00220）得到一定增強，能量代謝也較為旺盛，這可能與發酵環境中豐富的游離氨基酸有關，氨基酸可能作為此時菌群氮源兼部分碳源。卟啉與葉綠素代謝通路（ko00860）的中間產物是某些氨基酸次級代謝的輔酶，該通路的富集或許與發酵中卟啉環的合成底物（甘氨酸、乙酸等）較豐富有關[36]。此外，不難發現，脂肪代謝通路整體較匱乏，可能與缺少高產脂肪酶菌群（霉菌、酵母、芽孢桿菌屬、葡萄球菌屬等）有關[37]。與豆豉室溫發酵相比，高溫快速發酵可能僅含有少量制曲及發酵初期殘留的脂肪酶，脂肪的生物分解幾乎陷入停滯。

從有機酸合成關鍵酶的基因來看，產甲酸、丁酸的酶基因通常被認為與梭菌、真桿菌等屬有關，其中梭菌屬被報道具有較強的產氨能力[38]，這或將給該豆豉引入某些不良風味。另外，伴隨著乳酸發酵菌屬的消失，在發酵中后期一些參與形成豆豉獨特風味的關鍵物質如乙醇、乙酸、乙偶姻等關鍵酶基因下降明顯，其中乙偶姻回落最為明顯，暗示諸多以乳酸發酵為主的菌屬，可能在豆豉酯類、醛類、吡嗪類等風味的形成中扮演關鍵角色[5]。

3.3 物種與功能貢獻度分布

從菌屬的功能貢獻度分布看出，魏斯氏菌的氨基酸合成基因豐度較高，表明魏斯氏菌可能更偏向游離氨基酸的分解及利用。片球菌屬的碳代謝貢獻度更高，暗示其可能更善于利用環境中碳源；GT家族常被認為胞外多糖的合成存在較大聯系，而片球菌屬的GT家族占比較高的現象一定程度解釋了張國華等報道的其高產胞外多糖的性狀[39]。擬桿菌在Day 5組中氨基糖與核糖代謝的高豐度的現象，推測與此時其自身的大量增殖有關；另外，注意到擬桿菌屬雖然碳代謝、氨基酸合成貢獻度與物種豐度基本一致，但其在CAZy中GH、GT、PL等家族的功能貢獻度明顯更高。某些腸道菌群研究指出，擬桿菌群在很多代謝障礙的機體較多出現，其往往與能量代謝負相關[33]，另有研究報道該菌屬具有較好的多糖降解能力[40]。結合本研究，推測在豆豉發酵中，擬桿菌門（擬桿菌屬、普氏菌屬等）可能在中后期較多地分解了糖苷、多糖，自我代謝后仍有剩余，剩余部分或被厚壁菌門某些菌屬生長利用，可能是促使厚壁菌門豐度在豆豉發酵的后期快速回升的原因之一。更高層面上看，豆豉發酵即是一個菌群廣泛分享各類中間代謝產物的過程，這將加速各類原料基質的分解。

然而，需要指出的是，以宏基因組學為代表的基因測序技術，往往揭示的并不是菌群的存活現狀及代謝狀態，從樣品中提取的可能只是菌群死亡后尚未完全分解的DNA[9]。最新的發酵食品研究中表明[41-42]，若要以更深入、細致的動態視角剖析豆豉風味與菌群結構變化的內在機理，明確各菌種在營養風味物質演變中的功能貢獻，將代謝組學、宏轉錄組學等技術與宏基因組深度測序進行多組學聯合將很有必要，這也是本課題組下一步研究工作的重要探索方向。

4 結論

厚壁菌門、擬桿菌門是曲霉型豆豉高溫發酵的兩大優勢門，二者存在明顯的競爭關系，厚壁菌門乳酸發酵菌屬為前期主要發酵菌屬，相對豐度達69.91%，擬桿菌門擬桿菌屬為發酵中后期主要發酵菌屬，發酵過程中菌種資源豐富，高達5900種。發酵初期與中后期在菌種數量、組成以及代謝表現上皆存在較大波動，中后期反而有趨于平穩之勢。豆豉高溫發酵過程中菌群代謝旺盛，活躍基因中與代謝相關的基因相對豐度超過50%，其中菌群氨基酸合成代謝的基因有一定增加，與部分有機酸合成相關的酶基因大幅減少；主要菌屬對氨基酸的利用、糖苷的水解合成等功能的貢獻存在一定差異?？焖侔l酵豆豉雖利用高溫較大縮減了生產周期，但可能對產品風味有一定程度的影響，未來有必要聯合代謝組學、宏轉錄組學，對其發酵過程進行更為深入的研究。