基于宏轉錄組學技術對豆醬中活菌群落分析方法的建立

安飛宇，武俊瑞，尤升波，解夢汐，陳 旭，趙 越，包海燕，烏日娜,*

（1.沈陽農(nóng)業(yè)大學食品學院，遼寧 沈陽 110866；2.山東省農(nóng)業(yè)科學院生物技術研究中心，山東 濟南 250100）

豆醬是以大豆和面粉為主要原料，經(jīng)自然發(fā)酵而成的半流動狀態(tài)發(fā)酵食品，在我國北方地區(qū)稱為大醬[1]。傳統(tǒng)自然發(fā)酵豆醬所具有的獨特色、香、味、體是在各種微生物的共同作用下，通過諸多生化反應而形成的[2]，有研究表明，細菌、霉菌和酵母菌是豆醬發(fā)酵過程中的主要微生物[3]，這些微生物對于自然發(fā)酵豆醬品質的形成起到了極其重要的作用。

近年來，自然發(fā)酵豆醬中微生物群落結構在國內(nèi)外得到廣泛的關注和研究[4-6]，同時也發(fā)展了許多研究微生物群落結構和功能的技術[7]。研究方法先后經(jīng)歷了傳統(tǒng)的平板培養(yǎng)法[8]、基于聚合酶鏈式反應-變性梯度凝膠電泳（polymerase chain reaction-denaturing gradient gel electrophoresis，PCR-DGGE）技術[9-10]、構建克隆文庫和DNA測序等[11-12]。近年來，隨著測序技術的不斷進步及測序成本的逐漸降低，科學界開始越來越多地應用高通量測序技術研究發(fā)酵食品中的微生物群落[13]。但僅在DNA層面上的測序，如16S rRNA基因測序和宏基因組測序，其結果中有很大一部分是已經(jīng)死亡的微生物，測序結果不能反映不同發(fā)酵階段菌落的真實情況，對發(fā)酵過程中實際發(fā)生的代謝反應也只能是推測[14]，因此還需使用新技術對發(fā)酵體系中的活菌群落進行分析。

宏轉錄組學是后基因組時代最具代表性的新技術。宏轉錄組指的是在某個特定條件或時空，群落中所有微生物基因轉錄本的總和，可用于衡量微生物群落宏基因組的表達水平，篩選出高表達活性功能基因和微生物[15-16]。由于其是以微生物群落的總RNA為研究對象，因此宏轉錄組測序技術在一定程度上彌補了基因組測序結果存在死亡菌種干擾的這一缺陷，其結果只反映取樣時體系中活菌的群落結構和相對豐度。近年來，宏轉錄組開始應用于如酸面團、干酪等發(fā)酵食品[17-19]，用于探索發(fā)酵食品中活菌的群落結構及功能，并通過分析各活菌群落在發(fā)酵過程中的代謝途徑揭示其對風味的貢獻，而由于自然發(fā)酵豆醬體系復雜，其總RNA的提取存在一定困難，目前鮮見有將宏轉錄組學技術應用于發(fā)酵豆醬的報道。

本實驗基于宏轉錄組學技術對自然發(fā)酵豆醬初期和末期的兩個樣本進行活菌群落檢測，建立了適用于豆醬這一復雜體系的總RNA提取方法，并對比分析了自然發(fā)酵豆醬在不同發(fā)酵時期的活菌群落結構及功能，以期為揭示豆醬發(fā)酵過程中真正發(fā)揮作用的微生物和優(yōu)良發(fā)酵菌種的選育提供基本的技術支持，也可作為其他半固態(tài)發(fā)酵食品菌群分析的技術參考。

1 材料與方法

1.1 材料與試劑

RNA PowerSoil?Total RNA Isolation Kit美國Mobio公司；Ribo-Zero Magnetic Gold Kit 美國Epicentre公司；TruSeq Stranded mRNA LT Sample Prep Kit美國Illumina公司；High Sensitivity DNA Kit美國Agilent公司。

1.2 儀器與設備

NanoDrop分析儀 美國Thermo Scientific公司；瓊脂糖凝膠電泳儀 北京市六一儀器廠；Centrifuge 5804R離心機 德國Eppendorf公司；2100生物分析儀 美國Agilent公司；HiSeq測序系統(tǒng) 美國Illumina公司。

1.3 方法

1.3.1 樣品采集

樣品取自遼寧省沈陽市遼中區(qū)一制醬經(jīng)驗豐富的農(nóng)家，分別采集發(fā)酵第0天和第56天的豆醬樣品于無菌滅酶管中，后續(xù)采用兩種不同的樣品前處理方式：1）采樣后低溫運輸回實驗室，4 ℃、8 000 r/min離心5 min，棄上清液，重復處理一次后再液氮速凍并轉移至-80 ℃冰箱保存；2）采樣時直接液氮速凍并轉移至-80 ℃冰箱保存。

1.3.2 RNA提取與質檢

采用RNA PowerSoil?Total RNA Isolation Kit，抽提總RNA，對抽提完成的RNA樣品，進行1.5%瓊脂糖凝膠電泳進行質量判斷，要求RNA條帶完整，無降解。利用紫外分光光度計對RNA進行定量。

馬思特（上海）化學有限公司是總部在美國的Master Chemical Corporation設立的金屬切削液及相關設備專業(yè)制造商，于2002年在上海注冊成立的。馬思特（上海）化學有限公司是TRIM牌金屬切削液產(chǎn)品在亞洲地區(qū)的生產(chǎn)、研發(fā)和技術服務中心，也是為機加工行業(yè)提供化學品專業(yè)管理的服務商。

1.3.3 宏轉錄組文庫構建及測序

首先將rRNA序列特異性的探針與總RNA雜交，然后利用磁珠去除rRNA探針復合物，最后再用乙醇沉淀法進一步純化mRNA。使用Illumina公司的TruSeq Stranded mRNA LT Sample Prep Kit試劑盒進行文庫的構建。包括RNA片段化、cDNA合成、cDNA文庫構建、PCR擴增已經(jīng)加上接頭的DNA片段、文庫片段選擇與純化等步驟，各步驟遵照試劑盒生產(chǎn)商的操作指導完成。

用2100生物分析儀對文庫進行質檢，合格的文庫應該有單一的峰，無接頭。然后，在Promega QuantiFluor上對文庫進行定量，合格的文庫計算后濃度應在2 nmol/L以上。最后對合格的文庫在Illumina HiSeq X-ten平臺上進行2h150 bp的雙端測序。將需要上機的文庫梯度稀釋，并按所需數(shù)據(jù)量比例混樣?；旌玫奈膸熳冃猿蓡捂溸M行上機測序。

1.4 數(shù)據(jù)處理

首先對測序下機的原始數(shù)據(jù)進行篩查和過濾。質控后對每個樣本的mRNA序列，分別采用Trinity進行從頭組裝拼接[20]。將樣本序列與NCBI-NT數(shù)據(jù)庫中的細菌、古菌、真菌和病毒序列進行BLASTN比對。運用MEGAN軟件進行物種分類注釋[21]。最后結合序列在各樣本中的豐度數(shù)據(jù)，使用QIIME軟件獲得各樣本在各個分類等級上的相對豐度分布表并繪制柱狀圖[22]。對各樣本組裝拼接后得到的所有轉錄本以0.95的相似度和0.9的最低覆蓋度進行歸并去冗余，并以最長的序列作為該Unigene的代表序列。將Unigene集與KEGG及CAZy數(shù)據(jù)庫比對，從而對轉錄本功能進行注釋分析。

2 結果與分析

2.1 總RNA瓊脂糖凝膠電泳分析

由于豆醬樣品中的大豆經(jīng)過高溫高壓蒸煮并長時間發(fā)酵，因此提取物中不含植物源RNA。此外，自然發(fā)酵豆醬體系復雜且處于開放式的發(fā)酵環(huán)境，豆醬中存在大量的核糖核酸酶，其會對總RNA的提取造成影響，使得從豆醬中提取RNA存在一定難度。因此與傳統(tǒng)動植物組織及發(fā)酵食品總RNA提取方式[23-26]不同，本實驗嘗試不同的前處理方式，通過對總RNA抽提效果的比較建立最適用于豆醬這一復雜體系的總RNA提取方法。瓊脂糖凝膠電泳檢測結果如圖1所示，不同前處理方式所獲得的RNA質量有明顯的差異。同時，不同發(fā)酵時期的總RNA提取效果也存在一定差異。

圖1 兩種不同前處理方式提取豆醬總RNA的電泳圖Fig. 1 Electropherograms of total RNA extracted from soybean paste by two different pretreatment methods

如圖1所示，經(jīng)過離心處理后的樣品能得到23S、16S條帶，條帶清晰，無明顯拖尾現(xiàn)象，且23S的亮度高于16S條帶，說明經(jīng)離心處理后提取的總RNA完整性較好，純度較高，可用于后續(xù)實驗；而未經(jīng)離心處理的樣品RNA條帶相對較暗，尤其是發(fā)酵第0天條帶拖尾嚴重，完整性較差，無法滿足宏轉錄組建庫要求。同時實驗表明，發(fā)酵第0天的豆醬電泳條帶亮度明顯弱于發(fā)酵第56天，說明發(fā)酵初期豆醬所抽提的RNA濃度可能低于發(fā)酵末期。

2.2 總RNA質量分析

表1 RNA檢測結果Table 1 Quality evaluation of total RNA

使用2100生物分析儀檢測RNA樣品的純度、質量濃度和完整性（表1）。除有明顯拖尾的H0（b）外，樣品通過離心處理所提取的RNA質量濃度較高，分別為30.45、44.68 ng/μL。該結果與瓊脂糖凝膠電泳檢測結果一致，說明發(fā)酵初期豆醬所抽提的RNA質量濃度低于發(fā)酵末期。其可能的原因為，在發(fā)酵第0天豆醬體系剛剛混入鹽水，水分含量較高，并且發(fā)酵環(huán)境的巨大改變也可能會造成RNA降解。到發(fā)酵后期，豆醬水分含量降低且豆醬體系穩(wěn)定，因此抽提相對容易。

純度較好的RNA其OD260nm/OD280nm值應介于1. 8～2.0。若OD260nm/OD280nm值小于1.8，則表明多糖、蛋白和酚類物質較多，若OD260nm/OD280nm值大于2.0，則表明RNA有部分降解[23]。本實驗各樣品所抽提RNA的OD260nm/OD280nm值在1.90～2.11之間，其中通過離心處理的樣品所抽提的RNA降解較少；而OD260nm/OD280nm值均低于1.0，說明可能仍有部分雜質殘留，但不影響后續(xù)文庫的建立。除未經(jīng)離心處理的H0樣品外，RNA完整數(shù)目（RNA integrity number，RIN）值均大于7，所提取的RNA完整性較好。檢測結果表明，經(jīng)離心處理提取的RNA質量濃度純度和總量均符合轉錄組測序要求。

2.3 測序數(shù)據(jù)統(tǒng)計及序列組裝

宏轉錄組測序涉及到擴增環(huán)節(jié)的有兩個步驟，一個是RNA反轉錄cDNA時的逆轉錄-PCR，一個是測序過程中的橋式PCR，均不具明顯的擴增偏好性。同時，為保證后續(xù)分析質量可靠，對測序下機的原始數(shù)據(jù)進行篩查和過濾，獲得可用于后續(xù)分析的有效序列集（Clean Data），并統(tǒng)計其占測序原始數(shù)據(jù)的比例。其結果如表2所示。去除低質序列后，各個樣品的Clean Data獲得率均在96%以上，Clean Reads獲得率達到99%以上，數(shù)據(jù)質量足夠滿足分析需要。

表2 測序數(shù)據(jù)統(tǒng)計Table 2 Statistical analysis of sequencing data

在得到高質序列后由Trinity軟件進行拼接組裝[19]，相應統(tǒng)計結果如表3所示，不同發(fā)酵時期豆醬分別獲得154 922 條和22 552 條序列，且拼接效果良好。

表3 組裝結果統(tǒng)計Table 3 Statistics of assembly results

2.4 宏轉錄組物種注釋和豐度分析

將樣本序列與NCBINT數(shù)據(jù)庫進行比對。使用Krona軟件進行菌群分類學組成的交互展示[27]（圖2），圖中圓圈從內(nèi)到外依次代表界、門、綱、目、科、屬6 個分類水平，扇形的大小反映了不同分類單元的相對豐度高低。在每個分類水平，各單元以不同的顏色加以區(qū)分。

如圖2所示，相對豐度最高的門類為厚壁菌門（Firmicutes）。發(fā)酵初期（H0）豆醬的活菌群落結構較為復雜，除厚壁菌門（81.6%）外，Mucoromycota（3.8%）和子囊菌門（Ascomycota，3.8%）為主要菌門。在發(fā)酵末期（H56）厚壁菌門（Firmicutes）占絕對優(yōu)勢，其相對豐度為99.0%。在屬水平上，自然發(fā)酵豆醬的主要微生物菌屬為四聯(lián)球菌屬（Tetragenococcus）、乳桿菌屬（Lactobacillus）、葡萄球菌屬（Staphylococcus）、太平洋海洋桿菌屬（Oceanobacillus）和腸球菌屬（Enterococcus），均為細菌菌屬。在發(fā)酵初期，其活菌群落結構相對復雜，除相對豐度最高的四聯(lián)球菌（Tetragenococcus，34.6%）外，還鑒定出大量的葡萄球菌屬（Staphylococcus，24.9%），但隨著發(fā)酵的進行，葡萄球菌屬相對含量迅速降低，說明其可能為制醬初期所引入的雜菌，對豆醬的發(fā)酵作用影響不大。此外，在發(fā)酵初期有一定數(shù)量的真菌被檢出，如根霉菌屬（Rhizopus，3.06%）、曲霉菌屬（Aspergilus，2.71%），但到了發(fā)酵末期，真菌菌屬豐度大幅降低，自然發(fā)酵豆醬中活菌群落結構趨于穩(wěn)定，主要的菌屬為四聯(lián)球菌屬（Tetragenococcus，70.9%）和乳桿菌屬（Lactobacillus，20.0%）。四聯(lián)球菌屬是促進健康的益生菌，它廣泛存在于豆醬、醬油、魚醬等發(fā)酵食品中，在改善豆醬風味方面具有重要作用[28]。有研究表明[29]，豆醬中的乳桿菌主要是植物乳桿菌，其不僅有利于食品的發(fā)酵，還可以改善食物的風味，且具有一定益生特性，包括維持腸道內(nèi)菌群平衡、抑制腫瘤細胞的形成、降低血清膽固醇等[30]。此外，有部分序列被注釋為病毒，但其相對含量極低，一些潛在致病菌如葡萄球菌屬（Staphylococcus）也僅在發(fā)酵初期才有較高豐度，其含量會隨著發(fā)酵時間的延長而降低，因此認為豆醬不具消費風險。本實驗證明宏轉錄技術完全可以用于自然發(fā)酵豆醬活菌群落結構分析，而且此技術具有一定的敏感性，可以顯示出豆醬發(fā)酵初期與末期間活菌群落的差異，同時也證明豆醬在自然發(fā)酵過程中微生物群落演替是一個逐漸趨于穩(wěn)定的過程。

圖2 基于Krona的分類學組成信息交互展示圖Fig. 2 Krona-based taxonomic composition information interactive display

2.5 宏轉錄組功能注釋和豐度分析

2.5.1 KEGG功能注釋

基于各Unigene在KEGG數(shù)據(jù)庫中的功能注釋結果，使用QIIME軟件，獲取各樣本對應于KEGG數(shù)據(jù)庫第1和第2功能等級的相對表達量分布，其結果如圖3所示。在第1等級的代謝通路分類中（圖3A），代謝表達量占比最高，達到總體表達量的35.9%。發(fā)酵0 d時，代謝功能占其總表達量的27.5%，后隨著發(fā)酵的進行迅速升高，發(fā)酵第56天時，其相對表達量為44.2%。

圖3 KEGG數(shù)據(jù)庫相對表達量分布統(tǒng)計圖Fig. 3 KEGG database relative expression distribution statistics at (A)fi rst and (B) second levels

進一步對代謝通路子功能進行分析，如圖3B所示，在第2等級的代謝通路分類中，碳水化合物代謝功能表達量占比最高，達到總體表達量的18.4%。在發(fā)酵第0、56天的相對表達量分別為12.0%和24.8%，這表明自然發(fā)酵豆醬菌群對淀粉具有較高的利用率。其中糖酵解/糖異生和TCA循環(huán)所注釋到的表達量最高，分別為總體表達量的4.5%和3.7%。微生物可通過糖酵解途徑將糖類分解成乙醇和二氧化碳[31]，其有助于豆醬的風味產(chǎn)生，因此在發(fā)酵后期該功能表達量升高。此外，氨基酸代謝功能表達量也相對較高，達到總體表達量的4.7%，在發(fā)酵第0天及第56天的相對表達量分別為3.4%和6.0%。該結果與其他組學研究結果一致[32-34]，說明宏轉錄技術可以用于分析自然發(fā)酵豆醬活菌群落功能。

2.5.2 CAZy功能比較分析

根據(jù)KEGG數(shù)據(jù)庫的對比結果發(fā)現(xiàn)自然發(fā)酵豆醬中的微生物活動大多與碳水化合物代謝有關，因此后續(xù)進行CAZy注釋分析?；诓煌l(fā)酵階段樣本在CAZy功能數(shù)據(jù)庫中所注釋到的相對表達量，可以比較它們在兩樣本間的倍數(shù)差異關系，從而評估各樣本所分別富集的功能類群。如圖4所示，除糖苷水解酶和糖基轉移酶在H0中有略高表達量外，其他4 個功能類群包括碳水化合物結合模塊、輔助氧化還原酶、碳水化合物酯酶以及多糖裂解酶均在H56中有較高的活性表達。

圖4 CAZy蛋白模塊比較分析圖Fig. 4 Comparative analysis of CAZy protein modules

微生物產(chǎn)生的糖基轉移酶和糖苷水解酶能生成低聚糖。相關研究表明，發(fā)酵末期有機碳和可溶性糖含量都逐漸下降并保持恒定[35]，說明糖苷水解酶及糖基轉移酶相關代謝反應已經(jīng)基本完成。此外，由于在發(fā)酵前期鹽濃度較高且含氧量較低，霉菌已經(jīng)基本停止生長，但由霉菌分泌的各種酶類仍繼續(xù)發(fā)揮作用[36]，而到了發(fā)酵末期，霉菌豐度極低，無法分泌相關酶類，這可能是發(fā)酵末期糖苷水解酶和糖基轉移酶表達活性降低的主要原因。該結果也與宏轉錄組物種注釋結果一致，再次印證宏轉錄技術可以用于分析自然發(fā)酵豆醬活菌群落結構及功能，且分析結果準確。

3 結 論

近年來，自然發(fā)酵豆醬中微生物群落結構及功能在國內(nèi)外得到廣泛的關注[4-6]，為反映自然發(fā)酵豆醬菌落的真實情況，免除傳統(tǒng)分子生物學手段所帶來的死亡菌群干擾，故嘗試使用宏轉錄組學技術對豆醬發(fā)酵體系中的活菌群落進行分析。

本實驗首先對豆醬總RNA提取的前處理方法進行了對比研究，結果表明，與傳統(tǒng)植物組織或發(fā)酵食品樣品前處理方式不同[23-26]，對豆醬樣品進行離心處理可以很好地解決復雜發(fā)酵體系對RNA提取的影響，這種樣品前處理方式對發(fā)酵初期豆醬總RNA提取的幫助尤為明顯。然后對豆醬活菌菌群結構及功能進行注釋，根據(jù)NCBI的注釋結果，以相對豐度不低于0.01為閾值，共鑒定出13 個屬的15 種微生物，其中細菌占絕大多數(shù)，主要的細菌屬為四聯(lián)球菌屬（Tetragenococcus）和乳桿菌屬（Lactobacillus）。根霉菌屬（Rhizopus）、曲霉菌屬（Aspergilus）僅在發(fā)酵初期有較高豐度，其可能是由醬醅引入的，但隨著發(fā)酵的進行，真菌豐度大幅降低，說明在醬醪發(fā)酵階段，細菌為優(yōu)勢微生物，在發(fā)酵過程中起主導作用。其物種注釋結果與先前報道可以相互印證[4,6,9-10]，證明宏轉錄技術完全可以用于自然發(fā)酵豆醬活菌群落結構分析，在顯示出豆醬發(fā)酵初期與末期間活菌群落的差異的同時，也證明豆醬在自然發(fā)酵過程中微生物群落演替是一個趨于穩(wěn)定的過程。在功能注釋方面，通過KEGG功能注釋和CAZy功能比較分析，發(fā)現(xiàn)自然發(fā)酵豆醬中微生物功能主要集中在碳水化合物和氨基酸代謝，其中在真菌豐度較高的發(fā)酵初期，糖苷水解酶和糖基轉移酶相關代謝反應活性較高，該結果也與宏轉錄組物種注釋結果一致，再次印證宏轉錄技術可以用于分析自然發(fā)酵豆醬活菌群落結構及功能，且具有良好的靈敏度及穩(wěn)定性。將該技術與其他組學技術數(shù)據(jù)整合，將有利于挖掘出更多的有益菌株和潛在的功能基因應用于未來的工業(yè)生產(chǎn)，同時還可將該實驗體系應用于其他半固態(tài)發(fā)酵食品的活菌群落分析。