新鮮曲拉與干燥曲拉細菌群落結構的比較分析

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

（1.西藏農牧學院食品科學學院，西藏林芝860000；2.西藏自治區藏藥審評認證中心，西藏拉薩850000）

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

曲拉作為一種藏區典型的發酵型乳制品，其制作多采用傳統的手工開放式自作，因此不同的制作環境直接關系到曲拉的微生物組成[6-7]，進而影響曲拉的風味、安全及營養價值[8-10]。

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

1 材料與方法

1.1 材料與試劑

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

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

1.2 儀器

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

1.3 方法

1.3.1 樣品的采集及制備

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

圖1 曲拉樣品Fig.1 Qula samples

1.3.2 DNA的提取和PCR擴增

1.3.3 文庫建立

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

1.3.4 Illumina MiSeq

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

2 結果與分析

2.1 測序數據統計及OTU組成分析

樣品OUT聚類和多樣性指數見表1。

表1 不同樣本曲拉的測序數據統計Table 1 Sequence data statistics in different Qula samples

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

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

2.2 稀釋曲線

XQL和GQL稀釋曲線見圖2。

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

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

2.3 各分類水平上的OUT數

樣品各等級OUT物種統計表見表2。

表2 樣品各等級OUT物種統計表Table 2 Samples of various grades OUT species statistics

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

聚類分析到種水平的數據比在屬水平上的少，分析原因可能是16SrDNA對親緣關系較近的種屬無法進一步區分的造成的。16SrDNA中的保守區反映不同物種間親緣關系，而可變區則主要體現不同物質之間的差異。由于16SrDNA高度保守，對親緣關系較近的種屬分辨率不高，因此通過高通量測序技術基本上可以在屬水平上準確的區分細菌類群。

2.4 樣品在門水平上的分類比較

各樣品門水平菌群分布相對豐度見表3。

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

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

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

樣品細菌群落在屬水平的分類比較，在6份樣品中共檢出25個屬，見表4。

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

從屬水平分析，新鮮曲拉和干燥曲拉的優勢菌均為乳球菌屬、沙雷氏菌屬、不動桿菌屬和假單胞菌屬，豐度值分別為 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 物種豐度熱圖

根據OUT的豐度數據，熱圖可在屬水平上將不同的OUT分塊聚類，并根據顏色梯度反應不同樣品中細菌群落相似性、差異性及物種聚類關系。利用R語言繪制熱圖，可以直觀地顯示2種6個曲拉樣品中細菌的差異，見圖3。

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

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

2.7 樣本聚類分析

各樣品基于UniFrac距離的曲拉樣品的聚類圖見圖4。

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

2.8 樣品功能類群分布統計

通過將現有的16S rRNA基因測序數據與代謝功能已知的微生物參考基因組數據庫相對比，從而實現對細菌和古菌代謝功能的預測，見圖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 結論

本研究采用高通量測序技術對西藏林芝地區新鮮（軟）曲拉和干燥（硬）曲拉共2組樣品的細菌群落結構進行比較分析，結果從細菌多樣性角度分析，干（硬）曲拉的Simpson指數的平均值為0.90±0.01、Chao1指數平均值為760.21±133.80、ACE指數平均值為775.91±116.98、Shannon指數平均值為5.46±0.05均優于而新鮮曲拉的Simpson指數平均值為0.90±0.01、Chao1指數平均值為595.31±45.23、ACE指數平均值為602.15±49.33、Shannon指數平均值為5.39±0.17，這說明干（硬）曲拉的細菌多樣性更為豐富。干（硬）曲拉和新鮮（軟）曲拉在門水平上和屬水平上的比較發現二者優勢菌門和優勢菌屬一致均為厚壁菌門（Firmicutes）和變形菌門（Proteobacteria），乳球菌屬（Lactococcus）、沙雷氏菌屬（Serratia）、不動桿菌屬（Acinetobacter），但是曲拉在干燥過程中優勢菌門和菌屬的相對豐度發生了明顯的變化，其中厚壁菌門由0.53±0.023 7增至0.61±0.034 7；變形菌門由 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；不動桿菌屬（Acinetobacter）0.07±0.004 6降至0.04±0.004 4；由此可知曲拉的第一優勢菌乳球菌屬作為曲拉的主要功能菌[19-20]，其在干燥的過程中豐度發生了顯著的增加，而作為條件致病菌屬的沙雷氏菌屬在干燥的過程中豐度顯著降低[21]；由此可知曲拉的傳統形式干燥曲拉在食用功能性和安全性上占據一定的優勢。