基于16S rRNA高通量测序的西藏农、牧区牦牛酸奶菌群多样性分析

刘怡萱，许国琪，曹鹏熙，金彦龙，李小燕，刘 星,,

（1.西藏大学理学院，西藏 拉萨 850000；2.西藏大学生态学研究中心，极端环境生物资源与适应性进化实验室，西藏 拉萨 850000；3.武汉大学生命科学学院，武汉大学-西藏大学生态与保护联合研究中心，湖北 武汉 430072）

牦牛在青藏高原的驯化起始于一万年前，是极端环境条件下生存繁衍的优势畜种，牦牛奶是藏民族赖以生存的奶源和重要的经济来源之一[1-3]。除直接饮用外，经过加工产生的多种牦牛乳制品包括酸奶、曲拉、酥油、奶酪、乳清、奶茶等也是高原人民生产生活中不可或缺的食物，其中牦牛酸奶更是藏族人民最常食用的乳制品之一[4]。在西藏强紫外线、低压缺氧、昼夜温差大的极端环境下，对于发酵温度及时间没有精准控制的自然发酵牦牛酸奶的营养价值越来越受到关注，其中特殊的微生物资源也逐渐成为研究的热点[5]。传统发酵牦牛酸奶具有独特的风味和口感，参与发酵的各种微生物进入人体消化系统，部分菌群在人体的肠道菌群的建成中起重要作用，对人类健康产生影响[6-7]。

对发酵乳制品中微生物组成的研究始于19世纪末，早期研究中，采用传统分离培养的方法对其参与发酵微生物组成有了初步的认识[8-10]。然而自然界中可培养研究的微生物还不到整个微生物种群的20%，在各种环境样品中仅有小于1%的微生物能通过实验室培养技术获得纯培养，大部分微生物还属于不可培养或难以培养的微生物[11-12]。此外，传统的免培养技术，如变性梯度凝胶电泳[13-14]、限制性片段长度多态性[15]、荧光原位杂交[16]等所能鉴定到的微生物种类有限[6]。随着测序技术的发展，高通量测序已成为检测微生物多样性的重要手段，其读长高、精度高、通量高、无偏性的优势被广泛应用到微生物群落多样性的研究中，不仅可以检测到难培养和低丰度的物种以及曾经存活或是难以分离的部分微生物，还可以排除挑选菌株时所产生的偏见性，更能代表整体的微生物组成[5-6]。

极端环境下的微生物经过长期的自然选择，具备相对特殊的结构、机能以及遗传特性，以应答相应的胁迫因子[5]。目前，对于青藏高原传统发酵牦牛酸奶已有相关报道涉及加工工艺优化[17]、营养成分分析[18]、乳酸菌和酵母菌分离与鉴定[19]、乳酸菌益生作用[5]、挥发性风味成分分析[20]、微生物组成及其多样性分析[4,21-22]等，但进一步对西藏农、牧区牦牛酸奶中微生物差异及其原因的研究则鲜见报道。

本研究利用高通量测序方法对西藏不同生产方式下传统发酵的牦牛酸奶细菌群落组成和结构进行比较分析，并结合环境因子分析农、牧区牦牛酸奶间微生物多样性差异，以期为牦牛酸奶的进一步开发利用提供基础参考。

1 材料与方法

1.1 材料与试剂

共采集西藏牧民以传统方式发酵的牦牛酸奶样本19 个，包括在日喀则（RK）、拉萨（LA）等农区采集的样本5 个及在阿里（AL）、那曲（NQ）等牧区采集的样本14 个，样品信息见表1。样品收集在15 mL无菌离心管中，用-20 ℃车载冰箱运回实验室后，保存于-80 ℃冰箱备用。

聚合酶链式反应（plymerase chain reaction，PCR）Master Mix（M7505） 美国Promega公司；AxyPrep DNA凝胶回收试剂盒 美国Axygen公司；DL2000 DNA Marker 宝生物工程（大连）有限公司。

1.2 仪器与设备

C12S车载冰箱 美国Mobicool公司；DW86L-158超低温冰箱 中国江苏胜蓝公司；NanoDrop 2000超微量分光光度计 美国Thermo Scientific公司；5417R冷冻离心机德国Eppendorf公司；Veriti PCR仪 美国ABI公司；HiSeq 2500测序平台 美国Illumina公司。

1.3 方法

1.3.1 菌群DNA提取和PCR扩增测序

牦牛酸奶中菌群总DNA的提取参照Yu Jie等[23-24]方法。用提取的总DNA作为模板，以338F（5’-ACTCCTACGGGAGGCAGCAG-3’）、806R（5’-GGACTACHVGGGTWTCTAAT-3’）为引物，对细菌的16S rRNA V3-V4可变区进行特异性扩增。PCR扩增体系（50 μL）：10×Buffer 5 μL，4×dNTP混合物5 μL，MgCl23 μL，引物各2.6 μL，模板2 μL，Taq酶0.4 μL，加ddH2O补充至50 μL。反应程序为：95 ℃预变性2 min后，进行25 次循环（95 ℃变性30 s，55 ℃复性30 s，72 ℃延伸30 s），最后72 ℃延伸5 min。在超微量紫外分光光度计上检测PCR产物纯度和浓度，1%琼脂糖凝胶电泳检测DNA完整性，基于Illumina HiSeq 2500进行测序，PCR扩增及测序工作由广东仁秀测序公司完成。

基于上海美吉生物I-Sanger云平台（https://www.i-sanger.com/）、R软件（version 3.5.1）以及ITOL（https://itol.embl.de/）等软件平台进行生物信息学及数据统计分析。

1.3.2 环境因子的获取

利用长天HOLUX M-241无线GPS记录器对19 个酸奶样本进行定位并获取海拔等信息，利用CHELSA高精度数据库（Climatologies at high resolution for the earth’s land surface areas，http://chelsa-climate.org/）以及R软件geosphere软件包分别调取所在地年平均气温、地理距离信息。

1.3.3 菌群多样性分析

对测序得到各样本序列进行质量控制后，利用Usearch（vsesion 7.0）软件，按97%相似度的分类水平划分可操作分类单元（operation taxonomic units，OTU），并按最小样本序列数进行抽平，利用mothur（version v.1.30.1）计算α多样性指数，构建稀释曲线。

利用Quantitative Insights Into Microbial Ecology（QIIME，version 1.17）计算β多样性距离矩阵。基于bray_curtis距离进行相似性分析（analysis of similarities，ANOSIM），基于unweighted UniFrac距离对农、牧区牦牛酸奶进行主坐标分析（principal coordinates analysis，PCoA）。

进一步对菌群组成进行多级物种差异判别分析（linear discriminant analysis effect size，LEfSe），再根据分类学组成对样本按照不同分类水平进行线性判别分析（linear discriminant analysis，LDA，阈值为2），找出对样本划分产生显著性差异影响的菌群。进一步利用组间显著性差异检验，对不同区域之间的菌群进行假设检验，评估菌群丰度差异的显著性水平。

此外，基于地理距离、海拔差异以及年平均气温等环境因子数据，分析样本间差异原因。一方面，利用R软件geosphere程序包计算得到各样点间的地理距离矩阵，综合β多样性距离矩阵，对两者进行回归分析。另一方面，结合基于距离的冗余分析（distance-based redundancy analysis，db-RDA）研究海拔高度以及年平均气温与样本差异的关系。

2 结果与分析

2.1 农、牧区牦牛酸奶菌群多样性分析

2.1.1α多样性分析

对19 个样本进行测序共得到204 345 条有效序列，基于97%相似度的分类水平，共注释到48 个OTU，测序覆盖率大于99.97%。各样本群落丰富度Sobs指数稀释曲线和群落多样性Shannon指数稀释曲线趋于平缓，表明测序数据量合理，可以反映样本中绝大多数的菌群多样性信息（图1）。

α多样性分析表明（表2），Chao指数在农区样本和牧区样本间具有显著差异（P=0.009 2），农区样本群落丰富度大于牧区样本。Shannon指数（P=0.250 7）、Simpson指数（P=0.917 9）在农区样本与牧区样本间无显著差异，表明两组样本中菌群多样性类似。

表2 基于T检验的α多样性指数差异Table 2 Difference in α diversity index based on Student’s T test

2.1.2β多样性分析

ANOSIM（图2A）显示，农、牧区组间菌群结构差异显著大于组内差异（R=0.673 8，P=0.003），分组具有意义。PCoA（图2B）表明，来自农区和牧区的样本分别聚集在一起，其中，LA03样本介于两组之间，第1主成分对样本菌群结构差异的贡献率为42.08%，第2主成分为20.02%。β多样性综合分析显示，农、牧区牦牛酸奶样本中菌群结构存在明显差异。

图2 β多样性分析Fig.2 β-Diversity analysis

2.2 农、牧区牦牛酸奶菌群组成及差异分析

在所有样本中鉴定得到的48 个OTU，分属于4 个门，18 个属，在门水平上（图3A），厚壁菌门（Firmicutes）在农、牧区牦牛酸奶菌群中均占主要优势（96.589 5%、99.796 1%），而放线菌门（Actinobacteria）、变形菌门（Proteobacteria）、拟杆菌门（Bacteroidetes）在农区样本中的占比均高于牧区样本。从属分类水平看（图3B），农、牧区牦牛酸奶中乳杆菌属（Lactobacillus）（34.504 9%、70.379 9%）、乳球菌属（Lactococcus）（55.687 6%、0.116 9%）和链球菌属（Streptococcus）（3.548 1%、29.289 4%）为主要优势属，其他菌属相对丰度均小于1%。

图3 西藏农、牧区牦牛酸奶中菌群组成Fig.3 Bacterial composition of yak yogurt samples collected from agricultural and pastoral areas of Tibet

图4 西藏农、牧区牦牛酸奶中菌群差异分析Fig.4 Analysis of bacterial composition differences between yak yogurt samples collected from agricultural and pastoral areas of Tibet

LEfSe多级物种层级分析（图4A）显示，拟杆菌门、放线菌门在农区显著富集，而在牧区大量富集厚壁菌门；从属水平来看，乳球菌属等8 个菌属在农区显著富集，牧区仅有链球菌属、乳杆菌属显著富集。LDA（LDA阈值＞2.0）（图4B）与组间差异显著性检验（图5）共同表明，对农、牧区样本菌群结构具有显著影响的菌属包括乳球菌属、链球菌属、乳杆菌属等10 个菌属（P＜0.05）。

图5 西藏农、牧区牦牛酸奶中菌群组间差异的显著性检验Fig.5 Significant test of the differences in bacterial composition of yak yogurt samples between agricultural and pastoral areas of Tibet

2.3 环境因子对农、牧区牦牛酸奶菌群结构的影响

对地理距离矩阵和β多样性距离矩阵拟合得到方程为：Y=0.067 4lnx＋0.305 59（R2=0.028，P=0.028）（图6A），结果表明地理距离对菌群结构的差异解释度较低。基于unweighted UniFrac距离算法的db-RDA显示（图6B），海拔高度与年平均气温对农、牧区牦牛酸奶细菌组成均有显著性影响，且呈负相关，海拔高度（R2=0.730 1，P=0.001）对样本间差异的解释度大于年平均气温（R2=0.482 6，P=0.007）。

图6 环境因子对农、牧区牦牛酸奶菌群结构的影响Fig.6 Effects of environmental factors on the bacterial structures of yak yogurt samples collected from Tibetan agricultural and pastoral areas

3 讨 论

农区和牧区是西藏主要的两种传统生产生活方式，农区主要位于拉萨、山南、日喀则，集中分布在藏南河谷及河谷盆地中，以种植青稞、油菜等作物为主；牧区主要集中在阿里、那曲，天然草地面积位居全国第一[25]。牦牛作为西藏畜牧业的支柱，养殖主要集中在牧区，占据十分重要的地位[26]，农区由于地理环境原因，牦牛养殖量远小于牧区，养殖方式与牧区也有较大区别，人工添加饲料比重远大于牧区。然而，无论在农区还是牧区，传统发酵的牦牛酸奶都是农牧民发展经济、增加收入不可或缺的重要产品之一[27]。对自然发酵牦牛酸奶中菌群多样性的研究发现，厚壁菌门中乳杆菌属和乳球菌属为优势菌群[4,28]。本研究来自西藏的样品，菌群组成与之前研究结果基本一致，农区样本中菌群丰富度及多样性略高于牧区样本。

发酵食品中菌群多样性差异与地理环境密切相关[25,29-32]，酸奶中菌群种类在气候相对温和的地域比在相对寒冷和海拔较高的环境中更为丰富[5]，自然发酵酸奶中菌群的组成及数量还受发酵时间和制作方法等其他因素影响[22]。本研究除气候及海拔的因素外，还考虑了农区和牧区养殖方式的不同，农区由于受地理环境的限制，牦牛喂养过程中人工添加饲料比重大于牧区。研究发现，地理距离对酸奶中菌群结构影响较小，海拔影响大于温度影响。在海拔相对较低的农区牦牛酸奶中乳球菌属和乳杆菌属占优势，且乳球菌属含量高于链球菌属和乳杆菌属；而在海拔相对较高的牧区牦牛酸奶中乳杆菌属和链球菌属占优势，且乳杆菌属占比显著高于乳球菌属和链球菌属。乳杆菌属对酸的耐受力一直被认为高于链球菌属和乳球菌属，在发酵过程中进行微弱的水解作用，有助于维持较高的发酵活性[33]。在西藏牦牛酸奶的传统发酵方式中，对于温度及发酵时间没有精确的控制，由于牧区平均海拔高于农区，年平均气温相对较低，酸奶在牧区发酵时间略长于农区，因此导致牧区酸奶中耐酸并能提高发酵活性的乳杆菌含量高于农区酸奶。

生活在青藏高原的藏族人群，传统的饮食习惯以糌粑及肉类为主，果蔬类及其他蛋白质摄入量较少，饮食结构较为单一，牦牛酸奶作为其最常食用的乳制品之一[4]，对调节胆固醇含量、维持机体健康起到了重要作用[5,31-33]。由于青藏高原海拔差异较大，不同地域环境下的藏族人群饮食习惯也存在很大差异。人类饮食对肠道菌群有直接影响，而肠道菌群对人的健康更是发挥着至关重要的作用[34]。已有研究表明，藏族人群肠道菌群多样性在不同海拔[35]、不同生存环境[36]下有一定差异。西藏牦牛酸奶及其他发酵食品中的菌群对藏族人群肠道菌群及其对健康的影响机制尚不明确。因此，对青藏高原传统食品中菌群的进一步分析将有助于揭示自然环境与传统食品中菌群的关系以及环境、食品与不同生活方式下藏族人群肠道菌群三者之间的适应性关系。