基于Illumina MiSeq高通量技术比较甘肃藏区传统牦牛发酵乳制品细菌菌群多样性

朱 潇，梁 琪*，王湘竹，刘 瑛

（1 甘肃农业大学食品科学与工程学院 甘肃省功能乳品工程实验室 兰州730070 2 甘肃临夏州燎原乳业有限公司 甘肃临夏731100）

甘肃藏区分布于甘肃省西南部，平均海拔3 000 m，地形复杂多样，拥有独特的高原气候，强紫外线辐射（单日最高UVI 可达12.5），平均气温低（年平均气温低于5 ℃）且昼夜温差较大（近10年最高昼夜温差达到22 ℃），含氧量少（空气中的年平均氧含量仅为平原地区的50%），蕴含着丰富的微生物资源。牦牛是一种在青藏高原地区分布的珍稀物种，藏民们一般使用牦牛奶制作各种乳制品，这些乳制品是藏族家庭的主要生活来源[1]。酸牦牛乳是一种黏稠、芳香、营养丰富的发酵乳制品，藏民家庭采用的传统发酵工艺主要分为2 种：第1 种为自然发酵，不添加任何发酵剂[2]，这类酸乳可以很好地保存当地常年自然形成的优势菌种；第2 种是在少数藏区家庭，将藏灵菇粒作为发酵剂进行发酵。藏灵菇也被称为开菲尔，是一种颗粒状天然发酵剂[3-4]，含有稳定的微生物群落。此外，藏民将大量牦牛奶脱脂后经自然发酵并晾晒制成曲拉，主要提供日常消费，并且可以长期保存。制作曲拉时需要长时间将牦牛乳暴露于自然环境中，因此其含有高原环境中的多种微生物[5]。

国际上对于细菌群落组成的研究方法，主要有传统培养方法和基于分子生物学技术的非培养方法[6]。高通量技术属于非培养方法，能够在不偏向优势菌群的情况下检测总体菌群，可以更加真实、客观且全面地了解微生态系统中的细菌群落结构信息，目前已成功应用于食品微生物学研究[7-8]。近年来，随着高通量技术的发展，对于乳制品中微生物种类和丰度的研究成为国内外研究的热点。Ercolini 等[9]采用聚合酶链反应-变性梯度凝胶电泳分析蓝纹奶酪中细菌组成时发现，即使同一块奶酪，不同部位细菌种类也存在明显差异。Quigley 等[10]利用454 基因组测序平台对60 个爱尔兰软、半硬和硬奶酪的细菌菌群多样性进行分析，发现不同质地的奶酪中细菌群落组成存在显著差异。Alderet-Tapia 等[11]通过454 焦磷酸测序技术对原奶、发酵乳清、发酵乳和一种墨西哥手工自制成熟奶酪的细菌多样性进行分析，也发现制作奶酪过程中，不同类型乳制品的细菌群落存在差异性。目前，利用Illumina MiSeq 高通量测序平台分析甘肃藏区各类传统乳制品细菌菌群多样性的研究鲜见报道。

本研究对甘肃藏区酸牦牛乳、曲拉乳、藏灵菇酸奶的细菌菌群多样性进行表征。通过高通量测序技术和生物信息学比较分析，揭示细菌群落组成的差异性和相似性，对藏区益生性微生物资源的研究与发展奠定试验依据。

1 材料与方法

1.1 材料与试剂

其中藏灵菇粒在无菌脱脂牛奶中（28 ℃下培养24 h）被增活3 次，制成藏灵菇酸奶。

TruSeqTMDNA 试剂盒，美国Illumina 公司；凝胶提取试剂盒，美国Axygen Biosciences 公司；2%琼脂糖凝胶，西班牙biowest 公司；FastPfu 聚合酶，北京全式金生物；MetaVxTM文库构建试剂盒，美国GENEWIZ 公司；PCR 引物由金唯智公司设计。

1.2 仪器与设备

miniG 型离心机、漩涡混合器，德国IKA 公司；5424R 型高速台式冷冻离心机，德国eppendorf公司；Nanodrop 2000 型超微分光光度计，美国thermofisher 公司；Qubit 2.0 荧光计，德国Life 公司；tecan 型酶标仪，瑞士F200 公司；Agilent 2100生物分析仪，美国Agilent Technologies 公司；Illumina MiSeq 300PE，美国Illumina 公司；GeneAmpR9700 型PCR 仪，美国ABI 公司。

1.3 试验方法

1.3.1 DNA 提取 使用Qubit 2.0 荧光光度计得到DNA 样品的浓度，使用MetaVxTM文库构建试剂盒构建测序文库。

1.3.2 PCR 扩增及测序 以30～50 ng DNA 为模板，利用金唯智公司设计的一系列PCR 引物（包含“CCTACGGRRBGCASCAGKVRVGAAT”序列的上游引物和包含“GGACTACNVGGGTWTCT AATCC”序列的下游引物）扩增细菌的16S rDNA的V3-V4 的2 个高度可变区。PCR 扩增程序如下：95 ℃预变性2 min，95，55，72 ℃分别持续30 s共30 次循环，72 ℃延伸5 min。

使用Agilent 2100 生物分析仪检测文库质量，并且通过Qubit 2.0 荧光光度计检测文库浓度。DNA 文库混合后，按Illumina MiSeq 300PE仪器使用说明书进行PE250/300 双端测序，由MiSeq 自带的MiSeq Control Software（MCS）读取序列信息[12]。

1.3.3 分析项目 两两组装通过双端测序后的reads，保留序列长度大于200 bp 的序列，并过滤含有N 的序列。去除嵌合体后进行OUT 分析。使用Silva132[13]数据库比对16S rRNA，进行序列聚类（相似性97%）。在OUT 分析结果下，分别计算Shannon、Chao1 等α 多样性指数，并绘制稀释曲线[14]。基于unweighted unifrac 分析比较样本间是否有显著的细菌群落差异，并用PCoA（Principal co-ordinates analysis）及层次聚类显示出样本的β多样性。最后使用PICRUSt 进行功能预测分析。

为实施继续医疗教育，图书馆添购书籍增加藏书量是必不可少的举措。为了保障所添购书籍能够发挥实际作用，在添购之前就必须确定相关书籍的购买方向，避免花大量资金添购的书籍没有实质作用的尴尬局面。因此，在添购书籍之前有项必要的准备工作需要完成，就是统计医院医务人员所需要阅读、学习的医学知识，以此为基础来选择医学书籍的购买方向。随着社会科技的发展，网络已经完全融入到人们生活之中，目前年轻一代医务人员更多的倾向于网络书籍的阅读，因此图书馆也应改变以往常态，建立一个网络阅读平台，将书籍信息发布在平台上，更好的完善图书馆的供应需求，激发年轻一代医务人员学习兴趣，多多提倡网络交流，以此促进医疗事业发展[3]。

1.4 数据分析

利用cutadapt 1.9.1、Qiime 1.9.1、RDP Classifier 2.2、R 3.3.1、Cytoscape 3.5.1、STAMP 2.1.3、LEfSe 1.0 等软件进行数据处理分析。

2 结果与分析

2.1 测序数据的特点

图1是样品的PCR 扩增产物凝胶电泳图，大约在500 bp 处可见清晰、明亮的条带，且无明显的非特异扩增现象，说明PCR 扩增成功，结果满足后续测序要求。

图1 PCR 产物凝胶电泳Fig.1 Electrophoresis of PCR amplified product

通过高通量测序，9 份样品共获得1 389 098条原始序列；筛选后得到平均长度为470 bp 的511 070 条高质量16S rRNA 基因序列，表1显示了每个样本的分析读数。所有样本的覆盖估计量的观测值（高于0.999）进一步支持了以下结果。通过Shannon 和Simpson 指数对9 份样品中的细菌多样性进行评价。由表1看出，酸牦牛乳、曲拉乳、藏灵菇酸奶样品中细菌的平均Shannon 指数分别为0.171±0.100，0.014，0.337±0.105。Simpson 指数分别为：0.540±0.269，0.072，0.976±0.228。Shannon 指数越大，表示每个个体都属于不同的种属，从而多样性指数越高。因此，本研究中3 种类型传统乳制品多样性排序为：藏灵菇酸奶＞酸牦牛乳＞曲拉乳。

表1 甘肃藏区不同传统牦牛发酵乳制品α-多样性Table 1 α-Diversity of different traditional yak fermented dairy products in Tibetan areas of Gansu

2.2 稀释曲线

稀释曲线显示了在给定恒等百分比下，每个样本所获得的序列读数随识别的OTUs 数量的变化。理想情况下，通过曲线的平台期可以确定最佳覆盖范围，这表明增加读取次数并不会改变可以确定的OTUs 数量[15]。由图2可见，OTUs 的数量随着测序量的增加而不断上升，当到达一定测序深度时，各样品的稀疏曲线逐渐达到平稳，这说明随着测序量的增加不会再有较多的OTUs 和种系被发现，样品中细菌的多样性已经被充分展现。

2.3 样品细菌门水平组成

本研究采集的样品共鉴定出5 个细菌门（表2），分别为厚壁菌门（Firmicutes）、变形菌门（Proteobacteria）、拟杆菌门 （Bacteroidetes）、蓝菌门（Cyanobacteria）和绿弯菌门（Chloroflexi）。主要门为厚壁菌门，占细菌总序列的99.52%。一直以来，厚壁菌门是乳制品中研究的主要细菌[16]，该门分别占酸牦牛乳、曲拉乳和藏灵菇酸奶样本细菌序列的99.60%，99.97%和99.05%；变形杆菌是亚优势门，占总细菌序列的0.02%～0.49%。拟杆菌门、蓝菌门、绿弯菌门等小门占总细菌序列的0.08%。蓝菌门存在于样品HZ3、HZ4 及HZ6 中，而绿弯菌门仅存于HZ3 中。

图2 样品稀释曲线Fig.2 Sample dilution curve

2.4 样品细菌属水平组成

如表3所示，在属水平上分析了各样品细菌的相对丰度，鉴定出11 种细菌菌属。其中，乳杆菌属（Lactobacillus）和链球菌属（Streptococcus）来自于同一动物门（厚壁菌门）为高丰度，其次是球菌属（Macrococcus）、乳球菌属（Lactococcus）、大肠志贺菌属（Escherichia-Shigella）和拉乌尔菌属（Raoultella）。

表2 各样品门水平菌群分布相对丰度（%）Table 2 The relative abundance of horizontal flora distribution in phylum level of each sample （%）

在所有样品中均检测到了乳杆菌属和链球菌属。乳杆菌属在各样品中最为丰富，均达到90%以上，其中：酸牦牛乳为78.59%～96.19%，曲拉乳为99.37%，藏灵菇酸奶为94.71%～95.60%。在链球菌属水平下，酸牦牛乳中链球菌属含量最高（2.06%～16.74%），其次是藏灵菇酸奶（3.56%～4.22%），在曲拉乳中最低为0.58%；而对于球菌属水平，只有在酸牦牛乳中发现其存在，且在ZH5 中高达5.13%，其余2 类样品中尚未发现；未鉴定出属也具有较高的丰度，酸牦牛乳为0.01%～0.36%，曲拉乳为0.02%，藏灵菇酸奶为0.01%～0.19%。

此外，水栖菌属（Enhydrobacter）、狭义梭菌属（Clostridium_sensu_stricto_1）、尖叶菜豆属（Phaseolus_acutifolius）均在酸牦牛乳中被检测到，而在藏灵菇酸奶、曲拉乳中未被检测到或丰度值过低（＜0.05%）。

表3 各样品属水平菌群分布相对丰度（%）Table 3 The relative abundance of horizontal flora distribution in genus level of each sample （%）

2.5 样品细菌属水平丰度热图

图3是9 个样品在属水平下的物种丰度热图，图中颜色越深则说明所含的菌属越多。由图3可知，9 个样本中乳杆菌属含量最高，其次是球菌属。

其中，还有一些在数据库匹配中未鉴定出的属，也具有较高的丰度。在不同类型的传统牦牛发酵乳制品样品中，酸牦牛乳HZ4、HZ3、HZ6 及HZ1 所含菌属种类较多。曲拉乳菌属含量最少，说明每类样品中菌群的丰度存在差异。因此，不同类型的样品中菌群多样性存在差异性。

图3 属水平物种丰度热图Fig.3 The heat maps of species abundance on genus

2.6 样品细菌群落相似性分析

用R 语言分析3 类传统乳制品中的细菌群落，提取的前2 个主成分方差的累积贡献率共为97.76%，分别为95.75%及2.01%（图4）。在PC1轴上，酸牦牛乳存在部分交叠现象，说明第1 主成分是导致3 类样品细菌菌群差异的主要因素。而在PC2 维度上9 个样品均可以得到很好的区分，各样本所表示的点在空间上都有一定的距离，所产生的差异可能是由于原料、处理方式及环境因素不同所致[17]。

样本聚类分析是通过使用样本序列之间的进化信息，来比较环境样本是否在特定的进化谱系中具有显著的细菌群落差异，根据两样本间的分支长度判断，两样本越相似，样本间的分枝长度越短。从图5可知，9 份样品可分为两大类，其中6 份酸牦牛乳相似性较高，并与1 份藏灵菇酸奶聚为1 类，另1 份藏灵菇酸奶与曲拉乳聚为1类。在第1 类中，藏灵菇酸奶与酸牦牛乳产生了明显的分枝，而6 种酸牦牛乳也聚为4 组。第2类中藏灵菇酸奶与曲拉乳的分枝较长，相似性较低。以上聚类结果表明，不同种类的传统乳制品细菌群落多样性差异较大；而在不同来源的同种类乳制品中，细菌群落多样性也存在一定的差异性。这与PCoA 分析的结果一致。

2.7 样品PICRUSt 功能预测分析

PICRUSt 全称为 “Phylogenetic investigation of communities by reconstruction of unobserved states”。该软件通过被测微生物基因组全长的序列推断其共同祖先的基因功能，然后根据Greengenes 数据库推断该物种的基因功能，预测菌群代谢功能[17]。

图4 不同类型乳制品中菌群结构UniFrac非加权主坐标分析Fig.4 Unweighted UniFrac principal coordinate analysis of the bacterial communities in the different types of dairy products

图5 基于Unweighted UniFrac 距离矩阵的UPGMA 聚类分析图Fig.5 UPGMA cluster analysis graph based on Unweighted UniFrac distance matrix

如表4所示，除RNA 处理与修饰、染色质结构与动力学、核结构功能、真核细胞的胞外结构及细胞骨架缺失之外，另外20 个基因功能家族均具有数值。所有样品中，碳水化合物代谢、翻译、复制/重组/修复、基本功能预测及未知功能占比之和为49.84%～50.32%。基本功能预测占主导水平，占比为12.08%～12.18%。

表4 COG 功能分类丰度统计表（%）Table 4 Statistical table of abundance of COG functional classification （%）

3 讨论

系统掌握甘肃藏区牦牛发酵乳制品的细菌群落组成，具有一定的科学价值。首先，可以确认藏区牦牛发酵乳制品中的益生菌群，根据结果筛选合适的培养基，获得目标细菌的概率可以改善；其次，确定了不同细菌菌群在各类传统乳制品中的比例关系。

乳酸菌是来自不同动物乳中的主要菌群。采用高通量筛选方法，对牛、羊、水牛和人奶进行检测的报告，已经确定了乳酸菌reads，该序列占总序列的N40%[18]。Wu 等[19]采用传统的培养依赖法，报道发酵牦牛奶比新鲜牦牛奶含有更高水平的乳酸菌菌群，且乳酸菌为主要的细菌。牦牛由于常年生活在高原地区，它们对寒冷、低氧、强紫外线的耐受度很高，并且主要食用高原牧草，这些因素可能会影响其细菌菌群的组成。

在本研究的3 类传统牦牛发酵乳制品中，酸牦牛乳呈现出较高的物种丰度，6 种酸牦牛乳样品间也存在一定的细菌菌群差异，这是由于发酵过程中，原料、环境因素和其它外部环境条件变化导致[20]。在门水平上，酸牦牛乳的优势菌门为厚壁菌门和变形菌门；而在属水平上，乳杆菌属和链球菌属为酸牦牛乳中的高丰度，与国内外发酵动物酸奶（俄罗斯的酸牛奶[21]、意大利的酸驴乳[22]、新疆的酸羊奶[23]）的研究结果一致。其中，乳杆菌属为第1 优势菌属，丰度高达78.59%～96.19%，区别于其它动物生乳中乳球菌属为第1 优势菌[18]，这是由于乳杆菌属的耐酸性较强[23]。

Schornsteiner 等[24]发现奥地利手工硬奶酪的制作工艺和成熟条件影响了样品的微生物的分类。虽然在曲拉的生产过程中有一个加热环节，但仍有部分微生物存活了下来。因为乳清的酸化在制作曲拉的过程中尤为重要。Zhang 等[25]研究曲拉生产过程中微生物的变化，得出酸败乳清中各菌群的水平均高于乳清液，此时乳酸菌迅速成为优势菌群。这是由于酸化后的蛋白质中的营养物质比酸化前分布更加均匀，使得微生物更容易吸收。本文所采集的曲拉乳尚未进行自然环境的风干，可以将生产中的环境因素影响降低，更加客观地得出曲拉微生物的多样性。

目前，藏灵菇已被证明是一种具有功能性的食品微生物资源[26-28]，国内一些研究人员将重点放在了藏灵菇的微生物区系调查上，藏灵菇中的细菌和真菌以共生平衡的状态存在于多糖和蛋白质基质中[29-30]，在发酵过程中，藏灵菇颗粒将一定比例的菌群脱落到牛奶培养基中，随着颗粒的形成，生物量增加；发酵后，颗粒将从牛奶培养基中回收，可再次发酵。然而，大多数细菌（多达80%）属于乳酸菌属。Zhou 等[31]通过PCR-DGGE 确定了3个藏灵菇酸奶的优势菌种和酵母菌种，Gao 等[32]通过传统的培养法，在4 种藏灵菇酸奶中发现了可培养的细菌和酵母种类，随后又通过高通量测序技术揭示了这4 种藏灵菇酸奶的细菌多样性，共发现了11 种菌属，且优势菌属为乳酸杆菌属、乳球菌属及醋酸菌属，与本研究得到的优势菌属（乳酸杆菌属、链球菌属）略有不同，这是由于地域及原料奶不同所致。

通过对样品进行PICRUSt 功能预测分析发现，虽然藏区发酵乳的来源、类型和细菌多样性各不相同，但样品中各基因功能类似。藏区采用自然发酵，添加藏灵菇粒等沿用了数千年的方法制作酸牦牛乳及曲拉。这使产品具有良好的发酵特性，并极大的提高了其食用价值。

4 结论

采用高通量技术对采自甘肃藏区的传统发酵酸牦牛乳、曲拉乳、藏灵菇酸奶3 类传统乳制品进行细菌多样性分析。在所有样品中，共检测到5 个菌门、11 个细菌菌属，优势菌门均为厚壁菌门（Firmicutes）及变形菌门（Proteobacteria），优势菌属为乳杆菌属（Lactobacillus）、链球菌属（Streptococcus）及乳球菌属（Macrococcus），酸牦牛乳的细菌丰度高于曲拉乳及藏灵菇酸奶。此外基于PICRUSt 功能预测分析表明，基本功能占主导水平，这是甘肃藏区传统乳制品中的乳酸菌优良性能得以稳定延续数千年的原因之一。