高通量測序分析比較不同保存溫度下牦牛酸奶細菌多樣性

麻和平，張文齊，劉彩云，王 潔，彭章普，邵建寧*

（1.甘肅省科學(xué)院生物研究所甘肅省微生物資源開發(fā)利用重點實驗室，甘肅蘭州 730000；2.甘肅省科學(xué)院，甘肅蘭州 730000）

傳統(tǒng)發(fā)酵乳制品是經(jīng)過長期自然馴化優(yōu)良食品微生物的天然載體，發(fā)酵乳制品中的微生物決定了發(fā)酵乳制品的品質(zhì)[1-2]。近年來，有關(guān)傳統(tǒng)發(fā)酵乳制品中微生物多樣性的研究，既有傳統(tǒng)的分離培養(yǎng)技術(shù)研究，也有高通量測序、宏基因組分析等非培養(yǎng)技術(shù)研究[3-5]。利用傳統(tǒng)的分離技術(shù)雖然從樣品中獲得了一些微生物，但只能代表具體培養(yǎng)條件所能分離到的微生物，無法分離到樣品中所有微生物，不能完整地揭示樣品的菌群組成[6-9]。高通量測序[10-13]、宏基因組分析[14-16]等非培養(yǎng)技術(shù)研究無需分離培養(yǎng)微生物，能夠更好地揭示樣品中菌群組成及多樣性[17-18]。

目前對傳統(tǒng)發(fā)酵乳制品中參與發(fā)酵微生物的分離和菌群多樣性相關(guān)研究較多，孫思雨等[19]采用Illumina Miseq高通量測序?qū)Υㄎ鞲咴煌貐^(qū)的傳統(tǒng)自然發(fā)酵牦牛酸奶樣品進行細菌群落結(jié)構(gòu)分析；孫志宏等[20]采用454高通量測序技術(shù)從宏基因組角度對西藏地區(qū)自然發(fā)酵牦牛乳中微生物多樣性進行探索性研究；SUN Z等[21]對采集自中國內(nèi)蒙古、甘肅、四川以及蒙古國的17個自然發(fā)酵乳樣品進行多樣性分析；XU H Y等[22]對采集自中國新疆的22個自然發(fā)酵乳樣品進行多樣性研究；LIU W等[23]對采集自中國西藏的16個自然發(fā)酵牦牛乳樣品進行分析；丁武蓉[24]采用傳統(tǒng)培養(yǎng)法結(jié)合高通量測序?qū)η嗖馗咴瓊鹘y(tǒng)發(fā)酵牦牛酸奶中乳酸菌多樣性進行研究，然而，通過高通量測序分析比較不同保存溫度條件下傳統(tǒng)發(fā)酵乳制品微生物多樣性和優(yōu)化傳統(tǒng)發(fā)酵乳制品保存溫度卻鮮有報道。

甘南牧區(qū)位于甘肅省甘南藏族自治州境內(nèi)，地處青藏高原東北邊緣與黃土高原接壤帶，地理坐標位于東經(jīng)100°46′～104°44′，北緯33°06′～36°10′之間，地勢西北高東南低，地形狀況復(fù)雜，牧區(qū)平均海拔3 000 m以上，氣候主要以高原大陸性季風(fēng)氣候為主，地理季節(jié)分布有顯著差異，是我國重要的畜牧業(yè)生產(chǎn)基地。獨特的生態(tài)環(huán)境和自然長期馴化使得甘南牧區(qū)牧民沿用傳統(tǒng)方法制作的牦牛酸奶含有豐富且特有的食品微生物資源。本研究以不同保存條件下（常溫、4 ℃、-20 ℃和-80 ℃）甘南牧區(qū)牧民自制牦牛酸奶為研究對象，采用Illumina NovaSeq高通量測序技術(shù)，對牦牛酸奶細菌的16S rRNA V3-V4區(qū)基因序列進行測序，通過操作分類單元（operational taxonomic unit，OTU）聚類、α多樣性和β多樣性等分析比較不同保存溫度下牦牛酸奶的細菌群落結(jié)構(gòu)和多樣性，探究不同保存溫度對傳統(tǒng)發(fā)酵乳制品微生物菌群多樣性的影響，為傳統(tǒng)發(fā)酵乳制品采樣條件優(yōu)化提供參考依據(jù)。

1 材料與方法

1.1 材料與試劑

1.1.1 樣品采集

牦牛酸奶：2019年11月于甘肅省甘南藏族自治州卓尼縣（Z）、碌曲縣（L）實地采集牧民自制的牦牛酸奶，每份樣品采集4管，每管30 mL，密封后標記樣品號，分別于常溫、冰袋、干冰和液氮保存條件下運回實驗室后，分別在常溫、4 ℃、-20 ℃和-80 ℃溫度條件下保存?zhèn)溆?，樣品的具體信息見表1。

表1 甘南牧區(qū)牦牛酸奶采樣信息Table 1 Information of yak yogurt samples collected from Gannan pastoral area

續(xù)表

1.1.2 主要試劑

QIAampFastDNAStoolMiniKit、QIAquickGelExtraction Kit：德國QIAGEN公司；PowerSoil DNA Isolation Kit：美國MOBIO公司；KAPA HiFi Hotstart Ready Mix Kit：美國KAPA Biosystems公司；2×Phanta Master Mix：南京諾唯贊生物科技股份有限公司；AxyPrep DNA凝膠回收試劑盒：美國AXYGEN公司。

1.2 儀器與設(shè)備

5424R低溫高速離心機：德國Eppendorf公司；VDRTEX-5渦旋混合器：海門市其林貝爾儀器制造有限公司；Tanon 2500R凝膠成像系統(tǒng)：上海天能科技有限公司；DYY-6C電泳儀：上海天能科技有限公司；ND-2000C微量紫外分光光度計：美國賽默飛世爾科技公司；MiniAmp聚合酶鏈式反應(yīng)（polymerase chain reaction，PCR）儀：美國Life Technologies公司。

1.3 方法

1.3.1 牦牛酸奶菌群總DNA的提取

取2 mL牦牛酸奶樣品，采用試劑盒QIAamp Fast DNA Stool Mini Kit，按照說明書的要求提取牦牛酸奶菌群的總DNA，利用微量紫外分光光度計和1%瓊脂糖凝膠電泳進行總DNA的完整性和濃度質(zhì)檢。

1.3.2 PCR擴增與測序

以提取的總DNA為模板，341F（5'-CCTACGGGRSGCAGCAG-3'）和806R（5'-GGACTACVVGGGTATCTAATC-3'）為引物（引物中已帶特定barcode），使用KAPA HiFi Hotstart Ready Mix PCR kit對細菌的16S rDNA V3-V4區(qū)基因序列進行PCR擴增，PCR擴增條件：95 ℃預(yù)變性3 min；98 ℃變性20 s，58 ℃退火15 s，72 ℃延伸20 s，共30個循環(huán)；72 ℃再延伸5 min。PCR擴增體系：2×KAPA Library Amplification Ready Mix 15 μL，引物（10 μmol/L）各1 μL，模板DNA 50 ng，加雙蒸水（ddH2O）補充至30 μL。PCR擴增產(chǎn)物用2%瓊脂糖凝膠電泳檢測，使用AxyPrep DNA凝膠提取試劑盒切膠回收PCR擴增產(chǎn)物，文庫質(zhì)檢合格后，使用Qubit 2.0進行文庫定量，并根據(jù)每個樣品的數(shù)據(jù)量要求，進行相應(yīng)比例的混合，采用Illumina NovaSeq平臺PE250策略進行測序，Illumina NovaSeq高通量測序由上海銳翌生物科技有限公司完成。

1.3.3 數(shù)據(jù)分析

16S rDNA V3-V4區(qū)基因序列被控制在220～500 bp之間，保證每條reads平均質(zhì)量值不低于20，且含N堿基數(shù)不超過3個。將序列完全一樣的Clean Reads根據(jù)其豐度大小進行排序，過濾掉其中的Singletons。使用UPARSE（http：//drive5.com/uparse/）根據(jù)97%的相似度進行OTU聚類，并使用Userach（版本7.0.1090）鑒定和移除嵌合體序列。每個OTU都有一個代表性的序列，將該代表性序列與已知物種的RDP數(shù)據(jù)庫（http：//rdp.cme.msu.edu/）進行比對，置信度閾值設(shè)置為0.8，利用RDP Classifer（http：//rdp.cme.msu.edu/）對每個代表性序列進行物種注釋。不同樣本對應(yīng)的Reads數(shù)量差距較大，為避免因樣品數(shù)據(jù)大小不同而造成分析時的偏差，在樣品達到足夠的測序深度的情況下，對每個樣品進行隨機抽平處理。測序深度用alpha多樣性指數(shù)來衡量，OTU profiling table和alpha多樣性指數(shù)通過QIIME（version 1.9.1）的python腳本實現(xiàn)。

2 結(jié)果與分析

2.1 不同保存溫度下牦牛酸奶樣品的細菌菌群多樣性分析

通過對牦牛酸奶樣品細菌的16S rDNA V3-V4區(qū)基因序列進行測序，20份牦牛酸奶樣品質(zhì)控過濾后得到723 306條有效序列。采用稀釋曲線評估各樣品測序量是否足夠，結(jié)果見圖1。由圖1可知，隨著序列數(shù)的增加，超1（Chao1）指數(shù)趨于平緩，說明本次樣本測序的深度足夠，按最小樣本序列數(shù)24 822抽平足以反映各樣品中細菌物種多樣性，測序量滿足后續(xù)生物信息學(xué)分析的要求。各個樣品及分組信息、質(zhì)控后序列數(shù)、Observed species、Chao1指數(shù)、香農(nóng)（Shannon）指數(shù)、辛普森（Simpson）指數(shù)和測序深度見表2。由表2可知，所有樣本的測序深度均＞0.99，說明樣本測序結(jié)果可以反映樣品的真實情況。所有序列按97%的相似度進行OTU聚類，從20份牦牛酸奶樣品中共檢測到1 289個OTU，隸屬于26門、46綱、61目、127科、248屬。對不同保存溫度樣品組間α多樣性指數(shù)進行kw檢驗分析，不同保存溫度下樣品組間α多樣性指數(shù)Observed species、Chao1指數(shù)、香農(nóng)（Shannon）指數(shù)、辛普森（Simpson）指數(shù)和測序深度的P值分別為0.56、0.59、0.94、0.88和0.61，說明不同保存溫度樣品組間多樣性無顯著差異（P＞0.05）。

圖1 不存保存溫度下牦牛酸奶樣品的稀釋曲線Fig.1 Rarefaction curve of yak yoghurt samples at different storage temperature

表2 牦牛酸奶樣品細菌的序列信息及α多樣性指數(shù)Table 2 Sequence information and α-diversity indexes of bacteria of yak yoghurt samples

2.2 不同保存溫度下牦牛酸奶樣品細菌群落結(jié)構(gòu)的比較

2.2.1 組間細菌群落結(jié)構(gòu)比較

主坐標分析（principal coordinate analysis，PCoA）可以直觀地將樣品間、樣品組間差異程度體現(xiàn)在二維坐標圖上，樣品間距離越近則表示其物種組成越相似[25]，對不同保存溫度下牦牛酸奶樣品的細菌群落進行PCoA，并對UniFrac距離進行Adonis分析，結(jié)果見圖2。

圖2 各組牦牛酸奶樣品的細菌群落主坐標分析Fig.2 Principal coordinate analysis of bacterial community of yak yoghurt samples of each group

由圖2可知，主坐標1和主坐標2分別解釋了不同保存溫度下牦牛酸奶樣品中細菌群落73.16%和12.27%的信息，兩主坐標貢獻率之和＞85.43%，表明兩個主坐標較好地代表了樣品中的細菌群落信息。圖2中4組樣品間距離較近，表明4組樣品的物種組成較為接近，Adonis分析表明4種不同保存溫度牦牛酸奶樣品細菌菌群多樣性沒有顯著差異（P＞0.05）。

2.2.2 不同保存溫度下牦牛酸奶樣品細菌門水平物種豐度分析

在細菌門水平，所有牦牛酸奶樣品中共鑒定出26個細菌門，將相對豐度排序靠后的6個細菌門并為other，結(jié)果見圖3。從常溫、4 ℃、-20 ℃和-80 ℃保存的牦牛酸奶樣品中分別共鑒定出13個、12個、20個、13個細菌門，其中厚壁菌門（Firmicutes）、變形菌門（Proteobacteria）和擬桿菌門（Bacteroidetes）為所有牦牛酸奶樣品中的優(yōu)勢菌門[26]，其相對豐度分別＞75.87%、7.20%和3.52%，未鑒定到門水平的細菌相對豐度均＜0.1%。通過門水平細菌群落結(jié)構(gòu)分析比較，不同保存溫度下耗牛酸奶樣品組間細菌群落結(jié)構(gòu)組成基本相同，僅組成比例存在一定差異。

圖3 門水平不同保存溫度下牦牛酸奶樣品細菌群落結(jié)構(gòu)分析結(jié)果Fig.3 Analysis results of bacterial community structure of yak yogurt samples at different storage temperature at phylum level

2.2.3 不同保存溫度下牦牛酸奶樣品細菌屬水平物種豐度分析

在細菌屬水平，所有耗牛酸奶樣品中共鑒定出248個細菌屬，對于平均豐度值大于1%的細菌屬，從常溫保存的牦牛酸奶中鑒定出乳桿菌屬（Lactobacillus）、鏈球菌屬（Streptococcus）、醋桿菌屬（Acetobacter）、金黃桿菌屬（Chryseobacterium）、假單胞菌屬（Pseudomonas）、叢毛單胞菌屬（Comamonas），相對豐度分別為75.77%、11.88%、2.61%、2.50%、1.27%、1.18%；從4 ℃保存的牦牛酸奶中鑒定出乳桿菌屬（Lactobacillus）、鏈球菌屬（Streptococcus）、假單胞菌屬（Pseudomonas）、金黃桿菌屬（Chryseobacterium），相對豐度分別為75.40%、13.47%、3.02%、2.82%；從-20 ℃保存的牦牛酸奶中鑒定出乳桿菌屬（Lactobacillus）、鏈球菌屬（Streptococcus）、金黃桿菌屬（Chryseobacterium）、假單胞菌屬（Pseu-domonas）、叢毛單胞菌屬（Comamonas），相對豐度分別為61.82%、14.97%、5.27%、5.19%、1.25%；從-80 ℃保存的牦牛酸奶中鑒定出乳桿菌屬（Lactobacillus）、鏈球菌屬（Streptococcus）、假單胞菌屬（Pseudomonas）、金黃桿菌屬（Chryseobacterium）、叢毛單胞菌屬（Comamonas）、氣單胞菌屬（Aeromonas），相對豐度分別為68.86%、14.59%、5.2%、3.91%、1.27%、1.00%。圖4顯示牦牛酸奶樣品在屬水平相對豐度前20的物種屬，其他合并為other。通過屬水平細菌群落結(jié)構(gòu)分析比較，不同保存溫度下耗牛酸奶樣品組間細菌群落結(jié)構(gòu)組成相似，僅組成比例存在一定差異，乳桿菌屬和鏈球菌屬為優(yōu)勢菌屬[26]，平均豐度值分別＞61.82%和11.88%。

圖4 屬水平不同保存溫度下牦牛酸奶樣品細菌群落結(jié)構(gòu)分析結(jié)果Fig.4 Analysis results of bacterial community structure of yak yogurt samples at different storage temperature at genus level

2.2.4 不同保存溫度下牦牛酸奶樣品組間物種差異分析

采用秩和檢驗得到不同保存溫度下牦牛酸奶樣品組間的差異細菌屬共2個（P＜0.05），分別為Enhydrobacter、Peptoniphilus。為了直觀展示各組及各樣品間細菌屬差異，對差異細菌屬進行熱圖（Heatmap）分析，結(jié)果見圖5。

圖5 不同保存溫度下牦牛酸奶樣品間差異細菌屬的Heatmap分析結(jié)果Fig.5 Heatmap analysis results of differential bacterial genera between yak yoghurt samples at different storage temperatures

由圖5可知，在各牦牛酸奶樣品中Enhydrobacter的相對豐度高于Peptoniphilus，且在Z-20組中Enhydrobacter相對豐度較高，即-20 ℃保存條件下的耗牛酸奶樣品中Enhydrobacter較高。為了直觀展示不同保存溫度下耗牛酸奶樣品分組之間細菌屬差異，對各組牦牛酸奶樣品中的差異細菌屬進行主成分分析（principal component analysis，PCA），結(jié)果見圖6。

由圖6可知，不同保存溫度樣品組Z-4、Z-20、Z-80組的距離接近，細菌群落組成較為接近，Z-0組的距離較遠，細菌群落組成差異較大，說明4 ℃、-20 ℃、-80 ℃保存組牦牛酸奶樣品的物種組成更為接近，不同保存溫度下耗牛酸奶樣品組間細菌多樣性無顯著差異。

3 結(jié)論

本研究采用Illumina NovaSeq高通量測序技術(shù)對4種不同保存溫度條件下甘南牧區(qū)牦牛酸奶細菌群落多樣性及結(jié)構(gòu)進行分析，結(jié)果表明，不同保存溫度下的牦牛酸奶細菌群落多樣性無顯著差異，細菌群落組成基本相同，僅組成比例存在一定差異，優(yōu)勢細菌門為厚壁菌門（Firmicutes）、變形菌門（Proteobacteria）和擬桿菌門（Bacteroidetes），其平均相對豐度分別＞75.87%、7.20%和3.52%，優(yōu)勢細菌屬為乳桿菌屬（Lactobacillus）和鏈球菌屬（Streptococcus），其平均相對豐度分別＞61.82%和11.88%。