甘南地區牦牛曲拉中細菌群落結構

曹 磊，梁春御，曹瑛瑛，文開勇，文鵬程，楊 敏，馮曉蘶，張忠明,*，張衛兵

（1.甘肅農業大學食品科學與工程學院，甘肅 蘭州 730070；2.蘭州資源環境職業技術學院氣象學院，甘肅 蘭州 730021；3.甘肅食品藥品監督管理局，甘肅 蘭州 730070；4.甘肅農業大學理學院，甘肅 蘭州 730070）

牦牛曲拉，又稱奶渣，是將牦牛乳脫脂后，在自然條件下進行發酵使酪蛋白凝結、干燥后制成的一種發酵乳制品[1]。與新鮮牦牛乳比較，曲拉具有蛋白質含量較高、易貯藏、方便運輸的特點[1-2]。曲拉不僅是藏區牧民貯藏食品蛋白的重要手段，還是制作酸奶的發酵劑。另外，曲拉還可作為生產干酪素的原料[3]。由于曲拉加工環境相對開放，因此環境中的多種微生物參與了發酵過程[4]。

目前，微生物多樣性的分析方法主要包括傳統純培養方法和基于分子生物學技術的非純培養的方法。純培養的方法只能在實驗室條件下對可培養的微生物進行分析，不能全面了解微生物多樣性；非純培養的方法包括變性梯度凝膠電泳[5]、單鏈構象多態性[6]、溫度梯度凝膠電泳[7]、末端限制性片段長度多態性分析[8-10]等，其早期也廣泛應用于乳制品中微生物多樣性的研究[11-13]。然而，近年來新興的Illumina MiSeq高通量測序技術能夠全面準確地對研究對象中微生物種類組成和結構進行分析，相較于傳統的純培養方法及以16S rRNA為基礎的非純培養方法，能夠較多的產生測序覆蓋深度更大的數據量，更全面的呈現其微生物多樣性，明確各類微生物的結構比例及優勢菌群。該技術已在原乳[14]、發酵酸牛奶[9]、發酵酸駝奶和發酵酸馬奶[15]的微生物多樣性研究中得到了成功應用。

本研究采用Illumina MiSeq高通量測序技術對采自于甘肅省甘南藏族自治州的曲拉樣品中細菌多樣性進行分析，以期全面地解析曲拉中的細菌多樣性，系統地了解曲拉中細菌組成及群落結構。

1 材料與方法

1.1 材料與試劑

曲拉樣品采自于甘肅省甘南藏族自治州合作市那吾鄉瑪崗村（MG1、MG2、MG3）、卡斯合村（KS1、KS2、KS3）、紹瑪村（SM1、SM2、SM3）3 個村莊的9 個不同牧民家庭。將所有樣品置于自封袋中編號，并置于冰盒中運輸至實驗室以備實驗。

E.Z.N.A.Soil DNA試劑盒 美國Omega公司；Qubit 2.0 DNA檢測試劑盒 美國Invitrogen公司；Q5高保真DNA聚合酶 美國New England Biolabs公司；凝膠回收試劑盒 美國Axygen公司；TruSeq Nano DNA LT Library Prep Kit 美國Illumina公司。

1.2 儀器與設備

Pico-21型臺式離心機 美國Thermo Fisher科技公司；DYY-6C型電泳儀、DYCZ-21型電泳槽 北京市六一儀器廠；凝膠成像系統 美國UVP公司；Q32866型Qubit 2.0分光光度計 美國Invitrogen公司；T100TMThermal Cyeler型聚合酶鏈式反應（polymerase chain reaction，PCR）儀 伯樂生命醫學產品（上海）有限公司；MiSeq System SY-410-1003高通量測序儀 美國Illumina公司。

1.3 方法

1.3.1 曲拉微生物總DNA的提取

取0.2 g曲拉樣品，采用E.Z.N.A.Soil DNA Kit D5625-01提取試劑盒，參照說明書從樣品中提取DNA；然后使用Qubit 2.0分光光度計對提取的DNA進行定量，并通過0.8%瓊脂糖凝膠電泳確定DNA提取的完整性。

1.3.2 PCR擴增及測序

以稀釋后的基因組總DNA為模板，取30 ng進行PCR靶向擴增16S rRNA基因。采用V3-V4區338F（5’-ACTCCTACGGGAGGCAGCAG-3’）和806R（5’-GGACTACHVGGGTWTCTAAT-3’）分別為特異性引物。16S rRNA基因PCR熱循環譜如下：95 ℃預變性5 min；95 ℃變性30 s，56 ℃退火30 s，72 ℃延伸30 s，25 個循環；72 ℃退火10 min，最終保持在4 ℃條件下。PCR產物使用2%瓊脂糖凝膠電泳進行檢測，切割后使用Axygen公司的凝膠回收試劑盒回收。將擴增產物純化后，利用分光光度計精確定量后制備測序文庫，最后在MiSeq測序平臺進行雙端測序。PCR擴增、測序文庫制備及測序工作由上海派森諾生物科技股份有限公司完成。

1.3.3 高通量測序數據處理

采用Mothur（version 1.31.2）和QIIME（version 1.7.0）軟件處理及分析高通量測序數據[16-17]。先對測序結果進行整合，使每次讀數的堿基平均質量不小于Q20。然后將整合后的序列使用FLASH（version 1.2.7）軟件根據重疊堿基進行配對連接，并丟棄無法配對的序列[18]。使用UCHIME（version 4.2）軟件鑒定和去除嵌合體序列[19-20]。將序列分配入對應樣本，使用UCLUST算法進行序列聚類，以97%的序列相似度作為劃分閾值，將樣品聚類成可操作分類單位（operational taxonomic units，OTU）[18]。

1.3.4 群落多樣性和統計分析

利用Mothur（version 1.31.2）軟件進行α多樣性分析，其中包括Simpson指數、Chao1指數和Shannon指數，并在不同的分類水平上對群落結構進行了統計分析。

使用QIIME軟件對基于Unweighted和Weighted的UniFrac距離矩陣進行UPGMA聚類分析，并使用R（version 2.15.3）軟件進行可視化。

通過PICRUSt（version 1.0）軟件對群落基因功能特征進行預測[21]。

1.4 數據分析

利用Excel 2007軟件進行數據處理分析并作圖。

2 結果與分析

2.1 PCR擴增、測序序列及OTU基本情況

由圖1可以看出，PCR擴增后得到清晰明顯的條帶，可用于后續的研究分析。通過細菌的16S rRNA基因V3-V4區測序，9 份樣品共產生216 630 條有效序列，將得到的有效序列經過質量控制，在剔除一些疑問序列之后得到高質量序列為185 910 條，高質量序列占有效序列的比例為85.82%。將所有序列按97%的相似度進行OTU聚類，得到4 754 個OTU。由序列數及OTU聚類可以看出，曲拉中細菌種類繁多，物種豐富，且不同來源的曲拉樣品中存在一定差異。

圖1 PCR擴增產物凝膠電泳圖Fig. 1 Electrophoresis profile of PCR amplified product

2.2 稀疏曲線

稀疏曲線用來評判每個樣品在當前的測序深度下是否足以反映該樣品中所包含的細菌群落多樣性。由圖2可以看出，當樣品測序量較低時，每個樣品的Shannon指數呈顯著上升的趨勢，說明在當前測序量上樣品物種的多樣性較高，樣品中還有較多的物種沒有被檢測到。但隨著測序數量逐漸增加，各曲拉樣品中Shannon指數呈緩慢增加的趨勢，稀疏曲線的斜率逐漸降低，并且隨著測序數量的繼續增加，Shannon曲線基本穩定。雖然在當前的測序深度下每個樣品的稀疏曲線未能完全進入平臺期，但與X軸幾乎接近平行，已經達到飽和。因此，盡管隨著測序量的增加新的物種可能會被發現，但在此測序水平下，樣品中細菌的群落多樣性已能夠充分的展現。

圖2 Shannon指數稀疏分析圖Fig. 2 Sparse analysis of Shannon index

2.3 α多樣性分析

曲拉樣品中微生物α多樣性指數如表1所示。對于微生物群落而言，可以用Chao1指數、Shannon指數、ACE指數和Simpson指數反映樣品中物種的豐富度和多樣性。由表1可以看出，曲拉樣本中細菌的Chao1指數平均值為307.06±124.62，ACE指數為473.93±180.14，Simpson指數為0.74±0.13，Shannon指數為3.29±1.35。從微生物多樣性指數結果分析可以得出，曲拉樣品中細菌群落的物種多樣性及總體豐度相對較高，曲拉中的微生物數量處于較高的水平。同時，不同村莊來源的樣品微生物多樣性指數存在差異，這表明樣品制備環境可能會影響樣品中的微生物物種豐富度及多樣性。Matsoni[22]和Tarag[23]對于發酵乳的研究中也發現了相同的結果。在該研究中所有樣品的覆蓋率范圍為99.91%～99.98%，均大于99%，這表明測序水平能夠達到鑒定樣品中的大多數微生物群落的多樣性。

表1 微生物多樣性指數Table 1 Microbial diversity indexes

2.4 樣品細菌群落在門水平的分類及比較

表2 各樣品門水平菌群分布相對豐度Table 2 Relative abundance of bacterial flora distribution in Qula samples at phylum level

對9 份曲拉樣品從門的分類水平進行鑒定，結果如表2所示。在曲拉樣品中共檢測出8 個門，分屬于厚壁菌門（Firmicutes）、變形菌門（Proteobacteria）、放線菌門（Actinobacteria）、酸桿菌門（Acidobacteria）、藍藻門（Cyanobacteria）、擬桿菌門（Bacteroidetes）、綠彎菌門（Chlorof l exi）和假單胞菌門（Gemmatimonadetes）。由表2可知，厚壁菌門為曲拉樣品中的第1優勢門，其相對豐度在44.515%～81.010%的范圍內，平均相對豐度為60.751%；其次為變形菌門，為樣品中的第2優勢門，其相對豐度介于18.973%～55.363%之間，平均相對豐度為39.138%；放線菌門和酸桿菌門豐度較低，平均相對豐度僅為0.064%和0.026%，為樣品中的非優勢門；藍藻門（平均相對豐度為0.012%）、擬桿菌門（平均相對豐度為0.003%）、綠彎菌門和假單胞菌門（平均相對豐度為0.001%）的豐度都很低，為曲拉樣品中的稀有門。厚壁菌門和變形菌門也是牛乳、馬乳、駝乳及其發酵酸乳中的優勢門[15]，與本實驗的研究結果一致。另外，在發酵牛乳[24]、發酵牦牛乳[25]、發酵馬乳[15,26]中均檢測到綠彎菌門、酸桿菌門、藍藻門，且均為非優勢門，而在先前的研究中沒有檢測到假單胞菌門。

2.5 樣品細菌群落在屬水平的分類及比較

表3 各樣品屬水平菌群分布相對豐度Table 3 Relative abundance of bacterial flora distribution in Qula sample at genus level

對9 份曲拉樣品從屬的分類水平進行鑒定，共檢測出55 個屬，結果如表3所示。從屬分類水平來看，在曲拉樣品中檢測出的乳桿菌屬（Lactobacillus）、醋酸桿菌屬（Acetobacter）、乳球菌屬（Lactococcus）、假單胞菌屬（Pseudomonas）和明串珠菌屬（Leuconostoc）為豐度相對較高的屬，其相對豐度均大于5%。厚壁菌門的乳桿菌屬的相對豐度在11.628%～80.946%的范圍內，平均相對豐度為43.027%，為曲拉樣品中的第1優勢屬；變形菌門的醋酸桿菌屬的相對豐度在5.057%～52.408%的范圍內，平均相對豐度為24.739%，為第2優勢屬；厚壁菌門的乳球菌屬的相對豐度相對較高，平均相對豐度為11.402%，為第3優勢屬。假單胞菌屬和明串珠菌屬的平均相對豐度分別為7.971%和5.088%。不同來源的曲拉樣品在屬水平上細菌組成存在差異。另外，樣品中還包含許多在屬水平上未鑒定的屬（unclassified），未鑒定出屬也具有較高的豐度，其相對豐度在1.942%～10.928%的范圍內，平均相對豐度為4.922%。

腸球菌屬（Enterococcus）和葡萄球菌屬（Staphylococcus）為樣品中的低豐度屬，其平均相對豐度分別為0.054%和0.001%。腸球菌屬及葡萄球菌屬屬于食源性病原菌屬，說明曲拉相對粗放的生產環境有一定的影響。

不同乳制品的微生物組成有著豐富的多樣性，如傳統酸奶[23]、發酵牛乳[24,26-27]、發酵馬乳[26]、發酵牦牛乳[25]和開非爾[28]的優勢屬多為乳桿菌屬，與本實驗的研究結果一致。奶酪中的優勢屬各不相同，愛爾蘭奶酪[29]為乳球菌屬，丹麥原奶酪[28]中優勢菌株歸屬于乳桿菌屬和鏈球菌屬。有些發酵乳制品樣品中也檢測到埃希氏桿菌屬（Escherichia）和沙門氏菌屬（Salmonella）[23]，葡萄球菌屬和志賀氏菌屬（Shellogella）[14]等也有檢測到。

2.6 物種豐度熱圖

豐度熱圖顏色越綠則說明所含的菌屬較少。由圖3可以看出，9 個樣本的前5 種菌屬的顯示顏色比較紅，含量較多，分別是乳桿菌屬、醋酸桿菌屬、乳球菌屬、假單胞菌屬、明串珠菌屬。另外，一些沒有與數據庫匹配的未鑒定出屬也具有較高的豐度。不同來源的曲拉樣品菌屬中有的樣品顏色比較紅，有的樣品顏色比較綠，說明每個樣品中菌群的豐度存在差異，因此，不同來源的樣品中菌群多樣性存在差異性。

圖3 屬水平物種豐度熱圖Fig. 3 Heatmap of species abundance at genus level

2.7 基于UniFrac距離的樣本聚類分析

圖4 曲拉中菌群結構UniFrac非加權主坐標分析Fig. 4 Unweighted UniFrac principal coordinate analysis of bacterial communities in Qula samples

采用基于UniFrac距離的Weighted和Unweighted主坐標分析對樣品中細菌群落結構進行比較，結果分別見圖4、5。圖4為基于非加權UniFrac距離的主坐標分析，其第1主成分、第2主成分和第3主成分的貢獻率分別為47.68%、28.60%和4.76%，不同來源的樣品在PC1和PC2維度上明顯的呈聚集和分離趨勢，個別樣品之間存在部分交疊，能夠將不同樣品較好的區分，而在PC3維度上所有樣品可以達到很好區分效果。

圖5為基于加權UniFrac距離的主坐標分析，其第1主成分、第2主成分和第3主成分的貢獻率分別是71.75%、27.93%和0.15%，不同樣品在PC1和PC2維度上有樣品之間存在交疊現象，而在PC3維度上所有樣品可以很好地區分。

由樣品主坐標分析可知，所有樣品沒有緊密地組合在一起，各樣品所代表的點在空間分布有一定的距離，說明不同位置來源的曲拉樣品的細菌群落組成具有差異性。本研究中使用的曲拉樣品是通過傳統方法生產的，因此，樣品間細菌群落結構的差異不同很可能歸因于原料、地理位置和其他環境因素。從圖6可以看出，不同采樣地的樣品聚集到不同的類別，說明不同來源樣品的微生物多樣性存在一定的差異性，這與主坐標分析的結果一致。

圖5 曲拉中菌群結構UniFrac加權主坐標分析Fig. 5 Weighted UniFrac principal coordinate analysis of bacterial communities in Qula samples

圖6 基于Unweighted UniFrac距離矩陣的UPGMA聚類分析圖Fig. 6 Dendrogram from UPGMA cluster analysis based on unweighted UniFrac distance matrix

2.8 樣品中與細菌的基因代謝有關的功能特征

表4 與樣品中細菌的基因代謝有關的功能特征Table 4 Functional features related to genes involved in metabolism from bacteria in Qula samples

表5 樣本中細菌的其他功能基因Table 5 Other functional genes from bacteria in Qula samples

PICRUSt軟件用于預測樣品中存在的細菌類群的功能基因及其代謝途徑[30]。對曲拉樣本中的功能基因及代謝途徑預測情況如表4所示。在樣本中的功能基因中，44.04%～50.37%的微生物基因與以下代謝途徑有關：氨基酸代謝、碳水化合物代謝、能量代謝、核苷酸代謝、輔助因子及維生素代謝、外源物質的生物降解與代謝、脂質代謝、酶家族、糖生物合成與代謝、其他氨基酸的代謝、萜類化合物與多酮類化合物的代謝以及其他次生代謝產物的生物合成。其中，碳水化合物代謝和氨基酸代謝的基因在曲拉樣品中占主導地位，分別占所有基因的10.40%和8.68%。

樣品中存在的細菌的其他功能基因與細胞處理、環境信息處理、遺傳信息處理、人類疾病、機體系統等有關，具體情況如表5所示。其中，涉及膜轉運的基因占主導地位，占所有基因的12.85%。

以不同樣品中前50 種功能基因繪制熱圖，結果如圖7所示。相同來源的曲拉樣品之間的功能基因相似程度較高，而不同來源的樣品的差異性較大，這一結果與前面細菌群落的聚類分析結果一致。

圖7 基于樣本前50 個功能基因的聚類分析熱圖Fig. 7 Heatmap from clustering analysis based on top 50 functional genes in Qula samples

3 結 論

本研究基于Illumina MiSeq高通量測序平臺分析曲拉樣品中細菌群落結構及多樣性。結果表明，不同來源的曲拉樣品中的微生物多樣性存在差異性。與傳統方法相比較，高通量測序的研究結果能夠更全面地反映曲拉中菌群分布多樣性。曲拉中細菌群落組成分析表明，曲拉樣品中優勢門為厚壁菌門和變形菌門；優勢屬為乳桿菌屬、醋酸桿菌屬、乳球菌屬，細菌群落組成表明在不同來源的樣品中群落組成存在差異。樣品中存在的細菌的功能基因預測表明，不同來源樣品的細菌群落之間存在差異。