基于宏基因組學測序技術分析3 個輪次高溫大曲微生物群落

王國崢，陳 筆，盧建軍，陳良強，張巧玲，羅 寒

（貴州茅臺股份有限公司，貴州省釀酒工程技術研究中心，貴州 遵義 564501）

曲是我國傳統白酒釀造的重要原料，曲研究對于釀酒乃至整個傳統釀造行業都有十分重要的意義[1-2]。醬香型白酒風味復雜，這與其發酵過程中曲用量高密切相關[2]，其生產可分為下沙、造沙、1～7輪次，每個輪次都需要添加大量曲粉作為糖化劑與發酵劑，但各輪次生產因所處的環境氣候不同可能導致大曲微生物群落結構與功能出現差異[3]。下沙、造沙輪次是醬香型白酒生產的投糧階段，高粱糊化程度較低[4]，1輪次在9 個輪次中環境均溫最低、環境濕度最低[5]，這些不利因素都會阻礙發酵進程，因此大曲的功能在前期輪次的生產中更為重要，不僅直接決定發酵的進程，也關系到后期輪次微生物代謝產物、前體物質及香味物質的積累與形成[6]。因此本研究選取這3 個生產輪次的生產用曲作為實驗對象，進行微生物群落結構及功能分析。

隨著現代分子生物學技術的發展，高通量測序技術逐漸被用于食品微生物的研究中[7-8]，旨在全面揭示微生物群落結構與功能信息。目前對釀造微生物菌群結構及功能研究大多采用擴增子測序技術[9-15]，例如Wang Li等[15]利用16S rRNA擴增子測解析茅臺酒發酵過程中細菌群落動態變化。然而擴增子測序得到的分子信息較為簡單[16]，難以精細解析大曲微生物群落特征，而宏基因組學測序技術則是通過獲取樣品中全部基因序列解析微生物群落信息方法，可以準確獲取大量生態信息數據。宏基因組測序已成為研究微生物多樣性和群落特征的重要手段[17]，相比于高通量測序技術，宏基因組測序技術可以更為全面挖掘出微生物基因組信息，便于研究者挖掘菌群的生物多樣性、群落結構、功能特性和相互關系[18]。現已有一些學者將宏基因組測序應用于釀造微生物群落的研究當中[18-23]，其中也可見一些高溫大曲的研究，例如Fan等[18]利用宏基因組測序技術對汾酒大曲進行全面解析，發現乳酸菌是群落中的主要物種，而釀酒酵母僅為群落總數的0.005%，Yang Yang等[19]通過宏基因組測序技術解析了中溫大曲發酵過程中微生物群的分類和功能譜。目前還鮮有學者利用宏基因組測序技術對不同輪次的高溫大曲群落進行比較研究，導致對高溫大曲微生物群落的認識的短缺。因此本研究基于宏基因組學技術對3 個輪次（下沙、造沙、1 輪次）高溫大曲群落進行精細解析，挖掘其群落結構與功能，并對不同輪次大曲間結構和功能差異進行分析，旨在為高溫大曲的品質提升及生產優化提供理論依據。

1 材料與方法

1.1 材料與試劑

大曲樣本于2021—2022年在貴州茅臺酒廠采集，參考并優化麻穎垚等[23-24]采樣方法，分別在下沙（XS）、造沙（ZS）、1輪次（YLC）取得25 個大曲樣本，其中下沙輪次7 個（采樣時間為10月5日—10月30日），造沙輪次9 個（采樣時間為11月10日—12月10日），1輪次9 個（采樣時間為12月25日—1月20日），裝入無菌自封袋中。樣品采集之后置于冰盒回實驗室后迅速放入-80 ℃冰箱保存。

氫氧化鈉、氯化鈉、鄰苯二甲酸氫鉀、可溶性淀粉、己酸、無水乙醇、葡萄糖 上海阿拉丁試劑有限公司；E.Z.N.A.?Soil DNA Kit 美國Omega BioTek公司；瓊脂糖 生工生物工程（上海）股份有限公司；FastPfuPolymerase 北京全式金生物技術有限公司；AxyPrep DNA Gel Extraction Kit 美國Axygen公司。

1.2 儀器與設備

NanoDropTM8000超微量分光光度計、VeritiPro聚合酶鏈式反應（polymerase chain reaction，PCR）儀美國Thermo Fisher Scientific公司；ExperionTM全自動電泳系統 美國Bio-Rad公司；高速臺式冷凍離心機 湖南湘儀儀器有限公司；MGISEQ-2000測序儀 深圳華大智造公司。

1.3 方法

1.3.1 大曲樣品DNA提取與檢測

取2.0 g大曲樣品置于50 mL離心管中，加入20 mL磷酸鹽緩沖液（phosphate buffered saline，PBS），渦旋混勻10 min，然后2000 r/min、4 ℃離心5 min，將上清液轉移至100 mL離心管中。收集到的上清液于12000 r/min、4 ℃離心10 min，沉淀用10 mL PBS洗滌一次，再2000 r/min、4 ℃離心10 min，收集沉淀，最后將收集的沉淀用1mL的PBS重懸備用。重懸的沉淀采用E.Z.N.A.?Soil DNA Kit進行提取，詳細的實驗步驟參照該試劑盒說明書，提取的DNA用50 μL RNA free water溶解，經1%瓊脂糖凝膠電泳和分光光度計檢測濃度及純度。按照上述步驟重復提取大曲樣品總DNA 3 次，將3 次提取的DNA樣品混合并混勻，于-80 ℃保存備用。

1.3.2 宏基因組測序

文庫構建：取1 μg大曲基因組DNA，使用超聲波儀打斷，打斷后的樣品用磁珠進行片段選擇，使得樣品條帶集中在200～400 bp左右；修復雙鏈cDNA末端，并在3’末端加上A堿基并對連接產物進行擴增，擴增產物用磁珠進行產物純化回收；將回收后的擴增產物變性為單鏈，接著在環化反應體系中充分混勻適溫反應一定時間，可得到單鏈環形產物，消化掉未被環化的線性DNA分子后，即得到最終的文庫并進行檢測。

DNBSEQ平臺測序：檢測合格文庫安排上機測序（DNBSEQ），單鏈環狀DNA分子通過滾環復制形成一個包含多個拷貝的DNA納米球，接著將得到的DNB采用高密度DNA納米芯片技術，加到芯片上的網狀小孔內，通過聯合探針錨定聚合技術進行測序（送深圳華大基因測序），測序數據已上傳至NCBI（項目編號：PRJNA884532）。

1.3.3 大曲酶活力測定

采取二硝基水楊酸法[25]，分別對大曲糖化酶與α-淀粉酶活力進行檢測。

1.3.4 生物信息分析

測序數據質控：剔除含有10%不確定堿基（N堿基）的reads；剔除含有測序接頭序列的reads（有15 個堿基或更長區域比對到接頭序列）；剔除含有50%以上低質量堿基（Q＜20堿基）的reads。

序列拼接組裝：使用Megahit拼接軟件對Clean Data進行拼接組裝，過濾掉300 bp以下的片段。

非冗余基因集的構建：將所有樣品的預測基因序列用CD-HIT軟件進行聚類，構建非冗余基因集。將樣品所有的高質量reads與非冗余基因集進行比對，統計基因在對應樣品中的豐度信息從而構建geneprofile。

基因預測：使用MetaGeneMark（http://topaz.gatech.edu/GeneMark/meta_gmhmmp.cgi）對拼接結果中的contig進行基因預測，并將其翻譯為氨基酸序列。基因預測是指通過對組裝的基因組序列進行分析，根據已知生物的基因結構知識或數據庫序列識別其所包含的基因等功能區域；編碼基因預測是識別基因組序列上所包含的蛋白質編碼區域，通過在基因組序列上尋找開放閱讀框實現。而高質量基因集Gene Catalogs則是在基因預測基礎上進行去冗余操作得到。

物種稀疏性曲線分析：為驗證本次大曲樣本數量具有足夠的物種覆蓋程度，從樣本中隨機抽取一定數量的個體，統計出這些個體所代表物種數目，并以個體數與物種數構建曲線。它可以用來比較測序數量不同的樣本物種的豐富度，也可以用來說明樣本的取樣大小是否合理。分析采用對樣本進行隨機抽樣的方法，以抽到的樣本的物種個數構建稀疏曲線，末端部分呈現上升的趨勢表明抽樣量不足，當曲線末端上升趨勢趨于平緩時則表明采樣量足夠。

物種注釋與分析：使用軟件Kraken2[26]對宏基因組數據物種注釋，獲得物種的表達豐度，分析不同樣品或者分組間的相似或差異性。宏基因組物種分析包括物種豐度的柱狀圖和熱圖分析、α多樣性指數（Shannon指數）、β多樣性分析（Bray-Curtis距離）和線性判別分析（linear discriminant analysis effect size，LEfSe）[19]。

基因功能注釋與分析：對非冗余基因使用軟件Diamond（https://github.com/bbuchfink/diamond）的BLASTP功能進行功能注釋[27]，主要注釋數據庫為京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）[28]及碳水化合物-活性酶（carbohydrate-active enzymes，CAZy）[29]，并對注釋結果進行功能差異分析，主要包括主坐標分析（principal co-ordinates analysis，PCoA）、Anosim組間（內）差異分析。

2 結果與分析

2.1 大曲宏基因組測序與組裝質量分析

本次宏基因組測序共25 個樣品，平均組裝長度為100.5 Mb，共檢測到的基因目錄為438720 個。表1是對大曲樣本測序數據進行組裝分析以及基因預測和去冗余后的結果，顯示隨著輪次的進行，大曲樣本的組裝總長、基因數及去冗余基因數呈現遞增趨勢。

表1 大曲測序數據組裝以及基因集結果統計Table 1 Statistics of sequencing data and gene assembly of microbial community in Daqu

本次宏基因測序樣本基因片段長度集中在0～2000 nt 范圍內，短片段長度基因占據大多數（圖1a）；Venn圖顯示各輪次樣本基因種類具有較大差別（圖1b）；物種稀釋曲線分析（圖1c）表明3 個輪次的樣本數量均具有足夠物種覆蓋度，說明其鑒定結果總體結果具備統計學意義，能夠代表每個輪次的物種組成。

圖1 大曲樣本基因預測分析Fig.1 Gene prediction analysis of microbial community in Daqu

2.2 大曲樣本多樣性分析

通過Kraken2軟件解析大曲在種水平上的物種結構組成，首次將多個輪次高溫大曲微生物結構進行解析，共鑒定出5104 個物種。圖2為大曲群落在屬水平上的豐度前30的堆積圖，可以發現大曲中主要以克羅彭斯特菌屬（Kroppenstedtia）和芽孢桿菌屬（Bacillus）為優勢菌種。圖3顯示豐度前30的種，3 個輪次大曲種組成基本相同，但豐度稍有差異，大曲樣本中優勢菌種主要為K.eburnea和芽孢桿菌（B.sonorensis、B.velezensis、B.licheniformis），其中豐度最高的為K.eburnea，該菌種在芝麻香型白酒高溫大曲中也有發現，在系統發育中屬于高溫放線菌科（Thermoactinomycetaceae），能夠較好地適應大曲高溫發酵[30]，在本研究中發現其為醬香型高溫大曲的主要菌種，且在樣本中占比的相似度最高。

圖2 大曲樣本微生物群落在屬水平上的組成Fig.2 Microbial composition at genus level in Daqu

圖3 大曲樣本微生物群落在種水平上的組成Fig.3 Microbial composition at species level in Daqu

通過計算Shannon指數得出各輪次α多樣性指數（圖4a），發現各輪次大曲樣品內多樣性水平基本一致。而根據豐度矩陣計算各輪次內樣本間Bray-Curtis距離衡量β多樣性（圖4b），發現造沙多樣性較其他2 個輪次高，推測該輪次作為下沙與1輪次之間的過渡所致，造沙階段（11—12月）氣候波動較大導致群落干擾程度增加，根據群落中度干擾理論[31]，大曲中不同樣本間物種多樣性會更豐富。LEfSe（圖5）可知，造沙差異微生物只有2 種，1輪次有17 種，說明造沙輪次與其余2 個輪次的群落物種差異較小，而1輪次與其余2 個輪次的群落差異較大，具有更多統計學差異的生物標識[13]。

圖4 大曲樣本多樣性Fig.4 Diversity of microbial community in Daqu samples

圖5 不同輪次生產用曲的標記微生物Fig.5 Marker microorganisms in Daqu samples

2.3 大曲宏基因功能注釋與分析

2.3.1 基因集KEGG注釋結果

將大曲基因集注釋到KEGG二級功能分類的19 條通路途徑中，全局和概述圖譜被注釋到的基因數最多，其次是碳水化合物代謝與氨基酸代謝，三者均屬于KEGG一級功能“代謝”分類下的二級代謝功能（圖6）。全局和概述圖譜的KEGG三級功能分類包括主要代謝通路，其次是碳代謝、翻譯、單向轉導、折疊、分選和降解以及其他次生代謝的生物合成途徑，包含生物生存所必須的主要代謝功能或過程，因此在大部分研究中其分析結果基因數目最多[25]。大曲原料主要為小麥，其含有較為豐富的碳水化合物與蛋白質，因此碳水化合物代謝與氨基酸代謝相關基因也會比較豐富。

圖6 KEGG二級分類注釋圖Fig.6 Secondary classification annotation in KEGG

2.3.2 基因集CAZy注釋及淀粉水解相關酶活力測定

CAZy數據庫中主要涵蓋6 大功能類：糖苷水解酶（glycoside hydrolases，GHs），糖基轉移酶（glycosyl transferases，GTs），多糖裂合酶（polysaccharide lyases，PLs），碳水化合物酯酶（carbohydrate esterases，CEs），輔助氧化還原酶（auxiliary activities，AAs）和碳水化合物結合模塊（carbohydrate-binding modules，CBMs）。如圖7a所示，在CAZy注釋6大功能中，GHs注釋結果最多，說明大曲在執行碳水化合物代謝功能中，GHs是促進多糖的水解為主要酶系。大曲中淀粉降解相關的關鍵酶主要為α-淀粉酶和糖化酶[32-33]，酶活力檢測結果顯示（圖7b），造沙大曲糖化酶和α-淀粉酶活力最高，下沙和造沙酶活力顯著高于1輪次（P＜0.05），造沙輪次GHs相關基因數目為2050，高于其余2 個輪次。

2.3.3 大曲基因功能差異分析

不同輪次大曲的主要功能基本一致（圖8），但在不同輪次間存在差異。從圖8a、d可以看出，大曲功能基因數組間組內差異基本一致，組間差異（藍色柱子）顯著大于下沙和1輪次組內差異，但與造沙相比無顯著差異。在圖8c、f顯示出相似規律，下沙和1輪次在功能聚類上有顯著差別，這可能是下沙輪次主要生產時間（10月份）與1輪次生產（1月份）處于溫濕度差異比較大的2 個階段所致。雖主要功能種類二者無異，但在功能基因數上還是有較為顯著差異（圖8b、e），造沙作為1輪次和下沙之間的過渡輪次，其組內與組間差異一致（圖8a、d），在PCoA聚類分析上造沙輪次包含其余2 個輪次，而其余2 個輪次區分度較高，在聚類上進一步反映了造沙作為過渡輪次用曲的功能特點。

圖8 功能差異分析Fig.8 Functional difference analysis

3 結論

本研究基于宏基因組技術分析高溫大曲的微生物群落結構與功能，挖掘不同輪次生產用曲群落多樣性信息與功能特征，得出如下主要結論：

1）本研究將宏基因組學技術應用于多輪次高溫大曲微生物結構及其功能解析，補充高溫大曲在種水平上的研究結果，明確高溫大曲中以K.eburnea為優勢菌種，多種芽孢桿屬細菌為亞優勢種共同構建醬香型大曲微生物群落；不同大曲物種α多樣性水平基本一致，造沙β多樣性略高于其余2 個輪次；1輪次差異物種最豐富，具有更多統計學差異的生物標識。

2）通過KEGG與CAZy 2 個通用數據庫對大曲基因功能進行注釋與分析，發現各輪次大曲主要執行代謝功能，碳水化合物代謝與氨基酸代謝為其中最主要的功能；淀粉水解相關酶系中，GHs被注釋到的基因數最多；不同輪次大曲功能種類基本一致，但功能基因數目存在差異，其中下沙與1輪次大曲存在差異，造沙與其余2 個輪次無明顯差異，可能由于造沙處于其余2 個輪次的過渡階段所致。

3）本研究只采用了生產過程中前3 個輪次的大曲樣本，發現不同輪次大曲微生物群落結構與功能組成基本一致，在實際使用過程中功能表現出現差異可能是受到生產環境干擾，該研究結果為后續大曲微生物資源保護、開發大曲判別標準及生產質控奠定了一定研究基礎。后續研究中可進一步對優勢菌株進行篩選應用，分析優勢菌株在主要代謝通路中的作用，并進一步驗證受到干擾條件下關鍵酶基因的應激表達情況，以期為大曲的品質提升、酒體品質的提升、微生物功能基因庫的挖掘等方面提供更加系統的理論依據。