基于16S rDNA高通量测序技术分析北京豆汁儿微生物多样性和功能预测的研究

,*

(1.北京农业职业学院,北京 102400; 2.广东轻工职业技术学院,广东高校特色调味品工程技术开发中心(GCZXB1103),广州 510300)

豆汁儿是由绿豆自然发酵而来的北京特有传统食品,具有降燥解毒、促进脾胃功能、消暑降温的功效[1-6]。豆汁儿闻之有酸臭味(与泔水味相似),入口以酸味为主伴有回甘,其独特的风味是由微生物发酵产生的。目前,豆汁儿中微生物的相关报道较少,有文献表明豆汁儿中存在乳酸乳球菌乳酸亚种[7]。通过宏基因组学对微生物的种属组成、丰度及其优势菌进行研究分析是目前微生物学研究的热点[8-10],被广泛运用于食品微生物研究领域[11-14]。通过微观层面的科学研究,人们对微生物菌落有了全面的认识[15-18]。以环境样品中总微生物基因组为研究对象,通过测序分析和功能基因筛选等对样本中的微生物进行较全面的分析,如:菌落多样性分析、样本与物种关系分析等[19-22]。其中,16S rDNA测序分析技术最为常用,在不对微生物进行纯化培养的情况下,直接从食品生态系统中获取基因信息对其进行测序分析[23-26]。

本文通过16S rDNA测序对豆汁儿中所含微生物的种属类型、种类的相对丰度以及样本与物种的关系等进行分析,从而了解豆汁儿中菌落组成及其优势菌种,以期为豆汁儿的科学安全生产奠定基础。

1 材料与方法

1.1 材料与仪器

豆汁儿 来自北京不同城区具有代表性的6家豆汁儿店,编号分别为A、B、C、D、E、F。

HR/T16M台式高速冷冻离心机 湖南赫西仪器装备有限公司;YXQ-LS-50SII压力蒸汽灭菌器 上海笃特科学仪器有限公司;JA2003 电子天平 上海精密科学仪器有限公司;高通量测序仪 Illumina(Miseq);中科美菱DW-HL290超低温储存箱 上海旦鼎国际贸易有限公司。

表1 Alpha多样性分析Table 1 Alpha diversity analysis

1.2 实验方法

1.2.1 微生物16S rDNA测序的样品前处理方法 分别取6种豆汁儿,立即用12000 r/min转速的高速冷冻离心机4 ℃条件下离心30 min,去除上清液后放入-80 ℃的超低温储存箱保存,同一家店分三批次取样混合均匀。

1.2.2 16S rDNA PCR扩增 选择799F和1193R作为引物,引物序列分别为AACMGGATTA GATACCCKG和ACGTCATCCCCACCTTCC。为了保证后续数据分析的准确性及可靠性,在可行范围内使用低循环数扩增并保证每个样品扩增的循环数统一。扩增体系为:TransStart Fastpfu DNA Polymerase,20 μL反应体系,其中包括:4 μL的5×FastPfu Buffer,42 μL的2.5 mmol dNTPs,0.8 μL的Forward Primer(5 μmol/L),0.8 μL的Reverse Primer(5 μmol/L),0.4 μL的FastPfu Polymerase,0.2 μL的BSA,10 ng的Template DNA,补ddH2O至20 μL。PCR反应参数:a. 95 ℃预变性3 min;b. 循环数×(95 ℃变性30 s;55 ℃退火30 s;72 ℃延伸45 s);c. 72 ℃延伸10 min,至10 ℃停止。

二轮扩增:一轮引物 799F,退火55 ℃,27个循环;二轮引物 799F-1193R,退火55 ℃,13个循环。测序工作由北京美吉桑格生物医药科技有限责任公司完成。

1.2.3 16S rDNA功能预测分析 16S rDNA功能预测是通过PICRUSt软件(PICRUSt软件[27]存储了Greengene ID对应的Clusters of Orthologous Groups of Proteins(COG)信息和KEGG Ortholog(KO)信息)对OTU丰度表进行标准化,即去除16S marker gene在物种基因组中的copy数目的影响;再通过每个OTU对应的Greengene ID,比对到COG和KEGG库,获得OTU对应的COG家族信息和KO信息;并计算各COG的丰度和KO丰度。根据COG数据库的信息,可以从EggNOG(evolutionary genealogy of genes:Non-supervised Orthologous Groups,http://eggnog.embl. de/)数据库中解析到各个COG的描述信息,以及其功能信息,从而得到功能丰度谱;根据KEGG数据库(Kyoto Encyclopedia of Genes and Genomes,京都基因和基因组百科全书,http://www.genome.jp/kegg/)的信息,可以获得KO、Pathway、EC信息,并能根据OTU丰度计算各功能类别的丰度。针对Pathway,运用PICRUSt软件可获得代谢通路的3个水平信息及各个水平的丰度表。

1.3 数据处理

采用Circos-0.67-7、PICRUSt软件、R语言工具统计和作图,Excel进行数据处理分析。

2 结果与分析

2.1 样本OTU及Alpha多样性分析

除去低于置信阈值的分类谱系外,6个豆汁儿样品总共获得有效序列204144条。按照97%相似水平对OTU进行统计分析,6个豆汁儿样品中,OTU50数量最多,共有53613个,其对应的菌种是食窦魏斯氏菌(Weissellacibaria);OTU35次之,共有50080个,其对应的菌种是沙克乳酸杆菌(Lactobacillussakei)。豆汁儿样本中包含的物种分类个数统计结果为:域1个;界1个;门5个;纲8个;目13个;科21个;属35个;种51个;OTU 70个。

各样本OTU数目、Ace丰富度指数、Chao1丰富度指数、测序深度指数、Shannon指数及Simpson指数见表1。由表1可知,样品E的Ace丰富度指数为76.02;Chao1丰富度指数为88.25,是所有样品中菌群丰度最高的。由Shannon指数可知,样品B多样性最高,Shannon指数最大为1.96。6个豆汁儿样品的Simpson指数值在0.21～0.27区间内,说明菌落多样性相似度较高。

2.2 物种组成分析

2.2.1 物种Venn图分析 图1为6个豆汁儿样品在种分类水平上的物种Venn图。如图1所示,6个样品共有物种数目22种;样品E和样品F共有物种数目最多,有42种。由此得出样品E与样品F物种相似度高。共有物种数目在35种以上的群组有:E&F、D&E、B&E、D&E&F、D&F。

图1 OTU交集韦恩图Fig.1 Intersected Venn diagram analysis of OTU

6个豆汁儿样品合并后微生物物种分布如图2所示。食窦魏斯氏菌(Weissellacibaria)占比最多为25.73%;次之为沙克乳酸杆菌(Lactobacillussakei),占比为24.63%。与样本OTU统计分析分析结果相吻合。

图4 菌落热图Fig.4 Community heatmap

图2 物种分布饼图Fig.2 Species distribution pie map

2.2.2 菌落组成分析 基于目分类水平柱形图3结果表明:不同样本在目水平优势物种相同,乳杆菌目(Lactobacillales)占比最大,6个样本中乳杆菌目都在95.8%以上,其他菌目占0.2%以下;B样品中有3.7%的双歧杆菌目(Bifidobacteriales)。由此看出,豆汁儿的发酵是以乳杆菌目为主的,不同来源的豆汁儿,其所含微生物相似度极高。

图3 基于目分类水平的菌落柱形图Fig.3 Community barplot analysis based on category level

2.2.3 群落Heatmap图 样品A～F这6个样本的群落属水平物种组成热图和样本聚类树分析如图4所示,在属分类水平上,根据豆汁儿样本间物种相对丰度的相似性进行聚类,对聚类后各样品中不同菌属所含序列的丰度作菌落热图,可使高丰度和低丰度的物种分块聚集,反映出在属分类水平上所有豆汁儿样品菌落结构的相似性和差异性。

热图共展示35个属,不同豆汁儿样本中TOP5门类组成存在差异。其中样品A、样品B、样品D、样品E、样品F占比最大的为乳杆菌属(Lactobacillus),是豆汁儿中的优势菌属;而样品C占比最大的是魏斯氏属(Weissella)。且样品A～F 6种样品中优势菌种的组成相似,菌种类别统一度高,均含有乳杆菌属(Lactobacillus)、乳球菌(Lactococcus)和魏斯氏属(Weissella);6种样品菌属的丰度指数存在差异,表明菌属的占比量不同。

2.2.4 样本与物种关系分析 在种分类水平上的Circos图分析结果如图5。A～F 6种样品中优势菌种的组成相似,但占比不同,其中B、D和E的优势菌种均为沙克乳酸杆菌(Lactobacillussakei),该菌种在样本中的占比分别为37%、32%和31%;样品C和F的优势菌种是食窦魏斯氏菌(Weissellacibaria),在样本中的占比分别是:44%和35%;样品A的优势菌种是一种未分类的乳杆菌,该菌种在样本的占比为37%。有研究表明:沙克乳酸杆菌能够抑制食物腐败菌和食源性病原体[28],本研究中,样品B、D、E中存在该菌种,因此推测对豆汁儿中的腐败菌具有抑制作用。

图5 Circos样本与物种关系图Fig.5 Circos sample and species diagram

样品A、D和E这3个样品中,均存在食窦魏斯氏菌(Weissellacibaria),分别占比均为22%、26%、28%,食窦魏斯氏菌不产生物胺类物质,耐低酸、胆盐,可促进发酵食品的风味形成[29-30]。

2.3 功能预测分析

16S功能预测是通过PICRUSt软件对OTU丰度表进行标准化后,通过每个OTU 对应的Greengene Id获得OTU对应的COG家族信息和KEGG Ortholog(KO)信息,并计算丰度。由于乳杆菌目(Lactobacillales)占比95.8%以上,是豆汁儿发酵过程中的优势菌目,因此以乳杆菌目(Lactobacillales)对应的OTU数据(筛选后的OTU命名a～f),进行功能预测分析。

2.3.1 COG数据库注释分析 COG是由NCBI创建并维护的蛋白数据库[31]。图6是COG功能分类统计图。由图6可知,样本a～f的COG功能组成相似度高,功能注释结果分为25 类。功能性强的前五类分别是:G:碳水化合物转运和代谢,均值963750个;R:一般功能预测,均值为794154个;E:氨基酸转运与代谢,均值743039个;K:转录,均值726421个;J:翻译,均值669363个;此外发现均值1005492个功能未知的基因,有待今后进一步研究。

图6 COG功能分类统计图Fig.6 COG functional classification statistical map

注:A. RNA加工与修饰;B.染色质结构与动力学;C.能量产生与转化;D.细胞周期调控、细胞分裂,染色体分隔;E.氨基酸转运与代谢;F.核苷酸转运和代谢;G.碳水化合物转运和代谢;H.辅酶转运与代谢;I.脂质转运与代谢;J.翻译,核糖体结构和生物合成;K.转录;L.复制,重组和修复;M.细胞壁/细胞膜的生成;N.细胞运动;O.翻译后修饰,蛋白质转换,伴侣;P.无机离子转运与代谢;Q.次级代谢产物的合成、转运和代谢;R.一般功能预测;S.未知功能;T.信号转导机制;U.细胞内的转运、分泌和囊泡运输;V.防御机制;W.胞外结构;Y.核结构;Z.细胞骨架。

2.3.2 KEGG数据库注释分析 KEGG(Kyoto Encyclopedia of Genes and Genomes)即京都基因和基因组百科全书,始于1955年,将基因、基因组信息以及更高层次的功能信息结合起来,对基因的功能进行系统化分析[32]。表2是KEGG代谢通路Level 2丰度分析结果。Level 2丰度分析共有41种代谢通路,表2展示a～f代谢通路丰度值前15的结果。主要代谢通路有膜运输、碳水化合物代谢、复制和修复、氨基酸代谢、翻译、特征差和核苷酸代谢等。

其中,a～f碳水化合物代谢路径丰度值占比分别为:14.97%、15.29%、14%、15.33%、15.19%、14.79%。有关研究表明:能够代谢碳水化合物及酶活性种类最多的菌种是鼠李糖乳杆菌和副干酪乳杆菌,而植物乳杆菌碳水化合物利用能力和酶活性却相对较少;乳酸杆菌对碳水化合物的代谢特征也可为特异性益生元的选择提供参考[33]。a～f氨基酸代谢丰度值占比分别为:9.14%、8.74%、9.94%、8.68%、8.77%、9.09%。氨基酸经过不同的代谢途径,在转氨酶、裂解酶和脱氢酶等的作用下产生各种挥发性香气成分[34]。

通过KEGG数据库分析,a～f中共有780余种酶,表3是丰度均值前50的酶类。由表3看出乙醛脱氢酶/乙醇脱氢酶、醇脱氢酶等丰度均值较高,由此推测豆汁儿中的风味物质主要来源于氨基酸代谢。乳酸杆菌α-半乳糖苷酶和β-半乳糖苷酶的转糖基作用可用于生产益生元α-低聚半乳糖和β-低聚半乳糖[33]。β-半乳糖苷酶属于糖苷水解酶,不仅能将乳糖水解为半乳糖和葡萄糖,而且能催化转糖苷反应合成半乳糖寡糖(Galactooligosaccharide,GOS),GOS是一种双歧因子,能促进肠道内双歧杆菌增殖,对人类的健康有许多益处[35]。由表3可知,a～f中β-半乳糖苷酶丰度均值较高,由此推测在豆汁儿发酵过程中,乳杆菌目(Lactobacillales)的酶活性特征可促使β-低聚半乳糖及GOS的产生,益于肠道菌群健康。

表2 KEGG数据库主要代谢通路丰度值Table 2 Abundance of main metabolic pathways in KEGG database

表3 KEGG数据库主要酶丰度值Table 3 Abundance of main enzyme in KEGG database

3 结论

通过微生物多样性检测分析,微生物16S rDNA测序可以初步断定:豆汁儿中乳杆菌目(Lactobacillales)占比95.8%以上,是优势菌目;进一步分析(图4所示)表明,乳杆菌属(Lactobacillus)是豆汁儿中的优势菌属;食窦魏斯氏菌(Weissellacibaria)和沙克乳酸杆菌(Lactobacillussakei)是优势菌种(表1和图5所示)。通过豆汁儿中乳杆菌目(Lactobacillales)的功能预测分析结果推测:豆汁儿中的风味物质来源于氨基酸代谢，发酵过程中产生β-低聚半乳糖及GOS,有益于肠道菌群健康。