汉中地区浆水发酵过程中微生物群落演替及其功能预测

刘 阳，王珊珊，2，3*，高 琪，甄子辰，王 令，3，4，路宏朝，2，5，张 涛，3，4

（1.陕西理工大学 生物科学与工程学院，陕西 汉中 723001；2.陕西省资源生物重点实验室，陕西 汉中 723001；3.陕南秦巴山区生物资源综合开发协同创新中心，陕西 汉中 723001；4.陕西理工大学 秦巴生物资源与生态环境省部共建国家重点实验室（培育），陕西 汉中 723001；5.秦巴特色肉制品质量提升与安全控制陕西省高校工程研究中心，陕西 汉中 723001）

浆水是一种在中国西北地区有两千多年历史的无盐发酵蔬菜[1-2]，富含益生菌、维生素C、有机酸等多种易被人体消化吸收的营养物质[3]，具有促进消化，降低胆固醇，改善肠道功能等益生作用[4-5]。研究发现，地域、原料和发酵阶段的不同是微生物结构存在差异的主要因素，这些差异同时与其安全性、益生特性和风味品质形成密切相关[6]。因此，浆水中的微生物菌群结构研究成为热点，魏本良等[6-7]的研究通过高通量测序技术揭示了陕西、甘肃多个地区来源的浆水中优势细菌为乳酸杆菌属（Lactobacillus）。LIU Z G等[8]比较了不同地区浆水和酸菜中细菌群落差异，发现酸菜中细菌多样性高于浆水，且菌群结构差异显著。周书楠等[9]对琚湾酸浆面浆水中的微生物群落分析发现，细菌以乳酸杆菌属为主；尚雪娇等[10]研究了琚湾酸浆面浆水中的真菌群落，其中子囊菌门（Ascomycota）平均相对含量高达99.68%，主要是假丝酵母属（Candida）和毕赤酵母属（Pichia）。

目前对浆水的研究多围绕陕西关中、甘肃等不同地区之间浆水中微生物结构的差异比较，但是关于浆水发酵过程中微生物演替的研究鲜有报道。传统浆水多为家庭作坊式生产，生产环境较为开放，发酵过程中有可能产生潜在致病菌或腐败菌，随着传统发酵蔬菜的生产标准逐步规范化，明确发酵过程中微生物组成及变化有利于其品质的保证和提升。此外，汉中地处陕西南部，属于北亚热带气候区，北有秦岭、南有大巴山脉两大屏障，气候温和湿润、干湿有度，与关中、甘肃等地区气候差异明显，这可能使汉中浆水具有独特的微生物群落[11]。

本研究采用Illumina MiSeq测序技术对发酵前期、中期和后期3个阶段的汉中浆水样品进行细菌16S rRNA V3-V4区与真菌内源转录间隔区1（internally transcribed spacer 1，ITS1）序列测序，解析发酵过程中微生物菌群结构动态变化并预测微生物功能分布，为深入挖掘功能微生物以及规模化生产高品质浆水提供理论依据与数据支撑。

1 材料与方法

1.1 材料与试剂

芥菜（Brassica juncea）：汉中市市售；发酵引子：2022年1月采集自陕西汉中市家庭自制（以芥菜和面汤为原料自然发酵10 d）的浆水；强力土壤脱氧核糖核酸（deoxyribonucleic acid，DNA）提取试剂盒：美国MOBioLaboratories公司；Monarch DNA胶回收试剂盒：美国NEB公司。

1.2 仪器与设备

Legend Micro 21高速离心机：德国EPPENDORF公司；G560E振荡器：美国Scientific Industries公司；9902型96 well聚合酶链式反应（polymerase chain reaction，PCR）仪：美国Applied Biosystems公司；1795型微型离心机：天根生化科技（北京）有限公司；Novaseq6000 PE250高通量测序平台：美国Illumina公司。

1.3 实验方法

1.3.1 浆水的制备

选择汉中当地新鲜芥菜5 kg，清水洗净，控干水分，开水烫菜后捞出晾凉控水，放入坛中，添加含0.5%面粉的面汤10L和200 mL发酵引子，密封后置于阴凉通风处自然发酵，每天打开瓦罐上下翻动一次，并观察发酵状态，分别在发酵第1、5、10天取样（编号为D1、D5、D10），每次取3份平行样（分别编号为D1S1、D1S2、D1S3，D5S1、D5S2、D5S3，D10S1、D10S2、D10S3）。浆水样品经纱布过滤后离心获得菌体用于测序。

1.3.2 浆水细菌和真菌Illumina MiSeq测序

使用强力土壤DNA提取试剂盒完成浆水发酵液中基因组DNA提取，并用1%琼脂糖凝胶电泳检测。用引物338F（5'-ACTCCTACGGGAGGCAGCA-3'）和806R（5'-GGACTACHVGGGTWTCTAAT-3'）、ITS1F（5'-CTTGGTCATTTAGAGGAAGTAA-3'）和ITS1R（5'-GCTGCGTTCTTCATCGATGC-3'）分别扩增细菌16S rRNA基因V3-V4区和真菌ITS1区基因。PCR扩增体系：上下游引物各0.9 μL，基因组DNA（50±10）ng，KOD FX Neo（TOYOBO）0.6 μL，KOD FX Neobuffer15μL，脱氧核糖核苷三磷酸（deoxy-ribonucleoside triphosphates，dNTP）（2 mmol/L each）6 μL，双蒸水ddH2O补至30 μL。PCR扩增条件：98 ℃预变性2 min；98 ℃变性30 s，50 ℃退火30s，72℃延伸60s，共30个循环；72 ℃再延伸5 min。将PCR产物用1.8%琼脂糖凝胶电泳检测，使用Monarch DNA胶回收试剂盒切胶回收PCR产物。高通量测序委托北京擎科生物科技有限公司完成。

1.3.3 高通量测序数据分析

将高通量测序获得的双端序列数据，经FLASH软件拼接，后用Trimmomatic软件进行质量过滤，再使用UCHIME软件去除嵌合体。通过USEARCH软件在97%的相似水平下将所有样本的有效序列聚类成操作分类单元（operational taxonomic unit，OTU）。基于Silva数据库（Relwese128）和Unite数据库（Release 7.2），采用核糖体数据库项目（ribosomal database project，RDP）Classifier（version 2.2）对OTU代表序列进行物种分类注释，并分别在各个分类水平统计样本的群落组成，利用QIIME软件生成不同分类水平上的物种分类图，再利用R语言工具绘制群落结构图。使用QIIME软件计算样本的超1（Chao 1）指数和ACE指数（菌群丰度）、香农（Shannon）指数和辛普森（Simpson）指数（菌群多样性）等α多样性指数；基于独立OTU使用Bray-Curtis加权算法进行主坐标分析（principal coordinates analysis，PCoA）比较不同样品细菌群落结构的相似性和差异程度。利用微科盟生科云平台（https：//www.bioincloud.tech）进行线性判别分析（linear discriminant analysis，LDA），研究了不同发酵时期样品中存在显著差异的类群，在属水平上采用t-test的方法两两进行比较。基于测序数据使用PICRUSt2、BugBase软件进行功能预测，Prism 8.0软件将结果可视化。实验重复3次。

2 结果与分析

2.1 样品16S rRNA和ITS基因序列质量评估及α-多样性分析

高通量测序结果显示，取自3个发酵阶段的9个浆水样品中，共产生650 688条细菌有效序列，685 317条真菌有效序列。根据97%相似性水平，共得到605个细菌OTUs和9061个真菌OTUs。稀释曲线反映了测序数据的合理性，汉中浆水发酵过程中样品的微生物稀释曲线见图1。对有效序列进行α-多样性分析，同时评估样品的物种丰度和多样性，结果见表1。由图1、表1可知，所有样品的稀释曲线逐渐平坦并且全部样品的OTU覆盖率为99.86%～100%，说明测序结果基本涵盖了样品中的微生物种类，能够准确反映全部样品中的微生物组成。OTU数量与样品中的微生物丰度呈正相关；Chao 1指数和ACE指数体现物种的丰富度信息，数值越大，多样性越高；Shannon指数和Simpson指数综合体现物种的丰富度和均匀度，数值越大，多样性越高。在整个发酵过程中细菌的丰度较为稳定，真菌在发酵后期多样性和丰富度指数达到最高值，表明通过发酵真菌群落多样性和丰富度得到提升。此外，整个发酵过程中的真菌Chao 1指数、Shannon指数和Simpson指数均大于细菌，说明浆水中真菌多样性与丰富度均大于细菌。

表1 汉中浆水发酵过程中细菌和真菌Alpha多样性分析结果Table 1 Alpha diversity analysis reesults of bacteria and fungi in Hanzhong Jiangshui during fermentation process

图1 汉中浆水样品中细菌（a）和真菌（b）的稀释曲线Fig.1 Dilution curves of bacteria and fungi in Hanzhong Jiangshui samples

2.2 不同发酵阶段菌群差异分析

对汉中浆水3个不同发酵阶段的微生物进行β-多样性分析，结果见图2。基于Bray Curtis的主坐标分析（PCoA），样本距离越接近表明样品的群落组成越相似。在细菌方面（图2a），PC1和PC2对细菌群落结构的方差贡献率分别为81.49%和10.27%，3组样本之间距离较大，而组内样本间距离较小，说明3个发酵阶段细菌群落组成存在一定差异。在真菌方面（图2b），PC1和PC2对真菌群落结构的方差贡献率分别为64.93%和14.62%。D1和D5组内样本距离较大，D10组内样本间距较小，不同发酵时间的样品存在一定差异。

图2 汉中浆水样品中细菌（a, c）和真菌（b, d）的主坐标分析和线性判别分析Fig.2 Principal coordinates analysis and linear discriminant analysis of bacteria(a,c)and fungi(b,d)in Hanzhong Jiangshui samples

线性判别分析效应大小（linear discriminant analysis effect size，LEfSe）用于发现和解释高维度生物标识，先检测并找到具有显著丰度差异特征的类群，然后采用线性判别分析（LDA）估算每个组分丰度对差异效果影响的大小。判别分析结果表明，在LDA 阈值为4 的条件下，不同阶段的样本中细菌（图2c）和真菌（图2d）分别有3个和16个分类水平存在显著性差异（P＜0.05）。D1阶段的差异显著菌属为腐败乳杆菌属（Loigolactobacillus）、明串珠菌属（Leuconostoc）；D5阶段的差异显著菌属为乳植杆菌属（Lactiplantibacillus）、被孢霉属（Mortierella）、镰刀菌属（Fusarium）；D10阶段的差异显著菌属为德巴利酵母属（Debaryomyces）。

2.3 微生物群落分布特征分析

2.3.1 细菌菌群结构特性分析

经注释，汉中浆水样品细菌分布在11个门、15个纲、24个目、35个科、54个属。在门水平，厚壁菌门（Firmicutes（99.32%～99.47%））在发酵各个阶段占据绝对优势地位（图3a）。在属水平，细菌优势微生物包括乳植杆菌属（Lactiplantbacillus）、腐败乳杆菌属（Loigolactobacillus）、魏斯氏菌属（Weissella）、寡碳乳杆菌属（Paucilactobacillus）、明串珠菌属（Leuconostoc）及发酵剂乳杆菌属（Levilactobacillus）等。乳植杆菌属相对丰度呈现先增加后降低的趋势，发酵前期（1 d）平均相对丰度为25.61%，中期（5 d）时增加至41.09%，后期（10 d）下降至34.33%；腐败乳杆菌属相对丰度从发酵前期的37.74%降至中期时的25.73%，后期时增加至30.69%。乳植杆菌属和腐败乳杆菌属在发酵过程中的相对丰度和在63.65%～66.82%呈较小波动。魏斯氏菌和寡碳乳杆菌属在整个发酵周期中含量呈现较小浮动，在各发酵阶段的平均相对丰度分别为18.44%和6.85%。明串珠菌属在整个发酵过程中呈下降趋势，最终相对丰度为3.47%。发酵剂乳杆菌属在发酵前期含量较低，随着发酵时间的延长其相对丰度从2.5%增加至4.45%（图3b）。

图3 汉中浆水样品中细菌在门（a）和属（b）水平上的相对丰度及属水平差异显著细菌的t-test检验（c）Fig.3 Relative abundance of bacteria at phylum (a) and genus (b)levels and t-test of bacteria with significant difference at genus level (c) in Hanzhong Jiangshui samples

对3个发酵时期优势菌属进一步进行t-test检验。乳植杆菌属和腐败乳杆菌属在3个阶段之间均存在极显著差异（P＜0.01）；发酵剂乳杆菌属发酵后期的相对丰度和另外两个时期之间差异极显著（P＜0.01），明串珠菌属初期时的相对丰度和另外两个时期之间差异极显著（P＜0.01）；乳球菌属（Lactococcus）在初期时显著高于后期（P＜0.05），极显著高于中期（P＜0.01）（图3c）。

2.3.2 真菌菌群结构特性分析

真菌分布在12个门、46个纲、100个目、191个科、367个属。优势真菌门为担子菌门（Basidiomycota）和子囊菌门。担子菌门在发酵前期的相对丰度为68.63%，随着发酵时间的延长其相对丰度降低，在发酵至中期时下降至31.03%，子囊菌门从前期的26.79%增加至58.21%。发酵至后期时组内样品间各菌门占比均趋于稳定，担子菌门平均相对丰度为41.02%，子囊菌门平均相对丰度51.22%。可以发现，虽然不同样品之间菌群结构存在差异，但优势菌门均为担子菌门和子囊菌门，说明担子菌门和子囊菌门是浆水发酵过程中的主要真菌门（图4a）。根据高通量测序结果，结合OTU聚类分析结果和研究需求，将真菌微生物在属水平上相对丰度最大的菌属进行统计分析。

图4 汉中浆水样品中真菌在门（a）和属（b）水平上的相对丰度及属水平差异显著真菌的t-test检验（c）Fig.4 Relative abundance of fungi at phylum (a) and genus (b)levels and t-test of fungi with significant difference at genus level (c) in Hanzhong Jiangshui samples

在属水平上，发酵前期双胞菌属（Cystofilobasidium）（29.28%）和白冬孢酵母属（Leucosporidium）（20.25%）处于优势地位。发酵中期双胞菌属（9.36%）和白冬孢酵母属（5.46%）的优势地位均下降，镰刀菌属（Fusarium）的相对丰度由初始的2.30%增加至5.99%。发酵后期德巴利酵母属（Debaryomyces）的相对丰度由发酵前期的1.01%增加至9.7%，而双胞菌属和白冬孢酵母属分别由29.28%和20.25%减少至14.98%和7.63%。样品中检测出的物种平均相对丰度＞1%的属有9个，分别是双胞菌属（14.9%）、白冬孢酵母属（7.63%）、镰刀菌属（4.26%）、德巴利酵母属（9.70%）、Holtermanniella（4.53%）、被孢霉属（Mortierella）（2.95%）、枝孢菌属（Cladosporium）（2.23%）、拟青霉属（Simplicillium）（2.16%）、曲菌属（Aspergillus）（1.93%）（图4b）。对3个发酵时期优势菌属进一步进行t-test检验，结果显示，德巴利酵母属在发酵后期时显著高于另外两个阶段（P＜0.05）（图4c）。

2.4 功能预测

在对不同发酵时期的汉中浆水核心微生物群进行解析的基础上，进一步预测菌群的潜在功能。使用BugBase基于16S rRNA序列OTU对3个阶段细菌的需氧情况、革兰氏阳性/阴性、潜在致病性、可移动原件、生物膜形成和耐受性等9个表型的生物水平覆盖情况进行了预测，结果见表2。由表2可知，细菌类群主要为革兰氏阳性菌，相对丰度为98.70～98.97%。随发酵进程好氧菌丰度上升，厌氧菌与兼性厌氧菌的丰度降低，这可能与浆水为半开放式发酵且过程中需要定期上下搅动，氧气充足有关。在微生物与宿主相互作用下，生物膜是由细菌及其所分泌的胞外多糖、蛋白质和DNA所组织的致密胞外结构，进而降低其生理活性，提高其抗生素耐药水平，因此生物膜常被认为在微生物致病中发挥重要作用[12-13]。在浆水中生物膜形成和潜在致病性的丰度均较低，说明浆水具有较好的食品安全性。

表2 汉中浆水样品发酵过程细菌表型预测的相对丰度Table 2 Predicted relative abundance of bacterial phenotypes during fermentation process of Hanzhong Jiangshui samples

运用PICRUSt2对不同发酵阶段的细菌进行功能预测。浆水发酵过程中，在京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes，KEGG）通路level1上的功能分为六大类：代谢、遗传信息处理、人类疾病、环境信息处理、细胞过程和有机系统，所占比例分别为（74.11±0.11）%、（12.07±0.19）%、（5.16±0.02）%、（3.71±0.1）%、（2.94±0.03）%、（2.03±0.05）%，其中代谢途径占比最高，有机系统最低。

不同阶段细菌在level2功能层上主要分布在44类通路，对D1与D5和D5与D10之间进行t-test检验，P阈值为0.05，差异分析见图5。D1与D5阶段间存在30个显著差异功能。D5与D10阶段间存在29个显著差异功能。随着发酵时间的延长，群落样本为适应环境变化发生了代谢功能改变，丰度变化规律主要分为两个类型：其一是在萜类化合物和聚酮类代谢、碳水化合物代谢、细胞群落-原核生物、信号转导、跨膜运输方面，相关细菌丰度在D1～D5阶段先显著增加，随后在D5～D10阶段又显著降低。其二在核苷酸代谢、折叠、分类和降解、翻译、能量代谢、复制和修复、辅因子和维生素代谢、脂质代谢代谢通路中，随着发酵时间的延长，相关细菌丰度先显著降低后显著增加。

图5 汉中浆水样品细菌功能通路KEGG功能预测分析Fig.5 KEGG function prediction analysis of bacterial functional pathway in Hanzhong Jiangshui samples

3 讨论

浆水在制作时会加入一定比例陈浆水作为引子，为浆水接入部分微生物和生物酶作为浆水发酵的启动助力。随着发酵的进行，浆水中微生物群落结构在不断变化，微生物群体出现更迭。本研究通过高通量测序技术对浆水发酵过程中微生物群落变化特征进行了分析。细菌方面，随着发酵的进行，不同发酵阶段细菌的丰度和多样性较为稳定。在门水平，厚壁菌门相对丰度均高于99%，占据绝对主导地位，在属水平各发酵阶段优势属均为乳植杆菌属、腐败乳杆菌属、魏斯氏菌属、寡碳乳杆菌属、明串珠菌属和发酵剂乳杆菌属。据此前研究报道，西北地区浆水主要菌群为乳酸杆菌属，且占据主导地位，占比达80%以上，魏斯氏菌属、乳球菌属和明串珠菌属在各地区的浆水中均有检出[6-8]。采集自西安的浆水中检测到独特的醋酸杆菌属（Acetobacter），采集自汉中的浆水含有较高丰度的片球菌属（Pediococcus，57.34%），这些未在本研究中检出或丰度较低。此前研究人员采集多地域来源的浆水，浆水样品均已发酵结束，本研究则采集了不同发酵时期的汉中浆水，检测到独特的寡碳乳杆菌属、发酵剂乳杆菌属，推测这些群落差异可能与浆水制作的地理环境和蔬菜原料不同有关。REHMAT I等[14]发现筛选自浆水的乳酸菌具有降胆固醇功能且具有较强耐酸和耐胆盐能力。研究表明乳酸菌产生多种生物活性代谢产物，包括脂肪酸、过氧化氢、有机酸、细菌素、蛋白质化合物和酚酸等，具有抗真菌和抗氧化等潜在的生物活性[15-18]。乳植杆菌属具有较强的还原能力、抑制脂质过氧化、螯合Fe2+和清除自由基活性，其细胞匀浆还具有谷胱甘肽过氧化物酶活性和超氧化物歧化酶活性[19]。魏斯氏菌主要代谢产物为乳酸和乙酸，其还可以产生寡糖和胞外多糖，一些特定的魏斯氏菌被认为是潜在的益生菌[20]。腐败乳杆菌属具有抗黄曲霉毒素产生的作用，对浆水的安全性具有重要作用[21]。

真菌方面，在门水平，以担子菌门和子囊菌门为主。双孢菌属和白冬孢酵母属在发酵前期是主要的优势真菌，所占比例分别为29.28%和20.25%，随发酵进行其丰度逐渐降低，但仍具有优势地位，同时真菌的多样性也随发酵进行而提高。镰刀菌属丰度呈现先增加后减少的趋势，可能在发酵过程中受到酸性环境或其他微生物的抑制[22-23]。德巴利酵母属在发酵后期大量出现，这是奶制品中主要的酵母菌，目前研究表明德巴利酵母属具有调节肠道微生物区系，促进细胞增殖和分化以及改善消化功能等益生功能[24]。ZHAO C等[25]对黑糯米酒发酵过程中微生物组与代谢组变化关联分析，发现双孢菌属与辛酸乙酯、癸酸乙酯、肉豆蔻酸乙酯和醋酸乙酯呈正相关，并确定了双孢菌属可以产生酯酶，催化酯类前体（醇和酸）反应产生酯，在黑糯米酒的风味形成中起着重要作用。Holtermanniella、白冬孢酵母属和被孢霉属可以代谢多种碳源积累脂类物质[26-28]。

微生物的生态功能与其群落结构息息相关。对微生物表型和基因功能进行了预测，使用BugBase对细菌表型预测，与细菌群落结构分析结果相吻合，优势细菌多为革兰氏阳性，潜在致病菌占比较低。运用PICRUSt2对细菌进行功能预测发现，与萜类化合物和聚酮类代谢、碳水化合物代谢、信号转导、跨膜运输等功能相关特异性微生物丰度在D1至D5阶段显著增加，在此期间发酵体系中以汉中特色芥菜和面汤作为底物，富含淀粉、纤维素、果胶等多糖。分析发现细菌中优势菌属多为乳酸菌，这类微生物可以产生多种水解酶，将这些多糖类底物代谢后形成成分更为丰富的发酵环境。聚酮类与微生物的防卫和细胞间沟通相关，推测与信号转导相关功能菌群的富集有一定相关性。萜类成分主要来源于植物，芥菜的加入促使了前期相关代谢菌群丰度的增加。D5至D10阶段与核苷酸代谢、能量代谢、辅因子和维生素的代谢、脂质代谢等功能相关特异性微生物丰度显著增加，与发酵底物成分的变化相关。核苷酸、维生素、脂质等物质的积累以及微生物适应环境过程中释放的信号分子促进了群落更迭。被孢霉属等在前期积累脂类过程丰度增加，到D10阶段丰度降低，此时能够代谢脂类的相关微生物则开始富集。

通过解析汉中浆水发酵过程中的菌群结构及微生物演替，发现其优势细菌多为潜在益生菌，但仍检测出部分致病菌。因此在浆水发酵过程中应注意控制发酵条件以减少致病菌污染。基于BugBase、PICRUSt2虽然可以预测浆水中已知功能细菌群落的相对丰度，但仍有大量细菌功能未知，以及对部分优势物种的研究较少，今后将加强对该方面的研究。

4 结论

本研究发现不同发酵时期差异显著的菌门主要为被孢霉门（Mortierellomycota），差异菌属主要为乳植杆菌属（Lactiplantibacillus）、明串珠菌属（Leuconostoc）、腐败乳杆菌属（Loigolactobacillus）、发酵剂乳杆菌属（Levilactobacillus）、乳球菌属（Lactococcus）、德巴利酵母属（Debaryomyces）、被孢霉属（Mortierella）和镰刀菌属（Fusarium）。功能预测显示微生物群落变化与发酵体系中底物成分息息相关，前期底物中主要是多糖等大分子，具有碳水化合物代谢功能的菌群富集，将底物代谢为核苷酸、维生素、脂质等小分子底物，同时细胞间信号传递加强；在发酵后期微生物群落为适应这些底物环境变化发生更迭。本研究通过Illumina Mi Seq测序技术对汉中浆水发酵过程中的前期、中期、后期共3个时期进行微生物多样性分析，并通过PICRUSt2预测相关功能分布解析了汉中浆水发酵过程中细菌和真菌菌群结构组成及变化规律，为汉中浆水中功能微生物挖掘奠定研究基础，也为汉中浆水的高品质生产提供理论参考。