米粉發酵過程中乳酸菌多樣性及功能分析

熊香元，張立釗，陳力力，龔慧可，周 玥，任佑華

（湖南農業大學食品科技學院，食品科學與生物技術湖南省重點實驗室，湖南 長沙 410128）

我國南部傳統米粉主要分為發酵法和非發酵法2 種方式生產，其中發酵法是在發酵池中加入適量的“續渣液”（前次浸泡大米的浸泡水）和自來水，將整粒秈稻米浸泡2～5 d，使其在自然狀態下進行大米浸泡發酵后，再經清洗、磨漿、蒸擠、保濕、成型、蒸粉、切斷和包裝等工序制成米粉。發酵米粉較非發酵米粉口感爽滑、勁道感十足，從而提高了米粉的食用品質[1-2]。人們采用微生物傳統培養方法，對發酵米粉中的微生物進行研究，認為在大米浸泡發酵過程中起主要作用的是不同種類的乳酸菌和少量的酵母菌[3-4]，然而自然發酵過程復雜，采用傳統培養方法不僅耗時耗力，還有大量難培養或不可培養的微生物無法分離且不能對分離物進行精確鑒定，無法全面評價米粉發酵液中微生物群落多樣性。

高通量測序技術具有高通量、高靈敏度、高準確性和低運行成本的特點[5-6]，該測序技術已被應用于各個領域中[7-9]。因此，本研究采用Illumina MiSeq高通量測序技術分析米粉發酵過程中發酵液乳酸菌多樣性菌群結構的演替規律和基因功能預測，旨在全面揭示米粉發酵過程中乳酸菌群落結構的動態變化，為進一步開展自然發酵過程中主要微生物乳酸菌及其發酵作用機制的研究，改良傳統發酵米粉的生產技術奠定基礎。

1 材料與方法

1.1 材料與試劑

從湖南一鮮濕米粉生產廠采集發酵0～4 d的大米發酵液樣品5 組各3 份，分別裝入滅菌的100 mL玻璃瓶中，并且編號0d1、0d2、0d3、1d1、1d2、1d3……4d3，帶回實驗室后置于4 ℃冰箱保存，20 h內檢測處理。

細菌基因組DNA提取試劑盒、凝膠回收試劑盒、聚合酶鏈式反應（polymerase chain reaction，PCR）擴增試劑盒 北京康為世紀生物科技有限公司。

1.2 儀器與設備

SW-CJ-1FD單人單面凈化工作臺 蘇州凈化設備有限公司；FE28-PH計 梅特勒-托利多儀器有限公司；瓊脂糖凝膠電泳儀 美國Bio-Rad公司；ABI GeneAmp9700型PCR儀 愛普拜斯公司；紫外-可見分光光度計 上海元析儀器有限公司；Quanti Fluor TM-ST藍色熒光定量系統 Promega公司；MiSeq高通量測序平臺 美國Illumina公司。

1.3 方法

1.3.1 樣品細菌DNA提取

每個樣品取5 mL采用Bacteria Genomic DNA kit試劑盒方法進行發酵液中細菌基因組總DNA提取，隨后取各樣品總DNA用Tris-HCl緩沖溶液進行適當的稀釋，用分光光度計測定OD260/280nm值下DNA的濃度，取樣品DNA 2 μL加樣1%的瓊脂糖凝膠電泳檢測DNA的完整性。

1.3.2 細菌16S rDNA序列擴增和Illumina MiSeq測序

選擇質量合格的DNA作為模板，根據Illumina MiSeq測序區域（V3-V4）要求，采用合成帶有barcode的特異引物Primer 338F（ACTCCTACGGGAGGCAGCAG）及Primer 806R（GGACTACHVGGGTWTCTAAT）進行細菌16S rDNA PCR擴增，PCR擴增為20 μL反應體系。將PCR擴增達到純化要求的產物割膠回收，送往上海美吉生物工程有限公司采用MiSeq平臺測序。

1.3.3 Illumina Miseq測序數據分析

將測得原始數據通過barcode分配樣品reads，得到每個樣本的有效序列，根據PE reads之間的overlap關系，采用flash軟件將成對的reads拼接成一條序列，同時對序列質量進行質控和過濾，根據barcode區分樣品后進行操作分類單元（operational taxonomic units，OTU）聚類分析和物種分類學分析。基于OTU聚類分析結果，采用mothur軟件及菌群豐度指數計算法對OTU進行群落的物種豐度和多樣性指數分析，用于指數評估的OTU相似水平97%；基于分類學信息，在不同水平上對各樣品進行群落結構的統計分析。使用PICRUSt軟件對各組樣品中微生物的基因功能進行預測[10]，并依照COG數據庫和KEGG數據庫進行功能注釋[11]。

2 結果與分析

2.1 樣品序列分析及取樣深度驗證分析

圖1 瓊脂糖凝膠電泳檢測圖Fig.1 Agarose gel electrophoresis of bacterial genomic DNA samples

各樣品細菌基因組DNA的OD260/280nm均屬正常范圍，基因組DNA及PCR擴增產物瓊脂糖凝膠電泳檢測結果如圖1，電泳目的條帶清晰，其濃度和大小均符合下一步檢測要求。Illumina MiSeq測序15 個樣品共獲469 955 條有效序列，平均每個樣本（31 330±462） 條（表1）。在OTU相似水平97%下，最終序列平均長度420～430 bp，且測序長度上下浮動不大，說明得到的序列可以用于后續分析。圖2以樣本中隨機抽取一定數量的序列與它們所代表的物種數目構建稀釋曲線，其稀釋曲線逐漸趨于平坦，說明數據量越多對于發現新OTU邊際會越小，可以進行下一步的數據分析。

圖2 稀釋曲線Fig.2 Rarefaction curves

2.2 OTU統計及多樣性指數分析

2.2.1 OTU統計分析

基于相似性大于97%的原則，將獲得的有效序列進行聚類，共獲得278 個OTU，各樣品OTU數目為10～31 個不等，隨著時間延長逐漸減少（表1）。Venn圖可用于統計多個樣本中所共有和獨有的物種數目（如OTU），可以比較直觀地表現樣本的物種數目組成相似性及重疊情況。由圖3可知，15 個樣品有10 個共有OTU，分別占各樣品OTU數目的32.26%～100%，重疊度較高，說明各個階段發酵液的乳酸菌種類具有較大的相似性。5 組樣品中僅0 d樣品有4 個獨有OTU，這可能是大米在自然發酵過程中，發酵初期由于環境、原料和人員等因素造成發酵液中微生物菌群復雜。發酵0 d樣品與發酵1～4 d樣品共有OTU個數分別為26、23、20、12，這說明了隨著發酵時間的延長，發酵液中的微生物菌群逐漸變得穩定，某些微生物由于不適應發酵液環境生長受到抑制。

圖3 OTU維恩圖Fig.3 Venn diagram of OTUs

2.2.2 多樣性指數分析

在生態學中Chao指數和ACE指數常用來估計物種總數和群落中OTU數目，指數越大菌群豐度越高[12-13]。由表1可知，0 d樣品菌群豐度高于1～4 d樣品。Shannon指數與Simpson指數常用于反映樣品中物種均勻度。Shannon指數值越大，Simpson指數值越小，說明群落多樣性越高。由表1可知，群落多樣性排序為：0 d＞1 d＞2 d＞4 d＞3 d，并且從OTU數目可知，隨著發酵時間延長，物種種類減少。

表1 樣本信息表和α多樣性指數Table 1 Information about the samples and analysis of alpha diversity estimator

結合樣本的OTU的種類及其相對豐度，計算樣本間加權遺傳距離矩陣，利用距離矩陣繪制主坐標分析（principal coordinate analysis，PCoA）聚類圖，分析所有樣品間的相似性。由圖4可知，第1主成分的貢獻率為51.54%，第2主成分的貢獻率為22.76%。15 個樣品中，D3聚類效果最好，3 個樣品相互聚類到一起，而其他組別均出現差異。整體來看，組間樣品有明顯差別，但組內樣品的相似性大于其差異性。

圖4 PCoA聚類分析Fig.4 Cluster analysis by principal coordinate analysis

2.3 門、目、屬水平群落結構分析

本研究共注釋得到5 個門，分別為藍細菌門（Cyanobacteria）、綠彎菌門（Chloroflexi）、棲熱菌門（Deinococcus-Thermus）、擬桿菌門（Bacteroidetes）、變形菌門（Proteobacteria），發酵初期（1 d）涵蓋了所有的5 個門，而發酵結束時，僅存在厚壁菌門；另外，0～4 d所有樣品中，厚壁菌門（Firmicutes）豐度均在99%以上（99.69%～100%）。在目分類水平上，發酵初期（0 d和1 d）包含了16 個目，發酵結束時，僅存在3 個目；而且各樣品均以乳桿菌目（Lactobacillales）豐度最大（99.9%±0.34%）。門、目水平群落結構分析結果與OTU一致，隨著發酵時間延長物種種群逐漸減少。

按照Berry細菌學手冊中的生化分類法，乳酸菌可分為乳桿菌屬、鏈球菌屬、明串珠菌屬、雙歧桿菌屬和片球菌屬共5 個屬。從表2可知，雖然樣品在屬分類水平上包含了27 個屬，但主要是乳酸菌，尤其是乳桿菌屬（Lactobacillus）占主導地位，整個發酵期間相對豐度均在99%以上。另外有豐度較小的乳球菌屬，明串珠菌屬保留，其他菌屬隨著發酵時間延長消失。

表2 屬水平物種組成Table 2 Microbial community analysis at the genus level%

2.4 種水平群落結構分析

如圖5所示，15 個發酵液樣品中共鑒定出37 種乳酸菌，相對豐度大于1%的為未分類的乳酸菌（unclassifiedLactobacillus，67.75%）、能動乳桿菌（L.agilis，11.34%）、唾液乳桿菌（L.salivarius，10.54%）、發酵乳桿菌（L.fermentum，10.12%）。在所有發酵液樣品中除了未分類的乳酸菌占優勢均在60%以上之外，0 d樣品中優勢乳酸菌為發酵乳桿菌（20.20%），1 d樣品中為發酵乳桿菌（15.00%）和能動乳桿菌（11.69%），2 d樣品中為唾液乳桿菌（13.94%）和能動乳桿菌（10.00%），3 d樣品中為能動乳桿菌（9.96%）和唾液乳桿菌（9.86%），4 d樣品中為能動乳桿菌（14.28%）和唾液乳桿菌（11.71%）。根據上述結果顯示，米粉發酵液發酵過程中主要的優勢菌為能動乳桿菌、發酵乳桿菌和唾液乳桿菌，隨著發酵時間的延長，其他種類菌逐漸消失，只剩下乳酸菌作為優勢菌。

圖5 種水平群落熱圖Fig.5 Heatmap of microbial community at the species level

2.5 乳酸菌種群差異分析

圖6 樣品差異檢驗柱形圖Fig.6 Histogram of difference test among samples

根據所有樣品在種水平下物種注釋及豐度信息，選取豐度排名前15的種，進行組間差異顯著性檢驗并繪制差異檢驗柱形圖（圖6）。圖中顯示在發酵過程中發酵液的物種相似，相對豐度最大的為未分類乳桿菌（unclassifiedLactobacillus）、發酵乳桿菌、能動乳桿菌和唾液乳桿菌。其中發酵乳桿菌相對豐度差異明顯，樣品在第0天時發酵乳桿菌與其他時期相對豐度最高，發酵液樣品在第0天時棒狀乳桿菌扭曲亞種（L.coryniformissubsp.torquens）相對豐度差異顯著，第0天時相對豐度最低；第4天時羅伊氏乳桿菌（L.reuteri）相對豐度差異顯著，第4天時相對豐度最高。

2.6 樣品乳酸菌基因PICRUSt功能預測分析

2.6.1 COG功能預測

利用軟件PICRUSt進行功能基因的同源蛋白簇預測（相對豐度＞0.1%），可以通過16S測序獲得物種構成，推測樣本中功能基因的構成，根據COG數據庫的信息，從eggNOG數據庫中解析到各個COG的描述信息及其功能信息，從而得到功能預測豐度譜，繪制COG功能分類統計圖（圖7）。結果顯示不同發酵時期樣品得到的COG功能預測信息組成具有很大的相似性，獲得的預測功能基因中主要是氨基酸轉運和代謝、碳水化合物轉運和代謝、脂質轉運和代謝、輔酶轉運和代謝、能源產生和轉換等與代謝相關的基因。此外還有與翻譯，核糖體結構和起源，轉錄、復制、重組和修復，細胞壁/膜/包膜起源，無機離子運輸和代謝以及信號轉導機制等相關的功能基因。

圖7 COG功能分類統計圖Fig.7 Statistical chart of COG functional classification

2.6.2 KEGG代謝通路功能預測

KEGG是將基因、基因組信息以及更高層次的功能信息結合對基因的功能進行全面分析的數據庫。根據KEGG數據庫中得到的代謝通路Level 2豐度分析結果，15 個樣品中共有39 種代謝通路，根據Level 1可以分為7 類，其中細胞生理過程3 種、環境信息處理3 種、遺傳信息處理4 種、人類疾病6 種、組織系統7 種、未分類4 種和新陳代謝12 種。新陳代謝通路是微生物獲得營養進行生長繁殖的主要代謝途徑，同時也是增加發酵食品品質和風味的主要途徑。由Level 2水平上新陳代謝通路的代謝途徑和豐度值（表3）可知，碳水化合物代謝途徑占新陳代謝通路豐度比值分別為0 d，23.78%；1 d，24.19%；2 d，24.51%；3 d，24.79%；4 d，25.05%，各時期豐度值無顯著差異，但都呈逐漸增加趨勢。三羧酸循環、丙酮酸代謝、磷酸戊糖和糖酵解途徑等是碳水化合物的主要代謝途徑，這些代謝會產生大量的檸檬酸、乙酸、丁酸等有機酸使米粉發酵液的pH值逐漸降低，從而抑制其他雜菌生長[14]，發酵液的酸性環境能增加大米淀粉的相對結晶度和溶解度，環境pH值的變化，有利于酶催化進行生物化學反應，改變米粉的化學組成。而與食品風味物質形成相關的代謝途徑[15]（糖代謝、氨基酸代謝、脂代謝、核酸代謝）各個時期占新陳代謝通路豐度比值分別為59.01%、59.17%、59.26%、59.37%、59.47%，在不同發酵時期稍有變化。乳酸菌在糖代謝過程中除了正型乳酸發酵產生乳酸外，還會異型乳酸發酵產生乙醛、丙酮、乙偶姻、雙乙酰等多種揮發性芳香化合物，能促進米粉香味的形成；在脂代謝過程中，乳酸菌產生的酯酶使大米中的脂肪形成游離脂肪酸，有研究表明，脂肪酸對大米的食味品質有顯著影響[16]。KEGG數據庫分析15 個樣品有1 400余種酶，按照國際系統分類法得到6大類酶，其中氧化還原酶類占12.31%、轉移酶類占41.07%、水解酶類占22.71%、裂合酶類占5.62%、異構酶類占6.63%、合成酶類占11.67%。本實驗中列出了豐度值前30的酶類，如表4所示，氧化還原酶類中的L-乳酸脫氫酶能催化微生物代謝產生的乳酸，變成丙酮酸，進而參與微生物的糖代謝、脂代謝等。轉移酶類中豐度值最高的酶類是易錯DNA聚合酶，平均占比達到34.43%，易錯DNA聚合酶具有跨損傷合成DNA的能力，增強細胞對DNA損傷的耐受性，但是當它們作用于正常的DNA模板時，易產生突變，不具有校正功能[17]，這為以后微生物育種技術在發酵米粉工業生產中的應用奠定了基礎。水解酶類是一種能將大分子底物水解成小分子物質的酶，在水解酶中除生物體自身生長遺傳有關的酶外，金屬肽酶豐度值相對較高。金屬肽酶是一種能催化肽類和蛋白質中肽鍵水解的一類酶，能降解大米中的蛋白質，改變大米的化學組成，從而提升米粉品質。異構酶中的磷酸甘油酸變位酶是一種參與糖代謝過程中的關鍵酶[18]，這與KEGG代謝通路中的結果一致。

表3 KEGG 數據庫新陳代謝通路豐度值Table 3 Abundances of metabolic pathways in KEGG database

表4 KEGG數據庫主要酶豐度值Table 4 Abundances of main enzymes in KEGG database

3 討 論

不同發酵時期米粉發酵液乳酸菌多樣性分析結果表明，乳桿菌屬為各個發酵時期米粉發酵液中絕對優勢菌群，相對豐度均在99%以上。這與前人報道結果相似，Yi Cuiping等[19]利用高通量測序技術，分析米粉發酵液的微生物組成，發現微生物主要屬是乳酸桿菌屬、乳球菌和明串珠菌。整個發酵過程中乳桿菌屬占絕對優勢，可能是與米粉自然發酵開始前向浸泡水中添加的“續渣液”有關，這說明“續渣液”中主要微生物是乳桿菌屬。近年來，乳酸菌作為益生菌廣泛應用于食品及藥品領域中，含有大量乳酸菌的食品是營養與保健功能兼備的現代人的理想食品之一[20]。國內外學者利用乳酸菌強化發酵發現，經過乳酸菌發酵后米粉的直鏈淀粉、γ-氨基丁酸等含量上升，蛋白質、脂肪含量等下降，并認為是發酵過程中乳酸菌產生的酸和酶等對大米淀粉起到了純化作用[21-23]。米粉在發酵過程中主要的乳酸菌是乳桿菌屬，但也存在一些乳球菌和明串珠菌屬，乳球菌屬和乳桿菌屬均有較強的耐酸能力[24]，大部分乳球菌和乳桿菌在缺氧環境下能夠生長繁殖，這與米粉自然發酵環境和發酵液中pH值變化、霉菌等數量變化結果一致。乳球菌通常具有提高食品營養和品質的作用，明串珠菌能夠代謝產生適量的CO2氣體及多種風味化合物[25]，這一定程度上說明了發酵米粉品質優良可能是多種乳酸菌共同作用的結果。史國英等[26]對廣西傳統發酵米粉中的乳酸菌進行分離篩選，得到了4 株乳酸菌屬于戊糖片球菌（Pediococcus pentosaceus），2 株植物乳桿菌（L.plantarum）。馬霞等[27]對大米發酵液中的優勢菌株分離純化，得到1 株優勢菌株乳酸乳球菌（L.lactis)。本研究結果與之有些不同，這與樣品的取樣地區不同有關，另外本研究存在大量未分類乳桿菌。

乳酸菌種群差異分析發現，雖然在整個發酵過程中乳酸菌組成相似，但是0 d樣品棒狀乳桿菌扭曲亞種，1 d樣品發酵乳桿菌，4 d樣品羅伊氏乳桿菌與其他時期相對豐度差異顯著。因為隨著發酵時間延長其他乳桿菌數量顯著增加，發酵乳桿菌生長速率相對放緩。發酵乳桿菌為異型發酵乳酸桿菌,能夠代謝乳糖、半乳糖等多種糖類產生乳酸、乙酸、琥珀酸、乙醇等代謝產物，在發酵米粉中對米粉品質有重要的影響，發酵乳桿菌在發酵米粉中對米粉品質的影響需要進一步研究。

本研究進一步使用PICRUSt軟件對發酵液中的乳酸菌的基因功能進行了預測，PICRUSt能通過微生物群落的豐富度與數據庫對比，從而在不可觀測的情況下推測出生物群落的功能信息[28]。其中預測功能基因中主要是氨基酸轉運和代謝、碳水化合物轉運和代謝、脂質轉運和代謝、輔酶轉運和代謝、能源產生和轉換等與代謝相關的基因。COG功能分類統計圖顯示，發酵各個時期樣品的相同功能顏色有細微差別，這說明各樣品具有相同的功能，這與測定樣品OTU的結果吻合。物質代謝通路主要是糖代謝、氨基酸代謝、脂代謝、核酸代謝通路，這說明在米粉發酵過程中大米中糖的含量、蛋白質和脂肪的含量會有變化，這與相關報道[29-30]結果一致。乳酸菌主要利用葡萄糖等單糖進行糖代謝[31]，而大米主要成分淀粉是一種多糖，因此通過碳水化合物轉運和代謝通路，可以為乳酸菌的生長繁殖提供碳源，同時可以產生乳酸，說明發酵期間乳酸菌利用碳水化合主要進行自身的生長繁殖。乳酸菌在氨基酸代謝途徑中，通過乳酸菌分泌的轉氨酶、裂解酶、脫氫酶等酶的催化作用，可以產生醋酸、丁酸等有機酸，這些有機酸能與醇類反應生成脂類，形成風味物質；脂肪代謝途徑中，在酯酶的作用下可以最終生成醛類、醇類、酮類等芳香族化合物；另外乳酸菌自身能進行核苷酸代謝產生腺苷、次黃嘌呤等。由此可以推測在大米發酵過程中風味物質的形成主要與氨基酸代謝途徑有關。KEGG二級代謝通路分析表明，基因主要富集在碳水化合物代謝、氨基酸代謝、脂質代謝等新陳代謝中，這與本研究COG數據庫預測的碳水化合物、氨基酸、脂質轉運和代謝結果一致。基因功能預測表明乳酸菌代謝途徑復雜，酶的種類豐富，若需要具體了解哪些基因功能對大米的理化特性產生特定影響還需要進一步研究。

4 結 論

以不同時期的米粉發酵液樣品為研究對象，通過Illumina MiSeq測序分析發酵液樣品中乳酸菌多樣性，結果表明發酵液樣品均具有較高的乳酸菌多樣性，且各個時期的種群組成相似，乳桿菌屬占主導地位，但還存在乳球菌屬和明串珠菌屬。隨著發酵時間的延長，乳酸菌的多樣性降低。通過功能預測分析表明優勢乳酸菌屬對蛋白質代謝、碳水化合物代謝、脂類代謝等具有顯著影響，乳酸菌可能在發酵過程中對大米理化特性有一定的影響。本研究揭示了米粉發酵過程中乳酸菌群落結構的動態演替，為今后進一步研制米粉發酵劑、改進米粉生產工藝提供一定的研究基礎。