临沂和周口地区中高温大曲细菌群落结构与基因功能差异性研究

马梦月，郭 壮，李学思，管桂坤，刘 宇，李富强，王玉荣,

（1.湖北文理学院湖北省食品配料工程技术研究中心，湖北襄阳 441053；2.河南牧业经济学院食品与生物工程学院，河南郑州 450046；3.山东兰陵美酒股份有限公司，山东兰陵 277731；4.河南省宋河酒业股份有限公司，河南周口 477262）

作为中国市场占有率最高的白酒，浓香型白酒因其独特的酿造工艺和绵柔醇厚的口感受到了广大消费者的喜爱[1-2]。浓香型白酒在酿造过程中使用大曲多为中温大曲或中高温大曲，其中，中高温大曲以软质小麦为原料，蕴含了丰富的微生物群和酶类，为白酒的发酵生香提供了动力[3]。王清龙等[4]采用MiSeq高通量测序技术对河南南部、东部、西部和中部地区浓香型白酒大曲微生物群落结构进行了差异研究，发现南部和西部大曲的细菌属以芽孢杆菌属（Bacillus）为主，东部和中部大曲的细菌属以青枯菌属（Ralstonia）和乳杆菌属（Lactobacillus）为主；MA 等[5]采用MiSeq 高通量测序技术对山东潍坊、江苏宿迁和四川绵阳地区浓香型大曲微生物群落结构进行了差异研究，发现山东潍坊大曲中细菌属以魏斯氏菌属（Weissella）和未分类的欧文菌科（unclassified Erwiniaceae）为主，江苏宿迁大曲中细菌属以糖多孢菌属（Saccharopolyspora）、Bacillus和Weissella为主，四川绵阳大曲细菌属以变形杆菌属（Proteus）、Lactobacillus和Weissella为主；吴树坤等[6]采用MiSeq 高通量测序技术对四川宜宾、泸州和遂宁地区浓香型大曲的微生物群落结构进行了差异研究，发现宜宾和泸州大曲中细菌属以Weisslla和高温放线菌属（Thermoactinomyces）为主，遂宁大曲中细菌属以Bacillus和Staphylococcus为主。由此可见，地域环境的不同会对大曲所含微生物的类群造成一定的影响。同处于黄淮流域超级产区的山东临沂和河南周口地区，孕育了兰陵和宋河等多数知名浓香型白酒。其中临沂市位于山东省东南部，地近黄海，该地区以沂、沭河为中心，西、北、东三面群山环抱，光照充足，雨量充沛，属温带季风气候[7]。周口市位于豫东平原，河南省东南部，该地区拥有沙颍河、涡河、西肥河和洪汝河四大水系，降水分布不均，日照、气温变幅大，属暖温带半湿润季风型气候[8]。由此可见，两地气候存在明显差异，且两地相差500 公里左右，其环境亦存在一定差异。因此，对临沂和周口地区的中高温大曲细菌类群进行比较分析，不仅可以加深研究人员对2 个地区中高温大曲中微生物资源进行全面了解，亦可以进一步论证地域环境的因素是否会对大曲细菌类群结构产生影响。

近年来，随着第二代高通量测序技术的迅速发展，众多学者不断地利用该项技术来解析酒醅[9]、酒曲[10]和窖泥[11]的微生物多样性，同时该方法也被广泛应用于解析不同香型白酒酿造过程中酒醅微生物群落结构的动态变化趋势。与传统的微生物培养方法相比，第二代测序技术具有无需培养、高通量、检测速度快和结果准确等优点，此外，该技术还能准确识别难以培养、丰度较低的微生物，全面和客观地揭示微生态系统中微生物群落结构[12]，为探索白酒大曲中微生物的多样性和功能提供了新的工具。

本研究采用Illumina MiSeq 高通量测序技术对临沂和周口地区共28 份中高温大曲细菌群落结构进行比较分析，同时基于直系同源蛋白数据库对2 个地区中高温大曲细菌潜在的功能差异进行解析，以期为黄淮流域超级产区酿酒微生物资源的挖掘提供数据支撑和理论指导。

1 材料与方法

1.1 材料与仪器

中高温大曲 分别采集自山东省临沂市兰陵县（E117°41′～118°18′，N34°37′～35°06′）和河南省周口市鹿邑县（E115°02′～115°37′，N33°43′～34°05′），样品随机选取，及时粉碎后备用。其中，临沂地区大曲样品采集13 份，周口地区大曲样品采集15 份。中高温大曲制作步骤基本如下：选取麦粒完整、饱满和无农药污染的小麦为生产原料，加入深井水对其进行润料再辊磨粉碎，粉碎后二次加入深井水进行拌料，拌料均匀后输送到压曲机制作曲坯，待曲坯成型后立即入房排列进行卧曲，再发酵培养28～30 d，曲坯培养至曲心水分基本全部蒸发，品温降至接近室温时，曲已培养成熟，即可将曲块出房入库，贮存3～6 个月。

QIAGEN DNeasy mericon Food Kit DNA 基因组提取试剂盒 德国QIAGEN 公司；5×TransStartTMFastPfu Buffer、dNTPs Mix、FastPfu Fly DNA Polymerase 北京全式金生物技术有限公司；引物338F/806R（正向引物序列338F：5’-ACTCCTACGGGAG GCAGCA-3’；反向引物序列806R：5’-GGACTACHV GGGTWTCTAAT-3’）上海桑尼生物科技有限公司；MiSeq 高通量测序配套试剂 美国Illumina 公司。

Veriti FAST PCR 仪 美国ABI 公司；DYY-12电泳仪 北京六一仪器厂；UVPCDS8000 凝胶成像分析系统 美国ProteinSimple 公司；ND-2000C 微量紫外分光光度计 美国Nano Drop 公司；PE300 MiSeq高通量测序平台 美国Illumina 公司；R930 机架式服务器 美国DELL 公司。

1.2 实验方法

1.2.1 宏基因组DNA 提取、PCR 扩增和高通量测序 按照试剂盒的操作步骤对中高温大曲中细菌的基因组DNA 进行提取，并以其为模板，使用已添加核苷酸标签（barcode）的引物338F/806R 对总DNA的16S rRNA V3～V4区域进行扩增，扩增体系与方法参照GUO 等[13]的进行，之后将检测合格的PCR 产物寄至上海美吉生物医药科技有限公司完成2×300 bp 的双末端测序。

1.2.2 生物信息学分析 本研究利用QIIME（v1.9.1）平台进行生物信息学分析。首先遵循重叠区的碱基数≥10 bp、重叠区错配率≤0.2、barcode 碱基错配或有效序列长度≥50 bp 的原则对测序返回的双端序列进行拼接和高质量序列筛选[14]，使用PyNAST进行校准并排齐序列[15]，按照100%和97%相似度进行序列划分构建分类操作单元（operational taxonomic units，OTU）[16]，使用ChimeraSlayer 去除含有嵌合体的OTU[17]，最后选取代表性序列在RDP[18]、Greengenes[19]和SILVA[20]数据库中进行同源性比对并确定其分类学地位；计算各样品中细菌类群的超1、发现物种数、香农和辛普森等α多样性指数；基于加权和非加权UniFrac 距离进行主坐标分析2 个地区大曲细菌类群的β多样性；采用PICRUSt（phylogenetic investigation of communities by reconstruction of unobserved states）软件对大曲中细菌基因功能进行预测[21]。

1.3 数据处理

使用两独立样本的非参数（Mann-Whitney）检验和多元方差分析（multivariate analysis of variance，MANOVA）解析2 个地区细菌类群的差异。使用R 软件（v4.1.3）绘制箱型图和PCoA 散点图；使用在线网站（http://jvenn.toulouse.inra.fr/app/example.html）绘制Venn 图；使用Origin 2017 软件绘制堆积柱状图；使用在线网站（http://huttenhower.sph.harvard.edu/galaxy/）绘制LEfSe 图；使用STAMP 软件绘制菌群功能差异图。

2 结果与分析

2.1 两个地区中高温大曲α 多样性比较分析

基于高通量测序技术，28 份中高温大曲共获得1477414 条高质量16S rRNA 基因序列，每份样品平均52764 条，采用100%和97%相似度进行序列划分共得到9330 个OTU。超1 指数、发现物种数和香农指数、辛普森指数反映了微生物物种的丰富度和多样性，指数越高代表物种的丰富度和多样性越高[22]。本研究首先对2 个地区中高温大曲的α多样性进行了比较分析，结果如图1 所示。

图1 两个地区中高温大曲α 多样性比较分析Fig.1 Comparative analysis of α diversity in medium-high temperature Daqu in two regions

由图1 可知，2 个地区中高温大曲的平均超1 指数分别为2310 和2625，平均发现物种数分别为1353 和1611，平均香农指数分别为5.35 和5.68，平均辛普森指数分别为0.85 和0.86，经Mann-Whitney检验分析发现，2 个地区大曲菌群在4 个α多样性指数上差异均不显著（P＞0.05）。

2.2 两个地区中高温大曲细菌群落结构分析

经序列比对发现，2 个地区中高温大曲共鉴定到40 个门和904 个属，其中，临沂地区大曲鉴定到22 个门和522 个属，周口地区大曲鉴定到39 个门和797 个属。本研究将平均相对丰度＞1.0%的细菌门和属定义为优势细菌门和属。基于门和属水平2 个地区中高温大曲细菌群落结构组成如图2 所示。

图2 基于门（a）和属（b）水平的两个地区中高温大曲中细菌群落组成分析Fig.2 Analysis of bacterial community composition in medium-high temperature Daqu in two regions based on phylum (a) and genus (b) levels

由图2a 可知，2 个地区中高温大曲中优势细菌门有4 个，分别为硬壁菌门（Firmicutes）、变形菌门（Proteobacteria）、放线菌门（Actinobacteria）和拟杆菌门（Bacteroidetes）。在临沂和周口地区大曲中，Actinobacteria 的平均相对丰度分别为1.05%和4.38%，Bacteroidetes 分别为0.81%和2.86%，经Mann-Whitney 检验分析发现两者差异显著（P＜0.05）。

由图2b 可知，2 个地区中高温大曲中平均相对含量大于1.0%的属有6 个，分别为芽孢杆菌属（Bacillus）、高温放线菌属（Thermoactinomyces）、不动杆菌属（Acinetobacter）、克罗彭斯特德菌属（Kroppenstedtia）、地杆菌属（Geobacter）和泛菌属（Pantoea）。在临沂和周口地区大曲中，Acinetobacter的平均相对丰度分别为8.08%和4.28%，Geobacter分别为0.002%和3.44%，Pantoea分别为0.43%和1.83%，经Mann-Whitney 检验分析发现差异均显著（P＜0.05）。此外，魏斯氏菌属（Weissella）和黏液乳杆菌属（Limosilactobacillus）在周口地区大曲中的平均相对丰度仅为0.34%和0.04%，而其作为优势细菌属在临沂地区大曲中平均相对丰度分别为1.74%和1.55%。有研究报道，在大曲发酵过程中Weissella和Limosilactobacillus等乳酸菌的代谢活动为发酵提供了酸性环境，有助于微生物群落的有效自我驯化，其中Weissella还参与了2,3-丁二醇的代谢，促进了大曲风味物质的积累[23]，除此之外，乳酸菌是大曲中的优势细菌[24]，在白酒发酵过程中具有促进美拉德反应和保持酿酒微生态环境等作用[25]，是白酒风味物质演替的推动者之一。短波单胞菌属（Brevundimonas）、葡萄球菌属（Staphylococcus）和金黄杆菌属（Chryseobacterium）虽为周口地区大曲中的优势细菌属，平均相对丰度分别为1.50%、1.44%和1.07%，但其在临沂地区大曲中的平均相对丰度仅为0.26%、0.18%和0.25%。ZHANG 等[26]发现Staphylococcus参与了2,3-丁二醇的代谢，且其代谢产生的水解长链脂肪酸和甘油的脂肪酶可生成游离脂肪酸和甘油，因而构成了白酒中香味物质本身或前体物质；姚亚林等[27]发现Brevundimonas具有较强水解淀粉和蛋白质的能力，在中温和高温大曲中与水分和淀粉的生成呈现正相关。亦有研究指出，Thermoactinomyces、Pantoea、Staphylococcus和Chryseobacterium在大曲开始发酵阶段只存在于空气中，之后随着发酵温度的降低才逐渐进入大曲中得到积累[28]。由此可见，2 个地区大曲细菌类群在属水平上存在明显差异，且在后续研究中积极探索周围环境因素对大曲发酵过程中细菌群落的累积与代谢是极为必要的。

本研究进一步通过LefSe 对2 个地区中高温大曲中的生物标志物进行了甄别，结果如图3 所示。图中柱长显示了线性判别值，值越大其富集程度越高[29]。

图3 LDA 线性判别分析图Fig.3 LDA linear discriminant analysis diagram

由图3 可知，在LDA 阈值得分为4.0 时，2 个地区中高温大曲中共有9 个细菌类群存在显著差异（P＜0.05），其中各有2 个、3 个、1 个、2 个和1 个分别体现在门、纲、目、科和属水平上。在属水平，Geobacter在周口地区大曲中富集，由此可见，Geobacter可作为周口地区大曲中的生物标志物；在属水平，没有细菌类群在临沂地区大曲中富集，因而临沂地区大曲中的生物标志物在属水平未能体现。

2.3 两个地区中高温大曲菌群分布及β 多样性分析

通过上述分析发现，2 个地区中高温大曲细菌类群在门属构成上存在一定的差异，因而为进一步直观地展现其群落结构异同，本研究在OTU 水平上绘制了Venn 图并基于UniFrac 距离进行了主坐标分析，结果如图4 所示。

图4 两个地区中高温大曲中细菌类群OTU 分布（a）及基于加权（b）和非加权（c）UniFrac 距离的PCoA 分析Fig.4 OTU distribution of bacterial taxa (a) and PCoA analysis based on weighted (b) and unweighted (c) UniFrac distance in medium-high temperature Daqu in two regions

由图4a 可知，本研究共划分了9330 个细菌OTU，其中2912 个OTU 为临沂特有，1212 个OTU 为周口特有，5206 个OTU 为2 个地区共有。由此可见，2 个地区大曲既具有各自独特的细菌类群，同时还含有共有的细菌群系。由图4b 可知，在考虑各OTU 序列相对丰度的情况下进行主坐标分析时，第一主成分和第二主成分的总贡献值为56.42%，2 个地区样品间虽存在明显的重叠现象，但亦具有一定的分离趋势，经多元方差分析其方差值为5.81，结果可信度大于99.99%；由图4c 可知，在仅考虑OTU 中是否含有序列的情况下进行主坐标分析时，临沂地区大曲位于X 轴的正方向，周口地区大曲位于X 轴负方向，2 个地区样品间存在明显的分离趋势，经多元方差分析其方差值为4.83，结果可信度大于99.99%。由此可见，2 个地区中高温大曲中细菌类群在整体上存在明显差异。

2.4 两个地区中高温大曲发酵菌群功能差异分析

在比较了2 个地区中高温大曲细菌菌群结构的差异后，本研究对其菌群的基因功能亦进行了差异分析，将序列上传至PICRUSt 软件进行比对后得到具有差异的功能类群如图5 所示。

图5 两个地区中高温大曲菌群功能差异分析Fig.5 Analysis of microflora functional differences in medium-high temperature Daqu in two regions

由图5 可知，相较于周口地区大曲，临沂地区大曲菌群基因功能在未知功能、转录功能、氨基酸转运与代谢功能、一般功能预测和核苷酸转运与代谢功能上的表达均显著偏高（P＜0.05），在细胞内运输功能、翻译后修饰，蛋白质周转，伴侣功能、信号转导机制功能和能量生产和转换功能上的表达均显著偏低（P＜0.05）。焦晶凯[30]的研究发现，乳酸菌的蛋白水解系统可为其生长提供所需氨基酸，结合大曲中细菌属的分析结果来看，Weissella和Limosilactobacillus在临沂地区大曲中平均相对丰度均＞1.0%，而其在周口地区大曲中平均相对丰度均不足0.4%，这可能是造成临沂地区大曲中细菌氨基酸转运与代谢功能较强的原因。由此可见，2 个地区中高温大曲的细菌菌群基因功能亦存在明显差异。

3 结论

本研究通过MiSeq 高通量测序技术对黄淮流域超级产区临沂和周口地区中高温大曲细菌群落结构和菌群基因功能进行了比较分析，结果发现，相较于周口地区大曲，临沂地区大曲中Weissella和Limosilactobacillus平均相对丰度明显偏高，Brevundimonas、Staphylococcus和Chryseobacterium平均相对丰度明显偏低，且其在转录、氨基酸和核苷转运与代谢功能表达上明显偏高；Geobacter为周口地区大曲中的生物标志物，而临沂地区大曲中的生物标志物在属水平未能体现。由此可见，临沂和周口地区中高温大曲的细菌类群结构和菌群基因功能均存在明显差异，并进一步证实了地域因素会对大曲细菌类群结构产生影响。