宏基因组学探究原料乳冷藏过程菌群变化规律

王媛媛，剧 柠，苟 萌，罗玉龙，李璞钰

（宁夏大学食品与葡萄酒学院，银川 750021）

0 引 言

原料乳因其营养全面、接近中性的pH值、较高的水分活性为微生物的生长提供了良好的环境。近年来，原料乳微生物多样性的研究逐步深入。研究显示，原料乳中的乳酸杆菌、嗜热链球菌能够促进乳制品的发酵；假单胞菌、梭菌、芽孢杆菌是导致乳品腐败的主要原因；李斯特菌、沙门氏菌、大肠杆菌等可引起食源性疾病；此外乳中还含有能够产生真菌毒素的真菌[1]。原料乳中微生物多样性受生理状态（如挤奶次数、内源性抗菌酶系统等）以及环境因素的影响有较大差异。研究表明健康奶牛与乳腺炎奶牛相比，其生产的原料乳具有更高的微生物多样性和丰富度[2]。有报道还说明原料乳微生物随季节及牧场维度的变化存在显著差异[3-4]。冷藏是抑制原料乳微生物繁殖的重要方法之一，但无法抑制其中的嗜冷微生物。嗜冷微生物通过合成磷脂来适应冷藏温度，同时水解蛋白质与脂肪使原料乳变色、产生异味、黏度增加、凝胶化、酸败等[5]。有研究利用遗传多样性分析探究了冷藏原料乳中与腐败相关的嗜冷菌，87.5%的菌株表现出蛋白水解和脂解活性[6]。现阶段对原料乳的研究主要集中在嗜冷菌多样性、菌体生物膜形成能力和产生的酶的特性等方面。有关原料乳嗜冷菌腐败潜力的研究中指出嗜冷菌的产酶、细胞生物膜形成是菌群数量增加的主要原因[7]。然而目前对原料乳冷藏期间的菌群变化规律鲜有报道。

自然界中 99%的微生物无法培养，随着测序技术的发展，依托16S rRNA高变区的高通量测序技术在乳品微生物领域得到普遍应用，但其测序深度和准确性决定了该技术对微生物的注释仅在属水平范围，无法满足对细菌更深层次的分析研究。宏基因组测序技术将微生物基因组随机打断重新组装成较长序列来探究微生物进化，能够将物种注释到种水平乃至菌株水平，并解析基因功能，成为更为有效的研究手段[8]。了解原料乳微生物及其随冷藏时间的变化情况，对生鲜乳保藏、液态乳灭菌控制、奶酪和乳粉加工及乳制品安全具有重要意义[9]。因此，本试验对 4 ℃冷藏原料乳中的微生物，利用 Illumina Hiseq宏基因组测序平台，研究贮藏期间细菌群落的变化规律，从微观揭示其品质变化的机理，从而为冷藏乳及乳制品加工提供理论依据。

1 材料与方法

1.1 原料乳收集及细菌计数

乳罐原料乳样品于2019年冬季12月自宁夏银川市当地牧场采样，该牧场原料乳除牧场自产外，还收集当地农户优质原料乳。牧场原料乳平均体细胞数＜20万个/mL，菌落总数＜10 万个/mL，蛋白质质量浓度＞3.0 g/100 mL。采样时无菌取样并保持4 ℃冷藏，2 h内运回实验室无菌分装至10 mL分装瓶继续4 ℃下冷藏。根据国家食品安全卫生标准对原料乳进行菌落总数测定、大肠菌群计数、嗜冷菌计数和乳酸菌计数。细菌计数冷藏时间点为 2、24、48、72 h，每个时间点取3个乳罐原料乳样品，共计12个样品。

1.2 DNA提取

宏基因组学通常会从 16S rRNA分析结果中选取适宜样本点进一步分析。前期16S rRNA预试验中，对物种组间差异进行2组比较，发现差异均不显著（P＞0.05）。于是结合细菌计数的方法筛选样本点。细菌计数结果显示组间两两比较具有差异性，个别24 h样本与48 h样本差异不显著，可舍去一个时间点。通过16S rRNA物种柱形图（图1），发现24 h样本为门水平物种丰度的转折点，因此宏基因测序时间点选择2、24、72 h。每个时间点取3个乳罐原料乳DNA样，共计9个样品，分组编号A0(A01/A02/A03)、A1(A11/A12/A13)、A3(A31/A32/A33)。采用试剂盒改良法进行乳样DNA的提取，具体如下：将10 mL均质原料乳样在5 000 ×g，4 ℃下离心10 min，去除上清液和脂肪，再加入3 mL无菌水，使用涡旋混匀器将无菌水和菌体沉淀混匀后，5 000 ×g，4 ℃下离心5 min，去除上清液和脂肪，加入2 mL磷酸缓冲盐溶液（Phosphate Buffer Saline，PBS），混匀后吸取2 mL菌液。使用 QIAGEN DNeasy PowerFood MicrobialKit（100）试剂盒，完成DNA抽提后，运用核酸浓度测试仪和1% 琼脂糖凝胶电泳检测DNA完整性。DNA检测合格后保存在-20 ℃以便后续分析。

图1 基于16S rRNA分析门水平细菌组成Fig.1 16S rRNA analysis of bacteria at phylum level

1.3 宏基因组测序

1.4 基因集构建与丰度计算

使用CD-HIT软件对预测基因序列进行聚类（参数设置为：95% identity、90% coverage），每类最长的基因序列用于构建非冗余基因集。使用SOAPaligner软件将每个样品的高质量reads 与非冗余基因集对比（95% identity），统计基因在对应样品中的丰度信息。丰度计算方法选用RPKM（Reads Per Kilobase Million）法。

1.5 物种与功能注释

使用BLASTP（BLAST Version 2.2.28+）将非冗余基因集序列与非冗余蛋白库（Non-Redundant Protein Sequence Database，NR）、直系同源蛋白分组比对数据库（evolutionary genealogy of genes：Non-supervised Orthologous Groups，eggNOG）（Version 4.5）、京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）进行比对（BLAST比对参数设置期望值e-value为1e-5），获得对应的物种注释与功能注释后计算丰度。

1.6 统计分析

原料乳细菌计数按照平均值±标准差计算。物种与功能分析将分组样本按照求均值进行计算，将丰度小于0.01的物种或功能合并为 others。使用 R语言工具进行Venn图绘制、群落柱形图绘制、箱式图绘制以及基于Bray-Curtis 距离的主坐标分析（Principal Co-ordinates Analysis，PCoA）。组间显著性差异检验运用单因素方差分析（One-way ANOVA）。LEfSe分析选用all-against-all多组比较策略，对样本按照不同的分组进行线性判别分析（Linear Discriminant Analysis，LDA），找出对样本划分产生显著性差异影响的功能。

2 结果与分析

2.1 原料乳细菌计数

原料乳中主要可培养微生物随冷藏时间的变化见表1，结果显示：原料乳冷藏72 h期间，菌落总数、嗜冷菌数、乳酸菌数、大肠菌群数均显著变化（P＜0.001），呈上升趋势。单因素方差分析显示嗜冷菌数、乳酸菌数、大肠菌群数在冷藏24 h与48 h无显著差异，菌落总数在冷藏24 h与48 h虽呈现显著差异，但数量相近。说明24 h与48 h两样本点结果差异较小。另外，各细菌计数冷藏48 h至72 h增幅较大。72 h可做为冷藏期间微生物控制的重要时间点。

表1 原料乳4 ℃冷藏下细菌计数结果Table 1 Bacterial count of raw milk under refrigeration at 4 ℃

2.2 宏基因组测序和质量控制

本研究对3个时期（2，24，72h）的样本A0、A1、A3（n=3）共计 9个样品进行了宏基因组测序，测序与组装信息见表2，结果显示平均每个样本产生了43 291 358个总长度为6.5×109bp的原始序列，经过质控处理后，干净序列（clean reads）占原始序列的98%以上。拼接组装得到contigs数平均为1 003 897，通过基因预测，平均ORF序列数为398 999。以上数据保证了结果的完整性与有效性。

表2 宏基因测序与组装信息Table 2 Metagene sequencing and assembly information

2.3 物种群落多样性分析

原料乳4 ℃冷藏样本物种分布Venn图（图2）显示：属水平上共注释到589个菌属，种水平上共注释到1 668个菌种。表3为基于种水平的Alpha多样性分析结果，冷藏期间Chao1指数与ace指数呈上升趋势，说明物种随冷藏时间的增加而增加；Shannon指数与Simpson指数综合菌群组成丰富度与均匀度反映了样本的多样性，数据显示样本多样性A0大于A3大于A1，这可能因为原料乳冷藏24 h后不耐冷微生物的丰富度大幅降低，冷藏72 h后嗜冷微生物的大量繁殖使物种多样性增加。另外，各样本基于种水平的PCoA分析结果显示（图3），组内样本点较聚集，差异较小，各样本组间区分较明显，组间差异较大。说明分组有意义，物种在原料乳冷藏期间有显著变化。

表3 基于菌种水平的Alpha多样性Table 3 Alpha diversity based on species level

2.4 物种分类学注释

2.4.1 属水平菌群组成分析

本研究注释到相对丰度大于 0.01的菌属有 11个（图4a）。A0样本中优势菌属有不动杆菌属（Acinetobacter）（55.61%）、链球菌属（Streptococcus）（11.53%）、无浆体属（Anaplasma）（8.59%）、梭菌属（Clostridium）（2.40%），总占比达78.13%。A1样本中相对丰度最高的是不动杆菌属（93.97%）；其次为假单胞菌属（Pseudomonas）（1.18%）、黄杆菌属（Flavobacterium）（0.54%）、链球菌属（0.41%），总占比达 96.10%。A3样本中的优势菌属为不动杆菌属（54.21%）、黄杆菌属（18.26%）、假单胞菌属（13.57%）、乳球菌属（Lactococcus）（4.54%），总占比达90.58%。单因素方差分析结果显示（图4b）以上菌属皆为引起组间差异的差异菌属（P＜0.05）。随着冷藏时间的延长，不动杆菌属相对丰度先上升后下降，以冷藏24 h后最高；链球菌属、无浆体属、梭菌属相对丰度大幅下降；黄杆菌属、假单胞菌属、乳球菌属相对丰度均大幅上升。原料乳冷藏72 h内的菌属由不动杆菌属、链球菌属、无浆体属和梭菌属向黄杆菌属、假单胞菌属和乳球菌属逐渐演变。

图2 原料乳4 ℃冷藏样本物种分布Venn图Fig.2 Venn diagram of bacteria distribution of raw milk refrigerated samples at 4 ℃

图3 种水平细菌主成分分析Fig.3 PCoA of bacteria at species level

有关原料乳微生物多样性的报道中优势菌属及组成不尽一致。Hahne等[7]使用16S rRNA测序技术对冷藏原料乳进行物种多样性分析，其结果与本研究较一致，即不动杆菌属、链球菌属、假单胞菌属为原料乳主要菌属，其中不动杆菌属相对丰度最高。Porcellato等[3]指出寒冷季节的原料乳中假单胞菌属含量最多，其次为不动杆菌属、链球菌属、乳球菌属和黄杆菌属。这与本研究注释到的菌属种类较一致，但菌属占比有所差异。Tilocca等[9]指出原料乳中典型的菌属组成为：乳球菌属、链球菌属、乳杆菌属、明串珠菌属、肠球菌属、黄杆菌属、假单胞菌属、不动杆菌属、气单胞菌属。本研究中，注释到的明串珠菌属、气单胞菌属相对丰度远小于1%。可见，原料乳中的优势菌属组成大致相同，相对丰度却有所差异，产生差异的原因可能与原料乳挤出时的初始环境、乳本身营养成分及冷藏时间有关。

冷藏 72 h期间原料乳中的微生物菌属显著改变（P＜0.05）。链球菌属和梭菌属可能对原料乳体系的抗菌肽的作用比较敏感，如乳铁蛋白、溶菌酶、过氧化物等[10-12]，其相对丰度在冷藏24 h后大幅度降低，即链球菌属相对丰度从11.53%降至0.41%，梭菌属相对丰度从2.4%降至0.03%。无浆体属是革兰氏阴性严格寄生性致病菌，在牛的红细胞中生长繁殖[13]，Zhang等[14]从健康山羊分泌的乳汁中检测到该菌属。本研究中，无浆体属由于对环境的适应性较差，在冷藏24 h后相对丰度大幅降低（8.59%～0.32%）。与此同时，不动杆菌属由于其耐受低温的特性，大量的繁殖，成为相对丰度最高的菌属。这与Quigley等[1]指出不动杆菌在原料乳冷藏期间能够快速生长的结果相吻合。原料乳中的营养物质被微生物不断消耗，同时产生一系列代谢产物，这些代谢产物同样可以提供给微生物生长所用[9]。本研究中原料乳冷藏72 h，假单胞菌属、黄杆菌属、乳球菌属得到了足够的生长所需物质，成为优势菌属。

图4 属水平细菌组成分析Fig.4 Analysis of bacteria composition at genus level

2.4.2 种水平菌群组成分析

本研究注释到丰度大于0.01的菌种有16个（图5a）。A0样本中优势菌种有不动杆菌TTH0-4（Acinetobactersp.TTH0-4）（41.98%）、嗜吞噬细胞无浆体（Anaplasma phagocytophilum）（8.57%）、肺炎链球菌（Streptococcus pneumoniae）（6.06%）、变形链球菌（Streptococcus mutans）（5.45%），总占比达62.06%。A1样本中丰度最高的为不动杆菌 TTH0-4（73.02%），其次为鲍曼不动杆菌（Acinetobacter baumannii）（2.54%）、约氏不动杆菌（Acinetobacter johnsonii）（1.49%）、波希不动杆菌（Acinetobacter bohemicus）（1.36%），总占比达78.41%。A3样本中丰度最高的为不动杆菌 TTH0-4（41.63%），其次为寒冷黄杆菌（Flavobacterium frigidarium）（15.24%）、荧光假单胞菌（Pseudomonas fluorescens）（9.77%）、鱼乳球菌（Lactococcus piscium）（3.29%），总占比达69.93%。单因素方差分析结果（图5b）显示以上菌种皆为引起组间差异的菌种（P＜0.05）。所检测到的优势菌种变化趋势与相应优势菌属对应。随冷藏时间的延长，乳中菌种由不动杆菌TTH0-4、嗜吞噬细胞无浆体、肺炎链球菌和变形链球菌向寒冷黄杆菌、荧光假单胞菌和鱼乳球菌逐渐演变。

图5 种水平细菌组成分析Fig.5 Analysis of bacteria composition at species level

值得注意的是，除此前研究报道的原料乳中主要致腐菌荧光假单胞菌外，本研究检测到的不动杆菌TTH0-4、寒冷黄杆菌和鱼乳球菌鲜有报道。不动杆菌属携带多种抗生素耐药因子，在食品安全中受到了较高的重视[3,15]。同时不动杆菌属具有脂解或蛋白水解能力，也能产生生物膜[6]，具有一定的原料乳腐败潜力。本研究中，种水平上丰度最高的不动杆菌TTH0-4在原料乳中尚未见报道，此前曾被 Zhang等[16]分离自青藏高原的多年冻土区，并说明其能够很好适应寒冷的环境。该菌的检出可能与牧场的土壤环境有关。本研究中另一优势菌种，寒冷黄杆菌属于黄杆菌属，其对乳及乳制品的影响除了产生蛋白酶外，还产生磷脂酶C，可降解乳脂球膜的磷脂，继而改变原料乳苦味成分，对干酪的产量也具有一定影响[17]。该菌过氧化氢酶试验呈阳性[18]，其分解过氧化氢的能力可能使原料乳中乳糖氧化酶系统因过氧化氢含量的降低而受到影响甚至不被激活[12]。原料乳中的乳球菌属能将乳糖分解产生乳酸，能够参与蛋白质水解将氨基酸转化为风味化合物，也有助于脂肪代谢[19]，与乳品品质息息相关。本研究中注释到的鱼乳球菌曾被Carraro等[20]在原料乳中检出，该菌更倾向于在较低温度下生长，但是它如何参与原料乳的变化过程还未见明确的报道[21]。本研究也检测到了原料乳中普遍存在的荧光假单胞菌，它能够产生APRX碱性金属蛋白酶，该酶不仅可以降解具有类凝乳酶活性的 k-酪蛋白导致乳品在储存期间发生凝胶化，而且参与营养物质的利用，促进荧光假单胞菌的生长[22-24]。荧光假单胞菌还具有生物膜形成能力，增加该菌灭菌难度的同时，生物膜胞外聚合物（Extracellular Polymeric Substances，EPS）的保护机制会进一步增强酶的热稳定性[25]。因此，荧光假单胞菌是控制原料乳品质的关键菌种。

2.5 功能注释

为深入了解微生物演替产生的原因及对冷藏原料乳的影响，进行了微生物群落的COG功能注释（图6）。注释结果中丰度排名第一位的是S（功能未知）。LDA判别结果（图7）显示，S在A0样本中丰度存在显著差异，lg（LDA值）为 4.99，造成该结果的原因可能是原料乳冷藏初期未知微生物种类丰富，也可能由于冷藏初期微生物的功能较复杂，COG注释不足以涵盖所有功能基因。丰度排名其次为 L（复制重组修复）、J（翻译、核糖体结构与生物发生）、M（细胞壁/膜/包膜生物发生）。这说明乳环境中微生物群落为了能够快速适应变换的环境，抗逆性能增加[26]。众多微生物中发现的冷休克蛋白能够感知低温环境，作为RNA的分子伴侣可促进转录与翻译[27]。双组分系统同样能够使微生物快速适应多变的环境，其组成部分应答调控蛋白能够与DNA、RNA等结合以对转录或代谢过程进行调控，与组氨酸激酶的协同能够调控细菌生物膜的形成[28]。另外，微生物细胞壁、膜的产生与变化也能够感知冷刺激[29]。L、J、与M在A1样本中丰度存在显著差异，lg（LDA值）分别为 3.81、4.16、4.09。说明原料乳在冷藏24 h后微生物抗逆相关功能丰度较高。不动杆菌属在此阶段为最优势物种，表明其对多变环境的适应能力相对较强。邵毅等[30]研究也指出不动杆菌耐环境胁迫的能力较强。

原料乳中微生物进行着复杂的新陈代谢活动。与代谢相关功能丰度最高的为E（氨基酸转运与代谢），此外I（脂质转运与代谢）与 G（碳水化合物转运与代谢）也有注释。其他被注释的还有无机离子、辅酶、核苷酸、次级产物的转运与代谢等。体现出了较为全面的代谢系统。LDA判别结果显示，I在A1样本中丰度存在显著差异，lg（LDA值）为 3.98。不动杆菌属在此期间相对丰度较高，研究表明不动杆菌属能够利用碳源合成脂质，并使三脂酰甘油、游离脂肪酸等积累[31]。因此推测不动杆菌属对原料乳碳源及脂质的组成影响较大。A3样本中丰度存在显著差异的功能为E和G，lg（LDA值）分别为4.10、3.90。此时原料乳菌落总数、嗜冷菌数、乳酸菌数、大肠菌群数均显著增加，黄杆菌属、假单胞菌属、乳球菌属丰度占比上升。可能原因是这些微生物在此阶段充分利用原料乳中的碳源与氮源生长繁殖，同时，黄杆菌属、假单胞菌属产生的蛋白酶将蛋白质分解为氨基酸。有研究表明半胱氨酸、蛋氨酸中含硫，分解产生的硫化物浓度过高会使原料乳产生臭味[32]。亮氨酸能够脱氨脱羧生成异戊醛使原料乳产生麦芽气味[33]。说明黄杆菌属与假单胞菌属可能与原料乳风味的形成密切相关。另外，原料乳中乳糖经乳球菌属利用分解产生葡萄糖与半乳糖进而经糖酵解途径转换为乳酸或其他有机酸[34]，当原料乳酸度过高时，不耐酸微生物生长受到抑制，同时酪蛋白、乳球蛋白等可能发生沉淀，使原料乳出现凝块甚至乳清析出的现象。

图6 COG功能注释Circos图Fig.6 Circos diagram of COG function annotation

通过KEGG功能注释验证上述分析。结果与COG功能注释结果较一致（图8）。即遗传信息处理和碳水化合物、氨基酸、脂质代谢为主要功能。将氨基酸代谢、脂质代谢、碳水化合物代谢与主要优势微生物进行相关性分析（图9），结果表明，脂质代谢与不动杆菌的相关性更显著（P＜0.001），氨基酸代谢和碳水化合物代谢与假单胞菌的相关性更显著（P＜0.01）。该结果为上述COG的功能分析提供有效支撑。总而言之，原料乳在冷藏过程中抑菌物质减少，微生物利用蛋白质、脂肪、碳水化合物等营养物质得以大量繁殖。原料乳中的蛋白质与脂肪被分解产生醛、酮、醇、硫化物、苦味肽等代谢副产物[35-37]，导致乳成分发生改变，乳成分的改变又影响着微生物生长。可见微生物演替与乳环境营养物质的损耗以及代谢物的积累密切相关。

图7 COG功能差异判别分析Fig.7 Linear Discriminant Analysis of COG function

图8 KEGG功能注释热图Fig.8 Heat map of KEGG function annotation

图9 主要代谢功能与优势微生物相关性分析Fig.9 Correlation analysis between main metabolic functions and dominant microorganisms

3 结 论

利用宏基因组测序技术对原料乳 4 ℃冷藏过程中微生物群落变化及其功能特征进行研究。结果表明：

1）冷藏过程中微生物有显著的演替现象。其中不动杆菌TTH0-4、寒冷黄杆菌、鱼乳球菌在原料乳中的报道较少，需进一步研究。

2）原料乳冷藏前期，微生物通过L（复制重组修复）、J（翻译、核糖体结构与生物发生）、M（细胞壁/膜/包膜生物发生）等功能增加对多变乳环境的抗逆性。

3）不动杆菌属与脂质代谢显著相关（P＜0.001），对原料乳冷藏前期脂质的组成影响较大。假单胞菌属与碳水化合物代谢、氨基酸代谢显著相关（P＜0.01），对原料乳冷藏后期蛋白分解及风味形成影响较大。

下一步可通过代谢组学与宏基因组学进行关联分析，进一步发现原料乳冷藏过程菌群变化的机理。