高通量测序分析不同地区红腐乳细菌多样性

徐 琼，刘 洋,，曲勤凤，窦同海，陈羽菲，张娜娜，赵 磊，钟 江，翁史昱，杨捷琳，赵国屏

（1.上海市质量监督检验技术研究院，上海 200233；2.复旦大学生命科学学院，上海 200433；3.上海海关动植物与食品检验检疫技术中心，上海 200135）

腐乳是我国传统的发酵豆制品，至今已有一千多年的生产历史。它是用豆浆的凝乳状物经微生物发酵所制成的一种干酪型产品，故又被誉为“中国干酪”[1]。经微生物发酵后的腐乳除大豆本身所具有较高蛋白质、脂肪和大豆异黄酮[2]等营养外，还具有抗氧化[3-4]、抗高血压[5]、抗乙酰胆碱酯酶活性[6]等生理功能，深受人们的喜爱。

由于腐乳是半开放式生产，且大多为传统的手工操作，未有加热杀菌的过程，因此有可能被环境中的微生物杂菌污染[7]，从而影响产品的品质，甚至给消费者的安全带来隐患。如腐乳中常见的生物胺，主要是由微生物对游离氨基酸的脱羧作用产生[8]，诸如芽孢杆菌属（Bacillus）、梭菌属（Clostridium）、乳酸杆菌属（Lactobacillus）和假单胞菌属（Pseudomonas）均能产生脱羧酶[9-10]。此外，氨基甲酸乙酯[11]、肠毒素[12]、丙烯酰胺[13]也是常见的发酵过程中的微生物产生的副产物。

廖新浴等[14]利用纯培养和16S rDNA测序相结合对不同类型腐乳优势菌群进行了研究，发现不同腐乳中的菌群结构各不相同，同时存在优势的食源性致病菌如蜡状芽孢杆菌（Bacillus cereus）、粪肠球菌（Enterococcus faecalis）和恶臭假单胞菌（Pseudomonas putida）。由于该结果依赖于纯培养，对于含量较低和非可培养微生物可能有所遗漏，无法对样品中的所有微生物群落进行全面分析。

近年来，随着测序技术的不断发展，高通量测序技术也被称为下一代测序技术，由于成本低、可读量和质量高，是目前评估微生物多样性的有力工具[15-16]。许多研究已经将该技术用于发酵食品中菌群多样性的分析，如豆豉[17]、豆汁[18]、大酱[19]、曲拉[20]等，但运用高通量测序技术分析华东、华北和东北地区腐乳菌群多样性的研究鲜见报道。

本研究以我国不同地区红腐乳为研究对象，采用高通量测序技术对细菌16S rDNA V1-V3区进行测序，以期更加全面地了解不同地域红腐乳中细菌菌群多样性，明确各自的优势菌群，为建立地域性特色发酵食品提供理论依据。

1 材料与方法

1.1 材料与试剂

红腐乳样品从我国6 个省市当地超市购买，每例样品来源于不同生产厂家，共计9 例，具体地区与样品情况见表1。

表1 样品信息Table 1 Information about the samples tested in this study

聚合酶链式以应（polymerase chain reaction，PCR）Master Mix（M7505） 美国Promega公司；胶回收试剂盒（Code No. 9762）、引物 宝生物工程（大连）有限公司。

1.2 仪器与设备

SW22振荡水浴槽 德国Julabo公司；Centrifuge 5417R冷冻离心机 德国Eppendorf公司；Nanodrop2000超微量分光光度计 美国Thermo Scientific公司；VERITI梯度PCR仪 美国ABI公司；Sub-Cell电泳仪、Geldoc XR+型凝胶成像仪 美国Bio-Rad公司；MiSeq高通量测序平台 美国Illumina公司。

1.3 方法

1.3.1 核酸提取

采用十六烷基三甲基溴化铵（hexadecyltrimethylammonium bromide，CTAB）法[21]并稍作修改。取同批次样品3 份，从每份样品中倒取50 mL腐乳汤汁，共计150 mL。混合后，取混合样品汤汁5 mL，9 000 r/min离心10 min，去上清液，沉淀中加入10 mL的无菌水，在振荡器上振荡混匀，9 000 r/min离心10 min，去上清液，沉淀加入800 μL溶菌酶（10 mg/mL），37 ℃水浴1 h；然后加入40 μL的蛋白酶K和CTAB裂解液，56 ℃水浴2 h。每管裂解液中加入等量Tris饱和酚-三氯甲烷-异戊醇（25∶24∶1，V/V），混匀，13 000 r/min离心10 min。取上清液至新离心管中，加入等体积三氯甲烷-异戊醇（24∶1，V/V），混匀，13 000 r/min离心10 min。上清液移至新离心管中，加入2 倍体积无水冰乙醇，混匀，静置30 min左右，13 000 r/min离心2 min。弃上清液，用70%乙醇溶液洗涤沉淀2～3 次，得到的沉淀放置于室温条件下进行自然风干。风干后加入50 μL Tris-乙二胺四乙酸缓冲液，溶解沉淀。

1.3.2 PCR扩增及16S rRNA测序

以提取的总细菌基因组D N A为模板，引物针对16S rRNA基因V1-V3区合成特异引物，利用27F（5’-AGAGTTTGATCCTGGCTCAG-3’）和534R（5’-ATTACCGCGGCTGCTGG-3’）引物进行PCR扩增。PCR体系：基因组模板DNA 2 μL，2×PCR Master Mix 12.5 μL，正以引物各0.5 μL，ddH2O补足至25 μL。PCR扩增条件：94 ℃预变性2 min；94 ℃变性30 s，60 ℃退火50 s；72 ℃延伸45 s，扩增40 个循环；4 ℃保存。对所得的PCR产物用1%的琼脂糖凝胶电泳进行检测，对符合要求的PCR产物进行后续实验。扩增后的PCR产物经纯化回收后，使用Illumina MiSeq平台进行高通量测序及生物信息学分析。

1.4 数据处理

高通量测序得到的双端序列数据使用FASTQC（http://www.bioinformatics.babraham.ac.uk/projects/fastqc）进行质量控制。Raw data上的接头使用CUTADAPT（http://code.google.com/p/cutadapt/）去除。然后使用FLASH（https://sourceforge.net/projects/fl ashpage/ fi les/）连接Pair-end两端Reads。使用QIIME[22]筛选序列和嵌合体，嵌合体去除使用QIIME的Usearch[23]算法，参数默认。

基于可操作分类单元（operational taxonomic units，OTU）聚类分析结果，研究环境中微生物的多样性，使用QIIME分析流程中pick_open_reference方式进行OTU聚类。聚类算法使用Uclust[23]，数据库使用Greengenes 2013-08 release（http://greengenes.lbl.gov/Download/）版本，相似度阈值为97%，对所有序列进行OTU划分并进行生物信息统计分析，其余参数为QIIME默认参数。根据样本分析所得OTU情况使用QIIMEα_diversity分析模块进行α多样性分析，其中菌群丰度可由Chao 1指数进行评估，其数值越高表明群落物种丰富度越高；菌群多样性由Shannon指数、Simpson指数进行评估，Shannon指数与群落多样性呈正相关，与Simpson指数呈负相关。Coverage指数是测序深度指标，代表每个样品文库的覆盖率。

使用β多样性分析进行不同样品在微生物群落构成上的比较。使用软件QIIMEβ_diversity_through_plots,jackknifed_β_diversity模块进行β多样性分析。在之前分析步骤中，已经将序列按照其碱基组成的相似性，划归到各个OTU中。在进行分类学分析时，首先，将每一条典型序列都与Greengenes 2013-08 release（http://greengenes.lbl.gov/Download/）数据库进行比对，找出其最相近且可信度达80%以上的种属信息，以获得每个OTU的分类学信息。根据物种注释，统计每个样品在分类水平上的序列数目，分类水平分别是门水平（将相对丰度＜0.1%及unidentified归入others）和属水平（相对丰度＜5%及unidentified归入others）。使用软件QIIME，算法为rdp_classifier[24]v2.2。热图使用R语言的“pheatmap”包绘制。基于菌群之间Unifrac[25]算法加权的主坐标分析（principal co-ordinates analysis，PCoA）可以实现多个样品的分类。

2 结果与分析

2.1 细菌α多样性分析

α多样性以映样品内部微生物群落的物种丰富度，根据97%相似性水平下的OTU信息，采用α多样性指标的Shannon指数、Simpson指数和Chao 1指数对样品微生物物种的丰富度和多样性进行评估，结果见表2。所有样本的测序Coverage指数大于0.9，说明腐乳样品文库种序列基本都被测出，未被测到的概率较低，测序结果可以以映样品真实情况。由表2可知，东北地区样品（R7和R8）的Chao 1指数和Simpson指数均低于其他样品，说明该地区红腐乳样品中细菌群落丰度和群落多样性较低，与其他地区存在一定差异。

表2 红腐乳样品中细菌α多样性Table 2 α-Diversity of bacteria present in red sufu samples

采用对测序序列进行随机抽样的方法，以抽到的序列数与它们所能代表OTU数目构建曲线，即稀释性曲线。当曲线趋于平坦时，说明测序数据量合理，更多的数据量对发现新OTU的边际贡献很小；以之则表明继续测序还可能产生较多新的OTU。因此，通过作稀释性曲线，可以得出样品的测序深度情况。由图1可知，多数样本在序列数在接近10 000时，OTU已经接近饱和，再增加测序量对于发现新OTU的边际贡献很小。说明测序数据量合理，所含菌种数较多。现有测序量已经可以以映样品中细菌丰度信息。

图1 各样品中细菌稀释曲线Fig. 1 Rarefaction curves of bacteria in samples

2.2 细菌群落结构差异性分析

图2 门水平下红腐乳样品中细菌群落结构分布Fig. 2 Distribution of bacterial communities in red sufu samples at phylum level

由图2可知，红腐乳中主要细菌优势门分别为：平均相对丰度为73.93%的厚壁菌门（Firmicutes）、15.12%的放线菌门（Actinobacteria）以及9.63%的变形菌门（Proteobacteria）。不同地区红腐乳中细菌主要优势菌群基本相同，但是相对丰度差异较大，即样品R2和R8中放线菌门的相对丰度分别为42.21%和44.96%，而样品R9的放线菌门相对丰度只有2.34%；R1中变形菌门的相对丰度是所有红腐乳中最高的，达到31.63%。说明在门水平下，不同地区红腐乳中细菌群落结构组成相同，仅组成比例存在一定差异。其中厚壁菌门是不同地区红腐乳中的绝对优势菌群。

在属水平下，9 个红腐乳样品共检测出296 个细菌属，图3显示相对丰度前12的物种。其中共有的菌属分别为：棒状杆菌属（Corynebacterium）、四联球菌属（Tetragenococcus）、明串珠菌属（Leuconostoc）、Rummeliibacillus、乳球菌属（Lactococcus）、魏斯氏菌属（Weissella）、链球菌属（Streptococcus）、丛毛单胞菌属（Comamonas）、芽孢杆菌属（Bacillus）以及Trabulsiella。

图3 属水平下红腐乳样品中细菌群落结构分布Fig. 3 Distribution of bacterial communities in red sufu samples at genus level

从群落结构分布图可以发现不同样品中优势菌在种类及相对分布上有较大的差异。同一地区的样品间细菌群落组成相似：R3、R5和R6样品中优势属为明串珠菌属（29.45%、36.86%、30.72%）、乳球菌属（9.59%、12.55%、20.72%）以及四联球菌属（8.75%、11.36%、12.05%）；R4样品在丰度上与上述样品有所差异，分别对应明串珠菌属（13.33%）、乳球菌属（10.81%）、四联球菌属（47.11%）；而同样来自华东地区的样品R1，则展现了不同的细菌群落组成，相对丰度为22.83%的丛毛单胞菌属是该样品的优势菌属，这个属在其他样品中的相对丰度均小于2%。丛毛单胞菌属在1985年由De Vos等[26]建立，该属中的菌株更多地与环境污染物的降解有关[27-28]，在食品中鲜有报道。因此在食品发酵中对于该菌属的认识值得深入研究。

此外，样品R5、R6中的魏斯氏菌属（5.55%、4.03%）是相比于其他样品所特有的菌属。芽孢杆菌属和盐厌氧菌属（Halanaerobium）主要出现在R7、R8和R9样品中，芽孢杆菌属为R7和R9的优势菌，相对丰度分别为58.08%和52.24%。华东地区的优势菌属明串珠菌属在R7、R8和R9中相对丰度较低，均低于5%。R8中棒状杆菌属（44.08%）是该样品的优势菌属。相比于其他样品，R8中梭菌属（Clostridium）相对丰度最高，达到7.12%。梭菌属是能形成芽孢、厌氧生长的革兰氏染色阳性大杆菌，可以分解糖、蛋白质，通常可以产生混合的有机酸和醇类，不还原硫酸盐。它们的代谢物极富多样性，广泛分布在环境中，许多种可产生外毒素危害人类和牲畜生命健康[29]。

2.3 细菌β多样性分析

根据OTU的丰度数据，热图可在属水平上将不同的OTU分块聚类，并根据颜色梯度以映不同样品中细菌群落相似性、差异性及物种聚类关系。利用R语言绘制热图，可以直观地显示不同地区样品中微生物的差异。由图4可知，上侧为样品聚类树，左侧为OTU聚类树，不同地区的9 种红腐乳可分为3 个聚类：R7、R8和R9聚为一类，R3、R4、R5和R6聚为一类，R1和R2聚为一类，说明不同样品间细菌群落存在一定的差异性。

图4 样品中不同属组成和相对丰度的热图Fig. 4 Heatmap of composition and relative abundance of different genera

图5 不同地区腐乳样品中细菌群落主成分分析Fig. 5 Principal coordinate analysis of bacterial communities in red sufu samples

基于UniFrac距离的加权主坐标法对9 种红腐乳细菌群落结构进行分析。如图5所示，细菌群落第1主成分（PC1）贡献率为52.36%，第2主成分（PC2）贡献率为22.77%。两主成分之和大于70%，表明2 个成分较好地代表了样品中的细菌群落信息。从图5可知，不同腐乳样本中的菌群与地理位置存在一定关系，同一地区的红腐乳样品具有一定的聚类趋势，这与热图样品聚类结果相似：R3、R4、R5和R6聚为一类，均来源于华东地区；样品R7（东北）和R9（华北）聚为一类；R1（华东）、R2（华北）、R8（东北）3 个样品较为离散，各自独立。

3 讨 论

腐乳是中国独具特色的大豆发酵食品之一，大豆经微生物发酵后所制成的一种干酪型产品[30]。通常根据其色泽风味可分为白腐乳、红腐乳、青（臭）腐乳三大类[31]。红腐乳添加红曲浸润数月后酿制而成，故表面呈自然红色。腐乳风味的形成主要取决于微生物的发酵，腐乳中微生物群落特征的差异是不同地区风味不同的重要原因之一。有学者对红腐乳中的红曲霉进行遗传学比较，发现不同地区的红曲霉存在一定的遗传差异性[32]。由于腐乳特殊的生产工艺以及未有加热杀菌的过程，使得腐乳中微生物种类和数量受到广大学者的关注。因此对腐乳中微生物菌群的分析研究有助于为产品质量控制、挖掘潜在功能性微生物及风味改良提供理论支持。目前对于腐乳中微生物菌群的研究主要分为培养法和非培养法两种。利用传统微生物分离培养法，红腐乳中主要的优势菌种为蜡状芽孢杆菌和Virgibacillus senegalensis[14]；植物乳杆菌、嗜酸乳杆菌、乳杆菌Akporo7和嗜盐乳杆菌[33]也被发现存在于腐乳中。周婷婷等[34]分析了不同品牌白腐乳样品的蜡样芽孢杆菌总数，其中有2 个样品的蜡样芽孢杆菌数达到106CFU/g。Han等[31]对3 种类型腐乳中的不同微生物进行了鉴定和计数，发现腐乳中致病菌有蜡样芽孢杆菌、产气荚膜梭菌。值得注意的是，由于样品中存在大量的芽孢菌，干扰了金黄色葡萄球菌在选择性平板上的生长，样品中未能检测到金黄色葡萄球菌，但在白腐乳和臭腐乳样品中检测到金黄色葡萄球菌肠毒素A。可见，运用纯培养法对于腐乳中微生物菌群的认识有一定的局限性。

利用PCR-变性梯度凝胶电泳技术，比较不同品牌腐乳细菌多样性，发现乳酸杆菌属是腐乳中细菌的优势菌群[35]。耿予欢等[36]建立了不依赖纯培养的肠杆菌基因间重复一致序列-PCR技术，监测了腐乳发酵过程中微生物群落结构变化。上述两种方法虽然弥补了传统纯培养法的缺陷，但是均依赖于电泳技术，结果重现性相对较低。

高通量测序技术不仅能够快速地分析复杂微生物群落多样性，还能检测到低存在率以及不可培养的微生物。刘亚栋[37]利用高通量测序分析菌属发现同一品牌的红、白和青腐乳之间细菌组成具有较大的差异性。Huang Xiaoning等[38]利用Illumina MiSeq高通量测序平台发现肠杆菌属、不动杆菌属、乳球菌属是红腐乳发酵生产过程中的主要优势菌属。

本研究利用高通量测序技术对不同地区红腐乳中细菌多样性进行分析。从属水平上菌群结构分布可知，不同地区的腐乳中细菌的多样性存在差异。华东地区的腐乳主要优势菌属为明串珠菌属、乳球菌属以及四联球菌属；芽孢杆菌属则是东北地区红腐乳的主要优势菌属；此外，魏斯氏菌属、梭菌属则分别是华东地区（江苏省）、东北地区（黑龙江省）所特有的菌属。

对9 种红腐乳细菌群落结构进行主成分分析，发现同一地区的红腐乳样品中菌群存在一定的聚类关系。目前，我国腐乳生产发酵依据豆腐坯培菌的菌种不同，可分为毛霉型、根霉型和细菌型[39]，如R8样品是采用藤黄微球菌（Micrococcus luteus）进行发酵的细菌型腐乳，使得菌群结构不同于其他东北地区的样品，两个东北地区样品菌群结构距离较远。由于腐乳的生产过程多为开放式生产，且多为传统的手工操作，生产车间空气、主要接触面以及生产过程中的原辅料中的微生物均会对最终的产品产生影响[7]，使得每个厂家的腐乳中微生物菌群具有一定的独特性，这可能是R1样品与其他华东地区样品离散的原因；而R2样品中还添加了玫瑰花作为辅料，玫瑰花中的微生物在后发酵过程中影响了该样品的菌群结构。

4 结 论

本研究通过高通量测序技术对不同地区红腐乳中细菌多样性进行分析，对其中优势菌群组成进行了初步分析。结果发现，本次测序深度有效地覆盖了样品中微生物种类，不同地区红腐乳样品中细菌群落结构具有丰富多样性。在门水平上，厚壁菌门、放线菌门和变形菌门为红腐乳中主要的门；而在属水平上，不同地区红腐乳样品中优势菌在种类及相对分布上有较大的差异，明串珠菌属、乳球菌属以及四联球菌属为华北地区样品中主要的优势菌属，而东北地区腐乳样品中的芽孢杆菌属和盐厌氧菌属为主要的优势菌属。