基于高通量測(cè)序的液態(tài)奶中微生物多樣性研究

屠大偉,張清平,劉美艷,張巨鳳,翁盈秋

1.重慶工商大學(xué) 環(huán)境與資源學(xué)院，重慶 400067 2.重慶市天友乳業(yè)股份有限公司，重慶 400015 3.重慶萬(wàn)標(biāo)檢測(cè)技術(shù)有限公司，重慶 400714

1 引 言

牛乳因其營(yíng)養(yǎng)豐富且易被人體吸收而成為人民日常生活中重要的營(yíng)養(yǎng)來(lái)源。牛乳中含有豐富的營(yíng)養(yǎng)物質(zhì)，但在擠奶、原料乳的預(yù)處理、貯存以及運(yùn)輸?shù)冗^(guò)程中處理不當(dāng)會(huì)導(dǎo)致奶源中有害微生物的生長(zhǎng)和繁殖，從而造成牛乳中營(yíng)養(yǎng)物質(zhì)的破壞，影響原料乳的品質(zhì)以及食源性疾病的傳播[1]。

目前，常用于研究乳制品中微生物多樣性的方法主要有常規(guī)細(xì)菌鑒定技術(shù)、免疫學(xué)檢測(cè)技術(shù)、基因檢測(cè)技術(shù)和全基因組測(cè)序技術(shù)[2]等，傳統(tǒng)的檢測(cè)方法存在操作復(fù)雜，檢測(cè)靈敏度低的現(xiàn)象[3]。與傳統(tǒng)方法相比，高通量測(cè)序技術(shù)可以高效、準(zhǔn)確、全面地分析出不同被測(cè)樣品的微生物菌落組成結(jié)構(gòu)，同時(shí)可以分析一些不可培養(yǎng)且豐度低的微生物[4]?；诖?，本文采用高通量測(cè)序技術(shù)。目前，16S rDNA測(cè)序技術(shù)已成為研究微生物菌群結(jié)構(gòu)及其多樣性的重要方式[4-8]。張敏[33]等采用16S rDNA高通量測(cè)序方法比較了新疆西北部地區(qū)乳制品中微生物的多樣性，通過(guò)對(duì)新疆2個(gè)地區(qū)的7種乳制品進(jìn)行高通量測(cè)序分析發(fā)現(xiàn)，7種乳制品在門水平的優(yōu)勢(shì)菌門是厚壁菌門和變形菌門，而酸奶和原奶中的豐度不同，屬水平上的差異較大。同時(shí)還發(fā)現(xiàn)不同動(dòng)物來(lái)源的原奶和酸奶的微生物多樣性也存在著顯著的差異。姚宇秀[3]等采用16S rDNA高通量測(cè)序技術(shù)對(duì)原料乳及其采集過(guò)程中的相關(guān)設(shè)備以及奶牛的不同部位進(jìn)行樣本的分析，結(jié)果發(fā)現(xiàn)，輸入管路外側(cè)和盛裝原料奶的空貯罐中有金黃色葡萄球菌的存在。布仁其其格[9]等通過(guò)16S rRNA基因序列對(duì)酸馬奶傳統(tǒng)發(fā)酵過(guò)程中細(xì)菌群落結(jié)構(gòu)演替變化，結(jié)果表明，在發(fā)酵初期細(xì)菌多樣性最高，且在72 h時(shí)細(xì)菌豐度最高，門水平的優(yōu)勢(shì)菌門為厚壁菌門和變形菌門，屬水平的優(yōu)勢(shì)細(xì)菌屬為乳桿菌屬。

本研究的目的是利用高通量測(cè)序方法對(duì)重慶地區(qū)的不同液態(tài)奶的細(xì)菌菌群結(jié)構(gòu)和組成進(jìn)行分析比較，以期探索液態(tài)奶中的腐敗菌，實(shí)現(xiàn)對(duì)液態(tài)奶中有害菌的防治，從而實(shí)現(xiàn)提高液態(tài)奶的品質(zhì)。

2 材料與方法

2.1 材料與儀器

液態(tài)奶，重慶市天友乳液股份有限公司的4個(gè)牧場(chǎng)無(wú)菌生理鹽水，北京陸橋技術(shù)股份有限公司；植物基因組DNA提取試劑盒(TSP101-200)，TSINGKE；通用引物16S，TSINGKE；I-5TM2×High-Fidelity Master Mix(TP001)，TSINGKE；高純度低電滲瓊脂糖(TSJ001)，TSINGKE；DNA凝膠回收試劑盒(GE0101-200)，TSINGKE；DL2000 DNA Marker(TSJ011-500)，TSINGKE。

VITEK2 COMPACT全自動(dòng)細(xì)菌鑒定及藥敏分析系統(tǒng)，法國(guó)生物梅里埃公司； XH-B旋渦混合器，常州翔天實(shí)驗(yàn)儀器廠；BSC-1100A2-X生物安全柜，濟(jì)南鑫貝西生物技術(shù)有限公司；DK-98-IIA恒溫水浴鍋，北京市泰和潤(rùn)儀器有限公司；DW-86W50超低溫冰柜，浙江捷盛制冷科技有限公司；PD100A高壓滅菌鍋，致微(廈門)儀器有限公司； DHP-9272電熱恒溫培養(yǎng)箱，上海一恒科學(xué)儀器有限公司；高通量測(cè)序儀Illumina Hiseq 2500，美國(guó)因美納公司；離心機(jī)Legend Micro17，賽默飛世爾；電泳槽JYDF，北京君意東方電泳設(shè)備有限公司；PCR儀2720 thermal cycler，美國(guó)應(yīng)用生物系統(tǒng)公司；電泳儀JY300C，北京君意東方電泳設(shè)備有限公司；水浴鍋DFD-700，北京中興偉業(yè)；凝膠成像儀JY04S-3C，北京君意東方電泳設(shè)備有限公司；板式離心機(jī)L550，cence湘儀。

2.2 方法

2.2.1樣品采集

采集來(lái)自天友乳業(yè)重慶周邊4個(gè)不同牧場(chǎng)的原料乳(Sample 1、Sample 2、Sample 3、Sample 4)。所有的樣品經(jīng)無(wú)菌容器密封好后放入車載冰箱運(yùn)回實(shí)驗(yàn)室，樣品保存在-80℃超低溫冰箱，待試驗(yàn)。

2.2.2DNA的提取、擴(kuò)增和測(cè)序

取40 mL 原料乳樣本，6 000 r/min下4 ℃離心20 min，用滅菌脫脂棉去掉乳脂層，棄去上清液，加入無(wú)菌生理鹽水至10 mL，10 000 r/min高速冷凍離心10 min，用1 mL蔗糖緩沖液沖洗兩次，后懸浮400 μL蔗糖緩沖液及2U變?nèi)芫睾?00 μg溶菌酶，37 ℃培養(yǎng)1 h。選用806R(5’-GGACTACHVGGGTWTCTAAT-3’)和338F(5’-ACTCCTACGGGAGGCAGCAG-3’)對(duì)細(xì)菌16S rRNA V3-V4區(qū)進(jìn)行PCR擴(kuò)增。

PCR反應(yīng)體系為：12.5 μL I-5TM2×High-Fidelity Master Mix，1 μL 10 mmol/L Primer A，1 μL 10 mmol/L Primer B，加入模板DNA，補(bǔ)ddH2O至25 μL。用瓊脂糖凝膠電泳檢測(cè)，利用凝膠回收試劑盒切膠回收PCR產(chǎn)物。

PCR擴(kuò)增程序?yàn)椋?8 ℃預(yù)變性2 min；循環(huán)數(shù)(25～35次)×(95 ℃變性10 s，45～68 ℃退火10～15 s，72 ℃延伸5～15 s)；72 ℃終延伸1～5 min；4 ℃冷卻。

使用Illumina Hiseq 2500測(cè)序，由北京擎科生物科技有限公司重慶分公司完成。

2.3 數(shù)據(jù)處理

基于Illumina Hiseq 2500測(cè)序平臺(tái)，使用Flash 軟件對(duì)每個(gè)樣品的PE Reads進(jìn)行拼接，得到的序列即為原始序列，再通過(guò)Trimmomatic軟件對(duì)原始序列進(jìn)行優(yōu)化處理，得到較高質(zhì)量的優(yōu)化序列，最后通過(guò)UCHIME軟件鑒定并去除嵌合體序列，得到最終的有效序列。

MH7A細(xì)胞為永生化的RA關(guān)節(jié)滑膜細(xì)胞，購(gòu)自廣州吉妮歐公司。用高糖DMEM+10%胎牛血清于5%CO2，37℃的培養(yǎng)箱培養(yǎng)。購(gòu)買時(shí)為第3代，培養(yǎng)3代后用于后續(xù)實(shí)驗(yàn)研究。

3 結(jié)果與分析

3.1 測(cè)序數(shù)據(jù)統(tǒng)計(jì)及結(jié)果分析

為了獲得準(zhǔn)確、高質(zhì)量的生物信息分析結(jié)果，通過(guò)測(cè)序獲得的原始數(shù)據(jù)需要進(jìn)行優(yōu)化處理[10]。如表1所示。

表1 樣品測(cè)序數(shù)據(jù)統(tǒng)計(jì)Table 1 Sample sequencing data statistics

從表1中可以看出，4個(gè)樣品總計(jì)測(cè)得原始序列條數(shù)為242 576條，將低質(zhì)量的序列過(guò)濾后得到的優(yōu)化序列數(shù)總計(jì)233 818條，再將上述序列進(jìn)行冗余處理，最終得到的有效序列數(shù)為212 260條。

圖1為有效序列的長(zhǎng)度分布圖，從圖1中可以看出，大部分的序列均分布在400～450 bp之間，且430～440 bp之間的序列最多。從序列長(zhǎng)度的分布可以看出，與16S rDNAV3-V4區(qū)序列長(zhǎng)度大致吻合[11]。

圖1 有效序列長(zhǎng)度分布Fig.1 Distribution of effective sequence length

從樣本中隨機(jī)抽取一定數(shù)量的序列后，對(duì)獲得的序列所代表的物種數(shù)目進(jìn)行統(tǒng)計(jì)，以物種數(shù)目和序列數(shù)構(gòu)建的一條曲線為樣本的稀釋曲線[12]。稀釋曲線可以驗(yàn)證測(cè)序數(shù)據(jù)量是否可以反映樣品中的物種多樣性，并間接反映樣品的測(cè)序深度和物種的豐富度[12-13]。若在一定范圍內(nèi)，隨著序列數(shù)的增加，曲線表現(xiàn)出急劇上升的趨勢(shì)則反映樣本中仍有大量物種被發(fā)現(xiàn)；當(dāng)稀釋曲線趨于平緩則表示樣本量足夠充分[14]。圖2(a)是樣品中稀釋曲線圖，從該圖中可以看出，隨著樣本序列數(shù)的增加，每個(gè)樣品中OTU數(shù)目增加，當(dāng)OTU數(shù)目到達(dá)一定數(shù)量后，曲線趨于平緩，這表明實(shí)驗(yàn)中所得序列數(shù)基本可以反映出每組樣品中的菌群結(jié)構(gòu)。同時(shí)，可以看出Sample 4中OTU數(shù)目最多，Sample 1次之，Sample 2和Sample 3中OTU數(shù)目相對(duì)較少。

香農(nóng)指數(shù)曲線可以體現(xiàn)出各樣品在不同測(cè)序數(shù)量時(shí)樣本中微生物多樣性的指數(shù)，隨著樣本序列數(shù)的增加，曲線呈上升趨勢(shì)[30]；當(dāng)曲線趨于平坦時(shí)，說(shuō)明測(cè)序數(shù)據(jù)量足夠大，OTU種類不會(huì)在隨著序列量的增長(zhǎng)而增長(zhǎng)[31]。圖2(b)是樣品中香農(nóng)指數(shù)曲線，從該圖中可以看出，隨著測(cè)序樣本數(shù)量的增加曲線趨于平緩，這說(shuō)明本研究的測(cè)序深度是可以滿足分析要求[15]。

(a) 樣品稀釋曲線

3.2 樣品OTU統(tǒng)計(jì)分析

一般情況下，將相似性高于97%的不同16S rDNA序列定義為一個(gè)OTU，每一個(gè)OTU通常被視為一個(gè)微生物物種[32]。圖3展示了不同樣品的OTU數(shù)。從圖3中可以看出，各樣品的OTU數(shù)相差不大，均在70～85個(gè)之間，Sample 1～Sample 4各樣品中OTU個(gè)數(shù)分別為80，73，79，81；四個(gè)樣本中總共的OTU數(shù)目為119個(gè)。

圖3 樣品中OTU數(shù)目的分布圖Fig.3 Distribution of the number of OTUs in different samples

Venn圖體現(xiàn)了不同樣品之間微生物群落的相似性和差異性[4]，如圖4所示。從圖4中可以看出，4個(gè)樣品中共有OTU數(shù)目為43個(gè)，每個(gè)樣品特有的OTU個(gè)數(shù)分別為：12，5，7，10，分別占各樣本中OTU數(shù)量的15%、6.85%、8.86%和12.35%；Sample 1與Sample 2、Sample 1與Sample 3、Sample 2與Sample 4、Sample 3與Sample 4樣品之間共有，且其余樣品不具有的OTU數(shù)目較少，分別為2,3,2,3。因此，可以看出各樣品之間微生物群落沒(méi)有顯著的差異性。

圖4 樣品中OTU數(shù)目的Venn圖Fig.4 Venn diagram of the number of OTUs in different samples

3.3 樣品間α多樣性的分析

α多樣性分析體現(xiàn)了單個(gè)樣品內(nèi)部物種的豐富度和多樣性[16]，通常可以通過(guò)ACE指數(shù)、Chao1指數(shù)、Simpson指數(shù)以及Shannon指數(shù)等反映[17,18]。其中，ACE指數(shù)和Chao1指數(shù)用來(lái)反映微生物群落的豐富度[18]，ACE指數(shù)和Chao1指數(shù)越大表示樣品中微生物群落的豐富度越大[32]；Shannon指數(shù)和Simpson指數(shù)都是用來(lái)反映微生物群落的多樣性[30]，微生物群落的多樣性與Shannon指數(shù)呈正相關(guān)，與Simpson指數(shù)呈負(fù)相關(guān)[31]。

對(duì)不同樣本在97%一致性閾值下的α多樣性指數(shù)進(jìn)行統(tǒng)計(jì)和分析[4,19-20]，結(jié)果見(jiàn)表2。從表2中可以看出，所有樣品的覆蓋率(Coverage)均大于0.99，說(shuō)明本次研究所建立的數(shù)據(jù)庫(kù)可信度較高，數(shù)據(jù)的真實(shí)性較好，能夠有效地反映樣品中微生物菌群的多樣性。Sample 2樣品中ACE指數(shù)和Chao1指數(shù)最高，表明該樣品中菌群的豐富度較高；而Sample 3中ACE指數(shù)和Chao1指數(shù)最低，說(shuō)明該樣品中菌群的豐富度較低。同時(shí)，Sample 4樣品中其Simpson指數(shù)最小，Shannon指數(shù)最大說(shuō)明該樣品中菌群的多樣性最大；Sample 1次之，Sample 2樣品中菌群的多樣性最低。

表2 樣品間α多樣性指數(shù)統(tǒng)計(jì)Table 2 Statistics of α diversity index among samples

3.4 液態(tài)奶樣品間基于門水平的菌群結(jié)構(gòu)分析

圖5是液態(tài)奶各樣品中基于門水平的菌群結(jié)構(gòu)分布比例圖。從圖中可以看出，液態(tài)奶樣品中的10個(gè)細(xì)菌門被鑒定出，分別是異常球菌-棲熱菌門(Deinococcus-Thermus)、厚壁菌門(Firmicutes)、硝化螺旋菌門(Nitrospirae)、變形菌門(Proteobacteria)、藍(lán)藻細(xì)菌門(Cyanobacteria)、擬桿菌門(Bacteroidetes)、放線細(xì)菌門(Actinobacteria)、Patescibacteria、綠彎菌門(Chloroflexi)和Epsilonbacteraeota。從圖5中可以看出各樣品之間的優(yōu)勢(shì)菌門基本相同，但相對(duì)豐度差異較大。Sample 1中的主要菌群是厚壁菌門、變形菌門和擬桿菌門分別占38.9%、46.9%和14.0%；Sample 2中的主要菌群是厚壁菌門、變形菌門和放線細(xì)菌門分別占73.8%、10.0%和13.1%；Sample 3中的主要菌群是厚壁菌門、變形菌門和擬桿菌門分別占20.9%、28.0%和50.1%；Sample 4中的主要菌群是放線細(xì)菌門、厚壁菌門和異常球菌-棲熱菌門分別占30.0%、29.6%和22.8%。

圖5 基于門水平各樣品的菌群結(jié)構(gòu)分析Fig.5 Analysis of flora structure of each sample based on phylum level

由此可見(jiàn)，每組液態(tài)奶中門水平的菌群結(jié)構(gòu)差異性不大，厚壁菌門和變形菌門是四組液態(tài)奶的共有菌門，這與張敏等[33]的研究中乳制品的菌群結(jié)論相同。

物種豐度聚類熱圖分析是通過(guò)顏色變化與相似程度來(lái)反應(yīng)二維矩陣或表格中的數(shù)據(jù)信息，并呈現(xiàn)出群落物種的組成信息，表明了液態(tài)奶不同樣品間不同細(xì)菌門的相對(duì)豐度及細(xì)菌組成的差異性和樣品間的相似性[10]。采用物種豐度聚類熱圖分析液態(tài)奶樣品中含量前10個(gè)菌門和4個(gè)樣品之間的交互關(guān)系，見(jiàn)圖6。

圖6 基于門水平各樣品的物種豐度聚類熱圖Fig.6 Species abundance clustering heatmap for each sample at phylum level

從圖6中可知，不同樣品間門水平具有一定的豐度，且各樣品的菌群組成具有一定的差異性和相似性。同時(shí)，根據(jù)豐度聚類熱圖可以看出樣品間門水平的聚類關(guān)系，即4組樣品均聚在一起，說(shuō)明四組樣品間的豐度相似，其中Sample 2和Sample 4的相似度較高些，Sample 1與其余三組在門水平上的相似度較低。

3.5 液態(tài)奶樣品間基于屬水平的菌群結(jié)構(gòu)分析

圖7是液態(tài)奶各樣品中基于屬水平的聚類樹(shù)柱狀圖。從圖中可以看出，液態(tài)奶中的10個(gè)菌屬信息被鑒定出，分別是Olivibacter、克雷伯氏菌屬(Klebsiella)、微小桿菌屬(Exiguobacterium)、不動(dòng)桿菌屬(Acinetobacter)、異常球菌屬(Deinococcus)以及短波單胞菌屬(Brevunmdimonas)，寡養(yǎng)單胞菌屬(Stenotrophnomonas)、微桿菌屬(Microbacterium)、Cloacibacterium、芽孢桿菌屬(Bacillus)和其他未分類菌。

圖7 基于屬水平各樣品的聚類樹(shù)柱狀圖Fig.7 Cluster tree histogram of each sample based on genus level

各樣品在屬水平上優(yōu)勢(shì)菌屬分別是：Sample 1中的優(yōu)勢(shì)菌屬是厚壁菌門的芽孢桿菌屬(占35.9%)及變形菌門的短波單胞菌屬(占26.2%)；Sample 2中的優(yōu)勢(shì)菌屬是厚壁菌門的芽孢桿菌屬(占73.7%)；Sample 3中的優(yōu)勢(shì)均屬是擬桿菌門的Cloacibacterium(占49.9%)、變形菌門的寡養(yǎng)單胞菌屬(占27.6%)；Sample 4中的優(yōu)勢(shì)菌屬是厚壁菌門的芽孢桿菌屬(占29.2%)和微桿菌屬(占28.3%)、異常球菌-棲熱菌門的異常球菌屬(占22.8%)。從圖7中可以看出，Sample 1和Sample 4的細(xì)菌群落均勻性相對(duì)較高，Sample 2的細(xì)菌群落組成高度集中在芽孢桿菌屬中。

根據(jù)已有的研究報(bào)道表明，不動(dòng)桿菌屬、芽孢桿菌屬及克雷伯氏菌屬屬于嗜冷菌[21-24]。嗜冷菌在低溫存儲(chǔ)過(guò)程中，其依舊能夠繁殖并能夠產(chǎn)生影響原料乳品質(zhì)的腐敗酶：蛋白酶和脂肪酶等，破壞原料乳中的營(yíng)養(yǎng)物質(zhì)，影響原料乳的品質(zhì)[21,25-28]。研究表明：生鮮乳本身不存在大量的嗜冷菌，其主要來(lái)源于生產(chǎn)環(huán)節(jié)的污染，導(dǎo)致原料乳在運(yùn)送過(guò)程中出現(xiàn)大量嗜冷菌繁殖從而導(dǎo)致原料乳的品質(zhì)遭到破壞[22,29]。

4 結(jié)論

本實(shí)驗(yàn)通過(guò)16S rRNA基因的高通量測(cè)序技術(shù)對(duì)來(lái)自4個(gè)牧場(chǎng)液態(tài)奶的細(xì)菌群落結(jié)構(gòu)進(jìn)行研究。通過(guò)對(duì)比發(fā)現(xiàn)，在門水平上，4組液態(tài)奶的菌群結(jié)構(gòu)基本相似，但相對(duì)豐度差異較大，4組液態(tài)奶的優(yōu)勢(shì)菌門主要是厚壁菌門和變形菌門。在屬水平上，4組液態(tài)奶的優(yōu)勢(shì)菌屬和相對(duì)豐度都不相同，存在顯著差異。同時(shí)，還發(fā)現(xiàn)屬水平上芽孢桿菌屬、不動(dòng)桿菌屬屬于嗜冷菌，在液態(tài)奶的貯存中容易分解其中的蛋白質(zhì)、脂肪或碳水化合物，從而導(dǎo)致液態(tài)奶出現(xiàn)腐敗變質(zhì)等現(xiàn)象。這為后續(xù)對(duì)探索液態(tài)奶中腐敗菌的來(lái)源及防控提供了有力的研究依據(jù)。