Illumina MiSeq平臺高覆蓋率測定干酪中的細菌微生物多樣性

焦晶凱，莫蓓紅

（乳業生物技術國家重點實驗室，光明乳業股份有限公司乳業研究院，上海 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.