基于高通量测序的液态奶中微生物多样性研究

屠大伟,张清平,刘美艳,张巨凤,翁盈秋

1.重庆工商大学 环境与资源学院，重庆 400067 2.重庆市天友乳业股份有限公司，重庆 400015 3.重庆万标检测技术有限公司，重庆 400714

1 引 言

牛乳因其营养丰富且易被人体吸收而成为人民日常生活中重要的营养来源。牛乳中含有丰富的营养物质，但在挤奶、原料乳的预处理、贮存以及运输等过程中处理不当会导致奶源中有害微生物的生长和繁殖，从而造成牛乳中营养物质的破坏，影响原料乳的品质以及食源性疾病的传播[1]。

目前，常用于研究乳制品中微生物多样性的方法主要有常规细菌鉴定技术、免疫学检测技术、基因检测技术和全基因组测序技术[2]等，传统的检测方法存在操作复杂，检测灵敏度低的现象[3]。与传统方法相比，高通量测序技术可以高效、准确、全面地分析出不同被测样品的微生物菌落组成结构，同时可以分析一些不可培养且丰度低的微生物[4]。基于此，本文采用高通量测序技术。目前，16S rDNA测序技术已成为研究微生物菌群结构及其多样性的重要方式[4-8]。张敏[33]等采用16S rDNA高通量测序方法比较了新疆西北部地区乳制品中微生物的多样性，通过对新疆2个地区的7种乳制品进行高通量测序分析发现，7种乳制品在门水平的优势菌门是厚壁菌门和变形菌门，而酸奶和原奶中的丰度不同，属水平上的差异较大。同时还发现不同动物来源的原奶和酸奶的微生物多样性也存在着显著的差异。姚宇秀[3]等采用16S rDNA高通量测序技术对原料乳及其采集过程中的相关设备以及奶牛的不同部位进行样本的分析，结果发现，输入管路外侧和盛装原料奶的空贮罐中有金黄色葡萄球菌的存在。布仁其其格[9]等通过16S rRNA基因序列对酸马奶传统发酵过程中细菌群落结构演替变化，结果表明，在发酵初期细菌多样性最高，且在72 h时细菌丰度最高，门水平的优势菌门为厚壁菌门和变形菌门，属水平的优势细菌属为乳杆菌属。

本研究的目的是利用高通量测序方法对重庆地区的不同液态奶的细菌菌群结构和组成进行分析比较，以期探索液态奶中的腐败菌，实现对液态奶中有害菌的防治，从而实现提高液态奶的品质。

2 材料与方法

2.1 材料与仪器

液态奶，重庆市天友乳液股份有限公司的4个牧场无菌生理盐水，北京陆桥技术股份有限公司；植物基因组DNA提取试剂盒(TSP101-200)，TSINGKE；通用引物16S，TSINGKE；I-5TM2×High-Fidelity Master Mix(TP001)，TSINGKE；高纯度低电渗琼脂糖(TSJ001)，TSINGKE；DNA凝胶回收试剂盒(GE0101-200)，TSINGKE；DL2000 DNA Marker(TSJ011-500)，TSINGKE。

VITEK2 COMPACT全自动细菌鉴定及药敏分析系统，法国生物梅里埃公司； XH-B旋涡混合器，常州翔天实验仪器厂；BSC-1100A2-X生物安全柜，济南鑫贝西生物技术有限公司；DK-98-IIA恒温水浴锅，北京市泰和润仪器有限公司；DW-86W50超低温冰柜，浙江捷盛制冷科技有限公司；PD100A高压灭菌锅，致微(厦门)仪器有限公司； DHP-9272电热恒温培养箱，上海一恒科学仪器有限公司；高通量测序仪Illumina Hiseq 2500，美国因美纳公司；离心机Legend Micro17，赛默飞世尔；电泳槽JYDF，北京君意东方电泳设备有限公司；PCR仪2720 thermal cycler，美国应用生物系统公司；电泳仪JY300C，北京君意东方电泳设备有限公司；水浴锅DFD-700，北京中兴伟业；凝胶成像仪JY04S-3C，北京君意东方电泳设备有限公司；板式离心机L550，cence湘仪。

2.2 方法

2.2.1样品采集

采集来自天友乳业重庆周边4个不同牧场的原料乳(Sample 1、Sample 2、Sample 3、Sample 4)。所有的样品经无菌容器密封好后放入车载冰箱运回实验室，样品保存在-80℃超低温冰箱，待试验。

2.2.2DNA的提取、扩增和测序

取40 mL 原料乳样本，6 000 r/min下4 ℃离心20 min，用灭菌脱脂棉去掉乳脂层，弃去上清液，加入无菌生理盐水至10 mL，10 000 r/min高速冷冻离心10 min，用1 mL蔗糖缓冲液冲洗两次，后悬浮400 μL蔗糖缓冲液及2U变溶菌素和800 μg溶菌酶，37 ℃培养1 h。选用806R(5’-GGACTACHVGGGTWTCTAAT-3’)和338F(5’-ACTCCTACGGGAGGCAGCAG-3’)对细菌16S rRNA V3-V4区进行PCR扩增。

PCR反应体系为：12.5 μL I-5TM2×High-Fidelity Master Mix，1 μL 10 mmol/L Primer A，1 μL 10 mmol/L Primer B，加入模板DNA，补ddH2O至25 μL。用琼脂糖凝胶电泳检测，利用凝胶回收试剂盒切胶回收PCR产物。

PCR扩增程序为：98 ℃预变性2 min；循环数(25～35次)×(95 ℃变性10 s，45～68 ℃退火10～15 s，72 ℃延伸5～15 s)；72 ℃终延伸1～5 min；4 ℃冷却。

使用Illumina Hiseq 2500测序，由北京擎科生物科技有限公司重庆分公司完成。

2.3 数据处理

基于Illumina Hiseq 2500测序平台，使用Flash 软件对每个样品的PE Reads进行拼接，得到的序列即为原始序列，再通过Trimmomatic软件对原始序列进行优化处理，得到较高质量的优化序列，最后通过UCHIME软件鉴定并去除嵌合体序列，得到最终的有效序列。

MH7A细胞为永生化的RA关节滑膜细胞，购自广州吉妮欧公司。用高糖DMEM+10%胎牛血清于5%CO2，37℃的培养箱培养。购买时为第3代，培养3代后用于后续实验研究。

3 结果与分析

3.1 测序数据统计及结果分析

为了获得准确、高质量的生物信息分析结果，通过测序获得的原始数据需要进行优化处理[10]。如表1所示。

表1 样品测序数据统计Table 1 Sample sequencing data statistics

从表1中可以看出，4个样品总计测得原始序列条数为242 576条，将低质量的序列过滤后得到的优化序列数总计233 818条，再将上述序列进行冗余处理，最终得到的有效序列数为212 260条。

图1为有效序列的长度分布图，从图1中可以看出，大部分的序列均分布在400～450 bp之间，且430～440 bp之间的序列最多。从序列长度的分布可以看出，与16S rDNAV3-V4区序列长度大致吻合[11]。

图1 有效序列长度分布Fig.1 Distribution of effective sequence length

从样本中随机抽取一定数量的序列后，对获得的序列所代表的物种数目进行统计，以物种数目和序列数构建的一条曲线为样本的稀释曲线[12]。稀释曲线可以验证测序数据量是否可以反映样品中的物种多样性，并间接反映样品的测序深度和物种的丰富度[12-13]。若在一定范围内，随着序列数的增加，曲线表现出急剧上升的趋势则反映样本中仍有大量物种被发现；当稀释曲线趋于平缓则表示样本量足够充分[14]。图2(a)是样品中稀释曲线图，从该图中可以看出，随着样本序列数的增加，每个样品中OTU数目增加，当OTU数目到达一定数量后，曲线趋于平缓，这表明实验中所得序列数基本可以反映出每组样品中的菌群结构。同时，可以看出Sample 4中OTU数目最多，Sample 1次之，Sample 2和Sample 3中OTU数目相对较少。

香农指数曲线可以体现出各样品在不同测序数量时样本中微生物多样性的指数，随着样本序列数的增加，曲线呈上升趋势[30]；当曲线趋于平坦时，说明测序数据量足够大，OTU种类不会在随着序列量的增长而增长[31]。图2(b)是样品中香农指数曲线，从该图中可以看出，随着测序样本数量的增加曲线趋于平缓，这说明本研究的测序深度是可以满足分析要求[15]。

(a) 样品稀释曲线

3.2 样品OTU统计分析

一般情况下，将相似性高于97%的不同16S rDNA序列定义为一个OTU，每一个OTU通常被视为一个微生物物种[32]。图3展示了不同样品的OTU数。从图3中可以看出，各样品的OTU数相差不大，均在70～85个之间，Sample 1～Sample 4各样品中OTU个数分别为80，73，79，81；四个样本中总共的OTU数目为119个。

图3 样品中OTU数目的分布图Fig.3 Distribution of the number of OTUs in different samples

Venn图体现了不同样品之间微生物群落的相似性和差异性[4]，如图4所示。从图4中可以看出，4个样品中共有OTU数目为43个，每个样品特有的OTU个数分别为：12，5，7，10，分别占各样本中OTU数量的15%、6.85%、8.86%和12.35%；Sample 1与Sample 2、Sample 1与Sample 3、Sample 2与Sample 4、Sample 3与Sample 4样品之间共有，且其余样品不具有的OTU数目较少，分别为2,3,2,3。因此，可以看出各样品之间微生物群落没有显著的差异性。

图4 样品中OTU数目的Venn图Fig.4 Venn diagram of the number of OTUs in different samples

3.3 样品间α多样性的分析

α多样性分析体现了单个样品内部物种的丰富度和多样性[16]，通常可以通过ACE指数、Chao1指数、Simpson指数以及Shannon指数等反映[17,18]。其中，ACE指数和Chao1指数用来反映微生物群落的丰富度[18]，ACE指数和Chao1指数越大表示样品中微生物群落的丰富度越大[32]；Shannon指数和Simpson指数都是用来反映微生物群落的多样性[30]，微生物群落的多样性与Shannon指数呈正相关，与Simpson指数呈负相关[31]。

对不同样本在97%一致性阈值下的α多样性指数进行统计和分析[4,19-20]，结果见表2。从表2中可以看出，所有样品的覆盖率(Coverage)均大于0.99，说明本次研究所建立的数据库可信度较高，数据的真实性较好，能够有效地反映样品中微生物菌群的多样性。Sample 2样品中ACE指数和Chao1指数最高，表明该样品中菌群的丰富度较高；而Sample 3中ACE指数和Chao1指数最低，说明该样品中菌群的丰富度较低。同时，Sample 4样品中其Simpson指数最小，Shannon指数最大说明该样品中菌群的多样性最大；Sample 1次之，Sample 2样品中菌群的多样性最低。

表2 样品间α多样性指数统计Table 2 Statistics of α diversity index among samples

3.4 液态奶样品间基于门水平的菌群结构分析

图5是液态奶各样品中基于门水平的菌群结构分布比例图。从图中可以看出，液态奶样品中的10个细菌门被鉴定出，分别是异常球菌-栖热菌门(Deinococcus-Thermus)、厚壁菌门(Firmicutes)、硝化螺旋菌门(Nitrospirae)、变形菌门(Proteobacteria)、蓝藻细菌门(Cyanobacteria)、拟杆菌门(Bacteroidetes)、放线细菌门(Actinobacteria)、Patescibacteria、绿弯菌门(Chloroflexi)和Epsilonbacteraeota。从图5中可以看出各样品之间的优势菌门基本相同，但相对丰度差异较大。Sample 1中的主要菌群是厚壁菌门、变形菌门和拟杆菌门分别占38.9%、46.9%和14.0%；Sample 2中的主要菌群是厚壁菌门、变形菌门和放线细菌门分别占73.8%、10.0%和13.1%；Sample 3中的主要菌群是厚壁菌门、变形菌门和拟杆菌门分别占20.9%、28.0%和50.1%；Sample 4中的主要菌群是放线细菌门、厚壁菌门和异常球菌-栖热菌门分别占30.0%、29.6%和22.8%。

图5 基于门水平各样品的菌群结构分析Fig.5 Analysis of flora structure of each sample based on phylum level

由此可见，每组液态奶中门水平的菌群结构差异性不大，厚壁菌门和变形菌门是四组液态奶的共有菌门，这与张敏等[33]的研究中乳制品的菌群结论相同。

物种丰度聚类热图分析是通过颜色变化与相似程度来反应二维矩阵或表格中的数据信息，并呈现出群落物种的组成信息，表明了液态奶不同样品间不同细菌门的相对丰度及细菌组成的差异性和样品间的相似性[10]。采用物种丰度聚类热图分析液态奶样品中含量前10个菌门和4个样品之间的交互关系，见图6。

图6 基于门水平各样品的物种丰度聚类热图Fig.6 Species abundance clustering heatmap for each sample at phylum level

从图6中可知，不同样品间门水平具有一定的丰度，且各样品的菌群组成具有一定的差异性和相似性。同时，根据丰度聚类热图可以看出样品间门水平的聚类关系，即4组样品均聚在一起，说明四组样品间的丰度相似，其中Sample 2和Sample 4的相似度较高些，Sample 1与其余三组在门水平上的相似度较低。

3.5 液态奶样品间基于属水平的菌群结构分析

图7是液态奶各样品中基于属水平的聚类树柱状图。从图中可以看出，液态奶中的10个菌属信息被鉴定出，分别是Olivibacter、克雷伯氏菌属(Klebsiella)、微小杆菌属(Exiguobacterium)、不动杆菌属(Acinetobacter)、异常球菌属(Deinococcus)以及短波单胞菌属(Brevunmdimonas)，寡养单胞菌属(Stenotrophnomonas)、微杆菌属(Microbacterium)、Cloacibacterium、芽孢杆菌属(Bacillus)和其他未分类菌。

图7 基于属水平各样品的聚类树柱状图Fig.7 Cluster tree histogram of each sample based on genus level

各样品在属水平上优势菌属分别是：Sample 1中的优势菌属是厚壁菌门的芽孢杆菌属(占35.9%)及变形菌门的短波单胞菌属(占26.2%)；Sample 2中的优势菌属是厚壁菌门的芽孢杆菌属(占73.7%)；Sample 3中的优势均属是拟杆菌门的Cloacibacterium(占49.9%)、变形菌门的寡养单胞菌属(占27.6%)；Sample 4中的优势菌属是厚壁菌门的芽孢杆菌属(占29.2%)和微杆菌属(占28.3%)、异常球菌-栖热菌门的异常球菌属(占22.8%)。从图7中可以看出，Sample 1和Sample 4的细菌群落均匀性相对较高，Sample 2的细菌群落组成高度集中在芽孢杆菌属中。

根据已有的研究报道表明，不动杆菌属、芽孢杆菌属及克雷伯氏菌属属于嗜冷菌[21-24]。嗜冷菌在低温存储过程中，其依旧能够繁殖并能够产生影响原料乳品质的腐败酶：蛋白酶和脂肪酶等，破坏原料乳中的营养物质，影响原料乳的品质[21,25-28]。研究表明：生鲜乳本身不存在大量的嗜冷菌，其主要来源于生产环节的污染，导致原料乳在运送过程中出现大量嗜冷菌繁殖从而导致原料乳的品质遭到破坏[22,29]。

4 结论

本实验通过16S rRNA基因的高通量测序技术对来自4个牧场液态奶的细菌群落结构进行研究。通过对比发现，在门水平上，4组液态奶的菌群结构基本相似，但相对丰度差异较大，4组液态奶的优势菌门主要是厚壁菌门和变形菌门。在属水平上，4组液态奶的优势菌属和相对丰度都不相同，存在显著差异。同时，还发现属水平上芽孢杆菌属、不动杆菌属属于嗜冷菌，在液态奶的贮存中容易分解其中的蛋白质、脂肪或碳水化合物，从而导致液态奶出现腐败变质等现象。这为后续对探索液态奶中腐败菌的来源及防控提供了有力的研究依据。