王玉荣，杨成聪，葛东颖，尚雪娇，张振东，郭 壮*

（湖北文理学院食品科学技术学院，鄂西北传统发酵食品研究所，湖北 襄阳 441053）

我国传统发酵食品的种类众多且历史悠久，受制作工艺、地域环境和原材料的影响，其风味和微生物群系极具多元化[1]。为实现传统发酵食品的现代化和标准化生产，保障产品的质量安全，对传统发酵食品生产过程中微生物群系变化进行解析则显得极为重要[2]。相对于传统微生物培养方法和指纹图谱技术，以MiSeq技术为代表的第2代高通量测序技术具有无需培养、灵敏度高和检测速度快等优点，逐渐成为解析我国传统发酵食品微生物多样性的常规途径[3]，目前在腌干鱼[4]、浓香型白酒窖泥[5]、发酵黑茶[6]、酸奶[7]和茯砖茶[8]中有广泛的应用。

由于MiSeq高通量测序技术的读长不能完全覆盖16S rRNA的全长基因，所以研究人员在对传统发酵食品微生物多样性进行解析时，常对其单可变区或多可变区进行测序[9]。16S rRNA基因扩增的准确性直接决定了高通量测序技术的可信度，然而研究发现DNA聚合酶[10]、扩增体系[11]和引物[12]均会对聚合酶链式反应（polymerase chain reaction，PCR）结果产生偏好性。引物的选择对保证PCR扩增产物的特异性和质量尤为重要，只有选择合适的引物才可避免异源双链和单链DNA的形成[13]。目前研究人员通常使用引物扩增16S rRNA基因的V1～V2区[14]、V1～V3区[15]、V3区[16]、V3～V4区[17]、V4区[18]、V4～V5区[19]、V5～V6区[20]和V6～V8区[21]进行发酵食品细菌多样性高通量测序分析，而Sun Donglei等[12]指出采用V4～V5区进行细菌多样性解析时结果最接近真实情况。我国发酵食品众多且不同发酵食品之间基质存在较大差异，因而引物的选择尤为重要，其对细菌多样性结果的解析可能存在较大影响，然而关于扩增区域对我国特色发酵食品细菌MiSeq测序结果影响的报道尚少。

鲊广椒是以玉米面和新鲜辣椒为原料固体发酵而成的食品，在我国湖北、云南、重庆、四川和贵州等地区有着广泛的食用人群[22]。本研究使用高通量测序公司常提供的3 对引物27F/338R、338F/806R和515F/907R分别对细菌16S rRNA的V1～V2区、V3～V4区和V4～V5区进行扩增并进行MiSeq高通量测序，同时采用生物信息学和多元统计学相结合的方法，对不同引物扩增PCR产物的细菌多样性进行了解析，以期为食品科学领域研究人员在进行传统发酵食品细菌多样性解析时提供参考。

1 材料与方法

1.1 材料与试剂

DNeasy mericon Food Kit 德国QIAGEN公司；AxyPrepDNA凝胶回收试剂盒 康宁生命科学吴江有限公司；dNTPs Mix、DNA聚合酶（5 U/μL）、5×TransStartTM、FastPfuBuffer和FastPfuFly DNA Polymerase北京全式金生物技术有限公司；引物27F/338R、515F/907R和338F/806R由武汉天一辉远生物科技有限公司合成。

1.2 仪器与设备

CT15RE台式冷冻离心机 日本Hitachi公司；2100芯片生物分析仪 美国Agilent公司；ND-2000C微量紫外分光光度计 美国Nano Drop公司；vetiri梯度基因扩增仪 美国AB公司；FC3化学发光凝胶成像系统 美国FluorChem公司；MiSeq PE300高通量测序平台 美国Illumina公司；R920机架式服务器 美国Dell公司。

1.3 方法

1.3.1 样本采集及DNA提取

5 份以玉米粉和鲜红辣椒为主要原料制作的鲊广椒样品采集于恩施市舞阳坝体育场菜市场（109.47°N，30.3°E），样品装入无菌采样袋中置于含有冰袋的采样箱中低温运送回实验室，48 h内使用QIAGEN DNeasy mericon Food Kit试剂盒按照约束步骤完成DNA提取，DNA样品置于－20 ℃备用。

1.3.2 PCR扩增及产物纯化

选取16S rRNA基因某一区域作为目片段进行扩增，为便于区分每个样品的序列，在引物的前段分别添加含有7 个核苷酸序列的标签（barcode）。测序所用引物序列信息如表1所示。

表1 引物序列信息Table 1 Primer sequences used for sequencing

将提取的D N A作为模板进行P C R扩增，反应体系（2 0 μ L）为：D N A模板1 0 n g，正 反 向 引 物 （ 5 μ m o l/L） 各 0.8 μ L ， D N A聚合酶0.4 μ L，d N T P s m i x（2.5 m m o l/L）2 μL，PCR缓冲液（10×）4 μL，体系用ddH2O补充至20 μL。扩增条件为：95 ℃预变性3 min，95 ℃变性30 s，55 ℃退火30 s，72 ℃延伸45 s，35 个循环，72 ℃末端延伸10 min。

采用2.0%的琼脂糖凝胶电泳对PCR产物进行检测，并使用AxyPrepDNA凝胶回收试剂盒对PCR产物切胶回收，回收的产物分别稀释至100 nmol/L后混匀进行文库构建。

1.3.3 文库构建及MiSeq高通量测序

按照TruSeq®DNA Sample Preparation G uide中的Low Sample (LS) Protocol进行文库制备，使用PE 300平台进行双端测序。

1.3.4 序列质控

数据下机后，首先根据双端序列的重叠关系将其拼接为一条序列，然后依据barcode信息将序列划分到各样品，在进行序列方向校正的基础上，切除barcode和引物序列得到有效序列。在拼接过程中若成对序列符合下述条件中的一条则定义为低质量序列，在质控过程中予以剔除：1）重叠区的碱基数小于10 bp；2）重叠区的最大错配比率大于0.2；3）barcode碱基错配；4）超过1 个引物碱基发生错配；5）有效序列长度小于50 bp。

1.3.5 生物信息学分析方法

使用QIIME平台（v1.70）[26]进行生物信息学分析：1）PyNAST校准并排齐序列[27]；2）按照100%和97%相似性进行两步UCLUST归并[28]，建立无重复的单一序列集和分类操作单元（operational taxonomic units，OTU）；3）ChimeraSlayer剔除含有嵌合体序列的OTU[29]；4）使用RDP（Ribosomal Database Project，Release 11.5）[30]和Greengenes（Release 13.8）[31]数据库对OTU中的代表性序列进行比对，在界、门、纲、目、科和属水平上对每个样品的细菌菌系进行解析；5）若某OTU无法鉴定到界水平或在属水平上将其鉴定为玉米属（Zea）则将其定义为非特异性扩增，使用自编脚本将该OTU中所包含的所有序列剔除，并依据barcode信息将序列划分到各样品，使用自编脚本对各序列的长度进行统计，并将所有样品的序列文件合并为1 个fna文件后重复1～4步骤；6）使用FastTree软件绘制系统发育进化树[32]，计算Chao 1指数和Shannon指数[33]；7）基于UniFrac距离[34]进行主坐标分析和聚类分析。

1.3.6 核酸登录号

含有非特异性扩增的序列数据已提交至MG-RAST数据库（http：//metagenomics.anl.gov/），登录号为mgp85472。

1.4 数据统计分析

使用配对Kruskal-wallis检验对各引物的非特异性扩增比例进行显著性分析，使用Venny 2.1（http：//bioinfogp.c nb.csic.es/tools/venny/index.html）绘制维恩图，进而对不同引物扩增产物中的核心OTU数目进行展示；使用多元方差分析（multivariate analysis of variance，MANOVA）对不同引物扩增产物中细菌群落结构的差异进行显著性分析；使用Origin 8.6软件绘制其他图。

2 结果与分析

2.1 基于非特异性扩增及序列长度的引物对MiSeq测序的影响

由图1可知，在使用引物27F/338R、338F/806R和515F/907R对5 个鲊广椒DNA样品16S rRNA基因的V1～V2区、V3～V4区和V4～V5区进行扩增时产生的非特异性条带分别占到序列总数的30.01%、22.04%和13.56%。经配对Kruskal-wallis检验发现，扩增16S rRNA基因V4～V5区的非特异性序列所占比例显著低于其他两个可变区（P＜0.05）。在对特异性扩增序列进行剔除后，本研究共产出444 543 条高质量16S rRNA序列，平均每个样品产出19 636 条，序列长度统计如图2所示。

图1 非特异性扩增序列的比例Fig. 1 Percentages of non-specific amplification sequences

图2 序列长度统计Fig. 2 Statistic sequence lengths

由图2可知，去除引物和barcode后，使用引物27F/338R扩增16S rRNA基因V1～V2区的序列长度分别有11.98%和70.68%集中在327～329 bp和351～353 bp，使用338F/806R扩增16S rRNA基因V3～V4区的序列长度95.50%集中在450 bp，而使用引物515F/907R扩增16S rRNA基因V4～V5区的序列长度94.11%集中在395～397 bp。

2.2 基于各分类学地位相对含量的引物对MiSeq测序的影响

使用PyNAST将序列对齐，共有3 520 条序列因比对失败而被剔除，因而共有441 023 条序列进行OTU划分。根据100%和97%的相似性进行UCLUST，共得到112 266 条代表性序列和7 484 个OTU，经ChimeraSlayer检测去除7 484 个OTU，每个样品平均588 个。采用RDP和Greengenes数据库，本研究将序列鉴定为19 个门、43 个纲、70 个目、152 个科和347 个属。样品的测序序列信息和各分类单元数量统计如表2所示。

表2 序列信息和各分类单元数量统计Table 2 Information abou t the sequences and taxonomy

由表2可知，经配对Kruskal-wallis检验发现，使用引物27F/338R扩增16S rRNA基因V1～V2区的样品细菌Chao 1指数显著高于其他两个可变区（P＜0.05），而不同引物扩增的同一样品间Shannon指数差异不显著（P＞0.05）。由此可见，使用不同引物分别对鲊广椒中细菌16S rRNA基因的V1～V2区、V3～V4区和V4～V5区进行扩增时可能会对样品的丰度产生影响。门水平相对含量对比分析如图3所示。

由图3可知，鲊广椒中平均相对含量大于1.0%的细菌门及其含量分别为：Firmicutes（硬壁菌门，75.60%）、Proteobacteria（变形菌门，17.57%）、Cyanobacteria（蓝细菌门，4.81%）和Bacteroidetes（拟杆菌门，1.02%）。经配对Kruskal-wallis检验发现，不同引物扩增的同一样品Cyanobacteria（蓝细菌门）含量存在显著差异（P=0.009），引物338F/806R扩增16S rRNA基因的V3～V4区产物中该细菌门相对含量显著偏高，而其他优势门差异均不显著（P＞0.05）。属水平相对含量对比分析如图4所示。

图3 门水平相对含量对比分析Fig. 3 Comparative analysis of relative abundances at phylum level

图4 属水平相对含量对比分析Fig. 4 Comparative analysis of relative abundances at genus level

由图4可知，鲊广椒中平均相对含量大于1.0%的细菌属及其含量分别为：隶属于Firmicutes（硬壁菌门）的Lactobacillus（乳酸杆菌，68.64%）、Weissella（魏斯氏菌，3.25%）、Staphylococcus（葡萄球菌，1.60%）和Pediococcus（小球菌，1.00%）；隶属于Bacteroidetes（拟杆菌门）的Carnimonas（肉胞菌属，3.79%）、Enterobacter（肠杆菌属，3.06%）和Kushneria（克锡勒氏菌，1.58%）。经配对Kruskal-wallis检验发现，不同引物扩增的同一样品Pediococcus（小球菌）含量存在显著差异（P=0.016），引物338F/806R扩增16S rRNA基因的V3～V4区产物中该菌相对含量显著偏高，而其他优势属差异均不显著（P＞0.05）。

图5 基于属水平的维恩图Fig. 5 Venn diagram based on genus level

由图5可知，使用3 种引物均可扩增出的细菌属有110 个，其包含的序列数占所有质控后合格序列数的87.89%。虽然使用引物27F/338R、338F/806R和515F/907R扩增16S rRNA基因的V1～V2、V3～V4和V4～V5区的样品中独有48、75 个和43 个属，但其累计包含的序列数仅占所有质控后合格序列数的0.22%。由此可见，虽然不同引物可能在某些相对含量较少的种系型上的扩增效率存在差异，但其均能对样品中的优势微生物菌群进行有效扩增。

2.3 基于群落结构的引物对MiSeq测序的影响

图6 不同引物扩增样品细菌多样性主坐标分析Fig. 6 Bacterial diversity of amplified products with different primer pairs based on principal coordinate analysis

由图6可知，在以2 个权重最高的主成分作图时，3 种引物扩增的样品空间排布无明显的区分趋势，其中第1和第2主成分分别占全部变量79.62%和13.70%的权重。在采用非监督多变量统计学方法进行分析的基础上，本研究选取了主坐标分析前80%的主成分进行基于马氏距离的聚类分析和多元方差分析，结果如图7所示。

由图7可知，较之使用引物27F/338R扩增16S rRNA基因的V1～V2区，由引物515F/907R和338F/806R扩增16S rRNA基因的V3～V4和V4～V5区的同一样品细菌多样性较为相似，但经MANOVA发现三者差异均不显著（P=0.520）。由此可知，使用不同引物扩增出的同一样品其微生物群落结构无显著差异。

图7 不同引物扩增样品细菌多样性马氏距离聚类Fig. 7 Bacterial diversity of amplified products with different primerpairs based on Mahalanobis distance cluster

3 结 论

采用MiSeq高通量测序技术对传统发酵食品的微生物多样性进行解析，要求研究人员需具备较为完善的分子生物学、生物信息学和多元统计学相关知识，而这些可能是食品科学相关领域研究人员欠缺的，因而目前相当一部分发酵食品微生物多样性解析项目的测序和数据分析工作委托高通量测序公司完成，且所使用的引物通常由测序公司提供。自2010年以来，笔者所在团队在肠道微生物和传统发酵食品微生物多样性解析 方面积累了一定的实验操作和数据处理经验，通过454焦磷酸测序技术研究地域[35]、食用益生菌[36]、饮食变化[37]和疾病[38]等因素对人体肠道菌群群落结构影响的研究工 作，通过MiSeq高通量测序技术完成了襄阳大头菜腌制液[39]和琚湾酸酱面浆水[40]等传统发酵食品的微生物解析。在项目开展的过程中，研究发现虽然2代测序技术的出现，对肠道菌群和发酵食品的研究起到了巨大的促进作用，但由于在目的片段扩增过程中使用的引物为通用引物，所以其更侧重于肠道和发酵食品中所有细 菌或真菌菌群群落结构的研究，在PCR扩增过程中常会对一些宿主或发酵基质的基因进行扩增。本研究发现使用引物27F/338R、338F/806R和515F/907R对5 个鲊广椒DNA样品16S rRNA基因的V1～V2区、V3～V4区和V4～V5区进行扩增时产生的非特异性条带分别占到序列总数的30.01%、22.04%和13.56%。由此可见，在使用上述3 种引物进行PCR扩增时均会产生非特异性扩增，在使用QIIME平台进行生物信息学分析时需根据鉴定结果对特异性序列进行删除。

由于测序长度的限制，研究人员建议在使用454焦磷酸测序技术进行微生物多样性解析时选择V6区为测序靶点，因为该可变区域较短且遗传信息较为丰富[41]。但也有研究人员建议根据样品的微生物构成进行引物的选择，例如土壤和粪便这类以Firmicutes和Proteobacteria细菌为主的样品进行MiSeq高通量测序时通常以V4～V5区为测序靶点，为避免微生物多样性被过高评估V1区和V6区则不建议选择[42]。本研究发现，鲊广椒中的细菌主要隶属于Firmicutes和Proteobacteria，其平均含量分别为75.60%和17.57%，使用引物27F/338R以V1～V2区为测序靶点进行细菌丰度评价时，其Chao 1指数显著大于以V4～V5区为测序靶点，由此可见，使用MiSeq高通量测序技术对鲊广椒微生物多样性进行解析时，使用V1～V2区为测序靶点易过高评估样本菌群的丰度。

本研究发现，在使用MiSeq高通量测序技术以16S rRNA为测序靶点进行鲊广椒细菌多样性研究时，使用引物515F/907R扩增16S rRNA基因V4～V5区的非特异性序列所占比例显著低于其他两个引物，且使用引物27F/338R扩增16S rRNA基因的V1～V2区易过高评估样本菌群的丰度，因而虽然使用不同引物扩增出的同一样品其微生物群落结构无显著差异，但建议选择引物515F/907R扩增16S rRNA基因的V4～V5区。由于缺少以16S rRNA基因全长为测序靶点的对照组，因而本研究在一定程度上存在不足之处。在后续研究中积极引入单分子实时测序技术[43]，以16S rRNA基因全长为测序靶点，进一步探讨扩增区域对鲊广椒细菌多样性解析的影响具有积极的意义。