青海玉树地区牦牛乳源微生物多样性研究

孙璐

(青海大学畜牧兽医科学院 青海省高原放牧家畜动物营养与饲料科学重点实验室，西宁810016)

玉树藏族自治州位于青海省西南部青藏高原腹地的三江源头，生活在这里的人们处于寒冷、缺氧、强紫外线等高海拔环境中，一年中大部分时间无法补充蔬菜和水果，但是他们能够在极端环境中健康生活，很大程度上取决于牦牛乳，其产品构成他们日常饮食的主要组成部分［1］。牦牛乳中含有多种对哺乳动物后代健康很重要的复杂成分，主要包括蛋白质、乳糖和脂肪，以及一些相对分子量较低，具有特殊生物学功能的低聚糖、磷脂和鞘脂等化合物［2］。牦牛乳的产量较荷斯坦牛极低，但其蛋白质、干物质、脂肪及乳糖等营养成分含量高于其他乳，富含共轭亚油酸(CLA)、免疫球蛋白等功能性成分，是一种天然、独特的浓缩乳，被认为是制造婴儿、老年人及特殊需求人群食品的优质原料乳［3］。牦牛乳营养含量高，含有丰富的微生物群系。牦牛乳中的微生物从各种来源获得，发挥着多种作用，如促进乳品发酵(乳球菌、乳杆菌、链球菌、丙酸杆菌和真菌种群)、导致腐败(假单胞菌、梭状芽孢杆菌、芽孢杆菌)、有益于健康(乳酸杆菌和双歧杆菌)或引起疾病(李斯特菌、沙门氏菌、大肠杆菌、弯曲杆菌和产生霉菌毒素的真菌)等［4］。关于微生物早期的研究大多基于培养的方法获得微生物的组成，对于无法培养的微生物则不能通过该方法解析。高通量测序技术的发展，克服了微生物培养检测的局限性，被广泛用来鉴定和发现微生物［5］。本文为了解牦牛液态乳的质量和安全，对青海省玉树地区的生鲜牦牛乳源微生物进行高通量测序，为牦牛乳的开发利用及优势菌群的发掘等提供参考。

1 材料和方法

1.1 样品采集

2019年7～8月，选取玉树地区自然放牧，体况、年龄、胎次(4～5胎次)基本一致的泌乳中期母牦牛，母牛和犊牛分开过夜，每天人工挤奶1～2次但不完全挤完。3个样点采样方法相同。

采样地点为青海省玉树州称多县歇武镇、曲麻莱长江村、玉树州果青牧场3个样点。每个样点采集6个样(yak1、yak2、yak3)。采样点基本信息见表1。每份乳约50 mL，取5 mL于无菌无酶冻存管(康宁/Corning公司)中，液氮速冻，运回实验室后储存在-80℃冰箱备测。

表1 采样点基本信息

1.2 仪器与设备

DCC体细胞检测仪(利拉伐天津有限公司);超低温冰箱(-80℃，美国Thermo公司);Fresco21型高速冷冻离心机(美国Thermo公司);Tprofessiona PCR仪(德国Biometral公司);Gel DOC XR凝胶成像系统;PowerPacUniversal水平电泳仪;SUBCELL GT电泳槽(20 cm×25 cm，美国Bio-Rad公司)。

1.3 试验方法

1.3.1 基因组总DNA提取及PCR扩增、测序参照YU［16，17］等 的 方 法 提取3组牦 牛 乳 样 本 中 的DNA，以16S V4区引物(515F和806R)进行PCR扩增。通过TruSeq DNA PCR-Free Sample Preparation Kit构建文库，文库合格后，使用NovaSeq6000进行上机测序。PCR的扩增、PCR产物的混样、纯化，文库的构建和上机测序流程均由天津诺禾致源生物信息科技有限公司完成。基于天津诺禾致源云平台(https://magic.novogene.com)、SPSS Statistics软件(version 19.0)、I-Sanger云平台(https://www.isanger.com/)、R软件(version 3.5.1)以及ITOL(https://itol.embl.de/)等软件平台进行生物信息学及数据统计分析。

1.3.2 环境因子的获取利用天宝juno sa无线GPS记录器定位获得3个采样地点的坐标、海拔等信息，在国家气象科学数据中心(http://data.cma.cn/site/index.ht mL)获取采样点年平均气温。

1.4 数据计算

1.4.1 测序数据处理测序所得原始序列截去接头和引物序列后，使用FLASH(version 1.2.7)进行拼接获得原始Tags数据，原始Tags数据经Qiime(version 1.9.1，http://qiime.org/scripts/split_libraries_fastq.ht mL)质量控制后获得高质量的Tags数据(Clean Tags)，通过(https://github.com/torognes/vsearch/)与物种注释数据库进行比对检测嵌合体序列，并最终去除其中的嵌合体序列，得到最终的有效数据(Effective Tags)。

1.4.2 微生物多样性分析Shannon，Observed-OTUs，Chao及Simpson指数使用Qiime软件(Version 1.9.1)计算，Alpha多样性的指数组间差异分析使用R软件(Version 2.15.3)完成。利用QIIME软件(version 1.17)计算Beta多样性距离矩阵。相似性分析(Analysis of similarities，ANOSIM)基于bray_curtis距离分析，基于unweighted UniFrac距离对不同海拔牦牛乳微生物进行主坐标分析(PCoA)。对牦牛乳微生物组成进行差异判别分析(LEfSe)，再根据分类学组成对样本按照不同的分类水平进行线性判别分析(LDA，阈值为4)，找出对样本划分产生显著性差异影响的菌群。

2 结果

2.1 玉树地区牦牛乳微生物测序结果及多样性指数

质控后按照97%相似性进行可操作分类单元(operational taxonomic units，OTU)聚类，共得到12329个OTUs。绘制等级聚类曲线(Rank Abundance)和稀释曲线(Rarefaction curve)。稀释曲线是从样本中随机抽取一定测序量的数据，统计所代表物种数目(OTUs数目)，以抽取的测速数据量与对应的物种数来构建曲线，稀释曲线直接反映测序数据量的合理性，间接反映样本中物种的丰富度，图1曲线逐渐趋向平坦，表明本试验的测序数据量合理。

图1 Alpha多样性指数稀释曲线图

等级聚类曲线是将样本中OTUs按相对丰度由大到小排序得到对应的排序编号，再以OTUs的排序编号为横坐标，OTUs中的相对丰度为纵坐标，将这些点用折线连接，可直接反映样本中物种的丰富度和均匀度。如图2所示，在水平方向上，曲线横轴跨度逐渐变大，说明物种丰富度也随之变大;在垂直方向上，曲线逐渐平滑，说明物种分布均匀。

图2 等级聚类曲线

通过Shannon和Simpson指数评价微生物的多样性，Observed species和Chao1评价微生物菌群丰度。由表2可知，青海玉树地区不同样点牦牛乳样品中微生物的丰度存在显著差异，但3组样品的微生物多样性相似。

表2 Alpha多样性指数差异

2.2 玉树地区牦牛乳微生物组成

2.2.1 物种维恩图(Venn Diagram)分析对玉树地区不同样点的牦牛乳微生物制作物种Venn图，由图3可知，3组样品共有3104个OTUs。其中玉树州果青牧场(yak3)的牦牛乳样品中特有2697个OTUs，是3组中数量最多的;称多县歇武镇歇武村(yak1)牦牛乳样品中特有的OTUs个数为915;曲麻莱长江村(yak2)的牦牛乳样品特有OTUs 564个，是3组中最少的。

图3 韦恩图

2.2.2 玉树地区牦牛乳微生物门分类水平的组成对3组样品中的微生物进行门水平前10种物种相对丰度分析，结果如图4。相对丰度排名前10的微生物分别是:变形菌门(Proteobacteria)，厚壁菌门(Firmicutes)，拟杆菌门(Bacteroidetes)，放线菌 门(Actinobacteria)，蓝藻门(Cyanobacteria)，浮霉菌门(Planctomycetes)，酸酐菌门(Acidobacteria)，疣微菌门(Verrucomicrobia)，绿弯菌门(Chloroflexi)，软壁菌门(Tenericutes)。其中，变形菌门在yak1和yak2组中的含量最高，分别占比50.40%和44.52%;yak3组中变形菌门和厚壁菌门含量较多，占比接近，分别为35.09%和36.70%;软壁菌门在yak1组中的含量最少(0.40%)，疣微菌门在yak2组中的占比最小(0.18%)，浮 霉 菌 门 在yak3中 的 含 量 最 少(0.15%)。

图4 门水平玉树地区微生物的相对丰度

2.2.3 玉树地区牦牛乳微生物属分类水平的组成玉树地区不同样点牦牛乳微生物属水平物种丰度结果表明(图5)，微生物组成占比前10的分别是unidentified_Burkholderiaceae、金黄杆菌属(chryseobacterium)、耶尔森菌属(Yersinia)、乳杆菌属(Lactobacillus)、unidentified_Cyanobacteria、不动杆菌属(Acinetobacter)、假单胞菌属(Pseudomonas)、乳酸链球菌属(Lactococcus)、链球菌属(Streptococcus)、unidentified_Enterobacteriaceae。Yak1组中占比最高的为不动杆菌属(7.20%);yak2组和yak3组占比最高的分别是金黄杆菌属和乳杆菌属。乳酸链球菌属在yak1组和yak2组中的含量最少，unidentified_Burkholderiaceae在yak3组中占比最低。

图5 属水平下玉树地区牦牛乳微生物的相对丰度

2.2.4 玉树地区牦牛乳微生物种水平的组成玉树地区牦牛乳样品前10个种水平微生物组成比较表明，唐菖蒲伯克霍尔德菌(Burkholderia gladioli)和Chryseobacterium_vrystaatense在yak2组的比例最高，分别为12.40%和17.17%;鲁氏不动杆菌(Acinetobacter lwoffii)占比较高，在yak1组中含量最高;Lactobacillus_iners在yak1，yak2组的含量较低，唐菖蒲伯克霍尔德菌在yak3组中所占比例最低。

表3 种水平下玉树地区牦牛乳微生物的相对丰度

2.2.5 玉树地区牦牛乳物种主坐标分析基于unweightedUniFrac距离的PCoA主坐标分析结果显示，第一主成分对样品差异的贡献值为14.38%，第二主成分对样品差异的贡献值为10.36%。yak1、yak2及yak3组各组间距离较近，可能是由于3组样本采样地点距离、海拔较近，经纬度较为接近，且yak2和yak3组的微生物结构相似度较高。

图6 基于unweightedUniFrac距离的PCoA主坐标分析

2.3 玉树地区牦牛乳微生物组成组间差异分析

LDA值分布柱状图(图7A)(LDA score＞4)体现了玉树地区不同样点牦牛乳间微生物具有统计学差异的物种，由图可知yak3组中差异物种影响较大的是厚壁菌门、乳杆菌目(Lactobacillales)、芽孢杆菌纲(Bacilli);yak2组中差异物种影响较大的是黄杆菌目(Flavobacteriales)、金黄杆菌属(Chryseobacterium)、unidentified_Flavobacteriales;yak3组中差异物种影响较大的是耶尔森菌属、红细菌目(Rhodobacterales)、红细菌科(Rhodobacteraceae)。

玉树不同地区牦牛乳进化分支图(图7B)由内而外辐射的圆圈代表了由门至属(或种)的分类级别，图中无显著差异的统一着色为黄色，差异显著的物种跟随分组进行着色。由图可知，yak1组中差异显著的物种为unidentified_Cyanobacterin、红细菌科(Rhodobacteraceae)、红细菌目(Rhodobacterales)、鞘脂单胞菌科(Sphingomonadaceae)、鞘脂单胞菌目(Sphingomonadales)，yak2组中差异显著的物种为unidentified_Flavobacteria、黄杆菌目(Flavobacteriales)，yak3组中差异物种为乳杆菌目(Lactobacillales)、链球菌科(Streptococcaceae)、乳杆菌科(lactobacillaceae)、芽孢杆菌纲(Bacilli)。

图7 LEfSe组间差异分析 A：LDA值分布柱状图；B：进化分支图

3 讨论

牦牛乳中含有多种微生物，包括有益菌及各种致病菌。致病菌是食源性疾病的主要来源，对犊牛及消费者的健康造成威胁。本研究结果发现，门水平玉树地区牦牛生鲜乳中主要的微生物为变形菌门和厚壁菌门;属水平中主要的优势菌为不动杆菌属、金黄杆菌属和乳杆菌属，这与现有的研究结果一致［6，7］。

牛乳中的微生物多样性对于预测的乳品质量具有重要意义。牛乳中的微生物受来源动物、环境、温度、海拔等多种因素的影响［8，9］。变形菌门和耶尔森氏菌是常见的粪便菌群，其中变形菌门(如大肠杆菌、克雷伯氏菌、肠杆菌、变形杆菌和沙雷氏菌)在健康机体的肠道中很少见，多见于早期肠道炎症的机体内［10，11］。耶尔森氏菌是一种在粪便中存在的与肠炎有关的致病菌，耶尔森氏菌属中的3个物种对人类具有致病性，即鼠疫耶尔森氏菌、假结核耶尔森氏菌和小肠结肠炎耶尔森氏菌［12］。金黄杆菌属属于黄杆菌科，包含革兰氏阴性、黄色色素、杆状和不形成孢子的细菌种类，它们可以是自由生活或寄生［13］。其以促进巴氏杀菌牛乳变质和再污染的能力而闻名，从而降低牛乳质量，也是引起奶牛乳腺炎的机会性病原体［14，15］。这些菌均为致病菌，一般来说，不应该存在于牛乳中。可能是在挤奶过程中与空气、粪便接触进入牛乳。牧区应尽可能规范牦牛饲养、挤奶过程，做好挤奶前牦牛乳房的检查和消毒工作，尽量减少相关微生物群对产品质量的影响。不动杆菌属为原乳中常见的嗜冷菌属，在天气寒冷和生产率降低的泌乳后期逐渐占据优势，但是其大量繁殖会影响乳成分，导致腐败变质［16-18］。在本论文中不动杆菌在年平均气温-2.6℃，-3.3℃的称多县歇武镇歇武村和曲麻莱长江村较低，而在年平均气温为3℃的玉树州果青牧场最高。

4 结论

综合来看，玉树地区不同样点的牦牛乳Alpha多样性指数observed_species差异显著(P＜0.05)，指数Chao1、Shannon和Simpson差异不显著(P＞0.05)，微生物多样性和丰度相似，微生物均以厚壁菌门和变形菌门为主。