豆瓣酱是我国传统发酵调味品，色泽鲜艳，香气浓郁，深受广大消费者的喜爱[1]。它的生产主要经过制曲、蚕豆发酵、辣椒发酵和混合发酵等过程[2]。发酵时，原料中的糖类和蛋白质等会被分解产生酸、醇和酯等物质，从而赋予豆瓣酱特殊的风味[3]。研究显示，豆瓣酱中近40%的微生物为细菌[4]，细菌对与食品的发酵有着重要的作用[5]，因此对于豆瓣酱样品细菌的研究就显得尤为重要。目前，对于豆瓣酱的研究主要集中在菌株的分离鉴定[6]和生产工艺的优化[7]，而对于豆瓣酱中细菌多样性和品质的研究相对较少。
随着测序技术的快速发展，以Illumina MiSeq为代表的第二代测序技术能够快速准确对发酵食品中的微生物进行解析[8-9]。目前，该技术已广泛应用于白酒[10]和酸菜[11]等食品中微生物的解析；电子舌技术可以实现同时对食品中基本味和回味的相对强度进行测定，已广泛应用于酸奶[12]和啤酒[13]等发酵食品的滋味品质评价；电子鼻技术可以对食品中特定的挥发性风味物质进行定性和定量分析，被广泛应用于酱油[14]和虾酱[15]等发酵食品中。
本研究从湖北省荆州市荆州区的农户家中采集了13 个豆瓣酱，采用Illumina MiSeq测序平台对农家豆瓣酱中细菌多样性进行分析，同时采用电子舌和电子鼻技术对相关评价指标进行测定，并将细菌多样性和豆瓣酱品质指标进行了相关性分析，最终对豆瓣酱中细菌的基因功能进行了预测，以期对后续农家豆瓣酱品质和食用安全性的提升提供数据支撑。
1 材料与方法
1.1 材料与试剂
农家豆瓣酱：采集至湖北省荆州市荆州区。
dNTPs Mix、FastPfu Fly DNA Polymerase、5×TransStartTM FastPfu Buffer，北京全式金生物技术有限公司；QIAGEN DNeasy mericon Food Kit DNA基因组提取试剂盒，德国QIAGEN公司；阴阳离子溶液、参比溶液和味觉标准溶液，日本Insent公司。
1.2 仪器与设备
Illumina MiSeq高通量测序平台，美国Illumina公司；ND-2000C微量紫外分光光度计，美国Nano Drop公司；SA 402B电子舌，日本Insent公司；PEN3电子鼻，德国Airsense公司；UVPCDS8000凝胶成像分析系统，美国Bio-Rad公司。
1.3 方法
1.3.1 样品采集
于2018年10月上旬从湖北省荆州市荆州区郢城镇和李埠镇(E112°07′～112°22′，N30°32′～30°36′)农户家中采集豆瓣酱样品13 份。其制作时间约在70～80 d左右。
1.3.2 样品微生物宏基因组DNA提取
取2.0 g豆瓣酱，按DNA基因组提取试剂盒说明书中进行DNA提取，并对其进行检验，将检验合格后的DNA样品置于-20 ℃暂存备用。
1.3.3 细菌16S rRNA V4-V5区PCR扩增及测序
正反向引物分别为338F(5′- ACTCCTACGGGAGGCAGCAG-3′)和806R(5′- GGACTACHVGGGTWTCTAAT-3′)，参照郭壮的PCR扩增参数进行扩增[16]。将扩增后的DNA产物送至上海美吉生物医药科技有限公司进行测序。测序主要流程如下：(1)文库制备，对DNA进行片段化处理，并向两端添加特定的接头来构建测序文库。接头含有的互补序列使DNA片段结合到流动槽上，进行再扩增和片段纯化。(2)测序，将文库上样到流动槽后置于测序仪中，扩增DNA片段簇。边合成边测序过程中，核苷酸通过天然的互补性与DNA模板链结合。每个核苷酸均含有荧光标记和可逆终止子，荧光信号可指示出加入的核苷酸种类，而终止子被切割后，下一个碱基继续结合，读取正向DNA链后，序列会脱落，随后重复该过程，读取反向链，最终完成测序过程。
1.3.4 生物信息学分析
使用QIIME(v1.70)平台对质控后的有效序列进行细菌物种分析和多样性评价[16]。使用GREENGENE数据库对序列进行同源性比对，并对其分类学地位进行注释，从而对豆瓣酱中细菌的基因功能进行PICRUSt预测[17]。
1.3.5 核酸登录号
本研究中所有序列数据已提交至MG-RAST数据库，登录号为mgp89184。
1.3.6 基于电子舌技术对农家豆瓣酱滋味品质的评价
参考王玉荣等方法并做适当优化[18]。样品处理：取10 g豆瓣酱样品和90 mL蒸馏水混匀后10 000 r/min离心15 min，取上清液抽滤后置于4 ℃冰箱中24 h，去除上层油层后待用。数据处理：每个样品重复测定4 次，取后3 次数据为有效数据。
1.3.7 基于电子鼻技术对农家豆瓣酱风味品质的评价
参考王玉荣等对鲊广椒风味测定方法并做适当优化[19]。样品处理：取10 g豆瓣酱样品于电子鼻样品瓶中，60 ℃保温15 min，室温静置10 min后插入电子鼻测试电极进行测定，平行测定3 次。数据处理：响应曲线在60 s后达到稳定，选取63、64和65 s的响应值，并计算其平均值为测试值。
1.3.8 多元统计学分析
采用相关性分析法对平均相对含量大于0.5%的细菌属与品质评价指标之间的相关性进行研究，选取相关系数绝对值大于0.5，且矫正后P<0.5的评价指标，采用Cytoscape软件(v3.5.1)进行相关性网络图绘制；使用Pearson相关性分析法对豆瓣酱样品中优势细菌属和COGs功能类别之间的相关性进行计算，并使用热图对结果进行可视化。使用R软件(v3.3.2)进行相关性计算和作图；使用Origin 8.5软件(OriginLab Corp，MA，USA)进行绘图。
2 结果与分析
2.1 序列丰富度和多样性分析
本研究采用Illumina MiSeq高通量测序技术对农家豆瓣酱样品中细菌多样性进行了解析。16S rRNA测序结果及各分类学地位数量信息如表1所示。
表1 样品16S rRNA测序情况及各分类地位数量
Table 1 16s rRNA read counts and number of identifiable units at different taxonomical levels
样品编号序列数/条OTU数/个门/个纲/个目/个科/个属/个超1指数香农指数DBJ167 2762 91310203154922 6887.03DBJ269 5051 4738152335471 1781.98DBJ361 8672 687132140811442 5605.48DBJ468 7743 245112239841583 2336.46DBJ550 0191 1081734501011368434.33DBJ659 9601 9712045651132201 6985.98DBJ746 3951 4652455821402271 2367.34DBJ837 6964 4801841641192573 2338.33DBJ962 1682 0006111841641 6356.24DBJ1053 4622 588102340791702 4215.96DBJ1151 9792 2872154691292502 1275.63DBJ1250 6953 9401952761403273 2237.79DBJ1353 3812 8796152649912 8626.84
由表1可知，高通量测序共产生了733 177 条高质量的16S rRNA序列，平均每个样品56 398 条。按照100%和97%相似度对序列进行划分且去除嵌合体后共得到了33 036 个OTU。DBJ8样品的超1指数和香农指数均最大，分别为3 233和8.33。由此可见，DBJ8样品细菌多样性丰度和多样性均最大。
2.2 农家豆瓣酱中细菌微生物的构成分析
在门水平上，样品中的细菌隶属于23 个门，其中硬壁菌门(Firmicutes)、变形菌门(Proteobacteria)和放线菌门(Actinobacteria)的平均相对含量分别为61.92%、21.18%和11.98%。在属水平上，共发现530 个细菌属，优势细菌属(相对含量大于1.0%的细菌属)相对含量比较分析如图1所示。
由图1可知，农家豆瓣酱中相对含量大于1.0%的细菌属共有12 个，其中隶属于硬壁菌门的有芽孢杆菌属(Bacillus，20.69%)、葡萄球菌属(Staphylococcus，20.50%)、片球菌属(Pediococcus，7.09%)、乳酸杆菌属(Lactobacillus，6.30%)和四联球菌属(Tetragenococcus，1.17%)；隶属于放线菌门的有棒状杆菌(Corynebacterium，7.42%)、短杆菌属(Brevibacterium，1.71%)和短状杆菌属(Brachybacterium，1.19%)，而隶属于变形菌门(Proteobacteria)的有肠杆菌属(Enterobacter，3.64%)、克雷伯氏菌属(Klebsiella，2.68%)、雷氏菌属(Ralstonia，1.47%)、和沙门氏菌属(Salmonella，1.03%)。由此可知，纳入本研究的农家豆瓣酱中的细菌主要是由隶属于硬壁菌门的优势细菌属构成，其平均含量累计为55.65%。值得一提的是，芽孢杆菌属、葡萄球菌属和乳酸杆菌属在13 个样本中均存在，其累计平均含量为47.49%。基于OTU的UpSet图如图2所示。
由图2可知，DBJ5与其他样品的独有OTU数量最多为631 个，其次是DBJ2和DBJ7分别为546和524，分别占其豆瓣酱样品OTU总数的56.95%、37.07%和35.77%。同时，DBJ1和DBJ13中有102个相同OTU，DBJ2和DBJ5中有58个相同OTU。由此可见，虽然不同农家豆瓣酱样品中细菌属的相对丰度存在着较大差异，但其细菌菌属的种类却有一定的相似趋势，这也说明同一地区的豆瓣酱其菌群组成的相似性[20]。
图1 农家豆瓣酱中优势细菌属相对含量的比较分析
Fig.1 Relative abundances of the major bacteria genera among home-made soy bean paste samples
图2 基于OTU的UpSet图
Fig.2 UpSet diagram based on OTU
2.3 农家豆瓣酱优势细菌属及品质的相关性分析
本研究首先使用Illumina MiSeq高通量测序技术对农家豆瓣酱中细菌多样性进行解析，再使用电子舌和电子鼻对其滋味和风味品质进行了评价，并构建了优势细菌属与品质指标相关性的网络图。农家豆瓣酱样品各滋味指标相对强度的箱型图如图3所示。
图3 农家豆瓣酱滋味指标相对强度值的箱形图(n=13)
Fig.3 The box plot of relative intensity of each taste index in home-made soy bean paste samples (n=13)
由图3可知，农家豆瓣酱在酸味上的差异最大，其次是丰度(鲜的回味)、鲜味和咸味；而在后味B(苦的回味)、苦味和后味A(涩的回味)上的差异相对较小。由此可知，酸味、鲜味和咸味可能是导致农家豆瓣酱样品滋味品质差异的主要指标。豆瓣酱风味指标强度如表2所示。
农家豆瓣酱样品的整体风味品质差异主要集中在W1C、W5S、W2S和W1S等风味指标上。由此可知，不同的农家豆瓣酱其风味品质存在较大差异，且其风味品质的差异要大于滋味品质。因豆瓣酱的滋味和风味品质均受发酵菌群的影响[21],本研究对样本中优势细菌属和感官指标间的相关性进行分析。
表2 农家豆瓣酱风味指标强度表(n=13)
Table 2 The intensity of each flavor index of home-made soy bean paste (n=13)
金属电极性能描述平均值最小值最大值中位数变异度/%W1C对芳香类物质灵敏0.120.010.310.04101.12W5S对氮氧化物灵敏49.924.17127.8454.1086.09W3C对氨气、芳香类物质灵敏0.240.070.490.1568.27W6S对氢气有选择性2.351.074.692.1554.82W5C对烷烃、芳香类物质灵敏0.320.090.640.2066.50W1S对甲烷灵敏121.6514.07306.83141.4776.01W1W对有机硫化物、萜类物质灵敏34.409.3257.0242.0249.87W2S对乙醇灵敏52.606.19128.8849.7280.37W2W对有机硫化物灵敏16.576.0729.2420.0750.04W3S对烷烃灵敏3.781.378.243.7164.17
由图4可知，10 个品质评价指标与14 个优势菌属之间存在显著相关性(|r|>0.5, P<0.05)。值得注意的是，仅Weissella(魏斯氏菌属)和Citrobacter(枸橼酸杆菌属)与风味指标呈显著正相关。鲜味、丰度、W1C、W3C和W5C均为豆瓣酱的优良型指标，而苦味、咸味、后味A和涩味则为缺陷型指标，由此可见，优势细菌属中魏斯氏菌属、枸橼酸杆菌属和假单胞菌属(Pseudomonas)的增加，雷氏菌属、盐单胞菌属(Halomonas)、棒状杆菌(Corynebactenum)、短杆菌属、肠杆菌属、克雷伯氏菌属和沙门氏菌属的降低，且根瘤菌属和草螺菌属等维持在一定比例均有助于提高农家豆瓣酱的产品品质。
图4 农家豆瓣酱优势细菌属和感官指标相关性的网络图
Fig.4 Network diagram of the correlation between the dominant genus and sensory indicators of home-made soy bean paste
注：实线表示正相关，虚线表示负相关
2.4 基于多元统计学分析农家豆瓣酱细菌构成
本研究就样品间的微生物群落结构进行研究。基于OTU水平进一步采用非加权的UniFrac距离主坐标分析和UPGMA聚类分析对13 个豆瓣酱样品的β多样性进行了研究，结果如图5所示。
图5 基于分类操作单元的非加权UniFrac距离的主坐标分析(A)和UPGMA聚类分析(B)
Fig.5 Principal coordinate analysis and UPGMA clustering analysis of OTUs based on weighted UniFrac distance
由图5-A可知，13 个农家豆瓣酱样品呈现出明显的分类趋势，编号为DBJ1、DBJ3、DBJ9和DBJ13的豆瓣酱分为一类，编号为DBJ4、DBJ6、DBJ8、DBJ10、DBJ11和DBJ12的豆瓣酱样品分为一类。由图5-B可知，13 个农家豆瓣酱样品可分为2个聚类，其中聚类一由DBJ1、DBJ2、DBJ3、DBJ9和DBJ13构成，聚类二由其他样品构成，两者结果基本一致。由此说明，尽管纳入本研究的农家豆瓣酱样品中均存在大量的核心细菌菌群，但在不同的样品间其微生物群落结构也存在较大的差异。以聚类分组为依据，对豆瓣酱样品优势细菌属、滋味指标和风味指标进行Mann-Whiney检验发现，滋味和风味指标在两者中均不存在显著性差异，这可能是由于非加权UniFrac聚类分析未将菌群含量考虑在内，但在优势细菌属上面，芽孢杆菌属、雷氏菌属、短杆菌属、棒状杆菌和短状杆菌属均具有显著性差异(P<0.05)，由此说明，样品间的微生物群落结构存在较大的差异。
2.5 PICRUSt基因预测
本研究共计注释到4195 COGs(蛋白质直系同源簇)，这些COGs分别属于23 个功能大类。以聚类分组为依据，对功能大类进行Mann-Whiney检验，其显著差异功能大类相对丰度的箱型图如图6所示。
A-RNA的加工和修饰；B-染色质结构和动力学；C-能量生产与转换；G-碳水化合物的运输和代谢；H-辅酶转运和代谢；I-脂质转运和代谢；M-细胞壁/膜/包膜生物生成；N-细胞运动；Q-次生代谢产物的生物合成，转运和分解代谢；S-功能未知；U-细胞内运输，分泌和囊泡运输(下同)
图6 显著差异功能分类相对丰度的箱型图
Fig.6 Prominent functional categories assignment
由图6可知，23 个功能大类中有10 个功能大类具有显著性差异(P<0.05)，分别为A、B、C、G、H、I、M、N、Q和U。类别C、G、H、I、Q、U和M在所用样品中占主导性地位。值得注意的是，聚类一中有更多的基因代表序列G和M，因此，聚类一中的细菌菌群对于碳水化合物的利用更加高效，细菌的生长也更加旺盛，更有利于豆瓣酱的发酵。COGs功能类群与优势细菌属之间的秩相关热图如图7所示。
由图7可知，纳入本研究的农家豆瓣酱样品中的优势细菌属与COGs功能类群之间存在着明显的相关关系。其中细胞壁/膜/包膜生物生成与芽孢杆菌属呈现极显著正相关，而与葡萄球菌属、短杆菌属、棒状杆菌属和短状杆菌属呈现显著负相关；碳水化合物的运输和代谢与克雷伯氏菌属、肠杆菌属、沙门氏菌属和片球菌属呈现显著正相关，而与短杆菌属和短状杆菌属呈现显著负相关。由此可见，克雷伯氏菌属、肠杆菌属、片球菌属和芽孢杆菌属对豆瓣酱发酵过程具有促进作用。但结合图4的结论可知，其相对丰度过高易对豆瓣酱滋味和风味品质造成影响。
图7 COGs功能类群与优势细菌属之间的秩相关热图
Fig.7 Rank correlation heat map between functional classification of COGs and dominant bacterial genera
注：黑色表示正相关，白色表示负相关；*表示显著性差异(P<0.05);**表示非常显著性差异(P<0.01)；***表示极显著性差异(P<0.001)
3 结论
本研究以荆州地区的农家豆瓣酱为研究对象，使用Illumina MiSeq高通量测序、电子舌和电子鼻技术相结合的方式对其细菌多样性进行解析，同时探讨细菌对豆瓣酱品质的影响。研究发现硬壁菌门、变形菌门和放线菌门的累计平均相对含量高达95.08%，而芽孢杆菌属、葡萄球菌属、棒状杆菌属、片球菌属和乳酸杆菌属为其优势细菌属，相对含量占细菌总数的62%。纳入本研究的农家豆瓣酱中雷氏菌属、盐单胞菌属、棒状杆菌、短杆菌属、肠杆菌属、克雷伯氏菌属和沙门氏菌属的降低可以显著改善豆瓣酱的品质。此外，优势细菌属在碳水化合物的运输与代谢、能量生产与转换、细胞壁/膜/包膜生物生成、辅酶转运和代谢和脂质转运和代谢等方面发挥着积极作用。
[1] 邓维琴,陈功,张其圣,等.不同发酵时长郫县豆瓣酱品质对比研究[J].中国调味品,2018,43(11):78-84.
[2] 张灵.不同发酵模式下豆瓣酱化学成分和酶活力研究[D].长沙：湖南农业大学,2016.
[3] 李治华,黄驰,王自鹏,等.不同后熟发酵时间郫县豆瓣酱挥发性成分分析[J].食品科学,2014,35(16):180-184.
[4] 周红丽,李莎,等.变温发酵模式下豆瓣酱自然发酵过程中细菌多样性研究[J].食品科学,2017,38(14):120-126.
[5] 董丹.郫县豆瓣发酵过程细菌多样性及动态演化研究[D].成都：西华大学,2016.
[6] NIU C,FAN Z,ZHENG F,et al.Isolation and identification of gas-producing spoilage microbes in fermented broad bean paste[J].Food Control,2018,84(2):8-16.
[7] 李红玫,刘忠华.豆瓣酱发酵工艺的优化[J].中国调味品,2012,37(6):65-70.
[8] SCHIRMER M,IJAZ U Z,D′AMORE R,et al.Insight into biases and sequencing errors for amplicon sequencing with the Illumina MiSeq platform[J].Nucleic Acids Research,2015,43(6):e37-e37.
[9] 尚雪娇,王玉荣,杨江,等.应用Illumina MiSeq高通量测序技术解析梅干菜中细菌多样性[J].中国酿造,2019,38(1):140-144.
[10] SUN W,XIAO H,PENG Q,et al.Analysis of bacterial diversity of Chinese Luzhou-flavor liquor brewed in different seasons by Illumina MiSeq sequencing[J].Annals of Microbiology,2016,66(3):1 293-1 301.
[11] DU R,GE J,ZHAO D,et al.Bacterial diversity and community structure during fermentation of Chinese sauerkraut with Lactobacillus casei 11MZ-5-1 by Illumina MiSeq sequencing[J].Letters in Applied Microbiology, 2018,66(1):55-62.
[12] WEI Z,ZHANG W,WANG Y,et al.Monitoring the fermentation, post-ripeness and storage processes of set yogurt using voltammetric electronic tongue[J].Journal of Food Engineering,2017,203(6):41-52.
[13] BLANCO C A,DE LA FUENTE R,CABALLERO I,et al.Beer discrimination using a portable electronic tongue based on screen-printed electrodes[J].Journal of Food Engineering,2015,157(7):57-62.
[14] 李德茂,陈利梅,马淑凤.电子鼻在酱油识别中的应用研究[J].中国酿造,2010,28(4):107-109.
[15] FAN Y,XUE Y,LI Z,et al.Analyzing the flavor compounds in Chinese traditional fermented shrimp pastes by HS-SPME-GC/MS and electronic nose[J].Journal of Ocean University of China,2017,16(2):311-318.
[16] 郭壮,葛东颖,尚雪娇,等.退化和正常窖泥微生物多样性的比较分析[J].食品工业科技,2018,39(22):93-98.
[17] 孙峰,田伟,张菲,等.丹江口库区库滨带植被土壤细菌群落多样性及PICRUSt功能预测分析[J].环境科学,2019,40(1):421-429.
[18] 王玉荣,张俊英,胡欣洁,等.湖北孝感和四川成都地区来源的酒曲对米酒滋味品质影响的评价[J].食品科学,2015,36(16):207-210.
[19] 王玉荣,沈馨,董蕴,等.鲊广椒细菌多样性评价及其对风味的影响[J].食品与机械,2018,34(4):25-30.
[20] JIN H,MO L,PAN L,et al.Using PacBio sequencing to investigate the bacterial microbiota of traditional Buryatian cottage cheese and comparison with Italian and Kazakhstan artisanal cheeses[J].Journal of Dairy Science,2018,101(8):6 885-6 896.
[21] KIM M J,KWAK H S,JUNG H Y,et al.Microbial communities related to sensory attributes in Korean fermented soy bean paste (doenjang)[J].Food Research International,2016,89(11):724-732.