基于高通量測序技術(shù)比較不同來源原料乳中的細(xì)菌特征

周繼福，王娉，郭佳,3，趙曉美，陳穎

（1.中國檢驗檢疫科學(xué)研究院，北京 100176）（2.南京財經(jīng)大學(xué)食品科學(xué)與工程學(xué)院，江蘇南京 210023）（3.天津科技大學(xué)食品科學(xué)與工程學(xué)院，天津 300457）

原料乳是乳腺直接分泌、不做加工處理且富含蛋白質(zhì)、碳水化合物、脂肪酸和多種微量營養(yǎng)素的生鮮乳，也是唯一能為新生哺乳動物提供足量營養(yǎng)物質(zhì)的食物[1]。研究表明，原料乳已成為病原微生物傳播的重要途徑，Mungai等[2]發(fā)現(xiàn)美國2007～2012年非巴氏殺菌乳中病原菌引起的疾病暴發(fā)率是1993～2006年的三倍。美國及歐洲多國已制定相關(guān)原料乳飲食指南[3,4]，我國也制定頒布了GB 19301-2010《生乳》等一系列國家標(biāo)準(zhǔn)保證原料乳的質(zhì)量安全，規(guī)定了生乳中微生物、真菌等的限量及檢測標(biāo)準(zhǔn)，原料乳中微生物的檢測對奶牛健康狀況監(jiān)測、規(guī)范擠乳操作和牛奶質(zhì)量監(jiān)管有積極作用。目前用于原料乳中微生物鑒定的技術(shù)和方法包括流式細(xì)胞術(shù)[5]，酶聯(lián)免疫吸附測定技術(shù)[6]，平板培養(yǎng)法，基于PCR技術(shù)的（GTG）5-PCR指紋圖譜[7]、結(jié)合激光誘導(dǎo)熒光的微芯片電泳（microchip electrophoresis coupled with laser-induced fluorescence，MCE-LIF）[8]和環(huán)介導(dǎo)等溫擴增[9]等方法以及生化傳感器[10]，但這些依賴于微生物培養(yǎng)和儀器的特異性、穩(wěn)定性。Carvalho等[11]認(rèn)為實驗室條件下可培養(yǎng)的細(xì)菌不足1%。因此亟需采用不依賴于培養(yǎng)的方法對原料乳中的微生物多樣性進(jìn)行探究。

基于高通量測序技術(shù)的宏基因組學(xué)技術(shù)可對生態(tài)位中提取的所有微生物DNA進(jìn)行測序和分析，鑒定細(xì)菌種類和表征群落豐度。同時，技術(shù)的發(fā)展和成熟也使得高通量測序更多的用于乳品微生物群落分析，Dalmasso等[12]利用高通量測序技術(shù)分析奶酪中的微生物群落動態(tài)，Zhang等[13]運用高通量測序?qū)Σ煌a(chǎn)地山羊奶微生物多樣性進(jìn)行了研究。

本研究采用高通量測序技術(shù)分析不同地區(qū)、不同季節(jié)原料乳中微生物多樣性，比較不同來源原料乳之間的細(xì)菌群落結(jié)構(gòu)和相對豐度的差異，明確優(yōu)勢菌群和顯著性差異種屬，并對細(xì)菌功能基因代謝途徑進(jìn)行預(yù)測，為原料乳質(zhì)量和安全性提供相關(guān)的微生物信息。

1 材料與方法

1.1 材料與試劑

原料乳樣品為本實驗室采集，分別來自北京、河北、內(nèi)蒙古、天津四地，原料乳奶槽車載量約為30 t。

QIAamp?Fast DNA Stool Mini kit，德國Qiagen公司；Qubit3.0 DNA檢測試劑盒，美國life公司；2×TaqMaster Mix，南京Vazyme公司；MagicPure?Size Selection DNA Beads，北京全式金生物技術(shù)（TransGen Biotech）有限公司。

1.2 儀器與設(shè)備

高速冷凍離心機，德國D-37520 Osterode；Qubit?3.0熒光計，美國Invitrogen公司；聚合酶鏈?zhǔn)椒磻?yīng)（polymerase chain reaction，PCR）儀，美國BIO-RAD公司；凝膠成像系統(tǒng)，美國UVP公司；IlluminaMiSeq測序系統(tǒng)，美國Illumina公司。

1.3 方法

1.3.1 樣本采集

表1 原料乳主要信息表Table 1 Sampling information of raw milk

23份原料乳樣品收集自北京（BJ，n=11）、河北（HE，n=6）、內(nèi)蒙古（NM，n=5）和天津（TJ，n=1）四地，采樣時間跨夏（S，n=4）、秋（F，n=12）、冬（W，n=7）三季（表1），采樣時嚴(yán)格按照工廠微生物檢測要求進(jìn)行，冷藏條件下送至實驗室在-20 ℃下儲存并盡快完成細(xì)菌基因組DNA提取。

1.3.2 基因組DNA提取

使用Fast DNA Stool Mini kit提取原料乳中細(xì)菌DNA，略作修改，將離心獲得的沉淀中加入緩沖液（InhibitEX buffer）的體積從1000 μL減少為200 μL，以保證最終提取的DNA濃度，其余操作步驟按照試劑盒說明書進(jìn)行。提取的DNA低溫（-20 ℃）儲存待用。

1.3.3 高通量測序

利用Qubit熒光計對基因組DNA定量以確定PCR擴增體系中模板DNA的量，擴增區(qū)域為細(xì)菌16S rRNA 基因的 “ V3-V4 ” 區(qū)，引物為 341F：5’-CCTACGGGNGGCWGCAG-3’ ， 805R ：5’-GACTACHVGGGTATCTAATCC-3’[14]。擴增產(chǎn)物經(jīng)過瓊脂糖凝膠電泳檢測后使用核酸純化試劑盒對DNA進(jìn)行純化回收，Illumina MiSeq測序平臺用于完成高通量測序。

1.3.4 生物信息和數(shù)據(jù)統(tǒng)計分析

對23份原料乳樣品進(jìn)行測序，獲得的原始圖像數(shù)據(jù)經(jīng)過堿基識別后形成FASTQ格式的測序序列。剔除引物、接頭序列、PCR過程產(chǎn)生的錯誤序列以及低閾值序列。依據(jù)Filtered序列間的距離和相似性，采用軟件Usearch（version 5.2.236）進(jìn)行聚類形成操作分類單元（operational taxonomic units，OTU），聚類后的OTU常在97%的相似度下用于生物信息統(tǒng)計分析。Alpha多樣性指數(shù)用以衡量微生物群落多樣性和分布豐度。使用RDP classifier（version 2.12，數(shù)據(jù)庫：http://rdp.cme.msu.edu/misc/reso urces.jsp）和Blast分類方法對OTU分類序列進(jìn)行物種豐度計算，獲得OTU分類序列矩陣并繪制群落結(jié)構(gòu)分布柱狀圖。根據(jù)矩陣信息利用統(tǒng)計檢驗篩選差異物種分類單元，結(jié)果以誤差線圖表示。此外，為比較多樣本間的差異，進(jìn)行Beta多樣性分析。最后，基于京都基因與基因組百科全書（kyoto encyclopedia of genes and genomes，KEGG）對原料乳的細(xì)菌功能基因進(jìn)行預(yù)測。統(tǒng)計軟件包括STAMP（version 2.1.3）和Excel，當(dāng)p＜0.05時認(rèn)為具有統(tǒng)計學(xué)意義。

2 結(jié)果與分析

2.1 高通量測序信息

對高通量測序數(shù)據(jù)處理后共獲得1 722 199條高質(zhì)量序列，質(zhì)控后的序列平均長度為421.17 bp。OTU微生物物種分類結(jié)果顯示，原料乳中微生物共分為33門，64綱，117目，254科，732屬。

2.2 OTU聚類和分類注釋

為獲得樣品中微生物群落結(jié)構(gòu)組成信息，對高通量測序序列進(jìn)行歸類。韋恩圖（圖1）直觀展示微生物共有或獨特的OTU數(shù)目，北京、內(nèi)蒙古、河北、天津四地的OTU數(shù)分別為6787、3856、3876、251；獨有OTU數(shù)目分別為3020、981、935、52，說明四地區(qū)中，北京地區(qū)乳樣中微生物群落豐度最高。四地乳樣共有OTU數(shù)為70，北京和天津共有OTU數(shù)為178，北京和內(nèi)蒙古共有OTU數(shù)為2368，北京和河北共有OTU數(shù)為2441，內(nèi)蒙古和河北共有OTU數(shù)為1597，說明四地區(qū)中，北京和河北的乳樣最相似。夏、秋、冬三季的OTU數(shù)分別為2193、5254、6132，三個季節(jié)樣品共有OTU數(shù)為707，表明冬季乳樣中微生物群落豐度高于夏季和秋季；夏季和秋季共有的OTU數(shù)為1096，夏季和冬季的共有OTU數(shù)為1142，秋季和冬季的共有OTU數(shù)為2780，說明秋季和冬季的乳樣更相似。

圖1 不同來源原料乳中細(xì)菌OTUs分布韋恩圖Fig.1 Venn diagram of bacterial OTUs distribution in raw milk from different sources

使用物種分類方法統(tǒng)計門（phylum）、屬（genus）層次上的原料乳微生物群落組成，繪制門、屬水平上的原料乳微生物群落結(jié)構(gòu)分布柱狀圖，繪圖時將門、屬分類水平上相對豐度＜0.50%和＜5.00%的細(xì)菌群落及unclassified劃分為其他菌屬（others）[15]，以直觀反映原料乳中主要微生物數(shù)量及群落構(gòu)成。

在門分類水平上的原料乳中微生物細(xì)菌菌門和相對豐度各不相同，圖2顯示所有乳樣中均含有變形菌門（Proteobacteria），厚壁菌門（Firmicutes），擬桿菌門（Bacteroidetes）和放線菌門（Actinobacteria），平均相對豐度分別為76.54%，11.17%，7.10%和4.56%。此外，所有原料乳中還含有 Candidatus Saccharibacteria，綠彎菌門（Chloroflexi），芽單胞菌門（Gemmatimonadetes），異常球菌屬-棲熱菌門（Deinococcus-Thermus），柔壁菌門（Tenericutes）和梭桿菌門（Fusobacteria）。Nazmul等[16]在乳腺炎乳樣和健康乳樣中均發(fā)現(xiàn)四種菌門的存在，Rozo等[17]在原料乳制成的奶酪中也得到同樣的結(jié)果。此外，Sadan等[18]發(fā)現(xiàn)牛瘤胃中主要的菌群為厚壁菌門（Firmicutes）和擬桿菌門（Bacteroidetes），而瘤胃是分解飼料的主要場所，這間接說明厚壁菌門和擬桿菌門為原料乳中的常見微生物群落。

圖2 基于門水平的原料乳細(xì)菌群落分布柱狀圖Fig.2 Barplot of bacterial community distribution in raw milk based on phylum level

從樣品地區(qū)來源看，四個地區(qū)原料乳樣品中細(xì)菌主要來自于變形菌門，其相對豐度在68.99%～84.08%之間，擬桿菌門和放線菌門的細(xì)菌相對較少，擬桿菌門的相對豐度在10.00%以下，放線菌門的相對豐度則在5.00%以下；河北地區(qū)原料乳樣品中的厚壁菌門細(xì)菌的相對豐度比其他三個地區(qū)的均高2倍以上，達(dá)到20.99%。河北地區(qū)和內(nèi)蒙古地區(qū)乳樣中擬桿菌門和放線菌門細(xì)菌相對含量相似。從原料乳采樣季節(jié)來看，三個季節(jié)原料乳中細(xì)菌主要為變形菌門，平均相對豐度在62.54%～83.54%之間，夏季原料乳中厚壁菌門細(xì)菌的相對豐度是秋、冬兩季的2倍以上，為23.23%，三個季節(jié)的原料乳中擬桿菌門和放線菌門細(xì)菌的平均相對豐度較低，放線菌門的平均相對豐度在10.00%以下。冬季原料乳中放線菌門的相對含量高于擬桿菌門。夏季原料乳中擬桿菌門的相對豐度是秋季的2倍左右。除此之外，原料乳中相對豐度較低的CandidatusSaccharibacteria，綠彎菌門，芽單胞菌門，異常球菌屬-棲熱菌門，柔壁菌門和梭桿菌門在內(nèi)蒙古地區(qū)的秋季乳樣中均未發(fā)現(xiàn)。柔壁菌門在河北的冬季乳樣中不存在，這表明各地區(qū)在不同季節(jié)原料乳微生物群落結(jié)構(gòu)組成上存在差異。

在屬分類層次上（圖3），原料乳主要的菌屬為假單胞菌屬（Pseudomonas），不動桿菌屬（Acinetobacter），黃桿菌屬（Flavobacterium），水棲菌屬（Enhydrobacter），乳球菌屬（Lactococcus），鏈球菌屬（Streptococcus），苯基桿菌屬（Phenylobacterium），耶爾森氏菌屬（Yersinia），氣單胞菌屬（Aeromonas），考克氏菌屬（Kocuria），羅斯氏菌屬（Rothia），檸檬酸桿菌屬（Citrobacter），葡萄球菌屬（Staphylococcus），Epilithonimonas和Rheinheimera。其中假單胞菌屬（Pseudomonas），不動桿菌屬（Acinetobacter），黃桿菌屬（Flavobacterium）和水棲菌屬（Enhydrobacter）為原料乳共有菌屬，平均相對豐度分別為32.02%，28.53%，5.35%和3.81%。假單胞菌屬（Pseudomonas）和不動桿菌屬（Acinetobacter）的平均相對豐度之和占總微生物數(shù)的60.55%，為原料乳中兩種優(yōu)勢菌屬，該結(jié)果與已有報道類似[13,19]，于國萍等[20]發(fā)現(xiàn)不同牧場原料乳的主要菌屬中包含有不動桿菌屬和乳球菌屬。

圖3 基于屬水平的原料乳細(xì)菌群落分布柱狀圖Fig.3 Barplot of bacterial community distribution in raw milk based on genus level

按照不同地區(qū)分析發(fā)現(xiàn)，除天津地區(qū)外，其余三地區(qū)原料乳中的優(yōu)勢菌屬均是假單胞菌屬，相對豐度在30.25%～41.42%之間。不動桿菌屬作為天津地區(qū)乳樣中的主要菌屬，相對豐度為43.62%，明顯高于其余三地區(qū)原料乳中不動桿菌的豐度。除上述兩種菌屬外，北京地區(qū)乳樣中還包含相對豐度為6.03%的黃桿菌，河北地區(qū)乳樣中包含相對豐度為12.99%和6.09%的乳球菌屬和鏈球菌屬，這些在內(nèi)蒙古地區(qū)的乳樣中均不存在。此外，葡萄球菌屬和Rheinheimera也只在天津地區(qū)乳樣中存在，相對豐度分別為6.84%和5.65%。河北和內(nèi)蒙古原料乳中均不存在檸檬酸桿菌屬，北京、河北兩地的冬季乳樣中也不存在氣單胞菌屬。從不同季節(jié)來看，三個季節(jié)的原料乳樣品中均含有假單胞菌屬和不動桿菌屬，除夏季乳樣中的假單胞菌屬為7.53%外，其余季節(jié)的兩種菌屬的相對豐度均在22.95%以上。此外，夏季原料乳中還存在乳球菌屬和黃桿菌屬，相對豐度分別為16.82%和9.89%。

圖4 不同季節(jié)原料乳中微生物豐度誤差線圖Fig.4 Extended error bar of microbial abundance in different seasons

表2 原料乳樣品主要測序信息及α多樣性指數(shù)Table 2 Main sequencing information and α-diversity index of raw milk samples

注：BJ：北京地區(qū)；HE：河北地區(qū)；NM：內(nèi)蒙古地區(qū)；TJ：天津地區(qū)；S：夏季； F：秋季；W：冬季。

有研究表明，原料乳在冷藏條件下儲存一天，嗜冷菌的數(shù)量可以急劇增長到總微生物數(shù)量的90.00%以上[21,22]，這或許是假單胞菌屬和不動桿菌屬能成為原料乳優(yōu)勢菌屬的原因之一，實驗中分析發(fā)現(xiàn)的主要嗜冷菌包括假單胞菌屬，不動桿菌屬，氣單胞菌屬，黃桿菌屬，乳球菌屬，葡萄球菌屬和鏈球菌屬[23,24]。嗜冷菌除了產(chǎn)生影響原料乳質(zhì)量的酶以外，還可能存在潛在的致病性和抗生素耐藥性[21，Danmallam等[25]在患乳腺炎奶牛的牛奶中發(fā)現(xiàn)葡萄球菌屬細(xì)菌的存在。Decimo等[26]發(fā)現(xiàn)散裝牛奶中分離的假單胞菌屬對氨曲南展現(xiàn)出耐藥性。

根據(jù)采樣地區(qū)和采樣季節(jié)對原料乳中微生物相對豐度進(jìn)行顯著性差異分析，不同地區(qū)微生物群落雖占比不同，但并未體現(xiàn)出顯著性差異。而對不同季節(jié)原料乳中微生物群落相對豐度差異比較發(fā)現(xiàn)（圖4），夏秋、夏冬季節(jié)的假單胞菌屬（Pseudomonas）相對豐度存在顯著性差異（p＜0.05）（圖4a，4b），且差異比例超過30.00%。秋冬季節(jié)的氣單胞菌屬（Aeromonas）相對豐度也存在顯著性差異（p＜0.05）（圖4c）。

2.3 Alpha多樣性分析

Alpha多樣性既可以表示每個樣品微生物群落豐度，也能表示微生物群落多樣性，其衡量方法往往以指數(shù)的形式呈現(xiàn)（表2）。利用軟件Mothur（version 1.30.1）建立的Coverage指數(shù)范圍在0.99～1.00之內(nèi)，表明此次高通量測序可以真實反映原料乳中的微生物群落分布信息。對不同季節(jié)原料乳群落分布豐度指數(shù)（Chao和ACE）進(jìn)行統(tǒng)計分析，發(fā)現(xiàn)冬季乳樣的豐度指數(shù)差異顯著高于夏季原料乳（p＜0.05）。BJ-S1，NM-F1和TJ-S乳樣群落分布豐度指數(shù)低于其他同地區(qū)同季節(jié)樣品，這說明微生物物種總數(shù)可能受牧場規(guī)模、奶牛自身狀況以及季節(jié)等影響。

微生物群落分布多樣性指數(shù)（Shannon和Simpson）中，Simpson指數(shù)越小代表物種多樣性越高。對比不同季節(jié)、不同產(chǎn)地原料乳的多樣性指數(shù)并未發(fā)現(xiàn)顯著性差異。收集自冬季的兩份北京地區(qū)乳樣（BJ-W1和BJ-W3）的Shannon指數(shù)出現(xiàn)極端的大小差異，分別為5.94和0.70，這說明同一地區(qū)相同季節(jié)的原料乳微生物群落多樣性可能和牧場衛(wèi)生管理、奶牛自身健康、泌乳后操作及菌落總數(shù)的高低相關(guān)[27]。

2.4 Beta多樣性分析

圖5 基于UniFrac的原料乳中細(xì)菌主坐標(biāo)圖Fig.5 Principal coordinates of bacteria in raw milk based on UniFrac

Beta多樣性可用于分析原料乳樣本在特定的進(jìn)化譜系中的微生物群落差異。使用UniFrac對Beta多樣性進(jìn)行評估分析，圖5顯示主坐標(biāo)成分1（PCoA1）和主坐標(biāo)成分2（PCoA2）的貢獻(xiàn)率分別為34.00%和17.00%。圖中不同地區(qū)和不同季節(jié)原料乳樣品的聚類趨勢較弱，這說明不同地區(qū)和不同季節(jié)乳樣中的細(xì)菌群落差異不明顯。

2.5 功能預(yù)測分析

基于KEGG數(shù)據(jù)庫，采用PICRUSt軟件（version 1.0.0）對原料乳中微生物的功能基因代謝途徑進(jìn)行對比預(yù)測。結(jié)果分析發(fā)現(xiàn)：主要的功能基因相關(guān)代謝途徑為氨基酸代謝（amino acid metabolism），碳水化合物代謝（carbohydrate metabolism）和膜運輸（membrane transport）（圖6）。氨基酸代謝及碳水化合物代謝的豐度較高可能同原料乳中氨基酸和碳水化合物的高含量[13]及其生物功能相關(guān)，而膜運輸則可能與微生物生長活動所需的酶等物質(zhì)的運輸相關(guān)[28]。氨基酸在原料乳的蛋白質(zhì)、脂肪和乳糖生物合成中發(fā)揮重要作用，這也是氨基酸代謝活躍的一個原因。包括乳球菌和鏈球菌在內(nèi)的微生物對乳糖的利用也促進(jìn)了原料乳中碳水化合物的代謝。Zhang等[13]認(rèn)為乳的組成成分和微生物群落存在某種潛在聯(lián)系。吳建民等[29]發(fā)現(xiàn)奶牛瘤胃中的微生物通過增加氨基酸代謝對乳蛋白進(jìn)行調(diào)節(jié)。有報道稱，細(xì)胞信號轉(zhuǎn)導(dǎo)（cell communication）[30]和心血管疾?。╟ardiovascular diseases）[31]等代謝途徑曾在牛奶微生物中發(fā)現(xiàn)。但本研究中并未發(fā)現(xiàn)感覺系統(tǒng)（sensory system），細(xì)胞信號轉(zhuǎn)導(dǎo)甚至心血管疾病等代謝途徑。

在對代謝途徑的豐度差異分析中發(fā)現(xiàn)了一些同人類疾病相關(guān)的代謝途徑，并且在部分產(chǎn)地和季節(jié)樣品中存在顯著性差異。其中，北京-內(nèi)蒙古、河北-內(nèi)蒙古乳樣中的微生物在神經(jīng)退行性疾病（neurodegenerative diseases），循環(huán)系統(tǒng)（circulatory system）和癌癥（cancers）的功能豐度上存在顯著性差異（p＜0.05），神經(jīng)退行性疾病在內(nèi)蒙古樣品中的平均差異比例略高于北京、河北兩地。夏秋、夏冬季節(jié)樣品在神經(jīng)系統(tǒng)（nervous system），內(nèi)分泌系統(tǒng)（endocrine system）的功能豐度上也存在顯著性差異（p＜0.05）。這可能是原料乳中某些微生物和病原微生物有同源基因所導(dǎo)致的[28]。另外，夏秋、夏冬季節(jié)樣品在核苷酸代謝（nucleotide metabolism），聚糖的生物合成與代謝（glycan biosynthesis and metabolism）和遺傳信息處理（genetic information processing）也存在顯著性差異（p＜0.05）。

圖6 不同來源原料乳中微生物功能豐度熱圖Fig.6 Heatmap of microbial functional abundance of raw milk from different sources

3 結(jié)論

本實驗采用高通量測序技術(shù)對原料乳中微生物群落進(jìn)行分析，通過探究不同采樣季節(jié)、不同生產(chǎn)地區(qū)原料乳的優(yōu)勢菌群及豐度差異，并對功能基因代謝途徑進(jìn)行分析，以期為原料乳品質(zhì)及后續(xù)加工提供相關(guān)微生物信息。結(jié)果顯示，23份原料乳的優(yōu)勢菌門為變形菌門和厚壁菌門，優(yōu)勢菌屬除假單胞菌屬和不動桿菌屬外，還包括乳球菌屬和鏈球菌屬等其他主要菌屬。在門水平上，不同季節(jié)和地區(qū)原料乳菌群相對豐度比例雖然不同，但是并未體現(xiàn)出顯著差異。在屬水平上，不同季節(jié)和地區(qū)原料乳中部分微生物相對豐度存在顯著性差異（p＜0.05）。在功能預(yù)測分析中，氨基酸代謝、碳水化合物代謝和膜運輸為原料乳細(xì)菌主要代謝途徑，而研究表明這與微生物群落有關(guān)。研究并掌握原料乳中微生物組成，對于原料乳的后續(xù)生產(chǎn)和品質(zhì)評價具有重要意義。