焦晶凯，莫蓓红

（乳业生物技术国家重点实验室，光明乳业股份有限公司乳业研究院，上海 200436）

内蒙古呼伦贝尔草原以畜牧业为主，是传统奶制品的发源地之一，发酵乳制品种类繁多，奶酪是奶制品中最普遍的食品，为蒙古族人所喜爱，奶酪多以手工家庭制作为主，不添加任何商业菌株，以自然发酵为主，发酵菌株来自原料乳以及周围的环境，这些菌株赋予了奶制品特殊的风味和质地。已有研究显示，传统手工干酪中菌种丰富，包含乳杆菌（Lactobacillussp.）、明串珠菌（Leuconostocssp.）、肠球菌（Enterococcisp.）、嗜热链球菌（Streptococcus thermophilus）和乳酸乳球菌（Lactococcus lactis）等[1-4]。然而此前研究微生物多样性通常使用Sanger第一代测序技术，主要以DNA单链为模板，加入脱氧核糖核苷三磷酸（deoxy-ribonucleoside triphosphate，dNTP）和荧光双脱氧核苷三磷酸（ddNTP）进行聚合酶链式反应（polymerase chain reaction，PCR）扩增，然后进行电泳分离和激光诱导荧光颜色区分，获得碱基组成序列，但其成本高、测序通量低、信息不精确[5]。

第二代测序Illumina MiSeq方法分析有效的避免了通量低、操作复杂和准确率低等缺陷[6-9]，相对于Roche 454焦磷酸测序，具有操作简单、成本较低的优势，并且采用边合成边测序原理，结果可信度高，MiSeq高通量测序平台集中了Roche 454[10-11]和Illumina HiSeq 2500的优点，不仅可实现对多样品的多个可变区同时测序，而且在测序速度和测序通量上都有进一步提升，目前此平台已在微生物多样性群落结构研究方面受到了广大学者的认可[12-13]。全基因组测序对全面了解一个样本的微生物组成、分子进化等有着非常重要的意义。SCHMIDT P A等[14]使用Illumina高通量测序测定了土壤中真菌微生物多样性，并得到了较为理想的结果。赵爽[15]使用Illumnina高通量技术成功分析了对喷施氯苯嘧啶醇后的草坪根际土壤中真菌群落的多样性。

本研究使用Illumina MiSeq第二代测序方法测定来自内蒙古呼伦贝尔的4个地区的不同干酪样品，对其中微生物组成进行详尽的分析，同时分析了不同样品间物种差异和进化关系等。

1 材料和方法

1.1 材料与试剂

干酪样品来自于内蒙古呼伦贝尔的4个地区，分别是哈吉（编号C）、陈巴尔虎旗（编号D）、东乌珠尔（编号E）和西乌珠尔（编号F），这些样品均为牧民手工制作、自然成熟干酪，保存在4℃环境下，移至-70℃长期保存。

1.2 仪器与设备

Pico Green：美国Life Technologies公司；DP328粪便基因组DNA提取试剂盒：天根生化科技（北京）有限公司；AP-GX-500 DNA胶回收试剂盒：美国Axygen公司；P7589 PicoGreendsDNA 定量试剂盒：美国Invitrogen公司；MiSeq测序试剂盒v2、Illumina MiSeq测序仪：美国Illumina公司；DYY-6C电泳仪：北京天诚沃德生物技术有限公司；凝胶成像系统：美国伯乐公司；Pico17离心机：赛默飞世尔科技（中国）有限公司；Aligent 2100电泳生物分析仪：美国安捷伦科技有限公司；BioTek酶标仪：美国伯腾仪器有限公司。

1.3 试验方法

1.3.1 DNA提取

干酪样品中的微生物总DNA是用粪便基因组DNA提取试剂盒，根据说明书提取DNA。

1.3.2 PCR扩增

扩增16S rDNA的V4区域，引物520F：5′AYTGGGYDTAAAGNG3′，引物802R：5′TACNVGGGTATCTAATCC3′。PCR条件预变性是94℃5min，然后执行94℃30s、50℃30s、72℃30s共27个循环，退火温度72℃，时间是5min，最后保持在4℃条件下。1%琼脂糖电泳检测，切割回收，使用DNA胶回收试剂盒回收。

1.3.3 各样本定量

采用BioTek酶标仪对各个样品定量，取1μL样品加入24μL TE和25μL PicoGreen（稀释200倍）置于costar 96孔半量酶标仪上，在激发光485/20nm，发射光528/20nm下测定荧光强度。

1.3.4 DNA序列修饰

通过3′-5′核酸外切酶及聚合酶的共同作用，修复带有突出末端的DNA片段。在修复平整的DNA片段3′端引入单碱基“A”，接头3′末端含有单碱基“T”，从而保证DNA片段和接头能够通过“A”“T”互补配对连接，并防止接头连接DNA片段的过程中，DNA插入片段彼此相连，在连接酶的作用下，孵育含有标签的接头与DNA片段，使其相连。利用PCR选择性的富集两端连有接头的DNA片段，同时扩增DNA文库。PCR应尽量使用较少的环数，避免PCR扩增中文库出现错误。

1.3.5 验证和混合文库

利用PicoGreen和荧光分光光度计方法定量文库，使用Agilent 2100对PCR富集片段进行质量控制，验证DNA文库的片段大小分布。

1.3.6 上机测序

将混合好的文库（10nmol/L）逐步稀释定量至4～5pmol/L后进行上机测序，使用合成测序法，测定长度2×250bp。

1.3.7 数据分析

1.3.7.1 原始数据处理与样品序列数目统计

采用双峰（pair-end）测序，首先对原始数据进行质量控制，舍弃低质量序列（50个连续碱基平均质量＞Q30，不允许有N）。用软件Flash连接通过质量控制的序列对应的两端序列进行。设计无法连接的序列。对连接上的序列进行过滤（连续相同碱基＜6；模糊碱基N＜1），获得最终用于分析的序列。

1.3.7.2 OTU列表生成

应用Qiime，根据序列的相似度，将序列归为多个操作分类单元（operational taxonomic unit，OTU）。Qiime调用uclust对序列进行聚类，选取每个类最长的序列为代表序列。采用RDP-classifier，以RDP数据库的序列为训练集，对OTU代表序列进行注释。得到每个OTU的分类学信息。

1.3.7.3 稀释曲线及丰度分布

根据获得的OTU数据，作出每个样品的稀释曲线，以该曲线表明样品的取样深度。对各样品的OTU丰度大小排序，对丰度值取log2的对数作丰度分布曲线图，可提现样品中物种分布的均匀性。

1.3.7.4 群落结构分析

对OTU列表中获得的分类信息与丰度进行整理，在属层次下对各样品进行物种丰度统计、聚类分析及PCoA分析，可得到样品中群落组成结构及它们的相似性。

2 结果与分析

C、D、E、F4种样品分别读取了26 776、25 003、17 342、24 709条有效序列和26 704、24 927、17 287、24 640条优质序列。

2.1 样品所含操作分类单元及所含序列的数目

图1 操作分类单元维恩图Fig.1 Venn diagram of operational taxonomic unit

如图1所示，C样品中含有细菌种类541个，D样品中460个，E样品中含有563个，F样品中含有593个。C和F有相同种类230个，C和D有相同种类216个，C和E有相同种类218个，F和D有相同种类222个，F和E有相同种类236个，D和E有相同种类214个。C、F和D有相同种类156个，C、F和E有相同种类153个，C、D和E有相同种类151个，F、D和E有相同种类167个，在相似度0.97基础上分类出1 322个OTU，门1 142，纲1 057，目952，科716，属432。

2.2 物种丰度及群落结构分析

2.2.1 稀释曲线

图2表明每个样品的取样深度，从4种样品中随机抽取的测序条数如图2横轴所示，纵轴则表示基于该测序条数能构建的OUT数量，当曲线趋于平坦时，说明测序数据量合理，更多的数据量对发现新OTU的边际贡献很小；反之则表明继续测序还可能产生较多新的OTU。该分析是基于OTU序列差异水平在0.03即相似度97%的水平上进行运算的。

图2 4种样品的稀释曲线Fig.2 Dilution curve of four cheese samples

2.2.2 丰度分布曲线

对丰度值取log2的对数值得到图3所示的丰度分布曲线图。横轴所示为OTU相对丰度含量等级降序排列，纵轴对应的是OTU（物种）所占相对丰度比例，由图3可以看出4种样品的曲线斜率较大。

图3 4种样品的风度分布曲线Fig.3 Rank-abundance curve of four cheese samples

2.2.3 干酪中群落组成分析

在属的层次下对各样品物种丰度进行统计，如图4所示，4个样品所含物种基本相似，但物种分布差异较大，对所含主要种属进行了分析和归类。由图4可知，主要物种包括3类，乳酸杆菌（Lactobacillus）、醋酸菌（Acetobacter）和乳酸球菌（Lactococcus）；C样品与其他3个样品的不同在于其含量最多物种为醋酸菌醋酸菌（Acetobacter）而其他3个样品为乳酸杆菌（Lactobacillus）；F和D样品较为类似；E样品相对于其他样品来说，含有较少的醋酸菌，其他物种含量相对丰富些；此外在4种样品中还发现有巨型球菌（Macrococcus）、明串珠菌（Leuconostoc）、链球菌（Streptococcus）、肠杆菌（Enterobacteriaceae）和假单胞菌（Pseudomonas）等。

图4 4种干酪样品物种分布图Fig.4 Species distribution in four cheese samples

2.2.4 聚类分析

本实验在属的水平上利用Heapmap图分析了4个样品菌群的相似性和多样性，如图5所示，所含菌落分为3个大的分支，其中最为优势的OTU单元为硬壁菌门（Firmicutes）；OTU925包含序列数最多，其34 227条序列均属于乳酸杆菌（Lactobacillus）。由图5可知，在OTU861和9254包含序列也较多，分别属于醋酸杆菌（Acetobacter）和乳酸球菌（Lactococcus），OTU58、OTU717、OTU73、OTU479和OTU 1178等也属于乳酸杆菌（Lactobacillus），此外序列较多的还包含假单胞菌（Pseudomonas）、肠杆菌（Enterobacteriaceae）、链球菌（Streptococcus）和明串珠菌（Leuconostoc）等。

图5 系统发育进化树Fig.5 Phylogenetic tree analysis base on 16S rDNA gene of four cheese samples

2.2.5 主坐标PCoA分析

如图6所示，利用各样品序列间的进化信息计算样品距离，观察4种不同干酪样品多样性差异。由图6可知，PCoA分析中这些站点没有聚在一起，距离相对较远，表明4个干酪样品群落存在一定的差异性，因此，地域对各个群落组成有重要影响。

图6 4种干酪样品群落主坐标分析Fig.6 Principal coordinates analysis of microbiota in four cheese samples

3 结论

采用Illumina MiSeq测序方法，直接从环境样品中扩增高变区进行测序，避免了不可培养菌株的不可测性，客观还原菌群结构及丰度比例，严格控制PCR循环数，充分发挥高通量测序的大数据量优势。

结果表明，4个干酪样品共分离了细菌门1 142个，纲1 057个，目952个，科716个，属432个，主要包括乳酸杆菌（Lactobacillus）、醋酸杆菌（Acetobacter）和乳酸球菌（Lactococcus），本实验的稀释曲线分析了4种样品取样深度，每个样品都呈现相同或相似的变化趋势，在测序量较少时，每个样品OTU数目呈显著上升的趋势，而随着测序数量的增加，每个样品的OTU数目增加趋势逐渐变缓，最后基本达到饱和。由本实验4种样品的丰度曲线可以看出，相对于环境、土壤等样品来说曲线较窄，说明干酪中的菌种相对来说较为单一，曲线斜率较大，由此推算4种干酪样品中存在优势菌株。通过菌落组成分析，F、D、E 3个样品的主要优势菌为乳酸杆菌，所占比例均超过了50%，而C样品中醋酸菌所占比例超过50%，另外一种优势菌是乳酸球菌。PCoA分析可以看出4个样品菌群分布大不相同，说明地域对其微生物组成影响较大。

[1]CALLON C,MILLET L,MONTEL M C.Diversity of lactic acid bacteria isolated from AOC Salers cheese[J].J Dairy Res,2004,71(2):231-244.

[2]CHEBE OVÁ-TURCOVSKÁ V,ŽENIŠOVÁ K,KUCHTA T,et al.Culture-independent detection of microorganisms in traditional Slovakian bryndza cheese[J].Int J Food Microbiol,2011,150(1):73-78.

[3]NIKOLIC M,TERZIC-VIDOJEVIC A,JOVCIC B,et al.Characterization of lactic acid bacteria isolated from Bukuljac,a homemade goat's milk cheese[J].Int J Food Microbiol,2008,122(1-2):162-170.

[4]TERZIC-VIDOJEVIC A,VUKASINOVIC M,VELJOVIC K,et al.Characterization of microflora in homemade semi-hard white Zlatar cheese[J].Int J Food Microbiol,2007,114(1):36-42.

[5]秦 楠，栗东芳，杨瑞馥.高通量测序技术及其在微生物学研究中的应用[J].微生物学报，2011，51（4）：445-457.

[6]RUTVISUTTINUNT W,CHINNAWIROTPISAN P,SIMASATHIEN S,et al.Simultaneous and complete genome sequencing of influenza A and B with high coverage by Illumina MiSeq Platform[J].J Virol Methods,2013,193(2):394-404.

[7]WILLIAMS S,FOSTER P,LITTLEWOOD D.The complete mitochondrial genome of a turbinid vetigastropod from MiSeq Illumina sequencing of genomic DNA and steps towards a resolved gastropod phylogeny[J].Gene,2014,533(1):38-47.

[8]蔡元锋，贾仲君.基于新一代高通量测序的环境微生物转录组学研究进展[J].生物多样性，2013（4）：402-411.

[9]王兴春，杨致荣，王 敏，等.高通量测序技术及其应用[J].中国生物工程杂志，2012，32（1）：109-114.

[10]段 曌，肖 炜，王永霞，等.454 测序技术在微生物生态学研究中的应用[J].微生物学杂志，2011，31（5）：76-81.

[11]徐晓宇，刘 和.454 测序法在环境微生物生态研究中的应用[J].生物技术通报，2010（1）：73-77.

[12]董 逸.我国黄、东海典型海域微生物群落结构及其与环境变化的关系[D].青岛：中国科学院研究生院（海洋研究所）博士论文，2013.

[13]周林文.MiSeq：新一代个人化测序仪[J].生物技术世界，2011（10）：14-16.

[14]SCHMIDT P A,BÁLINT M,GRESHAKE B,et al.Illumina metabarcoding of a soil fungal community[J].Soil Biol Biochem,2013,65:128-132.

[15]赵 爽.土壤病原真菌的检测及土壤中真菌群落多样性的研究[D].南京：南京农业大学博士论文，2012.