基于16S rRNA高通量測序的西藏農、牧區牦牛酸奶菌群多樣性分析

劉怡萱，許國琪，曹鵬熙，金彥龍，李小燕，劉 星,,

（1.西藏大學理學院，西藏 拉薩 850000；2.西藏大學生態學研究中心，極端環境生物資源與適應性進化實驗室，西藏 拉薩 850000；3.武漢大學生命科學學院，武漢大學-西藏大學生態與保護聯合研究中心，湖北 武漢 430072）

牦牛在青藏高原的馴化起始于一萬年前，是極端環境條件下生存繁衍的優勢畜種，牦牛奶是藏民族賴以生存的奶源和重要的經濟來源之一[1-3]。除直接飲用外，經過加工產生的多種牦牛乳制品包括酸奶、曲拉、酥油、奶酪、乳清、奶茶等也是高原人民生產生活中不可或缺的食物，其中牦牛酸奶更是藏族人民最常食用的乳制品之一[4]。在西藏強紫外線、低壓缺氧、晝夜溫差大的極端環境下，對于發酵溫度及時間沒有精準控制的自然發酵牦牛酸奶的營養價值越來越受到關注，其中特殊的微生物資源也逐漸成為研究的熱點[5]。傳統發酵牦牛酸奶具有獨特的風味和口感，參與發酵的各種微生物進入人體消化系統，部分菌群在人體的腸道菌群的建成中起重要作用，對人類健康產生影響[6-7]。

對發酵乳制品中微生物組成的研究始于19世紀末，早期研究中，采用傳統分離培養的方法對其參與發酵微生物組成有了初步的認識[8-10]。然而自然界中可培養研究的微生物還不到整個微生物種群的20%，在各種環境樣品中僅有小于1%的微生物能通過實驗室培養技術獲得純培養，大部分微生物還屬于不可培養或難以培養的微生物[11-12]。此外，傳統的免培養技術，如變性梯度凝膠電泳[13-14]、限制性片段長度多態性[15]、熒光原位雜交[16]等所能鑒定到的微生物種類有限[6]。隨著測序技術的發展，高通量測序已成為檢測微生物多樣性的重要手段，其讀長高、精度高、通量高、無偏性的優勢被廣泛應用到微生物群落多樣性的研究中，不僅可以檢測到難培養和低豐度的物種以及曾經存活或是難以分離的部分微生物，還可以排除挑選菌株時所產生的偏見性，更能代表整體的微生物組成[5-6]。

極端環境下的微生物經過長期的自然選擇，具備相對特殊的結構、機能以及遺傳特性，以應答相應的脅迫因子[5]。目前，對于青藏高原傳統發酵牦牛酸奶已有相關報道涉及加工工藝優化[17]、營養成分分析[18]、乳酸菌和酵母菌分離與鑒定[19]、乳酸菌益生作用[5]、揮發性風味成分分析[20]、微生物組成及其多樣性分析[4,21-22]等，但進一步對西藏農、牧區牦牛酸奶中微生物差異及其原因的研究則鮮見報道。

本研究利用高通量測序方法對西藏不同生產方式下傳統發酵的牦牛酸奶細菌群落組成和結構進行比較分析，并結合環境因子分析農、牧區牦牛酸奶間微生物多樣性差異，以期為牦牛酸奶的進一步開發利用提供基礎參考。

1 材料與方法

1.1 材料與試劑

共采集西藏牧民以傳統方式發酵的牦牛酸奶樣本19 個，包括在日喀則（RK）、拉薩（LA）等農區采集的樣本5 個及在阿里（AL）、那曲（NQ）等牧區采集的樣本14 個，樣品信息見表1。樣品收集在15 mL無菌離心管中，用-20 ℃車載冰箱運回實驗室后，保存于-80 ℃冰箱備用。

聚合酶鏈式反應（plymerase chain reaction，PCR）Master Mix（M7505） 美國Promega公司；AxyPrep DNA凝膠回收試劑盒 美國Axygen公司；DL2000 DNA Marker 寶生物工程（大連）有限公司。

1.2 儀器與設備

C12S車載冰箱 美國Mobicool公司；DW86L-158超低溫冰箱 中國江蘇勝藍公司；NanoDrop 2000超微量分光光度計 美國Thermo Scientific公司；5417R冷凍離心機德國Eppendorf公司；Veriti PCR儀 美國ABI公司；HiSeq 2500測序平臺 美國Illumina公司。

1.3 方法

1.3.1 菌群DNA提取和PCR擴增測序

牦牛酸奶中菌群總DNA的提取參照Yu Jie等[23-24]方法。用提取的總DNA作為模板，以338F（5’-ACTCCTACGGGAGGCAGCAG-3’）、806R（5’-GGACTACHVGGGTWTCTAAT-3’）為引物，對細菌的16S rRNA V3-V4可變區進行特異性擴增。PCR擴增體系（50 μL）：10×Buffer 5 μL，4×dNTP混合物5 μL，MgCl23 μL，引物各2.6 μL，模板2 μL，Taq酶0.4 μL，加ddH2O補充至50 μL。反應程序為：95 ℃預變性2 min后，進行25 次循環（95 ℃變性30 s，55 ℃復性30 s，72 ℃延伸30 s），最后72 ℃延伸5 min。在超微量紫外分光光度計上檢測PCR產物純度和濃度，1%瓊脂糖凝膠電泳檢測DNA完整性，基于Illumina HiSeq 2500進行測序，PCR擴增及測序工作由廣東仁秀測序公司完成。

基于上海美吉生物I-Sanger云平臺（https://www.i-sanger.com/）、R軟件（version 3.5.1）以及ITOL（https://itol.embl.de/）等軟件平臺進行生物信息學及數據統計分析。

1.3.2 環境因子的獲取

利用長天HOLUX M-241無線GPS記錄器對19 個酸奶樣本進行定位并獲取海拔等信息，利用CHELSA高精度數據庫（Climatologies at high resolution for the earth’s land surface areas，http://chelsa-climate.org/）以及R軟件geosphere軟件包分別調取所在地年平均氣溫、地理距離信息。

1.3.3 菌群多樣性分析

對測序得到各樣本序列進行質量控制后，利用Usearch（vsesion 7.0）軟件，按97%相似度的分類水平劃分可操作分類單元（operation taxonomic units，OTU），并按最小樣本序列數進行抽平，利用mothur（version v.1.30.1）計算α多樣性指數，構建稀釋曲線。

利用Quantitative Insights Into Microbial Ecology（QIIME，version 1.17）計算β多樣性距離矩陣。基于bray_curtis距離進行相似性分析（analysis of similarities，ANOSIM），基于unweighted UniFrac距離對農、牧區牦牛酸奶進行主坐標分析（principal coordinates analysis，PCoA）。

進一步對菌群組成進行多級物種差異判別分析（linear discriminant analysis effect size，LEfSe），再根據分類學組成對樣本按照不同分類水平進行線性判別分析（linear discriminant analysis，LDA，閾值為2），找出對樣本劃分產生顯著性差異影響的菌群。進一步利用組間顯著性差異檢驗，對不同區域之間的菌群進行假設檢驗，評估菌群豐度差異的顯著性水平。

此外，基于地理距離、海拔差異以及年平均氣溫等環境因子數據，分析樣本間差異原因。一方面，利用R軟件geosphere程序包計算得到各樣點間的地理距離矩陣，綜合β多樣性距離矩陣，對兩者進行回歸分析。另一方面，結合基于距離的冗余分析（distance-based redundancy analysis，db-RDA）研究海拔高度以及年平均氣溫與樣本差異的關系。

2 結果與分析

2.1 農、牧區牦牛酸奶菌群多樣性分析

2.1.1α多樣性分析

對19 個樣本進行測序共得到204 345 條有效序列，基于97%相似度的分類水平，共注釋到48 個OTU，測序覆蓋率大于99.97%。各樣本群落豐富度Sobs指數稀釋曲線和群落多樣性Shannon指數稀釋曲線趨于平緩，表明測序數據量合理，可以反映樣本中絕大多數的菌群多樣性信息（圖1）。

α多樣性分析表明（表2），Chao指數在農區樣本和牧區樣本間具有顯著差異（P=0.009 2），農區樣本群落豐富度大于牧區樣本。Shannon指數（P=0.250 7）、Simpson指數（P=0.917 9）在農區樣本與牧區樣本間無顯著差異，表明兩組樣本中菌群多樣性類似。

表2 基于T檢驗的α多樣性指數差異Table 2 Difference in α diversity index based on Student’s T test

2.1.2β多樣性分析

ANOSIM（圖2A）顯示，農、牧區組間菌群結構差異顯著大于組內差異（R=0.673 8，P=0.003），分組具有意義。PCoA（圖2B）表明，來自農區和牧區的樣本分別聚集在一起，其中，LA03樣本介于兩組之間，第1主成分對樣本菌群結構差異的貢獻率為42.08%，第2主成分為20.02%。β多樣性綜合分析顯示，農、牧區牦牛酸奶樣本中菌群結構存在明顯差異。

圖2 β多樣性分析Fig.2 β-Diversity analysis

2.2 農、牧區牦牛酸奶菌群組成及差異分析

在所有樣本中鑒定得到的48 個OTU，分屬于4 個門，18 個屬，在門水平上（圖3A），厚壁菌門（Firmicutes）在農、牧區牦牛酸奶菌群中均占主要優勢（96.589 5%、99.796 1%），而放線菌門（Actinobacteria）、變形菌門（Proteobacteria）、擬桿菌門（Bacteroidetes）在農區樣本中的占比均高于牧區樣本。從屬分類水平看（圖3B），農、牧區牦牛酸奶中乳桿菌屬（Lactobacillus）（34.504 9%、70.379 9%）、乳球菌屬（Lactococcus）（55.687 6%、0.116 9%）和鏈球菌屬（Streptococcus）（3.548 1%、29.289 4%）為主要優勢屬，其他菌屬相對豐度均小于1%。

圖3 西藏農、牧區牦牛酸奶中菌群組成Fig.3 Bacterial composition of yak yogurt samples collected from agricultural and pastoral areas of Tibet

圖4 西藏農、牧區牦牛酸奶中菌群差異分析Fig.4 Analysis of bacterial composition differences between yak yogurt samples collected from agricultural and pastoral areas of Tibet

LEfSe多級物種層級分析（圖4A）顯示，擬桿菌門、放線菌門在農區顯著富集，而在牧區大量富集厚壁菌門；從屬水平來看，乳球菌屬等8 個菌屬在農區顯著富集，牧區僅有鏈球菌屬、乳桿菌屬顯著富集。LDA（LDA閾值＞2.0）（圖4B）與組間差異顯著性檢驗（圖5）共同表明，對農、牧區樣本菌群結構具有顯著影響的菌屬包括乳球菌屬、鏈球菌屬、乳桿菌屬等10 個菌屬（P＜0.05）。

圖5 西藏農、牧區牦牛酸奶中菌群組間差異的顯著性檢驗Fig.5 Significant test of the differences in bacterial composition of yak yogurt samples between agricultural and pastoral areas of Tibet

2.3 環境因子對農、牧區牦牛酸奶菌群結構的影響

對地理距離矩陣和β多樣性距離矩陣擬合得到方程為：Y=0.067 4lnx＋0.305 59（R2=0.028，P=0.028）（圖6A），結果表明地理距離對菌群結構的差異解釋度較低。基于unweighted UniFrac距離算法的db-RDA顯示（圖6B），海拔高度與年平均氣溫對農、牧區牦牛酸奶細菌組成均有顯著性影響，且呈負相關，海拔高度（R2=0.730 1，P=0.001）對樣本間差異的解釋度大于年平均氣溫（R2=0.482 6，P=0.007）。

圖6 環境因子對農、牧區牦牛酸奶菌群結構的影響Fig.6 Effects of environmental factors on the bacterial structures of yak yogurt samples collected from Tibetan agricultural and pastoral areas

3 討 論

農區和牧區是西藏主要的兩種傳統生產生活方式，農區主要位于拉薩、山南、日喀則，集中分布在藏南河谷及河谷盆地中，以種植青稞、油菜等作物為主；牧區主要集中在阿里、那曲，天然草地面積位居全國第一[25]。牦牛作為西藏畜牧業的支柱，養殖主要集中在牧區，占據十分重要的地位[26]，農區由于地理環境原因，牦牛養殖量遠小于牧區，養殖方式與牧區也有較大區別，人工添加飼料比重遠大于牧區。然而，無論在農區還是牧區，傳統發酵的牦牛酸奶都是農牧民發展經濟、增加收入不可或缺的重要產品之一[27]。對自然發酵牦牛酸奶中菌群多樣性的研究發現，厚壁菌門中乳桿菌屬和乳球菌屬為優勢菌群[4,28]。本研究來自西藏的樣品，菌群組成與之前研究結果基本一致，農區樣本中菌群豐富度及多樣性略高于牧區樣本。

發酵食品中菌群多樣性差異與地理環境密切相關[25,29-32]，酸奶中菌群種類在氣候相對溫和的地域比在相對寒冷和海拔較高的環境中更為豐富[5]，自然發酵酸奶中菌群的組成及數量還受發酵時間和制作方法等其他因素影響[22]。本研究除氣候及海拔的因素外，還考慮了農區和牧區養殖方式的不同，農區由于受地理環境的限制，牦牛喂養過程中人工添加飼料比重大于牧區。研究發現，地理距離對酸奶中菌群結構影響較小，海拔影響大于溫度影響。在海拔相對較低的農區牦牛酸奶中乳球菌屬和乳桿菌屬占優勢，且乳球菌屬含量高于鏈球菌屬和乳桿菌屬；而在海拔相對較高的牧區牦牛酸奶中乳桿菌屬和鏈球菌屬占優勢，且乳桿菌屬占比顯著高于乳球菌屬和鏈球菌屬。乳桿菌屬對酸的耐受力一直被認為高于鏈球菌屬和乳球菌屬，在發酵過程中進行微弱的水解作用，有助于維持較高的發酵活性[33]。在西藏牦牛酸奶的傳統發酵方式中，對于溫度及發酵時間沒有精確的控制，由于牧區平均海拔高于農區，年平均氣溫相對較低，酸奶在牧區發酵時間略長于農區，因此導致牧區酸奶中耐酸并能提高發酵活性的乳桿菌含量高于農區酸奶。

生活在青藏高原的藏族人群，傳統的飲食習慣以糌粑及肉類為主，果蔬類及其他蛋白質攝入量較少，飲食結構較為單一，牦牛酸奶作為其最常食用的乳制品之一[4]，對調節膽固醇含量、維持機體健康起到了重要作用[5,31-33]。由于青藏高原海拔差異較大，不同地域環境下的藏族人群飲食習慣也存在很大差異。人類飲食對腸道菌群有直接影響，而腸道菌群對人的健康更是發揮著至關重要的作用[34]。已有研究表明，藏族人群腸道菌群多樣性在不同海拔[35]、不同生存環境[36]下有一定差異。西藏牦牛酸奶及其他發酵食品中的菌群對藏族人群腸道菌群及其對健康的影響機制尚不明確。因此，對青藏高原傳統食品中菌群的進一步分析將有助于揭示自然環境與傳統食品中菌群的關系以及環境、食品與不同生活方式下藏族人群腸道菌群三者之間的適應性關系。