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

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

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

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

1 材料與方法

1.1 材料與試劑

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

TruSeqTMDNA試劑盒，美國Illumina公司；AP-MN-P-50凝膠提取試劑盒，美國Axygen Biosciences公司； MetaVxTM文庫構建試劑盒，美國GENEWIZ公司。

1.2 儀器與設備

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

1.3 試驗方法

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

1.3.2 PCR擴增、文庫構建及上機測序 以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′序列的下游引物對細菌16S rRNA的2個高度可變區(qū)(V3和 V4區(qū))進行PCR擴增[14]。擴增體系(25 μL)為：TransStart Buffer 2.5 μL，dNTPs 2 μL，上下游引物各1 μL，模板 DNA 20 ng，用ddH2O補齊至25 μL。PCR 反應程序：95℃預變性2 min；95℃、55℃、72℃分別持續(xù)30 s，共30次循環(huán)；72℃延伸5 min，得到的產物用1.5%瓊脂糖凝膠電泳檢測完整性。PCR擴增完成后使用 MetaVxTM文庫構建試劑盒構建測序所需文庫。使用 Agilent 2100 生物分析儀檢測構建完成的文庫質量，再通過Qubit3.0熒光計檢測文庫的濃度[15]。將構建好的文庫定量到10 nmol·L-1，按照Illumina MiSeq測序儀使用說明書進行PE250/FE300雙端測序[16]。PCR 擴增、文庫構建及測序由蘇州金唯智生物科技有限公司完成。

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

1.4 數據分析

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

2 結果與分析

2.1 PCR擴增基本情況

由圖1可知，13份樣品的PCR 擴增產物凝膠電泳在500 bp左右出現了清晰明顯的條帶，說明目的片段擴增成功。

2.2 樣品序列數特點

通過細菌的16S rRNA基因V3-V4 區(qū)測序，13份樣品共產生720 221條高質量序列，平均長度為463 bp。高質量序列占有效序列的比例為94%。將所有的有效序列按照97%的相似度進行OTU 聚類，共聚成166個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擴增產物凝膠電泳圖Fig.1 Electrophoresis of PCR amplified product

2.3 稀釋曲線

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

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

圖2 各樣品稀釋曲線Fig.2 Rarefaction curves of each sample

2.4 Alpha多樣性指數分析

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

2.5 樣品細菌菌群在門水平的分類

由表2可知，13份樣品中共鑒定出4個細菌門，分別是厚壁菌門(Firmicutes)、變形菌門(Proteobacteria)、藍菌門(Cyanobacteria)和擬桿菌門(Bacteroidetes)。厚壁菌門(Firmicutes)為13份樣品中的優(yōu)勢門，13份樣品的相對豐度介于99.95%～99.99%。變形菌門(Proteobacteria)為亞優(yōu)勢門，藍菌門(Cyanobacteria)和擬桿菌門(Bacteroidetes)為低豐度門，這兩類門在樣品QL1中均未鑒定到。由門水平分類來看，不同樣品在門水平由不同的菌落組成，對各菌門表現出不同的豐度。

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

2.6 樣品細菌菌群在屬水平的分類

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

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

2.7 β多樣性分析

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

2.7.1 主坐標分析 主坐標分析(PCoA)是通過R語言來分析樣品間的相似性或者差異性[24]。13份樣品的PCoA分析如圖3所示，PC1-PC2為第一和第二主坐標得到的PCoA圖，PC1-PC3為第一和第三主坐標得到的PCoA圖。PC1、PC2、PC3對樣品差異的貢獻率分別為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相距較遠存在顯著差異。

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

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

圖4 各樣品UPGMA菌群結構聚類分析Fig.4 UPGMACluster analysis of the bacterial communities in each sample

2.8 樣品的PICRUSt功能預測分析

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

3 討論

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

表4 各樣品COG功能分類相對豐度統計表Table 4 Statistical table of relative abundance of COG functional classification of each sample

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

4 結論

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