应用Illumina Miseq测序分析川西高原传统牦牛发酵酸奶中细菌多样性

(西南民族大学生命科学与技术学院,四川成都 610041)

牦牛酸奶是高原地区游牧民族最常食用的乳制品[1],川西的高原地区空气稀薄,气温较低,含氧量少,太阳辖射强,降水量较少,牦牛奶资源丰富。长期以来,牧民沿用传统发酵方法制作牦牛酸奶[2-3],其富含营养,风味独特,由于地域、气候环境以及生活习惯不同,形成了特有的微生物菌群,具有广泛的生物多样性[4]。研究表明牦牛酸奶含有大量对人体有益的微生物,如乳酸菌、酵母菌和芽孢菌[5],其中以乳酸菌为主。乳酸菌是一种益生菌,可以改善人体内部的微生物环境,改善肠道,加速食物的消化与吸收,有效降低体内的胆固醇,抑制腐败物质在肠道中的繁殖等[6-8]。如今,乳酸菌已经在食品、医药等领域中得到了非常广泛的应用[9-12]。

第二代DNA测序(next-generation DNA sequencing,NGS)技术就是指能在单次生化反应中同时检测来自数千(或数百万)DNA模板上碱基序列的DNA测序技术,通过读取多个短DNA 片段,拼接成完整的序列信息,具有高数据通量、短序列读长(read)及低于Sanger测序法准确度的特点[13]。由Illumina开发的小型单道测序仪Miseq的误差率较低,能输出较长的read[14]。近年来,Miseq测序技术已经被广泛应用于测试一些复杂环境中的微生物多样性,包括污泥、土壤、海水和食物等[15-19]。目前Illumina Miseq已运用于牦牛乳制品中,对细菌群落组成进行分析,为开发优良微生物资源提供了重要的依据[20-21]。这项技术与传统培养法相比,既可以产生测序覆盖深度更大的数据量,还可以及时发现一些低丰度的细菌菌群,保证微生物菌落分析的准确性,同时直观反应微生物的群落类型。

本研究应用Illumina Miseq测序技术,对川西高原传统发酵的牦牛酸奶中细菌群落结构进行分析,以期深入了解川西高原牦牛酸奶中包含的细菌群落类型、结构和优势菌群,为保护和开发利用我国宝贵的高原乳酸菌资源、建立我国乳酸菌种质资源数据库,以及工业化生产牦牛酸奶提供理论指导。

1 材料与方法

1.1 材料与仪器

牦牛酸奶样品 采自甘孜藏族自治州和阿坝藏族羌族自治州4个县不同村落的8份牦牛酸奶样品,如表1,将牧民家用自然发酵方式制作的牦牛酸奶,在其发酵器中充分搅拌后,用干净的一次性勺装入已灭菌采样瓶和离心管中,将采集的样品放入冰盒内冷却,保持在低温状态下带回实验室放于-80 ℃冰箱;DNA抽提试剂盒 美国Omega Bio-Tek公司;2%Agarose gels 西班牙biowest公司;FastPfuPolymerase 中国TransGen公司;AxyPrep DNA Gel Extraction Kit 美国Axygen公司;Illumina MiSeq platform 美国Illumina公司;FastDNASPINKitforSoil 美国MP公司;AxyPrep DNA Gel Extraction Kit 美国Axygen Biosciences公司。

表1 传统牦牛酸奶来源Table 1 Source of traditional yok yogurt

Eppendorf 5424R高速台式冷冻离心机 德国Eppendorf公司;NanoDrop2000超微量分光光度计 美国ThermoFisherScientific公司;DYY-6C电泳仪 北京市六一仪器厂;Miseq测序仪 美国Illumina公司;PCR仪 ABI GeneAmp®9700型。

1.2 实验方法

1.2.1 总DNA的提取与检测 DNA的提取步骤参照试剂盒说明书进行操作。提取的DNA浓度以及纯度通过NanoDrop2000进行分析,并通过1%琼脂糖凝胶电泳对得到的DNA质量进行检验。

1.2.2 细菌16S rRNA V3-V4区PCR扩增 以样本总DNA为模板,用(5′-GGACTACHVGGGTWT CTAAT-3′)引物实现V3～V4[17-18]可变区的PCR扩增处理。扩增程序为:95 ℃预变性3 min,27个循环(95 ℃变性30 s,55 ℃退火30 s,72 ℃延伸30 s),最后72 ℃延伸10 min。扩增体系为20 μL,4 μL 5×FastPfu缓冲液,2 μL 2.5 mmol/L dNTPs,0.8 μL引物(5 μmol/L),0.4 μL FastPfu聚合酶;10 ng DNA模板。

1.2.3 Illumina Miseq测序 反应中的PCR产物通过2%琼脂糖凝胶进行回收,并通过AxyPrep DNA Gel Extraction Kit对得到的物质实施纯化处理,用Tris-Hcl清洗,最后将产物通过2%琼脂糖电泳法进行检测。产物的具体含量通过QuantiFluorTM-ST(Promega,USA)法测试得到。利用Illumina Miseq平台相关操作规范来将得到的纯化扩增片段建立PE 2×300的文库。并通过Illumina公司的Miseq PE300平台来测得其序列。

1.3 数据处理与分析

试验中得到的测试数据利用Trimmomatic软件进行质控,并应用FLASH软件实现拼接;通过UPARSE软件按照97%的相似度对序列展开OTU聚类,对其中存在的嵌合体主要是通过UCHIME软件实现剔除[22]。以RDP classifier来对获得的每条序列进行物种分类,并根据Silva数据库进行序列对比分析,本文中的比对阈值设置为70%。

在经过Miseq测序之后,首先通过overlap关系将得到的PE reads实施拼接,另外针对序列的质量展开质控和过滤,随后对样本展开OTU聚类分析和物种分类学分析。通过OTU分析得到物种的多样性指数,同时实现对测序深度的检测;根据分类学信息,能够实现在各个分类水平上的群落结构统计分析[23-25]。经过上述的分析,可以针对群落的组成以及发育信息进行多元化分析和差异性检验,实现更深层次的统计分析。

2 结果与分析

2.1 发酵牦牛酸奶测序信息

本试验测序样品数8个,共获得375236条优化序列,平均每个样品46905条。将所有样品的优化序列进行聚类,在97%相似水平下的OTU进行生物信息统计分析,4个县8个样品共有88条不同的OTU。

2.2 发酵牦牛酸奶稀释曲线

本研究用随机抽样的方法获取序列数以及其对应的OTU数目,并绘制相应的稀释性曲线[16]如图1所示,样品GC稀释曲线到10000条序列趋于平缓,样品BY3、NMJ稀释曲线到30000条序列趋于平坦,样品LV、WQ、SZLK、LAN、BY1稀释曲线到40000条序列时已经接近平缓,趋于平坦,意味着继续增加测序量对于发现新的OTU没有影响,样品中序列没有被测出的概率极低,测序深度基本覆盖样品中所有细菌。

图1 样品稀释曲线Fig.1 Sample dilution curve

2.3 发酵牦牛酸奶alpha多样性指数统计

如表2所示为不同牧民自制传统发酵牦牛酸奶中细菌菌群α-多样性分析。8个样品共产生88条不同的OTU,通过观测到的物种数(Sobs)、赵氏指数(Chao[26])和物种丰富度指数(ACE)可知,样品LV、LAN中物种数目较多,而样品BY1、BY3和NMJ中物种数较少,样品WQ物种数虽然只要25种,但其ACE指数和Chao指数均较高,说明其菌群丰度较高。香农指数(Shannon[27])和辛普森指数(Simpson[28])通常用来表征物种多样性,其中Shannon是估计样品中微生物多样性,Shannon值越大,说明群落多样性越高;Simpson主要是以定量的方式表征某区域物种的多样性,其数值越大,就意味着该区域物种多样性越低[29]。据表2可知,样品GC的Shannon指数最高为1.21,Simpson指数最低为0.38,说明样品GC细菌群落多样性最高且分布均匀,而样品BY1中包含的物种丰富度程度最低。Coverage是反映群落覆盖度(Community coverage)的指数,其数值越高,则样本中序列没有被检测出的概率越低,表中Coverage均为1,说明97%的物种类型都能在样品库得到,再次说明测序深度基本覆盖了所有细菌。

表2 alpha多样性指数统计Table 2 The alpha diversity index of statistics

2.4 发酵牦牛酸奶Venn图分析

Venn图可用于统计多组或多个样本中所共有和独有的物种(如OTU)数目,可以比较直观地表现环境样本的物种(如OTU)组成相似性及重叠情况[30]。由图2可知,8个样品按照不同地区分组,四个地区共同拥有10条OTU,占总样品OTU的11.4%。巴塘县、红原县和阿坝县共同拥有17条OTU;巴塘县、阿坝县和白玉县共同拥有12条OTU;红原县、阿坝县和白玉县共同拥有12条OTU;红原县、巴塘县和白玉县共同拥有10条OTU。巴塘县和红原县共同拥有28条OTU;巴塘和白玉共同拥有19条OTU;红原县和阿坝县共同拥有22条OTU;阿坝县和巴塘县共同拥有20条OTU;阿坝县和白玉县共同拥有14条OTU;红原县和白玉县共同拥有13条OTU。每个样品还有其仅有的OTU,巴塘县独有31条OTU,数量最多,白玉县仅独有1条OTU,数量最少,红原县独有9条,阿坝县独有3条。因此,从总体看来,四个县的牦牛酸奶样品中的微生物群落的结构具有一定的相似性,但由于环境和地理位置等差异,不同牦牛酸奶样品中具有其独有的微生物群落。

图2 OTU韦恩图Fig.2 Venn chart of OTU

2.5 发酵牦牛酸奶中细菌群落组成

为了实现对测试样品中包含物种的多样性信息进行表征,需要对所有的有效序列进行聚类,并以97%的序列相似性作为标准将序列聚类成为OTU,然后再深度分析[31]。

从表3细菌在门水平分布来看,8个样品共检测出8个细菌门,分别为硝化螺旋菌门(Nitrospirae)、Saccharibacteria、异常球菌-栖热菌门(Deinococcus-Thermus)、绿弯菌门(Chloroflexi)、放线菌门(Actinobacteria)、拟杆菌门(Bacteroidetes)、变形菌门(Proteobacteria)和厚壁菌门(Firmicutes)。所有样本均含有厚壁菌门,且含量占比高,占据总样本数量的98.33%～99.96%,这是因为牦牛酸奶中富含大量的乳酸菌,乳酸菌均为厚壁菌门,说明8份牦牛酸奶样品在属水平上群落结构多样性相近,优势菌门均为厚壁菌门。玛依乐·艾海提等[32]应用454焦磷酸测序技术对新疆阿图什及乌什传统发酵酸奶进行分析;智楠楠等[21]应用Illumina Miseq对市场上销售的酸奶进行分析;呼斯楞等[33]通过传统培养的方法对内蒙古地区酸奶进行了分析;以上实验得出的结果在门水平上均与本实验结果一致。硝化螺旋门主要是污水处理厂中的功能菌群[34],其仅存在于样品LV中,可能是牧民游牧时生活在污水处理厂附近,导致所采样品受到些许污染。放线菌门、变形菌门、拟杆菌门、硝化螺旋菌门和绿湾菌门均是植物土壤中常有的微生物菌门[35],牧民生活环境恶劣,常常以地为桌,做好的酸奶普遍放在地上,这可能是样品中均检测出少量这些微生物的原因。

表3 细菌在门水平群落结构Table 3 Bacterial community structure at phylum level

从图3细菌在属水平分布来看,8个样本共检测68个已知属。其中相对丰度占比排名前六的有乳杆菌属(Lactobacillus,86.99%),链球菌属(Streptococcus,9.90%),厌氧芽孢杆菌属(Anoxybacillus,1.75%),巨型球菌属(Macrococcus,0.31%),魏斯氏菌属(Weissella,0.25%),乳球菌属(Lactococcus,0.10%)。在所有样品中乳杆菌属占主导地位,相对丰度占比69.31%～99.92%,为优势菌属,其所占比例高低依次是BY1、SZLK、LAN、NMJ、LV、BY3、WQ、GC;在样品GC、BY3和NMJ中链球菌属相对丰度也较高,占比分别为26.39%、17.06%和11.12%,说明链球菌属在这三份牦牛酸奶发酵过程中可能起着较为重要的作用;在样品LV中厌氧芽孢杆菌属相对丰度也较高,占比为12.21%;在样品LAN巨型球菌属相对丰度也较高,占比为2.47%。8份牦牛酸奶的优势菌属均为乳杆菌属,说明其在种水平群落结构相近。

图3 细菌在属水平群落结构Fig.3 Bacterial community structure at genus level

2.6 发酵牦牛酸奶样品多样性种类异同分析

Heatmap图通过颜色梯度来表征二维矩阵中的数据大小,并呈现群落物种组成信息,通过颜色变化与相似程度来反映不同分组(或样本)在各分类水平上群落组成的相似性和差异性[36]。根据所有样品在属水平上的信息,挑选丰度前35的属,通过其在每个样品中的丰度信息从物种和样品两个方面进行聚类,得到heatmap图,如图4所示。由图4可知,乳杆菌属在所有样品中含量均为最高,链球菌属在样品GC、WQ、SZLK、BY3和NMJ中含量较高。除此之外,样品GC中相对于其他样品含量更高的菌属有魏斯氏菌属、乳球菌属(Lactococcus)、Escherichia-Shigella、金黄杆菌属(Chryseobacterium)、Subdoligranulum、奇异球菌属(Deinococcus);样品LAN中相对于其他样品含量更高的有巨型球菌属、Tepidiphilus、Paraclostridium和芽孢杆菌属(Bacillus);样品BY3中相对于其他样品含量更高的菌属有醋酸杆菌属(Acetobacter);样品BY1中相对于其他样品含量更高的菌属有明串珠菌属(Leuconostoc);样品LV中相对于其他样品含量更高的菌属有厌氧芽孢杆菌属。其中金黄杆菌属某些细菌是条件致病菌,在动物机体处于应激或免疫力低下时容易引起感染[37],并且金黄杆菌属均具有较高的耐药性[38]。样品GC中含有的金黄杆菌属,可能是母牛产奶时正处于免疫力低下的时候引起的感染。

图4 属水平物种丰富度聚类热图Fig.4 Heatmap at genus level

2.7 PCoA分析

PCoA分析主要是分析样本中菌落相似性或差异性,主成分以纵坐标表示,百分比表示该成分对样品差异的贡献值[32]。图中各点代表不同地区牦牛酸奶中菌群组成,当样品距离比较接近时表示两种物种的结构更加相似,所以一般将相似度较高的样品聚集在一起,两种菌落差异较大的样品则会存在很大的间距。由图5可知,样品BY1、NMJ、SZLK和BY3距离较近,其群落结构较为相似;样品SZLK、WQ、GC和LAN距离较近,其群落结构相似;样品GC、LAN和LV距离较近,其群落结构相似。样品BY1、NMJ、SZLK和BY3距离样品LV较远,说明其群落结构差异较大。样品之间群落结构存在交叉,表示不同地区的牦牛酸奶中的菌群结构不仅存在差异性,同时也具有相似性。

图5 牦牛酸奶样PCoA主坐标分析Fig.5 Yak yoghurt sample PCoA analysis

3 讨论

本研究利用Illumina Miseq测序技术对川西高原阿坝藏族羌族自治州和甘孜藏族自治州4个县8户牧民家经传统自然发酵的牦牛酸奶细菌群落结构进行分析。结果表明,在门水平上,细菌的优势菌门为厚壁菌门(Firmicutes),在属水平上,相对丰度占比排名前六的有乳杆菌属、链球菌属、厌氧芽孢杆菌属、巨型球菌属、魏斯氏菌属、乳球菌属,乳杆菌属相对丰度占86.99%,占主导地位。据包秋华[2]、呼斯楞等[33]、扎木苏[5]、丁武蓉[39]和孙志宏等[40]通过不同方法对牦牛酸奶中细菌多样性的分析,牦牛酸奶中主要菌门均为菌群均为厚壁菌门,主要菌属均为乳杆菌属,本研究中的优势菌属和其他学者研究的牦牛乳中乳酸菌的优势菌属一致,其他主要菌属也大致相似,但出现了其他研究中少有的厌氧芽孢杆菌属和巨型球菌属。厌氧芽孢杆菌在土壤以及人体肠胃中大量存在,绝大部分病菌属于腐物寄生菌,部分致病菌能分泌出外毒素和侵袭性酶类,可以治愈一些常规的疾病,例如肉毒梭菌(Clostridiumdifficile)等[41]。本试验中样品出现厌氧芽孢杆菌属,可能是牧民在挤牛奶的时候不慎沾染到少许牛粪等腐败物,而在后续杀菌处理中未将其完全消除,而残留至最后酸奶成品中。巨型球菌属是革兰氏阴性球菌,其特性与果胶杆菌类似,常见于麦汁和啤酒当中产生醋酸、丙酸和硫化氢,使啤酒出现浑浊和臭鸡蛋味[42]。埃氏巨型球菌(Megasphaeraelsdenii)是牛瘤胃正常的优势菌,经丙烯酰辅酶A生成丙酸,是瘤胃丙酸特殊生成途径之一[43]。本试验牦牛酸奶出现的巨型球菌属可能是来源自牛瘤胃中。

在所有样品中乳杆菌属占主导地位,相对丰度占比69.31%～99.92%,为优势菌属,说明8份样品在属水平上的细菌结构非常的相似。根据同期运用传统培养的方法对样品中乳酸菌进行分离鉴定的结果发现,乳杆菌属和乳球菌属为主要菌属,其中89%是乳杆菌属,说明Illumina Miseq测序技术能较全面地分析出样品中细菌的菌落结构,为实际运用提供理论依据。

4 结论

本研究利用Illumina Miseq测序技术对川西高原阿坝藏族羌族自治州和甘孜藏族自治州四个县8户牧民家经传统自然发酵的牦牛酸奶细菌群落结构进行分析。结果表明,测序样品数8个,共获得375236条优化序列,平均每个样品46905条,共有88条不同的OTU;样品GC的α-多样性指数最佳,细菌群落多样性最高且分布均匀;巴塘县的两个样品相对其它县占仅有的OTU数最多,为31;8个样品在门水平上,细菌的优势菌门为厚壁菌门(Firmicutes);在属水平上,主要细菌属为乳杆菌属(Lactobacillus)、链球菌属(Streptococcus)、厌氧芽孢杆菌属(Anoxybacillus)、巨型球菌属(Macrococcus)和魏斯氏菌属(Weissella),乳杆菌属(Lactobacillus)是所有样品的优势菌属;根据样品的Heatmap和PCoA分析可知,8份样品在属水平和种水平上的群落结构相似,具有细微的差异,说明不同地区的牦牛酸奶中的群落结构具有差异性和相似性。