米粉发酵过程中乳酸菌多样性及功能分析

熊香元，张立钊，陈力力，龚慧可，周 玥，任佑华

（湖南农业大学食品科技学院，食品科学与生物技术湖南省重点实验室，湖南 长沙 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测序分析发酵液样品中乳酸菌多样性，结果表明发酵液样品均具有较高的乳酸菌多样性，且各个时期的种群组成相似，乳杆菌属占主导地位，但还存在乳球菌属和明串珠菌属。随着发酵时间的延长，乳酸菌的多样性降低。通过功能预测分析表明优势乳酸菌属对蛋白质代谢、碳水化合物代谢、脂类代谢等具有显著影响，乳酸菌可能在发酵过程中对大米理化特性有一定的影响。本研究揭示了米粉发酵过程中乳酸菌群落结构的动态演替，为今后进一步研制米粉发酵剂、改进米粉生产工艺提供一定的研究基础。