甘肅甘南傳統(tǒng)發(fā)酵牦牛乳中細(xì)菌菌群的多樣性研究

馬彩霞 梁 琪,* 王湘竹 劉 瑛

(1甘肅農(nóng)業(yè)大學(xué)食品科學(xué)與工程學(xué)院/甘肅省功能乳品工程實(shí)驗(yàn)室，甘肅 蘭州 730070； 2甘肅臨夏州燎原乳業(yè)有限公司，甘肅 臨夏 731100)

甘肅甘南位于甘肅省西南部，地處甘、青、川三省交界地帶，是甘肅藏區(qū)的重要代表，具有獨(dú)特的地理環(huán)境和生物多樣性[1-2]。甘肅甘南總面積約4.5萬(wàn)平方米，境內(nèi)海拔1 100～4 900 m，紫外線(xiàn)強(qiáng)度高(單日紫外線(xiàn)指數(shù)可達(dá)到12.5)，年平均氣溫約1.7℃，年降雨量約400～800 mm，屬于典型的高原大陸季風(fēng)氣候[3-4]。畜牧業(yè)是甘肅甘南的支柱產(chǎn)業(yè)，現(xiàn)有草場(chǎng)4 000 余萬(wàn)畝，可利用面積約為3 800萬(wàn)畝[5]。牦牛屬于牛亞科動(dòng)物，對(duì)于高寒環(huán)境有著極強(qiáng)的耐受性，常生活在高海拔的青藏高原地區(qū)[6]。甘南牦牛生活在海拔2 800米以上的地區(qū)，主要以天然牧場(chǎng)牧草為食[7]。2015年統(tǒng)計(jì)數(shù)據(jù)顯示，甘肅甘南牦牛總數(shù)約138萬(wàn)頭，占全省牦?？偭康?0%[5]，為甘南牧民提供了豐富的牦牛乳資源。藏區(qū)牧民沿襲傳統(tǒng)的發(fā)酵方法制作了多種牦牛乳制品，其中最常食用的是傳統(tǒng)發(fā)酵牦牛乳[8]。甘肅甘南地區(qū)牧民制作的發(fā)酵牦牛乳歷經(jīng)數(shù)千年傳承，在獨(dú)特的自然環(huán)境下經(jīng)過(guò)長(zhǎng)期自然馴化使得其中蘊(yùn)含豐富的微生物群落，因此將其作為研究對(duì)象極具代表性[9]。

近年來(lái)，高通量測(cè)序技術(shù)因具有成本低、測(cè)序質(zhì)量高、能更客觀(guān)全面地了解菌群結(jié)構(gòu)等優(yōu)勢(shì)，可用于更好地評(píng)估微生物多樣性，目前已廣泛應(yīng)用于發(fā)酵乳制品中微生物菌群結(jié)構(gòu)及多樣性的研究[10]。Shangpliang等[11]采用Illumina MiSeq高通量測(cè)序技術(shù)對(duì)印度阿魯納恰爾邦和錫金州的四大類(lèi)共54份傳統(tǒng)發(fā)酵乳制品中的細(xì)菌群落結(jié)構(gòu)進(jìn)行分析發(fā)現(xiàn)，不同種類(lèi)的發(fā)酵乳制品細(xì)菌菌群多樣性存在差異。Liu等[12]采用焦磷酸測(cè)序研究了俄羅斯卡爾梅克和赤塔州的蒙古族家庭的19種自然發(fā)酵牛乳(naturally fermented cow milk, NFCM)中細(xì)菌和真菌的群落多樣性，結(jié)果表明卡爾梅克和赤塔州的樣品存在菌群結(jié)構(gòu)差異，地理環(huán)境的差異可能是影響兩地樣品中微生物多樣性的重要因素。Sessou等[13]采用高通量測(cè)序技術(shù)對(duì)采自貝寧的手工自制奶酪樣品，以及分別采自貝寧和尼日爾的20份發(fā)酵乳樣品中酵母菌的多樣性進(jìn)行了研究，發(fā)現(xiàn)兩國(guó)樣品間的酵母組成存在顯著差異。目前鮮有關(guān)于甘肅甘南地區(qū)傳統(tǒng)發(fā)酵牦牛乳中細(xì)菌菌群多樣性的研究報(bào)道，因此，全面了解該地區(qū)傳統(tǒng)發(fā)酵牦牛乳中細(xì)菌菌群的多樣性十分必要。本研究基于Illumina MiSeq高通量測(cè)序技術(shù)對(duì)采自甘肅甘南的傳統(tǒng)發(fā)酵牦牛乳中細(xì)菌菌群的多樣性進(jìn)行表征，以期全面解析發(fā)酵牦牛乳中的細(xì)菌菌群結(jié)構(gòu)，系統(tǒng)解析牦牛乳中的細(xì)菌組成及群落結(jié)構(gòu)，為后期優(yōu)良微生物資源的研究與開(kāi)發(fā)奠定基礎(chǔ)。

1 材料與方法

1.1 材料與試劑

樣品于2019年7月采自甘肅省甘南藏族自治州合作市，勒秀鄉(xiāng)麻拉行政村(海拔3 163.5 m)共采集4份(編號(hào)為M1、M2、M3、M4)、那吾鄉(xiāng)加拉行政村(海拔2 889.4 m)共采集4份(編號(hào)為W1、W2、W3、W4)、那吾鄉(xiāng)多河爾行政村(海拔3 000.5 m)共采集4份(編號(hào)為D1、D2、D3、D4)。樣品Q(chēng)L1于2019年8月采自青海省海北藏族自治州祁連縣(海拔3 015.4 m)，所有樣品均由牧民使用傳統(tǒng)方法發(fā)酵制成。

TruSeqTMDNA試劑盒，美國(guó)Illumina公司；AP-MN-P-50凝膠提取試劑盒，美國(guó)Axygen Biosciences公司； MetaVxTM文庫(kù)構(gòu)建試劑盒，美國(guó)GENEWIZ公司。

1.2 儀器與設(shè)備

Qubit 3.0熒光計(jì)，德國(guó)Life公司；Illumina MiSeq測(cè)序儀，英國(guó)Illumina公司；7900HT型熒光定量 PCR儀，美國(guó)Applied Biosystems公司；Nanodrop 2000型超微量分光光度計(jì)，美國(guó)Thermo fisher公司；Agilent 2100生物分析儀，美國(guó)AgilenTTechnologies公司；5424R型高速臺(tái)式冷凍離心機(jī)，德國(guó)Eppendorf公司。

1.3 試驗(yàn)方法

1.3.1 樣品總DNA提取 使用DNA提取試劑盒依照說(shuō)明書(shū)從13份樣品中提取 DNA，用 Qubit?dsDNA HS Assay Kit 檢測(cè)所提DNA濃度。

1.3.2 PCR擴(kuò)增、文庫(kù)構(gòu)建及上機(jī)測(cè)序 以20～30 ng稀釋后的樣品總DNA為模板，以包含5′-CCT ACGG R R B GCA S CAG K V R V GAAT-3′序列的上游引物和包含5′-GGAC T AC N V GG GT W TC T AATCC-3′序列的下游引物對(duì)細(xì)菌16S rRNA的2個(gè)高度可變區(qū)(V3和 V4區(qū))進(jìn)行PCR擴(kuò)增[14]。擴(kuò)增體系(25 μL)為：TransStart Buffer 2.5 μL，dNTPs 2 μL，上下游引物各1 μL，模板 DNA 20 ng，用ddH2O補(bǔ)齊至25 μL。PCR 反應(yīng)程序：95℃預(yù)變性2 min；95℃、55℃、72℃分別持續(xù)30 s，共30次循環(huán)；72℃延伸5 min，得到的產(chǎn)物用1.5%瓊脂糖凝膠電泳檢測(cè)完整性。PCR擴(kuò)增完成后使用 MetaVxTM文庫(kù)構(gòu)建試劑盒構(gòu)建測(cè)序所需文庫(kù)。使用 Agilent 2100 生物分析儀檢測(cè)構(gòu)建完成的文庫(kù)質(zhì)量，再通過(guò)Qubit3.0熒光計(jì)檢測(cè)文庫(kù)的濃度[15]。將構(gòu)建好的文庫(kù)定量到10 nmol·L-1，按照Illumina MiSeq測(cè)序儀使用說(shuō)明書(shū)進(jìn)行PE250/FE300雙端測(cè)序[16]。PCR 擴(kuò)增、文庫(kù)構(gòu)建及測(cè)序由蘇州金唯智生物科技有限公司完成。

1.3.3 測(cè)序相關(guān)分析 將雙端測(cè)序所得的正反向reads進(jìn)行兩兩組裝，過(guò)濾去除組裝結(jié)果中含有N的序列且保留序列長(zhǎng)度大于200 bp的序列。所得的序列再經(jīng)過(guò)質(zhì)量過(guò)濾，去除嵌合體序列。使用VSEARCH(1.9.6)進(jìn)行序列聚類(lèi)，將97%的序列相似水平作為劃分閾值，聚類(lèi)成可操作分類(lèi)單位(operational taxonomic units，OTU)[17]?；贠TU的分析結(jié)果，分別計(jì)算得到ACE指數(shù)、Chao1 指數(shù)、Shannon指數(shù)等α多樣性指數(shù)，并繪制稀釋曲線(xiàn)[18]。使用QIIME軟件對(duì)基于非加權(quán)距離矩陣進(jìn)行算術(shù)平均組對(duì)法(unweighted pair group method with artithe metic mean, UPGMA)聚類(lèi)分析，用主坐標(biāo)分析(principal co-ordinates analysis，PCoA)可視化圖展示出樣品的β多樣性[19]。最后采用PICRUSt(V1.0)軟件對(duì)細(xì)菌群落的基因功能特征進(jìn)行預(yù)測(cè)分析[20]。

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

利用Vsearch(1.9.6)、Qiime(1.9.1)、Metastats、 LEfSE(1.0)、Cutadapt(v1.9.1)、PICRUST(v1.0.0)等軟件進(jìn)行數(shù)據(jù)分析。

2 結(jié)果與分析

2.1 PCR擴(kuò)增基本情況

由圖1可知，13份樣品的PCR 擴(kuò)增產(chǎn)物凝膠電泳在500 bp左右出現(xiàn)了清晰明顯的條帶，說(shuō)明目的片段擴(kuò)增成功。

2.2 樣品序列數(shù)特點(diǎn)

通過(guò)細(xì)菌的16S rRNA基因V3-V4 區(qū)測(cè)序，13份樣品共產(chǎn)生720 221條高質(zhì)量序列，平均長(zhǎng)度為463 bp。高質(zhì)量序列占有效序列的比例為94%。將所有的有效序列按照97%的相似度進(jìn)行OTU 聚類(lèi)，共聚成166個(gè)OTU。

注：1：M1；2： M2；3：M3；4：M4；5：W1；6：W2；7：W3；8：W4；9：D1；10：D2；11：D3；12：D4；13：QL1.圖1 PCR擴(kuò)增產(chǎn)物凝膠電泳圖Fig.1 Electrophoresis of PCR amplified product

2.3 稀釋曲線(xiàn)

稀釋曲線(xiàn)表示隨著測(cè)序量的增加，可能會(huì)檢測(cè)到的物種種類(lèi)隨之增加的狀況，被廣泛用于判斷測(cè)序是否達(dá)到一定的深度[21]。由圖2可知，隨著測(cè)序量的增加，樣品M1、M2、M3、M4、W1、W2、W3、W4、D1、D2、D3、D4、QL1的OTUs逐漸增加后趨于平穩(wěn)，表明測(cè)序深度充分，基本覆蓋了樣品中所有的物種。

表1 各樣品Alpha多樣性指數(shù)Table 1 Alpha diversity index of each sample

圖2 各樣品稀釋曲線(xiàn)Fig.2 Rarefaction curves of each sample

2.4 Alpha多樣性指數(shù)分析

Alpha多樣性用一系列的統(tǒng)計(jì)學(xué)指數(shù)來(lái)評(píng)估樣本的物種豐富度及多樣性，主要包括Chao1指數(shù)、ACE指數(shù)、Simpson指數(shù)、Shannon指數(shù)等[22]。Shannon指數(shù)越高該樣品的多樣性越豐富。由表1可知，13份樣品的多樣性指數(shù)均不相同。樣品Shannon指數(shù)從高到低依次為QL1>M4>W2>M2>D4>M1>D3>D1>M3=D2>W1>W4>W3，說(shuō)明不同來(lái)源樣品的細(xì)菌菌群多樣性不同。

2.5 樣品細(xì)菌菌群在門(mén)水平的分類(lèi)

由表2可知，13份樣品中共鑒定出4個(gè)細(xì)菌門(mén)，分別是厚壁菌門(mén)(Firmicutes)、變形菌門(mén)(Proteobacteria)、藍(lán)菌門(mén)(Cyanobacteria)和擬桿菌門(mén)(Bacteroidetes)。厚壁菌門(mén)(Firmicutes)為13份樣品中的優(yōu)勢(shì)門(mén)，13份樣品的相對(duì)豐度介于99.95%～99.99%。變形菌門(mén)(Proteobacteria)為亞優(yōu)勢(shì)門(mén)，藍(lán)菌門(mén)(Cyanobacteria)和擬桿菌門(mén)(Bacteroidetes)為低豐度門(mén)，這兩類(lèi)門(mén)在樣品Q(chēng)L1中均未鑒定到。由門(mén)水平分類(lèi)來(lái)看，不同樣品在門(mén)水平由不同的菌落組成，對(duì)各菌門(mén)表現(xiàn)出不同的豐度。

表2 各樣品門(mén)水平菌群分布相對(duì)豐度Table 2 The relative abundance of horizontal flora distribution aTThe phylum level of each sample

2.6 樣品細(xì)菌菌群在屬水平的分類(lèi)

由表3可知，13份樣品共鑒定出了4種不同的細(xì)菌菌屬以及在屬水平上未鑒定的屬(Unclassified)。其中，乳桿菌屬(Lactobacillus)為高豐度菌屬，鏈球菌屬(Streptococcus)為第二豐度菌屬，在13份樣品中均檢測(cè)到。魯杰氏菌屬(Ruegeria)和普雷沃菌屬(Prevotella)為低豐度菌屬，僅在部分樣品(M4、D1、D4)中檢測(cè)到，而在QL1中兩種屬均未檢測(cè)到。表明不同樣品在屬水平的分類(lèi)不同，對(duì)各屬表現(xiàn)出不同的豐度值。

表3 各樣品屬水平菌群分布相對(duì)豐度Table 3 The relative abundance of horizontal flora distribution aTThe genus level of each sample

2.7 β多樣性分析

β多樣性是對(duì)不同樣品中微生物群落組成的比較分析，從而反映出不同樣本間的多樣性的差異[23]。

2.7.1 主坐標(biāo)分析 主坐標(biāo)分析(PCoA)是通過(guò)R語(yǔ)言來(lái)分析樣品間的相似性或者差異性[24]。13份樣品的PCoA分析如圖3所示，PC1-PC2為第一和第二主坐標(biāo)得到的PCoA圖，PC1-PC3為第一和第三主坐標(biāo)得到的PCoA圖。PC1、PC2、PC3對(duì)樣品差異的貢獻(xiàn)率分別為99.4%、0.38%、0.13%。樣本之間距離代表了樣本微生物群落的相似性，距離越近，相似度越高[25]。由圖3可知，在PC1和PC2維度上部分樣品交疊在一起，而在PC3維度上所有樣品可以很好地區(qū)分。由樣品PCoA分析可知，在PC2維度上W1、W3、W4、D1、D2、M3距離較近有著相似的群落， M1、M2、D3、D4距離較近有著相似的群落，W2和M4距離較近有著相似的群落，整體甘南地區(qū)的樣品群落相似性較高，與QL1相距較遠(yuǎn)存在顯著差異。

圖3 各樣品PCoA分析Fig.3 Principal co-ordinates analysis of each sample

2.7.2 樣品菌群結(jié)構(gòu)聚類(lèi)分析 樣品聚類(lèi)分析是指使用算術(shù)平均組對(duì)法(UPGMA)分析樣本在特定進(jìn)化譜系中是否有顯著微生物群落差異[26]。兩樣本越相似，樣本間的分枝長(zhǎng)度越短。由圖4可知，13份牦牛酸乳樣品聚為兩大類(lèi)，第一大類(lèi)中W1、W3、W4、D1、D2、M3分支較小，相似度較高；M1、M2、D3、D4分支較小，相似度較高，W2和M4分支較小，相似性較高聚為一類(lèi)。QL1單獨(dú)聚為一類(lèi)，表現(xiàn)出與甘肅甘南地區(qū)樣品相比不同的微生物群落特征。表明甘南地區(qū)12份樣品(W1、W3、W4、D1、D2、M1、M2、D3、D4)的菌群差異不顯著，有著相似的進(jìn)化譜系，而QL1單獨(dú)聚為一類(lèi)，與甘南地區(qū)樣品差異較大，這與PCoA分析結(jié)果一致。

圖4 各樣品UPGMA菌群結(jié)構(gòu)聚類(lèi)分析Fig.4 UPGMACluster analysis of the bacterial communities in each sample

2.8 樣品的PICRUSt功能預(yù)測(cè)分析

PICRUSt軟件(phylogenetic investigation of communities by reconstruction of unobserved states)，將測(cè)序得到的菌群組成“映射”錄入Greengenes數(shù)據(jù)庫(kù)中，實(shí)現(xiàn)對(duì)菌群代謝功能的預(yù)測(cè)。目前已經(jīng)廣泛用于微生物生態(tài)學(xué)[27]。由表4可知，在13份樣品的功能基因中，基本預(yù)測(cè)功能占主導(dǎo)地位，占所有基因的11.76%，其次占比較多的是復(fù)制/重組/修復(fù)、碳水化合物代謝相關(guān)的基因，分別占所有基因的9.54%、9.39%。

3 討論

本研究基于Illumina MiSeq高通量技術(shù)對(duì)甘肅甘南傳統(tǒng)發(fā)酵牦牛乳中細(xì)菌菌群的多樣性進(jìn)行分析，結(jié)果表明，不同來(lái)源的傳統(tǒng)發(fā)酵牦牛乳樣品中細(xì)菌菌群的結(jié)構(gòu)存在一定差異。本研究使用的發(fā)酵牦牛乳是通過(guò)傳統(tǒng)方法生產(chǎn)的，因此樣品間細(xì)菌群落結(jié)構(gòu)的差異可能與原料、地理位置和其他環(huán)境因素有關(guān)[28-29]。在門(mén)水平上，厚壁菌門(mén)(Firmicutes)是優(yōu)勢(shì)菌門(mén)，每個(gè)樣品表現(xiàn)出不同的豐度，此外還鑒定出低豐度的擬桿菌門(mén)和藍(lán)菌門(mén)，這與新疆塔城地區(qū)發(fā)酵駱駝乳[30]、新疆昭蘇縣和特克斯縣發(fā)酵酸牛奶[31]、云南地區(qū)發(fā)酵山羊奶制品[32]的研究結(jié)果一致。在屬水平上，乳桿菌屬(Lactobacillus)為高豐度菌屬, 鏈球菌屬(Streptococcus)為第二豐度屬，這兩個(gè)屬通常與天然發(fā)酵乳制品有關(guān)，在乳制品的生產(chǎn)發(fā)酵過(guò)程中發(fā)揮著重要作用[33-34]。據(jù)報(bào)道，乳桿菌屬的酸耐受能力高于鏈球菌屬和乳球菌屬[35]。在牦牛乳發(fā)酵過(guò)程中，隨著酸度的上升，高比例的乳酸桿菌被保留富集，因此乳桿菌屬成為高豐度的菌屬。此外，在部分樣品中還出現(xiàn)了低豐度的魯杰氏菌屬(Ruegeria)和普雷沃菌屬(Prevotella)。魯杰氏菌屬最早由Uchino等[36]于1988年引入，屬于有益細(xì)菌屬，具有較高的磷酸二酯酶和磷酸單酯酶活性，在水生生物體體內(nèi)大量存在，但在以往的發(fā)酵乳制品中未見(jiàn)該屬的相關(guān)報(bào)道[37]。普雷沃菌屬來(lái)自擬桿菌門(mén)，在厭氧環(huán)境中繁殖，是腸道菌群的主要組成部分，有助于降解食物中的碳水化合物和多糖，參與合成多種維生素等[38]。Quigley等[38]利用高通量測(cè)序技術(shù)揭示了62種愛(ài)爾蘭手工奶酪中的細(xì)菌群落，結(jié)果發(fā)現(xiàn)，奶酪在發(fā)酵過(guò)程中有幾個(gè)菌屬參與，其中就包括普雷沃菌屬，分析原因是原料乳中帶入了普雷沃菌屬。因此推斷本研究中甘南地區(qū)樣品檢測(cè)出的低豐度的魯杰氏菌屬和普雷沃菌屬可能來(lái)源于原料乳或發(fā)酵環(huán)境。

表4 各樣品COG功能分類(lèi)相對(duì)豐度統(tǒng)計(jì)表Table 4 Statistical table of relative abundance of COG functional classification of each sample

有學(xué)者基于PICRUSt對(duì)哺乳動(dòng)物腸道微生物菌群的功能進(jìn)行綜合評(píng)估，發(fā)現(xiàn)哺乳動(dòng)物體內(nèi)的部分代謝與腸道微生物有關(guān)[39-40]。在本研究中，基于每個(gè)樣品的16S rRNA測(cè)序數(shù)據(jù)進(jìn)行PICRUSt分析發(fā)現(xiàn)，基本預(yù)測(cè)功能占主導(dǎo)地位，占所有基因功能的11.76%。其次占比較多的是復(fù)制/重組/修復(fù)、碳水化合物代謝相關(guān)基因，分別占所有基因的9.54%、9.39%。正是存在于乳酸菌中的各個(gè)功能基因，乳酸菌表現(xiàn)出了優(yōu)良的特性。

4 結(jié)論

本研究基于Illumina MiSeq高通量技術(shù)對(duì)甘肅甘南地區(qū)傳統(tǒng)發(fā)酵牦牛乳進(jìn)行細(xì)菌菌群的多樣性分析。結(jié)果表明，甘肅甘南地區(qū)傳統(tǒng)發(fā)酵牦牛乳中具有豐富的細(xì)菌多樣性，優(yōu)勢(shì)菌門(mén)為厚壁菌門(mén)(Firmicutes)，優(yōu)勢(shì)菌屬為乳桿菌屬(Lactobacillus)，不同來(lái)源的發(fā)酵牦牛乳在門(mén)屬水平呈現(xiàn)不同的豐度。β多樣性分析結(jié)果表明其群落相似性較高，有相似的進(jìn)化譜系。通過(guò)預(yù)測(cè)功能基因，發(fā)現(xiàn)傳統(tǒng)發(fā)酵牦牛乳中的絕大部分細(xì)菌與基本預(yù)測(cè)功能有關(guān)。