新鲜曲拉与干燥曲拉细菌群落结构的比较分析

刘振东，刘怡，李哲，毕娜，李梁，罗章，薛蓓，*

（1.西藏农牧学院食品科学学院，西藏林芝860000；2.西藏自治区藏药审评认证中心，西藏拉萨850000）

牦牛曲拉，又名奶渣，是由牦牛奶脱脂后在自然条件下发酵凝固并风干制成的粗奶酪[1-3]。作为一种传统的藏族食品，由于蛋白质含量高、益生菌种类多、保存时间长，携带方便等特点，曲拉在整个藏区广泛食用[4-5]。调查发现，藏民除食用传统的干曲拉即硬曲拉外，还有食用新鲜曲拉的习惯。由于冰箱在藏区家庭的普及，这种在室温下很难保存的湿曲拉，由于口感软糯，易于消化，老少皆宜等特点受到越来越多藏族同胞的喜爱。而目前关于新鲜曲拉的研究尚未见相关研究报道。

曲拉作为一种藏区典型的发酵型乳制品，其制作多采用传统的手工开放式自作，因此不同的制作环境直接关系到曲拉的微生物组成[6-7]，进而影响曲拉的风味、安全及营养价值[8-10]。

传统基于微生物分离、纯化、培养及鉴定的研究方法难以获得难培养及痕量微生物，难以将环境中完成的微生物的种群结构及生态关系展现出来[11-12]，高通量测序技术（high-throughput sequencing，HTS）因其不仅可以检测到那些难培养和低丰度的物种以及曾经存活或是难以分离的部分微生物，还可以排除挑选菌株时所产生的偏见性，更能代表整体的微生物组成，因此已经逐渐成为研究微生物组成及结构的一种重要方法[13]，近年该技术已广泛用于原料乳和各类发酵乳等乳制品生态位的研究[14-15]。刘怡萱等[16]通过该技术对西藏农、牧区牦牛酸奶菌群比较分析，发现厚壁菌门中乳杆菌属和乳球菌属为优势菌群；赵顺先等[17]通过该技术对新疆传统干奶酪乳酸菌多样性进行分析，发现干奶酪中乳酸菌归为2个目7个科12个属。其中乳杆菌属、双歧杆菌属、魏斯氏菌属和片球菌属是优势菌属；张敏等[18]利用该技术对新疆不同动物来源的原奶和酸奶细菌多样进行了比较分析，发现原奶样品中细菌Shannon-Wiener指数明显高于酸奶样品，原奶样品主要以变形菌门为主，而酸奶样品主要以厚壁菌门为主。本研究分别以西藏林芝地区农牧民自制的新鲜（软）曲拉和干（硬）曲拉为研究对象，利用Illumina MiSeq高通量测序技术对二者的细菌群落结构进行全面的分析，明确新鲜（软）曲拉和干（硬）曲拉各自的细菌群落组成及差异，为二者的功能菌群的挖掘和食用安全性的指导奠定一定的基础。

1 材料与方法

1.1 材料与试剂

牦牛曲拉样品：采自林芝市米林县邦仲村。

E.Z.N.A.Soil DNA试剂盒：美国OMEGA公司；Qubit2.0 DNA检测试剂盒：美国Invitrogen公司；Q5DNA聚合酶：美国Thermo Fisher公司；凝胶回收试剂盒：德国Qiagen公司。

1.2 仪器

5424R高速冷冻离心机：德国Eppendorf公司；DYY-6C型电泳仪：北京六一仪器厂；T-100 PCR仪：美国Bio-rad公司；MiSeq高通量测序仪：美国Illumina公司。

1.3 方法

1.3.1 样品的采集及制备

新鲜（软）曲拉（XQL）样品于2018年6月采自西藏林芝米林县米林镇帮仲村。将一半XQL样品置于无菌离心管中，冰盒保存并立即运送到实验室，在超净工作台上均分为 3 份（XQL1、XQL2、XQL3），置于-80冰箱中备用；将剩余一部份XQL在采集点进行传统日晒干制得到干（硬）曲拉样品（GQL），同时均分为3份（GQL1、GQL2、GQL3），置于无菌离心管中备用。

图1 曲拉样品Fig.1 Qula samples

1.3.2 DNA的提取和PCR扩增

1.3.3 文库建立

采用三代测序技术制备测序文库。第一步是把通过PCR扩增得到的产物的末端修复将黏性末端修整为平末端；在DNA序列的3’端添加A碱基以防止DNA片段自连；在DNA序列5’端添加含有文库特异性标签的测序接头；对上述连上接头的DNA片段进行PCR扩增，从而富集测序文库模板，并采用BECKMAN AMPure XP Beads再次纯化文库富集产物；最后对文库做最终的片段选择与纯化。

1.3.4 Illumina MiSeq

测序结果分析将Illumina MiSeq测序得到的是双端序列数据，进行质量控制得到以97%相似性分类OTUs。然后Blast在Gen Bank数据库进行比对，并提交到美吉生物的云分析平台进行多样性分析（http：//www.i-sanger.com/）。

2 结果与分析

2.1 测序数据统计及OTU组成分析

样品OUT聚类和多样性指数见表1。

表1 不同样本曲拉的测序数据统计Table 1 Sequence data statistics in different Qula samples

由表1可以看出XQL的测序序列数29 566.67±1 703.40，GQL 的测序序列数 33 604.67±2 042.39；XQL的 OTUs数为 981.00±14.80，GQL 的 OTUs数为 841.00±15.87。由序列数和OTUs数可知XQL和GQL中细菌种类丰富，且二者存在着明显的差异。

由表1可知XQL的Simpson指数平均值为0.90±0.01、Chao1 指数平均值为 595.31±45.23、ACE 指数平均值为602.15±49.33、Shannon指数平均值为5.39±0.17。GQL的 Simpson指数的平均值为 0.90±0.01、Chao1指数平均值为760.21±133.80、ACE指数平均值为 775.91±116.98、Shannon 指数平均值为 5.46±0.05。从多样性指数结果分析可以得出GQL细菌多样性指数相关参数均高于XQL，这表明曲拉样品在自然干燥的过程，影响了曲拉细菌的群落结构。

2.2 稀释曲线

XQL和GQL稀释曲线见图2。

图2 各样品中细菌稀释曲线Fig.2 Rarefaction curves of bacterial communities in different Qula samples

由图2可知当测序量较低时每个样品的Shannon指数均随着测序量的增加呈显著上升趋势，而当测序量达到5 000以后Shannon指数随着测序量的增加呈现缓慢平稳的趋势，稀释曲线的斜率逐渐降低；当测序量达到15 000以后，稀释曲线与X轴几乎接近平行，Shannon指数几乎达到上限。结合表1，可知XQL和GQL样品的测序序列数均显著高于15 000。结合稀释曲线可知在当前测序水平下，样品中的细菌群落多样性已能够充分覆盖。

2.3 各分类水平上的OUT数

样品各等级OUT物种统计表见表2。

表2 样品各等级OUT物种统计表Table 2 Samples of various grades OUT species statistics

由表2可知，XQL在门水平上的OUT数目8.00±1.00，在纲水平上的OUT数目9.67±1.53，在目水平上的OUT数目25.33±1.15，在科水平上的OUT数目34.00±2.65，在属水平上 OUT 数目 72.33±0.58，在种水平上OUT数目62.67±2.52。GQL在门水平上的OUT数 11.00±3.61，在纲水平上的 OUT 数目 14.33±4.93，在目水平上的OUT数目33.00±6.08，在科水平上的OUT数目 48.67±10.26，在属水平上 OUT 数目 95.67±18.18，在种水平上OUT数目71.33±9.45。

聚类分析到种水平的数据比在属水平上的少，分析原因可能是16SrDNA对亲缘关系较近的种属无法进一步区分的造成的。16SrDNA中的保守区反映不同物种间亲缘关系，而可变区则主要体现不同物质之间的差异。由于16SrDNA高度保守，对亲缘关系较近的种属分辨率不高，因此通过高通量测序技术基本上可以在属水平上准确的区分细菌类群。

2.4 样品在门水平上的分类比较

各样品门水平菌群分布相对丰度见表3。

表3 各样品门水平菌群分布相对丰度Table 3 The relative abundance of microbial communities in different Qula samples at the phylum level

样品细菌群落在门水平的分类比较，结果显示，在2组样品中共检测出6个门，分属于厚壁菌门、变形菌门、拟杆菌门、酸杆菌门、栖热菌门和螺旋体门。其中厚壁菌门为新鲜和干燥曲拉的第一优势菌门占比分别为 0.53±0.023 7，0.61±0.034 7；第二优势菌为变形菌门占比分别为 0.42±0.019 4，0.35±0.028 7；第三优势菌为拟杆菌门占比为 0.04±0.004 2，0.03±0.005 0。而酸杆菌门、栖热菌门和螺旋体门由于丰度较低，为样品的稀有菌门。

2.5 样品在属水平上的分类比较

样品细菌群落在属水平的分类比较，在6份样品中共检出25个属，见表4。

表4 各样品在属水平上的分类比较Table 4 The relative abundance of microbial communities in different Qula samples at the genus level

从属水平分析，新鲜曲拉和干燥曲拉的优势菌均为乳球菌属、沙雷氏菌属、不动杆菌属和假单胞菌属，丰度值分别为 0.50±0.023 5 和 0.58±0.040 0；0.19±0.011 5 和 0.14±0.011 4；0.07±0.004 6 和 0.04±0.004 4；0.05±0.004 0 和 0.05±0.004 7。

2.6 物种丰度热图

根据OUT的丰度数据，热图可在属水平上将不同的OUT分块聚类，并根据颜色梯度反应不同样品中细菌群落相似性、差异性及物种聚类关系。利用R语言绘制热图，可以直观地显示2种6个曲拉样品中细菌的差异，见图3。

由图3可知，上侧为样品聚类树，左侧为OUT聚类树，6个曲拉样品可分为2个大的聚类：新鲜曲拉可聚为一类，干燥曲拉样品聚为一类；物种丰富度上颜色越亮代表该菌属在该样品中的相对含量高于横向的其他样品，通过图2可知干燥曲拉样品主要优势菌上占据显著优势。

图3 各样品不同属组成和相对丰度热图Fig.3 Heatmap of the bacterial genera in different samples

2.7 样本聚类分析

各样品基于UniFrac距离的曲拉样品的聚类图见图4。

从图4可知，新鲜曲拉和干燥曲拉两类样品基本聚集到不同的类别，说明两类样品中细菌多样性存在一定的差异。

2.8 样品功能类群分布统计

通过将现有的16S rRNA基因测序数据与代谢功能已知的微生物参考基因组数据库相对比，从而实现对细菌和古菌代谢功能的预测，见图5。

由图5可知，新鲜曲拉样品和干燥曲拉样品在代谢功能上碳水化合物代谢（carbohydrate metabolism）、能量代谢（energy metabolism）、辅酶和维生素代谢（metabolism of cofactors and vitamins）、核苷酸代谢（nucleotide metabolism）为新鲜曲拉和干燥曲拉的主要代谢功能，而新鲜曲拉与干燥曲拉间的各代谢功能差异并不显著。

图4 各样品基于Weighted UniFrac距离矩阵的UPGMA（unweighted pair group method using arithmetic average）聚类分析图Fig.4 Dendrogram of cluster analysis for the different Qula samples

图5 各样品功能类群分布图Fig.5 The distribution of functional groups map for the different Qula samples

3 结论

本研究采用高通量测序技术对西藏林芝地区新鲜（软）曲拉和干燥（硬）曲拉共2组样品的细菌群落结构进行比较分析，结果从细菌多样性角度分析，干（硬）曲拉的Simpson指数的平均值为0.90±0.01、Chao1指数平均值为760.21±133.80、ACE指数平均值为775.91±116.98、Shannon指数平均值为5.46±0.05均优于而新鲜曲拉的Simpson指数平均值为0.90±0.01、Chao1指数平均值为595.31±45.23、ACE指数平均值为602.15±49.33、Shannon指数平均值为5.39±0.17，这说明干（硬）曲拉的细菌多样性更为丰富。干（硬）曲拉和新鲜（软）曲拉在门水平上和属水平上的比较发现二者优势菌门和优势菌属一致均为厚壁菌门（Firmicutes）和变形菌门（Proteobacteria），乳球菌属（Lactococcus）、沙雷氏菌属（Serratia）、不动杆菌属（Acinetobacter），但是曲拉在干燥过程中优势菌门和菌属的相对丰度发生了明显的变化，其中厚壁菌门由0.53±0.023 7增至0.61±0.034 7；变形菌门由 0.42±0.019 4 降至 0.35±0.028 7；乳球菌属（Lactococcus）由 0.50±0.023 5 增至0.58±0.040 0；沙雷氏菌属（Serratia）由 0.19±0.011 5 降至0.14±0.011 4；不动杆菌属（Acinetobacter）0.07±0.004 6降至0.04±0.004 4；由此可知曲拉的第一优势菌乳球菌属作为曲拉的主要功能菌[19-20]，其在干燥的过程中丰度发生了显著的增加，而作为条件致病菌属的沙雷氏菌属在干燥的过程中丰度显著降低[21]；由此可知曲拉的传统形式干燥曲拉在食用功能性和安全性上占据一定的优势。