王玉榮，楊成聰，葛東穎，尚雪嬌，張振東，郭 壯*

（湖北文理學院食品科學技術學院，鄂西北傳統(tǒng)發(fā)酵食品研究所，湖北 襄陽 441053）

我國傳統(tǒng)發(fā)酵食品的種類眾多且歷史悠久，受制作工藝、地域環(huán)境和原材料的影響，其風味和微生物群系極具多元化[1]。為實現(xiàn)傳統(tǒng)發(fā)酵食品的現(xiàn)代化和標準化生產，保障產品的質量安全，對傳統(tǒng)發(fā)酵食品生產過程中微生物群系變化進行解析則顯得極為重要[2]。相對于傳統(tǒng)微生物培養(yǎng)方法和指紋圖譜技術，以MiSeq技術為代表的第2代高通量測序技術具有無需培養(yǎng)、靈敏度高和檢測速度快等優(yōu)點，逐漸成為解析我國傳統(tǒng)發(fā)酵食品微生物多樣性的常規(guī)途徑[3]，目前在腌干魚[4]、濃香型白酒窖泥[5]、發(fā)酵黑茶[6]、酸奶[7]和茯磚茶[8]中有廣泛的應用。

由于MiSeq高通量測序技術的讀長不能完全覆蓋16S rRNA的全長基因，所以研究人員在對傳統(tǒng)發(fā)酵食品微生物多樣性進行解析時，常對其單可變區(qū)或多可變區(qū)進行測序[9]。16S rRNA基因擴增的準確性直接決定了高通量測序技術的可信度，然而研究發(fā)現(xiàn)DNA聚合酶[10]、擴增體系[11]和引物[12]均會對聚合酶鏈式反應（polymerase chain reaction，PCR）結果產生偏好性。引物的選擇對保證PCR擴增產物的特異性和質量尤為重要，只有選擇合適的引物才可避免異源雙鏈和單鏈DNA的形成[13]。目前研究人員通常使用引物擴增16S rRNA基因的V1～V2區(qū)[14]、V1～V3區(qū)[15]、V3區(qū)[16]、V3～V4區(qū)[17]、V4區(qū)[18]、V4～V5區(qū)[19]、V5～V6區(qū)[20]和V6～V8區(qū)[21]進行發(fā)酵食品細菌多樣性高通量測序分析，而Sun Donglei等[12]指出采用V4～V5區(qū)進行細菌多樣性解析時結果最接近真實情況。我國發(fā)酵食品眾多且不同發(fā)酵食品之間基質存在較大差異，因而引物的選擇尤為重要，其對細菌多樣性結果的解析可能存在較大影響，然而關于擴增區(qū)域對我國特色發(fā)酵食品細菌MiSeq測序結果影響的報道尚少。

鲊廣椒是以玉米面和新鮮辣椒為原料固體發(fā)酵而成的食品，在我國湖北、云南、重慶、四川和貴州等地區(qū)有著廣泛的食用人群[22]。本研究使用高通量測序公司常提供的3 對引物27F/338R、338F/806R和515F/907R分別對細菌16S rRNA的V1～V2區(qū)、V3～V4區(qū)和V4～V5區(qū)進行擴增并進行MiSeq高通量測序，同時采用生物信息學和多元統(tǒng)計學相結合的方法，對不同引物擴增PCR產物的細菌多樣性進行了解析，以期為食品科學領域研究人員在進行傳統(tǒng)發(fā)酵食品細菌多樣性解析時提供參考。

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化學發(fā)光凝膠成像系統(tǒng) 美國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 ℃?zhèn)溆谩?/p>

1.3.2 PCR擴增及產物純化

選取16S rRNA基因某一區(qū)域作為目片段進行擴增，為便于區(qū)分每個樣品的序列，在引物的前段分別添加含有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 個循環(huán)，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 序列質控

數(shù)據(jù)下機后，首先根據(jù)雙端序列的重疊關系將其拼接為一條序列，然后依據(jù)barcode信息將序列劃分到各樣品，在進行序列方向校正的基礎上，切除barcode和引物序列得到有效序列。在拼接過程中若成對序列符合下述條件中的一條則定義為低質量序列，在質控過程中予以剔除：1）重疊區(qū)的堿基數(shù)小于10 bp；2）重疊區(qū)的最大錯配比率大于0.2；3）barcode堿基錯配；4）超過1 個引物堿基發(fā)生錯配；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]數(shù)據(jù)庫對OTU中的代表性序列進行比對，在界、門、綱、目、科和屬水平上對每個樣品的細菌菌系進行解析；5）若某OTU無法鑒定到界水平或在屬水平上將其鑒定為玉米屬（Zea）則將其定義為非特異性擴增，使用自編腳本將該OTU中所包含的所有序列剔除，并依據(jù)barcode信息將序列劃分到各樣品，使用自編腳本對各序列的長度進行統(tǒng)計，并將所有樣品的序列文件合并為1 個fna文件后重復1～4步驟；6）使用FastTree軟件繪制系統(tǒng)發(fā)育進化樹[32]，計算Chao 1指數(shù)和Shannon指數(shù)[33]；7）基于UniFrac距離[34]進行主坐標分析和聚類分析。

1.3.6 核酸登錄號

含有非特異性擴增的序列數(shù)據(jù)已提交至MG-RAST數(shù)據(jù)庫（http：//metagenomics.anl.gov/），登錄號為mgp85472。

1.4 數(shù)據(jù)統(tǒng)計分析

使用配對Kruskal-wallis檢驗對各引物的非特異性擴增比例進行顯著性分析，使用Venny 2.1（http：//bioinfogp.c nb.csic.es/tools/venny/index.html）繪制維恩圖，進而對不同引物擴增產物中的核心OTU數(shù)目進行展示；使用多元方差分析（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區(qū)、V3～V4區(qū)和V4～V5區(qū)進行擴增時產生的非特異性條帶分別占到序列總數(shù)的30.01%、22.04%和13.56%。經配對Kruskal-wallis檢驗發(fā)現(xiàn)，擴增16S rRNA基因V4～V5區(qū)的非特異性序列所占比例顯著低于其他兩個可變區(qū)（P＜0.05）。在對特異性擴增序列進行剔除后，本研究共產出444 543 條高質量16S rRNA序列，平均每個樣品產出19 636 條，序列長度統(tǒng)計如圖2所示。

圖1 非特異性擴增序列的比例Fig. 1 Percentages of non-specific amplification sequences

圖2 序列長度統(tǒng)計Fig. 2 Statistic sequence lengths

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

2.2 基于各分類學地位相對含量的引物對MiSeq測序的影響

使用PyNAST將序列對齊，共有3 520 條序列因比對失敗而被剔除，因而共有441 023 條序列進行OTU劃分。根據(jù)100%和97%的相似性進行UCLUST，共得到112 266 條代表性序列和7 484 個OTU，經ChimeraSlayer檢測去除7 484 個OTU，每個樣品平均588 個。采用RDP和Greengenes數(shù)據(jù)庫，本研究將序列鑒定為19 個門、43 個綱、70 個目、152 個科和347 個屬。樣品的測序序列信息和各分類單元數(shù)量統(tǒng)計如表2所示。

表2 序列信息和各分類單元數(shù)量統(tǒng)計Table 2 Information abou t the sequences and taxonomy

由表2可知，經配對Kruskal-wallis檢驗發(fā)現(xiàn)，使用引物27F/338R擴增16S rRNA基因V1～V2區(qū)的樣品細菌Chao 1指數(shù)顯著高于其他兩個可變區(qū)（P＜0.05），而不同引物擴增的同一樣品間Shannon指數(shù)差異不顯著（P＞0.05）。由此可見，使用不同引物分別對鲊廣椒中細菌16S rRNA基因的V1～V2區(qū)、V3～V4區(qū)和V4～V5區(qū)進行擴增時可能會對樣品的豐度產生影響。門水平相對含量對比分析如圖3所示。

由圖3可知，鲊廣椒中平均相對含量大于1.0%的細菌門及其含量分別為：Firmicutes（硬壁菌門，75.60%）、Proteobacteria（變形菌門，17.57%）、Cyanobacteria（藍細菌門，4.81%）和Bacteroidetes（擬桿菌門，1.02%）。經配對Kruskal-wallis檢驗發(fā)現(xiàn)，不同引物擴增的同一樣品Cyanobacteria（藍細菌門）含量存在顯著差異（P=0.009），引物338F/806R擴增16S rRNA基因的V3～V4區(qū)產物中該細菌門相對含量顯著偏高，而其他優(yōu)勢門差異均不顯著（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檢驗發(fā)現(xiàn)，不同引物擴增的同一樣品Pediococcus（小球菌）含量存在顯著差異（P=0.016），引物338F/806R擴增16S rRNA基因的V3～V4區(qū)產物中該菌相對含量顯著偏高，而其他優(yōu)勢屬差異均不顯著（P＞0.05）。

圖5 基于屬水平的維恩圖Fig. 5 Venn diagram based on genus level

由圖5可知，使用3 種引物均可擴增出的細菌屬有110 個，其包含的序列數(shù)占所有質控后合格序列數(shù)的87.89%。雖然使用引物27F/338R、338F/806R和515F/907R擴增16S rRNA基因的V1～V2、V3～V4和V4～V5區(qū)的樣品中獨有48、75 個和43 個屬，但其累計包含的序列數(shù)僅占所有質控后合格序列數(shù)的0.22%。由此可見，雖然不同引物可能在某些相對含量較少的種系型上的擴增效率存在差異，但其均能對樣品中的優(yōu)勢微生物菌群進行有效擴增。

2.3 基于群落結構的引物對MiSeq測序的影響

圖6 不同引物擴增樣品細菌多樣性主坐標分析Fig. 6 Bacterial diversity of amplified products with different primer pairs based on principal coordinate analysis

由圖6可知，在以2 個權重最高的主成分作圖時，3 種引物擴增的樣品空間排布無明顯的區(qū)分趨勢，其中第1和第2主成分分別占全部變量79.62%和13.70%的權重。在采用非監(jiān)督多變量統(tǒng)計學方法進行分析的基礎上，本研究選取了主坐標分析前80%的主成分進行基于馬氏距離的聚類分析和多元方差分析，結果如圖7所示。

由圖7可知，較之使用引物27F/338R擴增16S rRNA基因的V1～V2區(qū)，由引物515F/907R和338F/806R擴增16S rRNA基因的V3～V4和V4～V5區(qū)的同一樣品細菌多樣性較為相似，但經MANOVA發(fā)現(xiàn)三者差異均不顯著（P=0.520）。由此可知，使用不同引物擴增出的同一樣品其微生物群落結構無顯著差異。

圖7 不同引物擴增樣品細菌多樣性馬氏距離聚類Fig. 7 Bacterial diversity of amplified products with different primerpairs based on Mahalanobis distance cluster

3 結 論

采用MiSeq高通量測序技術對傳統(tǒng)發(fā)酵食品的微生物多樣性進行解析，要求研究人員需具備較為完善的分子生物學、生物信息學和多元統(tǒng)計學相關知識，而這些可能是食品科學相關領域研究人員欠缺的，因而目前相當一部分發(fā)酵食品微生物多樣性解析項目的測序和數(shù)據(jù)分析工作委托高通量測序公司完成，且所使用的引物通常由測序公司提供。自2010年以來，筆者所在團隊在腸道微生物和傳統(tǒng)發(fā)酵食品微生物多樣性解析 方面積累了一定的實驗操作和數(shù)據(jù)處理經驗，通過454焦磷酸測序技術研究地域[35]、食用益生菌[36]、飲食變化[37]和疾病[38]等因素對人體腸道菌群群落結構影響的研究工 作，通過MiSeq高通量測序技術完成了襄陽大頭菜腌制液[39]和琚灣酸醬面漿水[40]等傳統(tǒng)發(fā)酵食品的微生物解析。在項目開展的過程中，研究發(fā)現(xiàn)雖然2代測序技術的出現(xiàn)，對腸道菌群和發(fā)酵食品的研究起到了巨大的促進作用，但由于在目的片段擴增過程中使用的引物為通用引物，所以其更側重于腸道和發(fā)酵食品中所有細 菌或真菌菌群群落結構的研究，在PCR擴增過程中常會對一些宿主或發(fā)酵基質的基因進行擴增。本研究發(fā)現(xiàn)使用引物27F/338R、338F/806R和515F/907R對5 個鲊廣椒DNA樣品16S rRNA基因的V1～V2區(qū)、V3～V4區(qū)和V4～V5區(qū)進行擴增時產生的非特異性條帶分別占到序列總數(shù)的30.01%、22.04%和13.56%。由此可見，在使用上述3 種引物進行PCR擴增時均會產生非特異性擴增，在使用QIIME平臺進行生物信息學分析時需根據(jù)鑒定結果對特異性序列進行刪除。

由于測序長度的限制，研究人員建議在使用454焦磷酸測序技術進行微生物多樣性解析時選擇V6區(qū)為測序靶點，因為該可變區(qū)域較短且遺傳信息較為豐富[41]。但也有研究人員建議根據(jù)樣品的微生物構成進行引物的選擇，例如土壤和糞便這類以Firmicutes和Proteobacteria細菌為主的樣品進行MiSeq高通量測序時通常以V4～V5區(qū)為測序靶點，為避免微生物多樣性被過高評估V1區(qū)和V6區(qū)則不建議選擇[42]。本研究發(fā)現(xiàn)，鲊廣椒中的細菌主要隸屬于Firmicutes和Proteobacteria，其平均含量分別為75.60%和17.57%，使用引物27F/338R以V1～V2區(qū)為測序靶點進行細菌豐度評價時，其Chao 1指數(shù)顯著大于以V4～V5區(qū)為測序靶點，由此可見，使用MiSeq高通量測序技術對鲊廣椒微生物多樣性進行解析時，使用V1～V2區(qū)為測序靶點易過高評估樣本菌群的豐度。

本研究發(fā)現(xiàn)，在使用MiSeq高通量測序技術以16S rRNA為測序靶點進行鲊廣椒細菌多樣性研究時，使用引物515F/907R擴增16S rRNA基因V4～V5區(qū)的非特異性序列所占比例顯著低于其他兩個引物，且使用引物27F/338R擴增16S rRNA基因的V1～V2區(qū)易過高評估樣本菌群的豐度，因而雖然使用不同引物擴增出的同一樣品其微生物群落結構無顯著差異，但建議選擇引物515F/907R擴增16S rRNA基因的V4～V5區(qū)。由于缺少以16S rRNA基因全長為測序靶點的對照組，因而本研究在一定程度上存在不足之處。在后續(xù)研究中積極引入單分子實時測序技術[43]，以16S rRNA基因全長為測序靶點，進一步探討擴增區(qū)域對鲊廣椒細菌多樣性解析的影響具有積極的意義。