传统发酵乳制品是经过长期自然驯化优良食品微生物的天然载体，发酵乳制品中的微生物决定了发酵乳制品的品质[1-2]。近年来，有关传统发酵乳制品中微生物多样性的研究，既有传统的分离培养技术研究，也有高通量测序、宏基因组分析等非培养技术研究[3-5]。利用传统的分离技术虽然从样品中获得了一些微生物，但只能代表具体培养条件所能分离到的微生物，无法分离到样品中所有微生物，不能完整地揭示样品的菌群组成[6-9]。高通量测序[10-13]、宏基因组分析[14-16]等非培养技术研究无需分离培养微生物，能够更好地揭示样品中菌群组成及多样性[17-18]。
目前对传统发酵乳制品中参与发酵微生物的分离和菌群多样性相关研究较多，孙思雨等[19]采用Illumina Miseq高通量测序对川西高原不同地区的传统自然发酵牦牛酸奶样品进行细菌群落结构分析；孙志宏等[20]采用454高通量测序技术从宏基因组角度对西藏地区自然发酵牦牛乳中微生物多样性进行探索性研究；SUN Z等[21]对采集自中国内蒙古、甘肃、四川以及蒙古国的17个自然发酵乳样品进行多样性分析；XU H Y等[22]对采集自中国新疆的22个自然发酵乳样品进行多样性研究；LIU W等[23]对采集自中国西藏的16个自然发酵牦牛乳样品进行分析；丁武蓉[24]采用传统培养法结合高通量测序对青藏高原传统发酵牦牛酸奶中乳酸菌多样性进行研究，然而，通过高通量测序分析比较不同保存温度条件下传统发酵乳制品微生物多样性和优化传统发酵乳制品保存温度却鲜有报道。
甘南牧区位于甘肃省甘南藏族自治州境内，地处青藏高原东北边缘与黄土高原接壤带，地理坐标位于东经100°46′～104°44′，北纬33°06′～36°10′之间，地势西北高东南低，地形状况复杂，牧区平均海拔3 000 m以上，气候主要以高原大陆性季风气候为主，地理季节分布有显著差异，是我国重要的畜牧业生产基地。独特的生态环境和自然长期驯化使得甘南牧区牧民沿用传统方法制作的牦牛酸奶含有丰富且特有的食品微生物资源。本研究以不同保存条件下（常温、4 ℃、-20 ℃和-80 ℃）甘南牧区牧民自制牦牛酸奶为研究对象，采用Illumina NovaSeq高通量测序技术，对牦牛酸奶细菌的16S rRNA V3-V4区基因序列进行测序，通过操作分类单元（operational taxonomic unit，OTU）聚类、α多样性和β多样性等分析比较不同保存温度下牦牛酸奶的细菌群落结构和多样性，探究不同保存温度对传统发酵乳制品微生物菌群多样性的影响，为传统发酵乳制品采样条件优化提供参考依据。
1 材料与方法
1.1 材料与试剂
1.1.1 样品采集
牦牛酸奶：2019年11月于甘肃省甘南藏族自治州卓尼县（Z）、碌曲县（L）实地采集牧民自制的牦牛酸奶，每份样品采集4管，每管30 mL，密封后标记样品号，分别于常温、冰袋、干冰和液氮保存条件下运回实验室后，分别在常温、4 ℃、-20 ℃和-80 ℃温度条件下保存备用，样品的具体信息见表1。
表1 甘南牧区牦牛酸奶采样信息
Table 1 Information of yak yogurt samples collected from Gannan pastoral area
续表
1.1.2 主要试剂
QIAampFastDNAStoolMiniKit、QIAquickGelExtraction Kit：德国QIAGEN公司；PowerSoil DNA Isolation Kit：美国MOBIO公司；KAPA HiFi Hotstart Ready Mix Kit：美国KAPA Biosystems公司；2×Phanta Master Mix：南京诺唯赞生物科技股份有限公司；AxyPrep DNA凝胶回收试剂盒：美国AXYGEN公司。
1.2 仪器与设备
5424R低温高速离心机：德国Eppendorf公司；VDRTEX-5涡旋混合器：海门市其林贝尔仪器制造有限公司；Tanon 2500R凝胶成像系统：上海天能科技有限公司；DYY-6C电泳仪：上海天能科技有限公司；ND-2000C微量紫外分光光度计：美国赛默飞世尔科技公司；MiniAmp聚合酶链式反应（polymerase chain reaction，PCR）仪：美国Life Technologies公司。
1.3 方法
1.3.1 牦牛酸奶菌群总DNA的提取
取2 mL牦牛酸奶样品，采用试剂盒QIAamp Fast DNA Stool Mini Kit，按照说明书的要求提取牦牛酸奶菌群的总DNA，利用微量紫外分光光度计和1%琼脂糖凝胶电泳进行总DNA的完整性和浓度质检。
1.3.2 PCR扩增与测序
以提取的总DNA为模板，341F（5'-CCTACGGGRSGCAGCAG-3'）和806R（5'-GGACTACVVGGGTATCTAATC-3'）为引物（引物中已带特定barcode），使用KAPA HiFi Hotstart Ready Mix PCR kit对细菌的16S rDNA V3-V4区基因序列进行PCR扩增，PCR扩增条件：95 ℃预变性3 min；98 ℃变性20 s，58 ℃退火15 s，72 ℃延伸20 s，共30个循环；72 ℃再延伸5 min。PCR扩增体系：2×KAPA Library Amplification Ready Mix 15 μL，引物（10 μmol/L）各1 μL，模板DNA 50 ng，加双蒸水（ddH2O）补充至30 μL。PCR扩增产物用2%琼脂糖凝胶电泳检测，使用AxyPrep DNA凝胶提取试剂盒切胶回收PCR扩增产物，文库质检合格后，使用Qubit 2.0进行文库定量，并根据每个样品的数据量要求，进行相应比例的混合，采用Illumina NovaSeq平台PE250策略进行测序，Illumina NovaSeq高通量测序由上海锐翌生物科技有限公司完成。
1.3.3 数据分析
16S rDNA V3-V4区基因序列被控制在220～500 bp之间，保证每条reads平均质量值不低于20，且含N碱基数不超过3个。将序列完全一样的Clean Reads根据其丰度大小进行排序，过滤掉其中的Singletons。使用UPARSE（http：//drive5.com/uparse/）根据97%的相似度进行OTU聚类，并使用Userach（版本7.0.1090）鉴定和移除嵌合体序列。每个OTU都有一个代表性的序列，将该代表性序列与已知物种的RDP数据库（http：//rdp.cme.msu.edu/）进行比对，置信度阈值设置为0.8，利用RDP Classifer（http：//rdp.cme.msu.edu/）对每个代表性序列进行物种注释。不同样本对应的Reads数量差距较大，为避免因样品数据大小不同而造成分析时的偏差，在样品达到足够的测序深度的情况下，对每个样品进行随机抽平处理。测序深度用alpha多样性指数来衡量，OTU profiling table和alpha多样性指数通过QIIME（version 1.9.1）的python脚本实现。
2 结果与分析
2.1 不同保存温度下牦牛酸奶样品的细菌菌群多样性分析
通过对牦牛酸奶样品细菌的16S rDNA V3-V4区基因序列进行测序，20份牦牛酸奶样品质控过滤后得到723 306条有效序列。采用稀释曲线评估各样品测序量是否足够，结果见图1。由图1可知，随着序列数的增加，超1（Chao1）指数趋于平缓，说明本次样本测序的深度足够，按最小样本序列数24 822抽平足以反映各样品中细菌物种多样性，测序量满足后续生物信息学分析的要求。各个样品及分组信息、质控后序列数、Observed species、Chao1指数、香农（Shannon）指数、辛普森（Simpson）指数和测序深度见表2。由表2可知，所有样本的测序深度均＞0.99，说明样本测序结果可以反映样品的真实情况。所有序列按97%的相似度进行OTU聚类，从20份牦牛酸奶样品中共检测到1 289个OTU，隶属于26门、46纲、61目、127科、248属。对不同保存温度样品组间α多样性指数进行kw检验分析，不同保存温度下样品组间α多样性指数Observed species、Chao1指数、香农（Shannon）指数、辛普森（Simpson）指数和测序深度的P值分别为0.56、0.59、0.94、0.88和0.61，说明不同保存温度样品组间多样性无显著差异（P＞0.05）。
图1 不存保存温度下牦牛酸奶样品的稀释曲线
Fig.1 Rarefaction curve of yak yoghurt samples at different storage temperature
表2 牦牛酸奶样品细菌的序列信息及α多样性指数
Table 2 Sequence information and α-diversity indexes of bacteria of yak yoghurt samples
2.2 不同保存温度下牦牛酸奶样品细菌群落结构的比较
2.2.1 组间细菌群落结构比较
主坐标分析（principal coordinate analysis，PCoA）可以直观地将样品间、样品组间差异程度体现在二维坐标图上，样品间距离越近则表示其物种组成越相似[25]，对不同保存温度下牦牛酸奶样品的细菌群落进行PCoA，并对UniFrac距离进行Adonis分析，结果见图2。
图2 各组牦牛酸奶样品的细菌群落主坐标分析
Fig.2 Principal coordinate analysis of bacterial community of yak yoghurt samples of each group
由图2可知，主坐标1和主坐标2分别解释了不同保存温度下牦牛酸奶样品中细菌群落73.16%和12.27%的信息，两主坐标贡献率之和＞85.43%，表明两个主坐标较好地代表了样品中的细菌群落信息。图2中4组样品间距离较近，表明4组样品的物种组成较为接近，Adonis分析表明4种不同保存温度牦牛酸奶样品细菌菌群多样性没有显著差异（P＞0.05）。
2.2.2 不同保存温度下牦牛酸奶样品细菌门水平物种丰度分析
在细菌门水平，所有牦牛酸奶样品中共鉴定出26个细菌门，将相对丰度排序靠后的6个细菌门并为other，结果见图3。从常温、4 ℃、-20 ℃和-80 ℃保存的牦牛酸奶样品中分别共鉴定出13个、12个、20个、13个细菌门，其中厚壁菌门（Firmicutes）、变形菌门（Proteobacteria）和拟杆菌门（Bacteroidetes）为所有牦牛酸奶样品中的优势菌门[26]，其相对丰度分别＞75.87%、7.20%和3.52%，未鉴定到门水平的细菌相对丰度均＜0.1%。通过门水平细菌群落结构分析比较，不同保存温度下耗牛酸奶样品组间细菌群落结构组成基本相同，仅组成比例存在一定差异。
图3 门水平不同保存温度下牦牛酸奶样品细菌群落结构分析结果
Fig.3 Analysis results of bacterial community structure of yak yogurt samples at different storage temperature at phylum level
2.2.3 不同保存温度下牦牛酸奶样品细菌属水平物种丰度分析
在细菌属水平，所有耗牛酸奶样品中共鉴定出248个细菌属，对于平均丰度值大于1%的细菌属，从常温保存的牦牛酸奶中鉴定出乳杆菌属（Lactobacillus）、链球菌属（Streptococcus）、醋杆菌属（Acetobacter）、金黄杆菌属（Chryseobacterium）、假单胞菌属（Pseudomonas）、丛毛单胞菌属（Comamonas），相对丰度分别为75.77%、11.88%、2.61%、2.50%、1.27%、1.18%；从4 ℃保存的牦牛酸奶中鉴定出乳杆菌属（Lactobacillus）、链球菌属（Streptococcus）、假单胞菌属（Pseudomonas）、金黄杆菌属（Chryseobacterium），相对丰度分别为75.40%、13.47%、3.02%、2.82%；从-20 ℃保存的牦牛酸奶中鉴定出乳杆菌属（Lactobacillus）、链球菌属（Streptococcus）、金黄杆菌属（Chryseobacterium）、假单胞菌属（Pseu-domonas）、丛毛单胞菌属（Comamonas），相对丰度分别为61.82%、14.97%、5.27%、5.19%、1.25%；从-80 ℃保存的牦牛酸奶中鉴定出乳杆菌属（Lactobacillus）、链球菌属（Streptococcus）、假单胞菌属（Pseudomonas）、金黄杆菌属（Chryseobacterium）、丛毛单胞菌属（Comamonas）、气单胞菌属（Aeromonas），相对丰度分别为68.86%、14.59%、5.2%、3.91%、1.27%、1.00%。图4显示牦牛酸奶样品在属水平相对丰度前20的物种属，其他合并为other。通过属水平细菌群落结构分析比较，不同保存温度下耗牛酸奶样品组间细菌群落结构组成相似，仅组成比例存在一定差异，乳杆菌属和链球菌属为优势菌属[26]，平均丰度值分别＞61.82%和11.88%。
图4 属水平不同保存温度下牦牛酸奶样品细菌群落结构分析结果
Fig.4 Analysis results of bacterial community structure of yak yogurt samples at different storage temperature at genus level
2.2.4 不同保存温度下牦牛酸奶样品组间物种差异分析
采用秩和检验得到不同保存温度下牦牛酸奶样品组间的差异细菌属共2个（P＜0.05），分别为Enhydrobacter、Peptoniphilus。为了直观展示各组及各样品间细菌属差异，对差异细菌属进行热图（Heatmap）分析，结果见图5。
图5 不同保存温度下牦牛酸奶样品间差异细菌属的Heatmap分析结果
Fig.5 Heatmap analysis results of differential bacterial genera between yak yoghurt samples at different storage temperatures
由图5可知，在各牦牛酸奶样品中Enhydrobacter的相对丰度高于Peptoniphilus，且在Z-20组中Enhydrobacter相对丰度较高，即-20 ℃保存条件下的耗牛酸奶样品中Enhydrobacter较高。为了直观展示不同保存温度下耗牛酸奶样品分组之间细菌属差异，对各组牦牛酸奶样品中的差异细菌属进行主成分分析（principal component analysis，PCA），结果见图6。
图6 不同保存温度下牦牛酸奶样品组间差异细菌属的PCA结果
Fig.6 PCA results of differential bacterial genera among groups of yak yoghurt samples at different storage temperature
由图6可知，不同保存温度样品组Z-4、Z-20、Z-80组的距离接近，细菌群落组成较为接近，Z-0组的距离较远，细菌群落组成差异较大，说明4 ℃、-20 ℃、-80 ℃保存组牦牛酸奶样品的物种组成更为接近，不同保存温度下耗牛酸奶样品组间细菌多样性无显著差异。
3 结论
本研究采用Illumina NovaSeq高通量测序技术对4种不同保存温度条件下甘南牧区牦牛酸奶细菌群落多样性及结构进行分析，结果表明，不同保存温度下的牦牛酸奶细菌群落多样性无显著差异，细菌群落组成基本相同，仅组成比例存在一定差异，优势细菌门为厚壁菌门（Firmicutes）、变形菌门（Proteobacteria）和拟杆菌门（Bacteroidetes），其平均相对丰度分别＞75.87%、7.20%和3.52%，优势细菌属为乳杆菌属（Lactobacillus）和链球菌属（Streptococcus），其平均相对丰度分别＞61.82%和11.88%。
[1]张和平.自然发酵乳制品中乳酸菌的生物多样性[J].生命科学，2015，27（7）：837-846.
[2]CHEN H C,WANG S Y,CHEN M J.Microbiological study of lacticacid bacteria in kefir grains by culture-dependent andculture-independent methods[J].Food Microbiol,2008,25(3):492-501.
[3]林莹，代道芳，辛志平，等.传统培养结合分子生物学法分析传统发酵食品中乳酸菌的多样性[J].中国酿造，2011，30（2）：52-54.
[4]周丹燕，戴世鲲，王广华，等.宏基因组学技术的研究与挑战[J].微生物学通报，2011，38（4）：591-600.
[5]孟飞，俞春娜，王秋岩，等.宏基因组与宏基因组学[J].中国生物化学与分子生物学报，2010，26（2）：116-120.
[6]李宇辉，王俊钢，刘成江，等.新疆伊犁牧区发酵乳制品中酵母菌的分离和多样性分析[J].食品与发酵工业，2013，39（7）：98-103.
[7]BAO Q H,LIU W J,YU J,et al.Isolation and identification of cultivable lactic acid bacteria in traditional yak milk products of Gansu Province in China[J].J Gen Appl Microbiol,2012,58(2):95-105.
[8]张亚川，蔡静静，剡文莉，等.新疆伊犁地区乳品中发酵菌种的筛选及产酸性能研究[J].中国乳品工业，2019，47（7）：4-7，18.
[9]秦艳婷.新疆地区传统发酵乳制品中乳酸菌的分离鉴定及生物多样性分析[D].呼和浩特：内蒙古农业大学，2014.
[10]曹晨霞，韩琬，张和平.第三代测序技术在微生物研究中的应用[J].微生物学通报，2016，43（10）：2269-2276.
[11]李伟程，侯强川，于洁，等.传统发酵乳制品中微生物多样性研究[J].食品工业科技，2018，39（1）：131-136.
[12]姚国强.传统发酵乳中细菌多样性及其功能基因研究[D].呼和浩特：内蒙古农业大学，2017.
[13]ZHONG Z,HOU Q,KWOK L,et al.Bacterial microbiota compositions of naturally fermented milk are shaped by both geographic origin and sample type[J].J Dairy Sci,2016,99(10):7832-7841.
[14]SUN Z,LIU W,BAO Q,et al.Investigation of bacterial and fungal diversity in tarag using high-throughput sequencing[J].J Dairy Sci,2014,97(10):6085-6096.
[15]CARAFA I,NARDIN T,LARCHER R,et al.Identification and characterization of wild lactobacilli and pediococci from spontaneously fermented Mountain cheese[J].Food Microbiol,2015,48(6):125-132.
[16]米其利，李雪梅，管莹，等.高通量测序在食品微生物生态学研究中的应用[J].食品科学，2016，37（23）：302-308.
[17]SHANGPLIANG H N J,RAI R,KEISAM S,et al.Bacterial community in naturally fermented milk products of Arunachal Pradesh and Sikkim of India analysed by high-throughput amplicon sequencing[J].Sci Rep,2018,8(1):1532.
[18]张冬蕾.应用焦磷酸测序技术对新疆地区传统发酵乳制品中微生物多样性的研究[D].呼和浩特：内蒙古农业大学，2016.
[19]孙思雨，陈炼红，王琳琳.应用Illumina Miseq测序分析川西高原传统牦牛发酵酸奶中细菌多样性[J].食品工业科技，2019，40（23）：98-103.
[20]孙志宏，刘文俊，张和平.基于宏基因组方法对西藏传统发酵牦牛奶中微生物多样性的研究[J].北京工商大学学报（自然科学版），2012，30（4）：19-24.
[21]LIU W,ZHENG Y,KWOK L Y,et al.High-throughput sequencing for the detection of the bacterial and fungal diversity in Mongolian naturally fermented cow's milk in Russia[J].Bmc Microbiol,2015,15(1):45.
[22]XU H Y,LIU W J,GESUDU Q,et al.Assessment of the bacterial and fungal diversity in home-made yoghurts of Xinjiang,China by pyrosequencing[J].J Sci Food Agr,2015,95(10):2007-2015.
[23]LIU W J,XI X X,SUDU Q,et al.High-throughput sequencing reveals microbial community diversity of Tibetan naturally fermented yak milk[J].Ann Microbiol,2015,65(3):1741-1751.
[24]丁武蓉.青藏高原传统发酵牦牛奶中乳酸菌多样性及其益生功能研究[D].兰州：兰州大学，2014.
[25]刘振东，毕娜，李哲，等.西藏不同产区曲拉细菌群落结构的比较分析[J].食品与发酵工业，2020，46（6）：60-66.
[26]张敏，张艳，黄丽丽，等.基于16S rDNA高通量测序方法比较新疆西北部地区乳品中微生物的多样性[J].食品科学，2017，38（20）：27-33.