西藏不同产区曲拉细菌群落结构的比较分析

刘振东，毕娜，李哲，李梁，罗章，薛蓓*，汪雯翰

1(西藏农牧学院 食品科学学院，西藏 林芝， 860000)2(西藏自治区藏药审评认证中心，西藏 拉萨，850000)3(上海市农业科学院食用菌研究所,上海， 201403)

曲拉(藏语，指奶干渣)是藏区传统发酵乳制品之一,由牦牛奶脱脂后，在自然条件下发酵酪蛋白凝固并风干制成，属于粗奶酪的一种[1-2]。作为传统的藏族食品，曲拉具有蛋白质含量高、保存时间长、携带方便等优点，在整个藏区被广泛食用[3-4]。

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

基于微生物分离、纯化、培养及鉴定的传统研究方法难以获得培养困难及痕量的微生物，而且难以展现环境中完整的微生物种群结构及生态关系[8-9]。目前，高通量测序技术已经逐渐成为研究微生物组成及结构的一种重要方法，近年该技术已广泛用于原料乳和各类发酵乳等乳制品生态位的研究[10-11]。

本研究旨在利用高通量测序技术，全面了解西藏7个不同产区牦牛曲拉的细菌群落结构及相互差异，为西藏不同产区曲拉的安全性评价和功能微生物的挖掘奠定一定的理论基础。

1 材料与方法

1.1 材料与试剂

21份牦牛曲拉样品，分别来自西藏拉萨(LS)、日喀则(RKZ)、山南(SN)、昌都(CD)、林芝(LZ)、那曲(NQ)和阿里共7个地区；E.Z.N.A.Soil DNA试剂盒，美国OMEGA公司；Qubit2.0 DNA检测试剂盒，美国Invitrogen公司；Q5®High-Fidelity DNA Polymerase，美国NEB公司；凝胶回收试剂盒，德国Qiagen公司。

1.2 仪器与设备

5424R高速冷冻离心机，德国Eppendorf公司；DYY-6C型电泳仪，北京六一仪器厂；T-100 PCR仪，美国Bio-rad公司；MiSeq高通量测序仪，美国Illumina公司；BG-gdsAUTO(130) 凝胶成像系统，北京百晶生物技术有限公司。

1.3 实验方法

1.3.1 样品的采集

2018年6月采集拉萨市林周县、日喀则市萨迦县、山南市隆子县、昌都市八宿县、林芝市米林县、那曲市申扎县、阿里地区改则县7个地区的曲拉样品，分别编号为LS、RKZ、SN、CD、LZ、NQ、AL。将所有样品置于无菌离心管中，放入冰盒保存并立即运送到实验室，置于-80 ℃冰箱中备用。

1.3.2 DNA提取和PCR扩增

对于不同来源的微生物组样本，根据E.Z.N.A.®soil试剂盒说明书提取各牦牛曲拉样品中的总DNA，然后通过0.8%琼脂糖凝胶电泳的方法检测总DNA的提取质量。本研究采用NEB公司的Q5高保真DNA聚合酶保证结果的可靠性，采用标准细菌V3-V4引物338F(5’-ACTCCTACGGGAGGCAGCA-3’)和806R(5’-GGACTACHVGGGTWTCTAAT-3’)进行PCR扩增。

PCR体系(25 μL):基因组DNA模板2 μL，dNTP 2 μL，ddH2O8.75 μL，缓冲液5 μL，Q5 DNA 0.25 μL，DNA模板2 μL，5×GC缓冲液5 μL。

PCR扩增条件：94 ℃预变性2 min；94 ℃变性30 s，55 ℃退火30 s，72 ℃延伸40 s，扩增40个循环；72 ℃延伸1 min。尽量保证同一批样本的扩增条件一致。

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。然后在GenBank数据库进行Blast比对，并提交到美吉生物的云分析平台进行多样性分析(http://www.i-sanger.com/)。

2 结果与分析

2.1 测序数据统计及Alpha多样性分析

采用Illumina MiSeq平台对7个地区曲拉样品测序结果进行分析统计，结果如表1所示。由表1可知21个样品总计测得有效序列条数为(241 148±23 200)条，以97%的序列相似度作为OTU划分阈值，共产生(3 994±78)个OTUs，由序列数及OTUs可以看出，曲拉中细菌种类繁多、物种丰富，且不同来源的曲拉存在明显差异。不同区域曲拉样品间Alpha多样性统计中的Chao1指数、Ace指数差异显著(P<0.05)。

表1 不同样本的测序数据统计

注：“-”表示无需统计；不同小写字母表示差异显著(下同)

2.2 稀释曲线

由图1可知，当样品测序量较低时，每个样品的Shannon指数随着测序量增加呈显著上升趋势；但随着测序量逐渐增加，各曲拉样品的Shannon指数增加逐渐趋于缓慢；当样品测序量到达5 000后各样品的Shannon指数基本稳定，说明样品测序量达到5 000时已足以反映21个曲拉样品中细菌群落多样性。

图1 各样品测序的稀释曲线

2.3 在门水平上各样品菌群分类及比较

由表2可知，7个产地的曲拉样品共检测出24个门，其中LS检测出全部24个菌门，RKZ中检测出6个菌门，SN中检测出7个菌门，CD中检测出6个菌门，LZ中检测出12个菌门，NQ中检测出20个菌门，而AL最少仅检测出2个菌门。7个地区曲拉样品的主要优势菌门均为厚壁菌门和变形菌门，除LS和CD外，其他5个地区曲拉样品的第一优势菌均为Firmicutes,其丰度值在61.01%～96.54%,LS和CD的第一优势菌Proteobacteria,其丰度值在40.60%～75.09%。

2.4 在属水平上各样品菌群分类及比较

如表3所示，对7个地区21份曲拉样品从属水平进行鉴定共鉴定出550个属，其中LS中检测出409个属，RKZ中检测出47个属，SN中检测出69个属，CD中检测出36个属，LZ中检测出99个属，NQ中检测出330个属，AL中检测出15个属。各地区样品在属水平上检出的种类差异较大，LS>NQ>LZ>SN>RKZ>CD>AL。

表2 各样品门水平上菌群分布相对丰度 单位：%

表3 各样品属水平上菌群分布相对丰度 单位：%

如表4所示，从属水平分析，多数产地曲拉样品中优势菌属主要为乳球菌属Lactococcus、乳杆菌属Lactobacillus，其中RKZ中Lactococcus为(63.73±8.29)%,Lactobacillus为(4.15±1.38)%;SN中Lactococcus为(25.82±1.42)%;CD中Lactococcus为(46.62±4.60)%,Lactobacillus为(9.20±6.49)%;AL中Lactococcus为(85.13±2.99)%,Lactobacillus为(8.57±1.33)%，与曹磊等[12]研究的青海甘南地区牦牛曲拉群落结构存在共性的结果一致。乳杆菌属和乳球菌属均为乳酸菌，可在曲拉制作的初期产生大量的乳酸、乙酸和一定量的胞外多糖影响曲拉风味[13-14]，从优势菌角度分析，以上几个地区的曲拉风味特点存在一定的相似性。而LS中优势菌属与其他地区差异较为明显，其第一优势菌为雷尔氏菌属Ralstonia为(17.25±7.05)%,且致病菌假单胞菌属Pseudomonas相对丰度达(16.05±3.22)%、苍白杆菌Ochrobactrum达(9.52±3.19)%，均明显高于其他地区，从优势菌角度分析其风味及品质特点相对其他地区具有一定特殊性。

2.5 物种丰富度热图

根据OTU的丰度数据，热图可在属水平上将不同的OTU分块聚类，并根据颜色梯度反应不同样品中细菌群落相似性、差异性及物种聚类关系。利用R语言绘制热图，可直观地显示不同产区曲拉样品中细菌的差异。由图2可知，上侧为样品聚类树，左侧为OTU聚类树，西藏7个地市的21个曲拉样品可分为2个大聚类：LS1和LS3可聚为一类，其余19个样品聚为一类；物种丰富度上红色越亮代表该菌属在该样品中的相对含量高于横向的其他样品，绿色越亮则相反，不同地区主要优势菌属存在一定差异。

表4 不同地区曲拉优势菌分布表 单位：%

图2 各样品不同属组成和相对丰度热图

2.6 基于UniFrac距离的样本聚类分析

由图3可知，细菌群落第1主成分(PC1)的贡献率是64.87%，细菌群落第2主成分(PC2)的贡献率是21.35%，两者之和为86.22%，大于70%，表明两个成分较好地代表了曲拉样品中的细菌群落信息。不同曲拉样本的菌群与地理位置存在一定关系，同一个采样区的样本具有一定的聚类趋势，同时发现LZ、RKZ和CD可聚为一类,而其他4个地区的样品较为离散，各自独立。

图3 PCA分析的样本二维排序图

2.7 UPGMA聚类分析图

聚类分析主要指层次聚类(hierarchical clustering)的分析方法，以等级树的形式展示样品之间的相似度，通过聚类树的分枝长度衡量聚类效果的好坏。样本根据彼此之间的相似度聚类，样本间的分枝长度越短，两样本越相似。由图4可知，不同采样地的样品聚集到不同的类别，说明不同来源样品的微生物多样性存在一定的差异性，这与PCA分析的结果一致。

图4 各曲拉样品聚类分析图

2.8 样品功能类群分布统计

通过将现有的16S rRNA基因测序数据与代谢功能已知的微生物参考基因组数据库相对比，从而实现对细菌和古菌代谢功能的预测。由图5可知，7个样品中细菌群落的代谢功能主要集中在氨基酸代谢(amino acid metabolism)、碳水化合物代谢(carbohydrate metabolism)和能量代谢(energy metabolism)，其中AL和SN样品在氨基酸代谢丰度上明显低于其他5个地区样品，而在碳水化合物代谢丰度上LS样品显著地区其它6个地区样品，而在能量代谢丰度上7个地区样品差异不显著。

图5 各曲拉样品功能类群分布图

3 讨论

本研究采用高通量测序技术对西藏7个地市共计21份曲拉样品的细菌群落结构进行分析，从属水平上群落结构分布可知，不同地区曲拉细菌多样性存在一定差异且存在一定的地域独立性，如拉萨(LS)和那曲(NQ)样品的细菌种类远高于其他5个地区，细菌丰度上除拉萨地区革兰氏阴性植物病原细菌雷尔氏菌属(Ralstonia)和致病菌假单胞菌属(Pseudomonas)占主要优势菌群外，其他6个地市的优势菌均主要集中于发酵乳制品常见菌乳球菌属(Lactococcus)和乳杆菌属(Lactobacillus)，这6个地区曲拉的主要优势细菌群落结构与SUN等[15]对青海牦牛酸奶进行纯化培养和16S rRNA鉴定的结果一致，同时与其他传统发酵乳制品如酸牛奶、酸马奶、酸驼奶等细菌群落多样性类似[16-17]，在门水平上常见原料乳的细菌组成差异不明显[18-19]。而本研究发现拉萨地区样品存在一定的植物病原菌菌属和条件致病菌菌属。通过功能分析发现曲拉中细菌群落有助于人体消化过程中的氨基酸和碳水化合物代谢。

由于本研究采样点较单一及样品量相对较少，并不能直接反映出制备环境直接造成的影响，为进一步探讨拉萨地区牦牛曲拉的细菌多样性及安全性，应增加采样点和采样量，进一步开展细致研究。