乳制品是全世界共认的高营养食品，千百年来成为人民喜爱的食品，同时也为不同的外源微生物提供了理想的繁殖生境，形成了复杂多样的微生物群落，促进了发酵乳制品的形成[1-2]。酸马奶是其中一种非常著名的发酵乳制品，其以马奶为原材料与发酵剂混合，在开放的环境中自然发酵2 d～3 d即可，且具有略带酸味、沁人心脾、清爽适口等特性[3]。此外，由于酸马奶含有丰富的微生物资源及营养，使其对结核病为主的传染性疾病、消化系统疾病、心血管系统疾病等均有明显的辅助治疗作用[4-5]。但是，酸马奶的这些治疗作用与其化学成分和菌群结构密切相关。酸马奶的品质不仅取决于原料的差异，也与其微生物组成有密切关系。菌群结构容易受到地理和气候条件变化的影响，同时也受温度和发酵时间变化的影响。这些因素的差异导致了不同地区酸马奶的化学成分、菌群结构和风味的差异[6]。因此，了解传统酸马奶的微生物多样性至关重要。
以往对酸马奶微生物多样性的研究多采用传统分离鉴定方法，发现酸马奶中的主要微生物是乳酸菌和酵母菌，其对风味、酸度和口感特别重要[5，7]。鉴于一些微生物具有不可培养的特性，传统分离鉴定方法会严重低估微生物多样性。近年来，16S rRNA高通量测序技术被广泛应用于传统酿造食品中复杂微生物群落的研究，高通量测序不仅具有高通量、高覆盖、高准确率等特点，还可以全面地揭示样本微生物群落的组成及多样性等信息[8-9]。为了解传统酸马奶中细菌群落多样性，本研究利用Illumina Miseq高通量测序技术解析传统酸马奶的细菌群落结构及组成，为工业化生产提供理论依据及筛选潜在价值的益生菌提供数据参考。
1 材料与方法
1.1 材料与试剂
酸马奶：取自新疆塔城和丰县5个不同农场生产的样品（命名为 C1、C2、C3、C4、C5）和周边 5 户牧民手工酿造的样品（命名为 M1、M2、M3、M4、M5），取样后立即装入采样管中封口，用液氮冷冻后带回实验室置于-80℃储存，以备后续试验。
FastDNA ®Spin Kit for Soil试剂盒：美国 MP Biomedicals公司；Truseq®DNA PCR-Free Sample Preparation Kit试剂盒、Illumina Miseq PE250型测序仪：美国Illumina公司；Agencourt AMPure XP核酸纯化磁珠：美国Beckman Coulter公司；细菌16S rRNA V3-V4区 PCR 扩增引物（338F：5′-ACTCCTACGGGAGGCAGCA-3′和 806R：5′-GGACTACHVGGGTWTCTAAT-3′）：北京诺禾致源公司。
1.2 仪器与设备
Illumina Miseq PE250型测序仪：美国Illumina公司；S1000型聚合酶链式反应（polymerase chain reaction，PCR） 仪、CFX96 型荧光定量 PCR 仪、Gel DocTM XR+型凝胶成像系统、低温冷冻离心机：美国Bio Rad公司；DYY-8C型水平电泳槽：北京六一仪器厂；VORTEX-3型旋涡混匀器：德国IKA公司；NaneDrop 2000型超微量分光光度计：美国Thermo公司。
1.3 试验方法
1.3.1 酸马奶DNA提取及PCR建库
利用FastDNA®Spin Kit for Soil试剂盒提取酸马奶样本中微生物的基因组DNA。使用NaneDrop 2000型超微量分光光度计测定其溶度和纯度。将合格的DNA用引物338F和806R进行PCR扩增，PCR扩增体系（20 μL）：模板 DNA 10 ng，正向及反向引物各10 pmol，2×高保真 DNA聚合酶混合物10 μL。扩增条件为94℃预变性2 min；94℃变性20 s，55℃退火30 s，72℃延伸1 min，共进行25个循环；最后72℃延伸10 min[10]。扩增产物使用2%的琼脂凝胶进行电泳检测。加入等体积的Agencourt AMpure XP核酸纯化磁珠对产物进行纯化。利用Truseq®DNA PCR-Free建库试剂盒制备测序文库并添加接头序列，经Nane－Drop 2000确定文库质量，使用Illumina Miseq PE250测序平台进行高通量测序。
1.3.2 生物信息学分析
采用FLASH软件[11]对测序的原始数据进行拼接，利用QIIME软件[12]对低质量的拼接序列进行过滤，使用UCHIME软件[13]去除嵌合体，得到优质序列。根据UCLUST软件[14]在相似性97%的水平上对优质序列进行聚类，得到各可操作分类单元（operational taxonomic unit，OTU），以各OTU中丰度最高的序列作为代表序列。在80%置信水平下，通过Silva 132数据库[15]对OTU进行分类学注释，将蓝藻门和不能注释到门水平的OTU删除。基于OTU得到分析结果，采用随机抽平的方法，利用QIIME软件[12]对样品进行多样性分析。基于OTU的主坐标分析（principal co-ordinates analysis，PCoA）、样本聚类分析、箱线图绘制、热图绘制、火山图绘制均通过R软件实现。
2 结果与分析
2.1 样本测序数据处理结果
传统酸马奶样本测序数据处理结果统计见表1。
表1 传统酸马奶样本测序数据处理结果统计
Table 1 Statistics of sequencing data processing results in traditional koumiss
样品编号 原始序列条数有效序列条数平均序列碱基长度/bp样本GC含量/% Q30/%C1 41 635 32 647 465 54.63 92.64 C2 42 187 34 784 461 55.06 91.89 C3 40 685 29 975 448 53.14 92.55
续表1 传统酸马奶样本测序数据处理结果统计
Continue table 1 Statistics of sequencing data processing results in traditional koumiss
注：GC含量为鸟嘌呤和胞嘧啶所占的比率；Q30为质量值＞30的碱基占总碱基数的百分比。
样品编号 原始序列条数有效序列条数平均序列碱基长度/bp样本GC含量/% Q30/%C4 40 142 28 647 457 54.64 92.39 C5 42 659 31 567 456 54.24 92.78 M1 47 248 35 649 464 54.29 93.01 M2 48 362 36 438 459 55.07 92.68 M3 43 547 35 614 453 54.85 92.88 M4 42 374 34 547 461 54.62 92.76 M5 40 316 31 459 462 54.37 92.92
根据表1可知，10份酸马奶样品共检出429 155条原始序列，通过相关软件过滤后最终得到331 327条有效序列。去除条形码（barcode）和引物后，10份酸马奶样品的平均碱基长度均在450 bp以上，且GC含量为54.49%。此外，酸马奶样品的Q30值均高于90%，说明测序数据可靠且质量较好。
稀疏性曲线是判定样本的取样大小是否合理及确定测序深度是否可以覆盖整个微生物体系。酸马奶样品细菌群落的稀疏性曲线见图1。
图1 酸马奶样品细菌群落的稀疏性曲线
Fig.1 Rarefaction curves of bacterial communities in traditional koumiss
不同的稀疏性曲线代表不同的酸马奶样品。
由图1可知，随着测序深度的增加，被检测发现的OTU呈先迅速增加，随后平稳增加。当测序深度＞25000时，酸马奶样品中能检测到新OTU数量越来越少。这些结果表明虽然通过增加测序深度可以识别一些新的微生物，但在分析中已捕获了样品的大多数微生物，进一步说明Illumina Miseq高通量测序技术可有助于全面了解酸马奶的细菌多样性。
2.2 细菌群落的多样性分析
在97%的相似度阈值下，对331 327条有效序列进行聚类分析，根据QIIME软件进行OTU划分，并根据得到OTU进行α-多样性分析。酸马奶样品细菌Alpha多样性结果统计见表2。
表2 酸马奶样品细菌Alpha多样性结果统计
Table 2 Statistics of α-diversity results of bacteria in traditional koumiss
注：不同的小写字母代表相互之间差异显著（P＜0.05）。
样品编号 OTU数 Chao1指数 Shannon指数 Simpson指数C1 251 262 1.4 0.5 C2 190 209 1.52 0.51 C3 175 194 1.58 0.52 C4 226 234 1.46 0.49 C5 241 248 1.89 0.55 C（平均值） 217 229 1.57b 0.51b M1 178 199 1.99 0.6 M2 265 277 1.94 0.59 M3 232 246 1.8 0.57 M4 216 229 1.85 0.58 M5 232 248 1.34 0.5 M（平均值） 225 240 1.78a 0.57a
由表2可知，通过OTU聚类后可得到1 083个OTU，OTU平均为 175个/样本～265个/样本，其中农场生产样品的平均OTU数量为217个/样本，而牧民手工酿造样品的平均OTU数量为225个/样本，两者之间没有显著差异。Chao1指数是反映物种丰富度的指标，它通过观测到的结果推算出一个理论的丰富度，而这个丰富度更接近真实的丰富度，通常能观测到的物种丰富度肯定会比实际少。根据表2可知，10份酸马奶样品所检测到OTU数与Chao1指数的差异均不超过20，从而进一步说明当前测序深度是可以基本覆盖整个酸马奶样品的微生物体系。Shannon指数和Simpson指数是评价群落多样性和均匀度的常用指数。农场生产的5份酸马奶样品的Shannon指数和Simpson指数分别在1.4～1.89和0.49～0.55之间，其平均值分别为1.57和0.51；而牧民手工酿造的5份酸马奶样品的Shannon指数和Simpson指数分别在1.34～1.99和0.5～0.6之间，其平均值分别为1.78和0.57。牧民手工酿造的样品中Shannon指数和Simpson指数要显著高于农场生产的样品（P＜0.05），表明酸马奶的细菌多样性会受到地理和发酵器具的影响[16]。
为了更好地了解不同地域生产酸马奶的微生物结构的变化，采用主坐标分析法和非加权配对算术平均（unweighted pair group method with arithmetic mean，UPGMA）聚类法进行微生物群落多样性研究。酸马奶样品细菌群落结构的主坐标分析和非加权配对算术平均聚类分析见图2。
图2 酸马奶样品细菌群落结构的主坐标分析和非加权配对算术平均聚类分析
Fig.2 PCoA and UPGMA clustering analysis of bacterial communities in traditional koumiss
A为主坐标分析；B为非加权配对算术平均聚类分析；C代表农场生产；M代表牧民手工酿造。
如图2A所示，PCoA的第一主坐标和第二主坐标分别占总变异的93.01%和4.39%。农场生产的样品可明显聚集在一起，位于PCoA图的左上方，而牧民手工酿造的样品也可聚集在一起，位于PCoA图的右侧，表明两者的细菌群落结构存在明显差异。此结果与UPGMA聚类分析结果（图2B）相一致，表明不同来源的酸马奶的微生物组成易受到当地环境因素的影响。
2.3 细菌群落的结构分析
在划分OTU后，挑选不同OTU中丰度最高的序列作为该OTU的代表序列，随后与微生物参考数据库进行97%相似性比对分析，从而对酸马奶样品测序得到的每个OTU进行分类学鉴定。在分类学门水平揭示酸马奶样品细菌群落的结构差异。酸马奶样品细菌群落结构在门水平的分布见图3。
图3 酸马奶样品细菌群落结构在门水平的分布
Fig.3 Distribution of bacterial communities at phylum level in traditional koumiss
如图3所示，酸马奶样品的细菌群落分别来自厚壁菌门（Firmicutes）、变形菌门（Proteobacteria）、放线菌门（Actinobacteria）、异常球菌—栖热菌门（Deinococcus-Thermus）和拟杆菌门（Bacteroidetes）。其中厚壁菌门（Firmicutes）和变形菌门（Proteobacteria）占整个测序序列的99%左右。这与BAO等[17]的研究结构相一致，说明这些门水平的微生物决定酸马奶的成熟。此外，厚壁菌门在酸马奶中相对丰度可达80.19%～97.65%，说明这些来自厚壁菌门的微生物具有厌氧、耐酸和分解蛋白等特性。
为了进一步准确揭示酸马奶细菌群落的结构，故从分类学的属水平探究其组成，并揭晓不同来源样品的差异。酸马奶样品细菌群落结构在属水平的分布见图4。
图4 酸马奶样品细菌群落结构在属水平的组成和主要乳酸菌及醋酸菌的差异
Fig.4 Composition of bacterial communities at genus level in traditional koumiss and the difference of main lactic acid bacteria and acetic acid bacteria
A为属水平的分布；B为主要乳酸菌和醋酸菌的差异。C代表农场生产；M代表牧民手工酿造。箱线图是用矩形表示值分布的图形，矩形的末端标记第一个和第三个四分位数，最大值和最小值用触须标记。中间线标记中间值，异常值用点标记。**表示差异显著，P＜0.01；*** 表示差异极显著，P＜0.001。
如图4A所示，共鉴定出78个属，其中所有样品中平均相对丰度＞1%的属有6个，分别为乳杆菌属（Lactobacillus）（47.87%）、链球菌属（Streptococcus）（26.88%）、乳球菌属（Lactococcus）（4.16%）、醋杆菌属（Acetobacter）（3.7%）、明串珠菌属（Leuconostoc）（1.56%）和肠杆菌属（Enterobacter）（1.09%）。这些乳酸菌的丰富度促使厚壁菌门在酸马奶中具有显著优势。乳酸菌可能来源于发酵原料和自然环境，其可将原料中的小分子糖类物质转化为酸类化合物，从而赋予酸马奶柔和的酸感；部分糖类物质被转化为醇类化合物，这些醇类物质和酸类代谢物反应形成酯类物质，进而赋予产品具有酯香和醇香等特征[18]。此外，乳酸菌可通过蛋白水解系统分解发酵基质中的可溶性蛋白形成短肽和氨基酸等小分子代谢产物，赋予酸马奶独特的鲜味特征[18]。众所周知，乳杆菌属在发酵乳制品中起着重要作用[19]。根据图4A可知，乳杆菌属是绝对优势属，所有酸马奶样品均以乳杆菌属为主，其中占每个样品细菌总数的41.20%～55.82%。不同来源的酸马奶样品中乳杆菌属的相对丰度存在明显差异，其中农场生产酸马奶样品中乳杆菌属的平均相对丰度为42.2%，显著低于（P＜0.01）牧民手工酿造酸马奶样品中乳杆菌属的平均相对丰度（53.55%）（图4B）。上述结果与TANG等[20]的研究结果相一致，证实了乳杆菌属广泛存在酸马奶中，且在不同地域来源酸马奶样品之间存在差异。此外，链球菌属在酸马奶样品中相对丰度仅次于乳杆菌属，其占每个样品细菌总数的16.49%～39.76%，不同样品之间差异明显。农场生产酸马奶样品中链球菌属的平均相对丰度为36.07%，显著高于（P＜0.01）牧民手工酿造酸马奶样品中乳杆菌属的平均相对丰度（17.69%）（图4B），表明链球菌属和乳杆菌属在酸马奶发酵体系中可能存在互利共生的关系。这个与丁秀云[21]的研究结果相一致，此研究还发现嗜热链球菌与保加利亚乳杆菌组合发酵可增加胞外多糖的含量。如图4B所示，醋杆菌属和乳球菌属的相对丰度在不同地域来源酸马奶样品中差异较大，其在农场生产酸马奶样品的平均相对丰度分别为6.42%和0.42%，而在牧民手工酿造酸马奶样品的平均相对丰度分别0.98%和7.9%。相比牧民手工酿造酸马奶样品，农场生产酸马奶样品中醋杆菌属的平均相对丰度提高了5.5倍（P＜0.001），但乳球菌属的平均相对丰度却降低了95%（P＜0.01）。醋杆菌属是食醋和果醋中的一种重要微生物，由于其可产生乙醇脱氢酶，能够将酵母产生的酒精转化为乙酸[22]，也能在酸马奶中被检测到[17]。YIN等[23]发现添加谷氨酸可以增加醋杆菌属细胞生长及可培养细胞数。正是由于乳酸菌将原料中蛋白分解为短肽和氨基酸，从而促进醋杆菌的生长繁殖。乳球菌也被认为是一种重要的乳制品微生物，其可将多余的糖转化为甘露醇和葡聚糖，而这两种物质非常适合乳酸菌的发酵。上述结果表明不同地域采集的不同样品中细菌属及相对丰度含量存在显著差异，与GESUDU等[24]的研究结果相一致。根据传统酸马奶发酵工艺的总体情况，分析得出：（1）自制发酵剂、马奶原料、发酵环境和器具中的微生物是酸马奶微生物体系的来源；（2）酸马奶的品质和生产的工艺参数（如温度、发酵时间、搅拌强度等）可能会影响不同地区的酸马奶中细菌群落的变化。
2.4 差异微生物的分析
差异表达分析是比较两组不同样品基于表达异同的基本方法，可获得一组样品相对于另一组样品表达显著上调和下调，从而进一步研究这些差异表达基因的功能，多用芯片和转录组测序，但也可用于16S rRNA高通量测序技术探究不同类型样品之间微生物的差异，从而揭示其差异的原因。两组不同样品的差异微生物结果见图5。
图5 火山图显示了两个不同环境来源酸马奶样品出现不同程度的OTU富集
Fig.5 Volcano plot showing differentially abundant OTUs enriched in traditional koumiss from two different places
标记来自乳杆菌属、链球菌属、肠杆菌属和肠球菌属的差异微生物。
由图5可知，共发现16个差异显著的OTU，其中6个OTU可显著在农场生产酸马奶样品中富集 [2个来自链球菌属、2个来自醋杆菌属、1个来自乳杆菌属、1个来自罗斯氏菌属（Rothia）]，而剩下10个差异显著的OTU则显著在牧民手工酿造酸马奶样品中富集[1个来自肠杆菌属、1个来自肠球菌属（Enterococcus）、1个来自微小杆菌属（Exiguobacterium）、2个来自乳杆菌属、2个来自乳球菌属、1个来自明串珠菌属（Leuconostoc）、1个来自假单胞菌属、1个来自栖热菌属（Thermus）]。上述结果说明牧民手工酿造的酸马奶具有较高的细菌多样性，这也与表2的结果相一致。
3 结论
本研究利用Illumina Miseq高通量测序技术探究了两个不同来源酸马奶样品细菌群落的多样性及其组成。研究发现牧民手工酿造的样品中Shannon指数和Simpson指数要显著高于农场生产的样品（P＜0.05），且两者的细菌群落结构存在明显差异。厚壁菌门（Firmicutes）和变形菌门（Proteobacteria）是主要优势菌门，两者相对丰度之和可达98%。在属水平上，乳杆菌属、乳球菌属、醋杆菌属、明串珠菌属和肠杆菌属是酸马奶样品的优势菌属，尤其是乳杆菌属的相对丰度可达41.20%～55.82%。农场生产酸马奶样品中链球菌属和醋杆菌属的相对丰度含量要显著高于（P＜0.01）牧民手工酿造酸马奶样品。此外，通过差异分析发现两个不同来源酸马奶样品中有16个OTU存在显著差异，这些OTU中绝大部分都来自乳酸菌。这些结果有助于了解环境异质性对传统酸马奶细菌群落的影响，但本研究是在同一个城市采集的样品，未来可进一步探究酸马奶微生物多样性的生物地理分布模式。综上所述，本研究全面了解传统酸马奶的细菌种群分布，为工业化生产提供理论依据及筛选潜在价值的益生菌提供数据参考。
[1] QUIGLEY L,O′SULLIVAN O,BERESFORD T P,et al.Highthroughput sequencing for detection of subpopulations of bacteria not previously associated with artisanal cheeses[J].Applied and Environmental Microbiology,2012,78(16):5717-5723.
[2]QUIGLEY L,O′SULLIVAN O,STANTON C,et al.The complex microbiota of raw milk[J].FEMS Microbiology Reviews,2013,37(5):664-698.
[3] 布仁其其格,高雅罕,任秀娟,等.不同发酵时期酸马奶细菌群落结构[J].食品科学,2016,37(11):108-113.BURENQIQIGE,GAO Yahan,REN Xiujuan,et al.Dynamic changes of bacteria community structure during koumiss fermentation[J].Food Science,2016,37(11):108-113.
[4] 敖其,王丹,包英姝,等.酸马奶对慢性胃炎患者临床症状的改善效果研究[J].中国乳品工业,2018,46(5):14-17.AO Qi,WANG Dan,BAO Yingshu,et al.Effect of koumiss on relieving the chronic gastritis[J].China Dairy Industry,2018,46(5):14-17.
[5] 孔杭如,唐善虎,胡洋,等.酸马奶研究现状及进展[J].中国乳品工业,2016,44(5):32-35.KONG Hangru,TANG Shanhu,HU Yang,et al.Present situation and progress of koumiss[J].China Dairy Industry,2016,44(5):32-35.
[6] SUN T,WANG J.Analysis of chemical composition and microorganism flora of traditionally home-made koumiss in Xinjiang[J].China Dairy Industry,2005,33(10):9.
[7]SUN Z H,LIU W J,ZHANG J C,et al.Identification and characterization of the dominant lactobacilli isolated from koumiss in China[J].The Journal of General and Applied Microbiology,2010,56(3):257-265.
[8]SHOKRALLA S,SPALL J L,GIBSON J F,et al.Next-generation sequencing technologies for environmental DNA research[J].Molecular Ecology,2012,21(8):1794-1805.
[9] 王兴春,杨致荣,王敏,等.高通量测序技术及其应用[J].中国生物工程杂志,2012,32(1):109-114.WANG Xingchun,YANG Zhirong,WANG Min,et al.Highthroughput sequencing technology and its application[J].China Biotechnology,2012,32(1):109-114.
[10]ZHAO Z,MU T H,SUN H N.Microbial characterization of five Chinese traditional sourdoughs by high-throughput sequencing and their impact on the quality of potato steamed bread[J].Food Chemistry,2019,274:710-717.
[11]MAGO
T,SALZBERG S L.FLASH:fast length adjustment of short reads to improve genome assemblies[J].Bioinformatics,2011,27(21):2957-2963.
[12]CAPORASO J G,KUCZYNSKI J,STOMBAUGH J,et al.QIIME allows analysis of high-throughput community sequencing data[J].Nature Methods,2010,7(5):335-336.
[13]EDGAR R C,HAAS B J,CLEMENTE J C,et al.UCHIME improves sensitivity and speed of chimera detection[J].Bioinformatics,2011,27(16):2194-2200.
[14]EDGAR R C.Search and clustering orders of magnitude faster than BLAST[J].Bioinformatics,2010,26(19):2460-2461.
[15]QUAST C,PRUESSE E,YILMAZ P,et al.The SILVA ribosomal RNA gene database project:improved data processing and webbased tools[J].Nucleic Acids Research,2012,41(D1):D590-D596.
[16]GUO L,YA M,GUO Y S,et al.Study of bacterial and fungal community structures in traditional koumiss from Inner Mongolia[J].Journal of Dairy Science,2019,102(3):1972-1984.
[17]BAO L,BAO X,DAI Y,et al.Bacterial community succession and metabolite changes during the fermentation of koumiss,a traditional Mongolian fermented beverage[J].International Dairy Journal,2019,98:1-8.
[18]ABDEL-RAHMAN M A,TASHIEO Y,SONOMOTO K.Lactic acid production from lignocellulose-derived sugars using lactic acid bacteria:overview and limits[J].Journal of Biotechnology,2011,156(4):286-301.
[19]YAO G,HE Q,ZHANG W,et al.Single molecule,real-time sequencing technology improves the sensitivity in detecting bacteria from koumiss,a traditional fermented dairy product[J].Science Bulletin,2020,65(24):2065-2069.
[20]TANG H,MA H,HOU Q,et al.Profiling of koumiss microbiota and organic acids and their effects on koumiss taste[J].BMC Microbiology,2020,20(1):1-11.
[21]丁秀云.嗜热链球菌KLDS SM基因组学分析及比较基因组学分析[D].哈尔滨:东北农业大学,2017.DING Xiuyun.Genome and functional genomic analysis of Streptococcus thermophilus KLDS SM applied in dairy production[D].Harbin:Northeast Agricultural University,2017.
[22]陈洋,汪超,高冰,等.高耐受性醋酸菌的筛选及发酵特性研究[J].中国酿造,2015,34(12):28-33.CHEN Yang,WANG Chao,GAO Bing,et al.Screening and fermentation characteristics of acetic acid bacteria with ethanol and high temperature tolerance[J].China Brewing,2015,34(12):28-33.
[23]YIN H S,ZHANG R K,XIA M L,et al.Effect of aspartic acid and glutamate on metabolism and acid stress resistance of Acetobacter pasteurianus[J].Microbial Cell Factories,2017,16(1):1-14.
[24]GESUDU Q,ZHENG Y,XI X X,et al.Investigating bacterial population structure and dynamics in traditional koumiss from Inner Mongolia using single molecule real-time sequencing[J].Journal of Dairy Science,2016,99(10):7852-7863.