劉振東，劉怡，李哲，畢娜，李梁，羅章，薛蓓，*

（1.西藏農(nóng)牧學(xué)院食品科學(xué)學(xué)院，西藏林芝860000；2.西藏自治區(qū)藏藥審評(píng)認(rèn)證中心，西藏拉薩850000）

牦牛曲拉，又名奶渣，是由牦牛奶脫脂后在自然條件下發(fā)酵凝固并風(fēng)干制成的粗奶酪[1-3]。作為一種傳統(tǒng)的藏族食品，由于蛋白質(zhì)含量高、益生菌種類多、保存時(shí)間長(zhǎng)，攜帶方便等特點(diǎn)，曲拉在整個(gè)藏區(qū)廣泛食用[4-5]。調(diào)查發(fā)現(xiàn)，藏民除食用傳統(tǒng)的干曲拉即硬曲拉外，還有食用新鮮曲拉的習(xí)慣。由于冰箱在藏區(qū)家庭的普及，這種在室溫下很難保存的濕曲拉，由于口感軟糯，易于消化，老少皆宜等特點(diǎn)受到越來(lái)越多藏族同胞的喜愛(ài)。而目前關(guān)于新鮮曲拉的研究尚未見(jiàn)相關(guān)研究報(bào)道。

曲拉作為一種藏區(qū)典型的發(fā)酵型乳制品，其制作多采用傳統(tǒng)的手工開(kāi)放式自作，因此不同的制作環(huán)境直接關(guān)系到曲拉的微生物組成[6-7]，進(jìn)而影響曲拉的風(fēng)味、安全及營(yíng)養(yǎng)價(jià)值[8-10]。

傳統(tǒng)基于微生物分離、純化、培養(yǎng)及鑒定的研究方法難以獲得難培養(yǎng)及痕量微生物，難以將環(huán)境中完成的微生物的種群結(jié)構(gòu)及生態(tài)關(guān)系展現(xiàn)出來(lái)[11-12]，高通量測(cè)序技術(shù)（high-throughput sequencing，HTS）因其不僅可以檢測(cè)到那些難培養(yǎng)和低豐度的物種以及曾經(jīng)存活或是難以分離的部分微生物，還可以排除挑選菌株時(shí)所產(chǎn)生的偏見(jiàn)性，更能代表整體的微生物組成，因此已經(jīng)逐漸成為研究微生物組成及結(jié)構(gòu)的一種重要方法[13]，近年該技術(shù)已廣泛用于原料乳和各類發(fā)酵乳等乳制品生態(tài)位的研究[14-15]。劉怡萱等[16]通過(guò)該技術(shù)對(duì)西藏農(nóng)、牧區(qū)牦牛酸奶菌群比較分析，發(fā)現(xiàn)厚壁菌門(mén)中乳桿菌屬和乳球菌屬為優(yōu)勢(shì)菌群；趙順先等[17]通過(guò)該技術(shù)對(duì)新疆傳統(tǒng)干奶酪乳酸菌多樣性進(jìn)行分析，發(fā)現(xiàn)干奶酪中乳酸菌歸為2個(gè)目7個(gè)科12個(gè)屬。其中乳桿菌屬、雙歧桿菌屬、魏斯氏菌屬和片球菌屬是優(yōu)勢(shì)菌屬；張敏等[18]利用該技術(shù)對(duì)新疆不同動(dòng)物來(lái)源的原奶和酸奶細(xì)菌多樣進(jìn)行了比較分析，發(fā)現(xiàn)原奶樣品中細(xì)菌Shannon-Wiener指數(shù)明顯高于酸奶樣品，原奶樣品主要以變形菌門(mén)為主，而酸奶樣品主要以厚壁菌門(mén)為主。本研究分別以西藏林芝地區(qū)農(nóng)牧民自制的新鮮（軟）曲拉和干（硬）曲拉為研究對(duì)象，利用Illumina MiSeq高通量測(cè)序技術(shù)對(duì)二者的細(xì)菌群落結(jié)構(gòu)進(jìn)行全面的分析，明確新鮮（軟）曲拉和干（硬）曲拉各自的細(xì)菌群落組成及差異，為二者的功能菌群的挖掘和食用安全性的指導(dǎo)奠定一定的基礎(chǔ)。

1 材料與方法

1.1 材料與試劑

牦牛曲拉樣品：采自林芝市米林縣邦仲村。

E.Z.N.A.Soil DNA試劑盒：美國(guó)OMEGA公司；Qubit2.0 DNA檢測(cè)試劑盒：美國(guó)Invitrogen公司；Q5DNA聚合酶：美國(guó)Thermo Fisher公司；凝膠回收試劑盒：德國(guó)Qiagen公司。

1.2 儀器

5424R高速冷凍離心機(jī)：德國(guó)Eppendorf公司；DYY-6C型電泳儀：北京六一儀器廠；T-100 PCR儀：美國(guó)Bio-rad公司；MiSeq高通量測(cè)序儀：美國(guó)Illumina公司。

1.3 方法

1.3.1 樣品的采集及制備

新鮮（軟）曲拉（XQL）樣品于2018年6月采自西藏林芝米林縣米林鎮(zhèn)幫仲村。將一半XQL樣品置于無(wú)菌離心管中，冰盒保存并立即運(yùn)送到實(shí)驗(yàn)室，在超凈工作臺(tái)上均分為 3 份（XQL1、XQL2、XQL3），置于-80冰箱中備用；將剩余一部份XQL在采集點(diǎn)進(jìn)行傳統(tǒng)日曬干制得到干（硬）曲拉樣品（GQL），同時(shí)均分為3份（GQL1、GQL2、GQL3），置于無(wú)菌離心管中備用。

圖1 曲拉樣品Fig.1 Qula samples

1.3.2 DNA的提取和PCR擴(kuò)增

1.3.3 文庫(kù)建立

采用三代測(cè)序技術(shù)制備測(cè)序文庫(kù)。第一步是把通過(guò)PCR擴(kuò)增得到的產(chǎn)物的末端修復(fù)將黏性末端修整為平末端；在DNA序列的3’端添加A堿基以防止DNA片段自連；在DNA序列5’端添加含有文庫(kù)特異性標(biāo)簽的測(cè)序接頭；對(duì)上述連上接頭的DNA片段進(jìn)行PCR擴(kuò)增，從而富集測(cè)序文庫(kù)模板，并采用BECKMAN AMPure XP Beads再次純化文庫(kù)富集產(chǎn)物；最后對(duì)文庫(kù)做最終的片段選擇與純化。

1.3.4 Illumina MiSeq

測(cè)序結(jié)果分析將Illumina MiSeq測(cè)序得到的是雙端序列數(shù)據(jù)，進(jìn)行質(zhì)量控制得到以97%相似性分類OTUs。然后Blast在Gen Bank數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，并提交到美吉生物的云分析平臺(tái)進(jìn)行多樣性分析（http：//www.i-sanger.com/）。

2 結(jié)果與分析

2.1 測(cè)序數(shù)據(jù)統(tǒng)計(jì)及OTU組成分析

樣品OUT聚類和多樣性指數(shù)見(jiàn)表1。

表1 不同樣本曲拉的測(cè)序數(shù)據(jù)統(tǒng)計(jì)Table 1 Sequence data statistics in different Qula samples

由表1可以看出XQL的測(cè)序序列數(shù)29 566.67±1 703.40，GQL 的測(cè)序序列數(shù) 33 604.67±2 042.39；XQL的 OTUs數(shù)為 981.00±14.80，GQL 的 OTUs數(shù)為 841.00±15.87。由序列數(shù)和OTUs數(shù)可知XQL和GQL中細(xì)菌種類豐富，且二者存在著明顯的差異。

由表1可知XQL的Simpson指數(shù)平均值為0.90±0.01、Chao1 指數(shù)平均值為 595.31±45.23、ACE 指數(shù)平均值為602.15±49.33、Shannon指數(shù)平均值為5.39±0.17。GQL的 Simpson指數(shù)的平均值為 0.90±0.01、Chao1指數(shù)平均值為760.21±133.80、ACE指數(shù)平均值為 775.91±116.98、Shannon 指數(shù)平均值為 5.46±0.05。從多樣性指數(shù)結(jié)果分析可以得出GQL細(xì)菌多樣性指數(shù)相關(guān)參數(shù)均高于XQL，這表明曲拉樣品在自然干燥的過(guò)程，影響了曲拉細(xì)菌的群落結(jié)構(gòu)。

2.2 稀釋曲線

XQL和GQL稀釋曲線見(jiàn)圖2。

圖2 各樣品中細(xì)菌稀釋曲線Fig.2 Rarefaction curves of bacterial communities in different Qula samples

由圖2可知當(dāng)測(cè)序量較低時(shí)每個(gè)樣品的Shannon指數(shù)均隨著測(cè)序量的增加呈顯著上升趨勢(shì)，而當(dāng)測(cè)序量達(dá)到5 000以后Shannon指數(shù)隨著測(cè)序量的增加呈現(xiàn)緩慢平穩(wěn)的趨勢(shì)，稀釋曲線的斜率逐漸降低；當(dāng)測(cè)序量達(dá)到15 000以后，稀釋曲線與X軸幾乎接近平行，Shannon指數(shù)幾乎達(dá)到上限。結(jié)合表1，可知XQL和GQL樣品的測(cè)序序列數(shù)均顯著高于15 000。結(jié)合稀釋曲線可知在當(dāng)前測(cè)序水平下，樣品中的細(xì)菌群落多樣性已能夠充分覆蓋。

2.3 各分類水平上的OUT數(shù)

樣品各等級(jí)OUT物種統(tǒng)計(jì)表見(jiàn)表2。

表2 樣品各等級(jí)OUT物種統(tǒng)計(jì)表Table 2 Samples of various grades OUT species statistics

由表2可知，XQL在門(mén)水平上的OUT數(shù)目8.00±1.00，在綱水平上的OUT數(shù)目9.67±1.53，在目水平上的OUT數(shù)目25.33±1.15，在科水平上的OUT數(shù)目34.00±2.65，在屬水平上 OUT 數(shù)目 72.33±0.58，在種水平上OUT數(shù)目62.67±2.52。GQL在門(mén)水平上的OUT數(shù) 11.00±3.61，在綱水平上的 OUT 數(shù)目 14.33±4.93，在目水平上的OUT數(shù)目33.00±6.08，在科水平上的OUT數(shù)目 48.67±10.26，在屬水平上 OUT 數(shù)目 95.67±18.18，在種水平上OUT數(shù)目71.33±9.45。

聚類分析到種水平的數(shù)據(jù)比在屬水平上的少，分析原因可能是16SrDNA對(duì)親緣關(guān)系較近的種屬無(wú)法進(jìn)一步區(qū)分的造成的。16SrDNA中的保守區(qū)反映不同物種間親緣關(guān)系，而可變區(qū)則主要體現(xiàn)不同物質(zhì)之間的差異。由于16SrDNA高度保守，對(duì)親緣關(guān)系較近的種屬分辨率不高，因此通過(guò)高通量測(cè)序技術(shù)基本上可以在屬水平上準(zhǔn)確的區(qū)分細(xì)菌類群。

2.4 樣品在門(mén)水平上的分類比較

各樣品門(mén)水平菌群分布相對(duì)豐度見(jiàn)表3。

表3 各樣品門(mén)水平菌群分布相對(duì)豐度Table 3 The relative abundance of microbial communities in different Qula samples at the phylum level

樣品細(xì)菌群落在門(mén)水平的分類比較，結(jié)果顯示，在2組樣品中共檢測(cè)出6個(gè)門(mén)，分屬于厚壁菌門(mén)、變形菌門(mén)、擬桿菌門(mén)、酸桿菌門(mén)、棲熱菌門(mén)和螺旋體門(mén)。其中厚壁菌門(mén)為新鮮和干燥曲拉的第一優(yōu)勢(shì)菌門(mén)占比分別為 0.53±0.023 7，0.61±0.034 7；第二優(yōu)勢(shì)菌為變形菌門(mén)占比分別為 0.42±0.019 4，0.35±0.028 7；第三優(yōu)勢(shì)菌為擬桿菌門(mén)占比為 0.04±0.004 2，0.03±0.005 0。而酸桿菌門(mén)、棲熱菌門(mén)和螺旋體門(mén)由于豐度較低，為樣品的稀有菌門(mén)。

2.5 樣品在屬水平上的分類比較

樣品細(xì)菌群落在屬水平的分類比較，在6份樣品中共檢出25個(gè)屬，見(jiàn)表4。

表4 各樣品在屬水平上的分類比較Table 4 The relative abundance of microbial communities in different Qula samples at the genus level

從屬水平分析，新鮮曲拉和干燥曲拉的優(yōu)勢(shì)菌均為乳球菌屬、沙雷氏菌屬、不動(dòng)桿菌屬和假單胞菌屬，豐度值分別為 0.50±0.023 5 和 0.58±0.040 0；0.19±0.011 5 和 0.14±0.011 4；0.07±0.004 6 和 0.04±0.004 4；0.05±0.004 0 和 0.05±0.004 7。

2.6 物種豐度熱圖

根據(jù)OUT的豐度數(shù)據(jù)，熱圖可在屬水平上將不同的OUT分塊聚類，并根據(jù)顏色梯度反應(yīng)不同樣品中細(xì)菌群落相似性、差異性及物種聚類關(guān)系。利用R語(yǔ)言繪制熱圖，可以直觀地顯示2種6個(gè)曲拉樣品中細(xì)菌的差異，見(jiàn)圖3。

由圖3可知，上側(cè)為樣品聚類樹(shù)，左側(cè)為OUT聚類樹(shù)，6個(gè)曲拉樣品可分為2個(gè)大的聚類：新鮮曲拉可聚為一類，干燥曲拉樣品聚為一類；物種豐富度上顏色越亮代表該菌屬在該樣品中的相對(duì)含量高于橫向的其他樣品，通過(guò)圖2可知干燥曲拉樣品主要優(yōu)勢(shì)菌上占據(jù)顯著優(yōu)勢(shì)。

圖3 各樣品不同屬組成和相對(duì)豐度熱圖Fig.3 Heatmap of the bacterial genera in different samples

2.7 樣本聚類分析

各樣品基于UniFrac距離的曲拉樣品的聚類圖見(jiàn)圖4。

從圖4可知，新鮮曲拉和干燥曲拉兩類樣品基本聚集到不同的類別，說(shuō)明兩類樣品中細(xì)菌多樣性存在一定的差異。

2.8 樣品功能類群分布統(tǒng)計(jì)

通過(guò)將現(xiàn)有的16S rRNA基因測(cè)序數(shù)據(jù)與代謝功能已知的微生物參考基因組數(shù)據(jù)庫(kù)相對(duì)比，從而實(shí)現(xiàn)對(duì)細(xì)菌和古菌代謝功能的預(yù)測(cè)，見(jiàn)圖5。

由圖5可知，新鮮曲拉樣品和干燥曲拉樣品在代謝功能上碳水化合物代謝（carbohydrate metabolism）、能量代謝（energy metabolism）、輔酶和維生素代謝（metabolism of cofactors and vitamins）、核苷酸代謝（nucleotide metabolism）為新鮮曲拉和干燥曲拉的主要代謝功能，而新鮮曲拉與干燥曲拉間的各代謝功能差異并不顯著。

圖4 各樣品基于Weighted UniFrac距離矩陣的UPGMA（unweighted pair group method using arithmetic average）聚類分析圖Fig.4 Dendrogram of cluster analysis for the different Qula samples

圖5 各樣品功能類群分布圖Fig.5 The distribution of functional groups map for the different Qula samples

3 結(jié)論

本研究采用高通量測(cè)序技術(shù)對(duì)西藏林芝地區(qū)新鮮（軟）曲拉和干燥（硬）曲拉共2組樣品的細(xì)菌群落結(jié)構(gòu)進(jìn)行比較分析，結(jié)果從細(xì)菌多樣性角度分析，干（硬）曲拉的Simpson指數(shù)的平均值為0.90±0.01、Chao1指數(shù)平均值為760.21±133.80、ACE指數(shù)平均值為775.91±116.98、Shannon指數(shù)平均值為5.46±0.05均優(yōu)于而新鮮曲拉的Simpson指數(shù)平均值為0.90±0.01、Chao1指數(shù)平均值為595.31±45.23、ACE指數(shù)平均值為602.15±49.33、Shannon指數(shù)平均值為5.39±0.17，這說(shuō)明干（硬）曲拉的細(xì)菌多樣性更為豐富。干（硬）曲拉和新鮮（軟）曲拉在門(mén)水平上和屬水平上的比較發(fā)現(xiàn)二者優(yōu)勢(shì)菌門(mén)和優(yōu)勢(shì)菌屬一致均為厚壁菌門(mén)（Firmicutes）和變形菌門(mén)（Proteobacteria），乳球菌屬（Lactococcus）、沙雷氏菌屬（Serratia）、不動(dòng)桿菌屬（Acinetobacter），但是曲拉在干燥過(guò)程中優(yōu)勢(shì)菌門(mén)和菌屬的相對(duì)豐度發(fā)生了明顯的變化，其中厚壁菌門(mén)由0.53±0.023 7增至0.61±0.034 7；變形菌門(mén)由 0.42±0.019 4 降至 0.35±0.028 7；乳球菌屬（Lactococcus）由 0.50±0.023 5 增至0.58±0.040 0；沙雷氏菌屬（Serratia）由 0.19±0.011 5 降至0.14±0.011 4；不動(dòng)桿菌屬（Acinetobacter）0.07±0.004 6降至0.04±0.004 4；由此可知曲拉的第一優(yōu)勢(shì)菌乳球菌屬作為曲拉的主要功能菌[19-20]，其在干燥的過(guò)程中豐度發(fā)生了顯著的增加，而作為條件致病菌屬的沙雷氏菌屬在干燥的過(guò)程中豐度顯著降低[21]；由此可知曲拉的傳統(tǒng)形式干燥曲拉在食用功能性和安全性上占據(jù)一定的優(yōu)勢(shì)。