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

王國崢，陳 筆，盧建軍，陳良強(qiáng)，張巧玲，羅 寒

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

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

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

1 材料與方法

1.1 材料與試劑

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

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

1.2 儀器與設(shè)備

NanoDropTM8000超微量分光光度計(jì)、VeritiPro聚合酶鏈?zhǔn)椒磻?yīng)（polymerase chain reaction，PCR）儀美國Thermo Fisher Scientific公司；ExperionTM全自動電泳系統(tǒng) 美國Bio-Rad公司；高速臺式冷凍離心機(jī) 湖南湘儀儀器有限公司；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，將上清液轉(zhuǎn)移至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進(jìn)行提取，詳細(xì)的實(shí)驗(yàn)步驟參照該試劑盒說明書，提取的DNA用50 μL RNA free water溶解，經(jīng)1%瓊脂糖凝膠電泳和分光光度計(jì)檢測濃度及純度。按照上述步驟重復(fù)提取大曲樣品總DNA 3 次，將3 次提取的DNA樣品混合并混勻，于-80 ℃保存?zhèn)溆谩?/p>

1.3.2 宏基因組測序

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

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

1.3.3 大曲酶活力測定

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

1.3.4 生物信息分析

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

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

非冗余基因集的構(gòu)建：將所有樣品的預(yù)測基因序列用CD-HIT軟件進(jìn)行聚類，構(gòu)建非冗余基因集。將樣品所有的高質(zhì)量reads與非冗余基因集進(jìn)行比對，統(tǒng)計(jì)基因在對應(yīng)樣品中的豐度信息從而構(gòu)建geneprofile。

基因預(yù)測：使用MetaGeneMark（http://topaz.gatech.edu/GeneMark/meta_gmhmmp.cgi）對拼接結(jié)果中的contig進(jìn)行基因預(yù)測，并將其翻譯為氨基酸序列?；蝾A(yù)測是指通過對組裝的基因組序列進(jìn)行分析，根據(jù)已知生物的基因結(jié)構(gòu)知識或數(shù)據(jù)庫序列識別其所包含的基因等功能區(qū)域；編碼基因預(yù)測是識別基因組序列上所包含的蛋白質(zhì)編碼區(qū)域，通過在基因組序列上尋找開放閱讀框?qū)崿F(xiàn)。而高質(zhì)量基因集Gene Catalogs則是在基因預(yù)測基礎(chǔ)上進(jìn)行去冗余操作得到。

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

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

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

2 結(jié)果與分析

2.1 大曲宏基因組測序與組裝質(zhì)量分析

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

表1 大曲測序數(shù)據(jù)組裝以及基因集結(jié)果統(tǒng)計(jì)Table 1 Statistics of sequencing data and gene assembly of microbial community in Daqu

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

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

2.2 大曲樣本多樣性分析

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

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

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

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

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

圖5 不同輪次生產(chǎn)用曲的標(biāo)記微生物Fig.5 Marker microorganisms in Daqu samples

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

2.3.1 基因集KEGG注釋結(jié)果

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

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

2.3.2 基因集CAZy注釋及淀粉水解相關(guān)酶活力測定

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

2.3.3 大曲基因功能差異分析

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

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

3 結(jié)論

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

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

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

3）本研究只采用了生產(chǎn)過程中前3 個輪次的大曲樣本，發(fā)現(xiàn)不同輪次大曲微生物群落結(jié)構(gòu)與功能組成基本一致，在實(shí)際使用過程中功能表現(xiàn)出現(xiàn)差異可能是受到生產(chǎn)環(huán)境干擾，該研究結(jié)果為后續(xù)大曲微生物資源保護(hù)、開發(fā)大曲判別標(biāo)準(zhǔn)及生產(chǎn)質(zhì)控奠定了一定研究基礎(chǔ)。后續(xù)研究中可進(jìn)一步對優(yōu)勢菌株進(jìn)行篩選應(yīng)用，分析優(yōu)勢菌株在主要代謝通路中的作用，并進(jìn)一步驗(yàn)證受到干擾條件下關(guān)鍵酶基因的應(yīng)激表達(dá)情況，以期為大曲的品質(zhì)提升、酒體品質(zhì)的提升、微生物功能基因庫的挖掘等方面提供更加系統(tǒng)的理論依據(jù)。